###################################Header for General commands#######################################################

atom_style	granular		#Simulation of particles

boundary	f f f			#fixed boundaries -> particles will be deleted if leaving the simulation box

units		si

communicate	single vel yes	#default

newton		off				#default


###################################System variables##################################################################

#Definition of boundaries
variable xmin equal -2.55
variable xmax equal 2.55

variable ymin equal -2.55
variable ymax equal 2.55

variable zmin equal -0.05
variable zmax equal 5.0

#Definition of the timestep

variable dt equal 1e-4			#timestep = 0.0001 second; Each iteration step represents 0.0001 seconds.


#################################Specific variables for current simulation#########################################

variable natoms equal 2 	#1 -> particle #2->hopper,frame and ground, lid

####variable for material properties####

####Young Modulus####
variable	youngmodulus1	equal 	1e8		#N/mm²
variable	youngmodulus2	equal 	1e8		#N/mm²

####Poission ratio####
variable 	poission1	equal	0.3
variable 	poission2	equal	0.3

####variable for contact properties####

####coefficient of restitution####
variable	CoR11	equal	0.6
variable	CoR12	equal	0.6	
variable	CoR21	equal	0.6
variable	CoR22	equal	0.6

####sliding friction coefficient####
variable	sf11	equal	0.3
variable	sf12	equal	0.3
variable	sf21	equal	0.3
variable	sf22	equal	0.3

####rolling friction coefficient####
variable	rf11	equal	0.8
variable	rf12	equal	0.8
variable	rf21	equal	0.8
variable	rf22	equal	0.8

###################################Definition of cohesionEnergy Density#############

variable	CeD11	equal	100000
variable	CeD12	equal	100000
variable	CeD21	equal	100000
variable	CeD22	equal	100000





####variable for particle####

#Number of particle radius
variable	nradii	equal	2

variable	radius1	equal	0.02	#m
variable	radius2	equal	0.04	#m

variable	frac1	equal	0.5		#20%
variable	frac2	equal	0.5		#80%

variable	density	equal	2000	#kg/m³

####filling parameters####

variable	filltime	equal	5		#seconds

variable	fillmass	equal	5000	#kg

variable	fillmassrate	equal	${fillmass}/${filltime}	#kg/s

variable	fillsteps	equal	${filltime}/${dt}		#Transform time to iteration steps

####settle time####

variable	settletime	equal	1 #second

variable	settlesteps	equal	${settletime}/${dt}		#Transform time to iteration steps

####open parameter of the lid####

variable	openvel		equal	0.4	#m/s

variable	opentime	equal	1 #second

variable	opensteps	equal	${opentime}/${dt}		#Transform time to iteration steps


####discharge time#####

variable	dischargetime	equal	25 #seconds

variable	dischargesteps	equal	${dischargetime}/${dt}	#Transform time to iteration steps


###################################Definition of simulationbox#######################################################

region		reg block ${xmin} ${xmax} ${ymin} ${ymax} ${zmin} ${zmax} units box

create_box	2 reg

neighbor	0.004 bin		#default

neigh_modify	delay 0		#default

#################################Definition of the contact models####################################################
#pair_style	gran model hertz tangential history cohesion sjkr #Hertzian with cohesion 
pair_style gran model hertz tangential history rolling_friction epsd2	#contact model

pair_coeff	* *				#default

timestep ${dt}

fix		integrator all nve/sphere			#default

fix		gravi all gravity 9.81 vector 0.0 0.0 -1.0		#gravity of 9.81 m/s² in negative z direction

###################################Definition of Material properties#################################################


fix 		m1 all property/global youngsModulus peratomtype ${youngmodulus1} ${youngmodulus2}

fix 		m2 all property/global poissonsRatio peratomtype ${poission1} ${poission2}

fix 		m3 all property/global coefficientRestitution peratomtypepair ${natoms} 0.6 0.6 0.6 0.6

fix 		m4 all property/global coefficientFriction peratomtypepair ${natoms} ${sf11} ${sf12} ${sf21} ${sf22}

fix 		m5 all property/global coefficientRollingFriction peratomtypepair ${natoms} ${rf11} ${rf12} ${rf21} ${rf22}

#fix m6 all property/global characteristicVelocity scalar 2. # Required for hooke model
#fix m7 all property/global cohesionEnergyDensity peratomtypepair ${natoms} ${CeD11} ${CeD12} ${CeD21} ${CeD22} # <--- This is the density for cohesion energy 
###################################Generation and Loading of the Geometry .stl#######################################

fix			silo 	all mesh/surface file Silo.stl 	type 2  scale 0.001		#load mesh from STL file. Type 2 for geometry. Scale down to transform mm to meters
fix			lid   	all mesh/surface file Lid.stl 	type 2  scale 0.001
fix			ground 	all mesh/surface file Ground.stl type 2  scale 0.001 
fix			frame 	all mesh/surface file Frame.stl  type 2  scale 0.001


fix 		walls 	all wall/gran model hertz tangential history rolling_friction epsd2 mesh n_meshes 4 meshes silo lid ground frame

###################################Generation and Insertion of the particles#########################################

fix		pts1 all particletemplate/sphere 10487 atom_type 1 density constant ${density} radius constant ${radius1}	

fix		pts2 all particletemplate/sphere 11887 atom_type 2 density constant ${density} radius constant ${radius2}	

fix 	pdd1 all particledistribution/discrete 32452867 ${nradii} pts1 ${frac1} pts2 ${frac2} 

fix		ins_mesh all mesh/surface/planar file Insertionsface.stl type 1 scale 0.001 

fix 	ins all insert/stream seed 86028157 distributiontemplate pdd1 &
        mass ${fillmass} massrate ${fillmassrate} overlapcheck yes all_in yes vel constant 0 0 -0.75 &
        insertion_face ins_mesh extrude_length 0.25

		
########################################Compute Stress ######################################################

#compute	peratom all stress/atom
#compute 	all mobile stress/atom
#compute 	thermo_press all pressure thermo_temp

###################################Dumping of the data for post-processing to visualize############################

shell 	mkdir post

#Definition of the dumptime

variable dumptime equal 0.05	# Every 0.05 s 1 image 

variable dumpstep equal ${dumptime}/${dt}	#Transform to iteration steps


dump	dmpparticle all custom ${dumpstep} post/particles_*.liggghts id type x y z vx vy vz fx fy fz radius mass
dump	dmpground all mesh/stl ${dumpstep} post/Ground*.stl ground
dump	dmpsilo all mesh/stl ${dumpstep} post/Silo*.stl silo 
dump	dmpframe all mesh/stl ${dumpstep} post/Frame*.stl  frame
dump	dmplid all mesh/stl ${dumpstep} post/Lid*.stl lid
dump   dmp all custom/vtk 1000 post1/chute_*.vtk id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius 
dump  dumpstress all mesh/gran/VTK 1000 post1/mesh_*.vtk stress wear
dump 	dmp0 all custom 1000 all/dump*.data id type x y z vx vy vz fx fy fz


####################################RUN the simulation filling###########################################################

run	 ${fillsteps}

unfix	ins

##################################RUN the simulation settling#############################################################

run	 ${settlesteps}

#################################RUN the simulation Open Door#############################################################
fix		MoveLid	all move/mesh mesh lid	linear ${openvel} 0. 0.

run	 ${opensteps}

unfix	MoveLid

#################################RUN the simulation Outflow###############################################################


run	 ${dischargesteps}