echo            both
log             ../DEM/log.liggghts append
thermo_log      ../DEM/thermo.txt

#######################################################
# variables                                           #
#######################################################
variable r0 equal 0.00325 # 0.02 # 0.011
#######################################################

# read packed bed and calc convective heat transfer
atom_style      granular
atom_modify     map array
boundary        m m m
newton          off

communicate     single vel yes

units 			si
processors      2 8 2 #4 4 4 #2 4 2 # 1 1 2

# read the restart file
read_restart    ../DEM/restart/liggghts.restart

neighbor        ${r0} bin
neigh_modify    delay 0 binsize 0.022 # 0.01

## Variables 
variable 	theta equal 4 # angle in degrees
variable    IV    equal 1 # insertion velocity in m/s
variable    rpm   equal 3 # rpm of the revolution of rotary kiln
variable 	massrate_value equal 0.017 # massrate to be fed, 60 kg/hr = 10/3600 = 0.017 kg/sec
variable 	feeding_time equal 1 # 3 # time to feed or run the liggghts_init simulation
variable 	YM_P equal 5.e6 # Young's modulus
variable 	PR_P equal 0.45 # Poisson's ratio
variable 	CoR equal 0.3 # Coefficient of restituion
variable 	CoF equal 0.5 # Coefficienct of friction 
variable  	RF equal 0.3 # Rayleigh time factor
variable 	rho equal 4000 # particle density
variable 	mps equal 0.0003 # minimum particle radius

## Calculations
variable    cos_angle equal cos(${theta}*PI/180)
variable    sin_angle equal sin(${theta}*PI/180)
variable 	IV_z 	  equal -${IV}*${cos_angle}
variable 	IV_y 	  equal ${IV}*${sin_angle}
#variable 	mass_value equal 1

# Material properties required for granular pair styles
fix         m1 all property/global youngsModulus peratomtype ${YM_P}
fix         m2 all property/global poissonsRatio peratomtype ${PR_P}
fix         m3 all property/global coefficientRestitution peratomtypepair 1 ${CoR}
fix         m4 all property/global coefficientFriction peratomtypepair 1 ${CoF}

# pair style
pair_style  gran model hertz tangential history # Hertzian without cohesion
pair_coeff  * *

print "******************************************timestep order calculation purpose only****************************************************"
variable 	SM equal ${YM_P}/2.0/(1+${PR_P})	# Shear modulus
variable 	timestep_size equal ${RF}*PI*(${mps})*sqrt(${rho}/${SM})/(0.163*${PR_P}+0.8766)	# fraction of Rayleigh timestep
print 		"timestep = ${timestep_size}"

# timestep, gravity
variable 	timestep_size equal 0.00001
timestep	${timestep_size}
fix         gravi all gravity 9.81 vector 0.0 ${sin_angle} -${cos_angle}

# walls
fix 		ywalls1 all wall/gran model hertz tangential history primitive type 1 yplane 0.
#fix    	ycylwall all wall/gran model hertz tangential history primitive type 1 ycylinder 0.1 0.0 0.0 shear y 0.5
# Please see shear option to rotate cylindrical wall with given vshear magnitude (+ve vshear for clockwise rotation and otherwise)
# https://www.cfdem.com/media/DEM/docu/fix_wall_gran.html?highlight=xplane

#Mesh
fix  		rotaryKiln all mesh/surface file ../DEM/geometry/Rotary_Kiln.stl type 1 
#fix  		ins_mesh all mesh/surface/planar file geometry/ins_face.stl type 1 scale 1.0
fix  		granwalls all wall/gran model hertz tangential history mesh n_meshes 1 meshes rotaryKiln

# heat transfer
fix         ftco all property/global thermalConductivity peratomtype 0.89 	# lambda in [W/(K*m)] # (1.2+0.58)/2
fix         ftca all property/global thermalCapacity peratomtype 775   	 	# cp in [J/(kg*K)]
fix         heattransfer all heat/gran initial_temperature 673 # 600. # (980+570)/2
# Ref: 
# 1. Heat transfer analysis of iron ore spherical pellet in the direct reduction process
# 2. Effect of iron (III) oxide powder on thermal conductivity and diffusivity of lime mortar

# set particle temperature for the bed
run         0
region     total block INF INF INF INF INF INF units box
set        region total property/atom Temp 673.

# cfd coupling
fix         cfd all couple/cfd couple_every 100 mpi
fix         cfd2 all couple/cfd/force

# convection and initial temperature
fix         cfd3 all couple/cfd/convection T0 673 # 600

# Variable total mass
variable 		mass_total equal mass(all)

# apply nve integration to all particles that are inserted as single particles
fix         integr all nve/sphere

# output settings, include total thermal energy
compute         rke all erotate/sphere
thermo_style    custom step atoms v_mass_total ke c_rke f_heattransfer vol
thermo          1000
thermo_modify   lost ignore norm no
compute_modify  thermo_temp dynamic yes

# dump commands
dump        	dmp all custom 10000 ../DEM/post/dump*.liggghts_run id type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius f_Temp[0] f_heatFlux[0]
dump   			dumpKiln all mesh/stl 10000 ../DEM/post/kiln*.stl

# rotating the rotary kiln
variable 	time_period equal 60/${rpm}
fix			moveKiln all move/mesh mesh rotaryKiln rotate origin 0. 0. 0. axis 0. 1. 0. period ${time_period}

run         0