###################################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 -0.25
variable xmax equal 0.25

variable ymin equal -0.4
variable ymax equal 0.4

variable zmin equal 0
variable zmax equal 2.3

#Definition of the timestep

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


































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

variable natoms equal 2 	#1 -> particle #2->loputa, razvod

####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.36
variable	sf21	equal	0.36
variable	sf22	equal	0.0

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





####variable for particle####

#Number of particle radius
variable	nradii	equal	2

variable	radius1	equal	0.0030	#m
variable	radius2	equal	0.0030	#m

variable	frac1	equal	0.2		#20%
variable	frac2	equal	0.8		#80%

variable	density	equal	7800	#kg/m³









































####filling parameters####

variable	filltime	equal	3		#seconds

variable	fillmass	equal	80	#kg

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

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

####settle time####

variable	settletime	equal	2 #second

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

####open parameter of the loputa####

variable	openvel	equal	-0.4	#m/s

variable	opentime	equal	0.5 #second

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


####discharge time#####

variable	dischargetime	equal	12 #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 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}



















###################################Generation and Loading of the Geometry .stl#######################################

fix			razvod 	all mesh/surface file Razvod.stl 	type 2  scale 0.001		#load mesh from STL file. Type 2 for geometry. Scale down to transform mm to meters
fix			loputa   	all mesh/surface file Loputa.stl 	type 2  scale 0.001

fix	MesureMesh all mesh/surface/planar file mesureface.stl type 1  scale 0.001    	#load mesh for mesurment



fix 		walls 	all wall/gran model hertz tangential history rolling_friction epsd2 mesh n_meshes 2 meshes razvod loputa








































###################################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 1 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 NInsertionsface.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

		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		







###################################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/vtk ${dumpstep} post/particles_*.vtk id radius mass vx vy vz

dump	dmprazvod all mesh/stl ${dumpstep} post/Razvod*.stl razvod

dump	dmploputa all mesh/stl ${dumpstep} post/Loputa*.stl loputa






###################################Load Mesh and active mesh for the mesurment ############################

fix	MesureMesh all mesh/surface/planar file mesureface.stl type 1 scale 0.001    	#load mesh for mesurment

fix	mesurePipe all massflow/mesh mesh MesureMesh vec_side 0. 0. -1			#Activate mesurment mesh


###################################SAVE data of mesurePipe to variable ##################################

variable	MassTotalPipe		equal f_mesurePipe[1]	#field 1 ->sum total mass that has been passed [kg]

variable	NoPTotalPipe		equal f_mesurePipe[2]	#field 2 ->sum total NoP that has been passed [particles]

variable	MassflowRatePipe	equal f_mesurePipe[3]	#field 3 ->actual mass Flow rate [kg/s]

variable	NoPflowRatePipe	equal f_mesurePipe[4]	#field 4 ->actual particle flow rate  [particle/s]


###################################Write Data of Mesurment###############################################

variable outputtime equals 0.05 #s

variable outputsteps equal ${outputtime}/${dt}

variable t equal step*dt

fix output2 all print ${outputsteps} "$t,${MassTotalPipe},${NoPTotalPipe},${MassflowRatePipe},${NoPflowRatePipe}" screen no file PipeMesurments.csv title "Time,NoPTotalPipe,MassTotalPipe,NoPflowRatePipe,MassflowRatePipe" #


###################################Za tem je samo še RUN simulation######################################





































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

run	 ${fillsteps}

unfix	ins

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

run	 ${settlesteps}

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

run	 ${opensteps}

unfix	MoveLoputa


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


run	 ${dischargesteps}