DEM only simulation - Add dragforce - Recirculation

Hi Matt,
Thank you for your quick response and its really heartening to see senior members helping new learners.
I have included LIGGGHTS script below for some details. I would try to use "periodic" boundary condition, as you suggested for recirculation of particles.
For grid size, I would say I have no preference or pardon me if I am not able to find a link between porosity calculation and grid size value. I have attached picture of domain that is 1.5m in diameter and around 3 m in height, and particle diameter is around 60mm. I would appreciate if guidance may be provided on selecting grid size and finding porosity value which would help in finding realistic drag force to apply.
Thank you again.
Sohaib

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

variable radius1 equal 0.0300 #m

variable frac1 equal 1.0 #100%

variable density equal 4000 #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 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 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 pdd1 all particledistribution/discrete 32452867 ${nradii} pts1 ${frac1}

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

#Adding Fluid Drag
variable Alpha equal (1-0.66)^(2.0/3.0)
variable Uinf equal 1.03
variable U equal v_Uinf/v_Alpha
variable Rho equal 1.65
variable Mu equal 3.86e-5
variable Diameter equal 0.06
variable Across equal 1.0/6.0*v_Diameter^3*3.1415
variable Vrx atom -vx
variable Vry atom -vy
variable Vrz atom v_U-vz
variable Vr atom sqrt(v_Vrx^2+v_Vry^2+v_Vrz^2)
variable Re atom v_U*v_Rho*v_Diameter/v_Mu #variable Re atom v_Vr*v_Rho*v_Diameter/v_Mu
variable Cd atom 0.44 #variable CdV atom (1+0.15*v_Re^0.687)*24/v_Re
#variable Cd atom 0.44*(v_Cdv<0.44)+v_Cdv*(v_Cdv>=0.44)
#variable Dragx atom 0.5*v_Rho*v_Vr*v_Cd*v_Across*v_Vrx
#variable Dragy atom 0.5*v_Rho*v_Vr*v_Cd*v_Across*v_Vry
variable Dragz atom 0.5*v_Rho*abs(v_Vrz)*v_Cd*v_Across*v_Vrz*9.109
fix adddrag1 all addforce 0 0 v_Dragz
#fix DragCustom all print 1000 "${v_Dragx} ${v_Dragy} ${v_Dragz}" screen no title "Dragx Dragy Dragz" file DragCustomf.txt

fix surface_mf all mesh/surface file Insertionsface.stl type 1 scale 0.001 move 0.0 0.0 0.0
fix surface_mf2 all mesh/surface file Lid.stl type 2 scale 0.001 move 0.0 0.0 0.0
fix mymassflow all massflow/mesh mesh surface_mf count multiple vec_side 0.0 0.0 -1.0 writeTime file massflow_particle_data.txt
fix mymassflow2 all massflow/mesh mesh surface_mf2 count multiple vec_side 0.0 0.0 -1.0 writeTime file massflow_particle_data2.txt
variable time equal step*dt # current simulation time (step and dt are built-in variables)
variable m equal f_mymassflow[1] # total mass that has crossed the mesh since simulation start
variable np equal f_mymassflow[2] # total number of particles that have crossed the mesh since simulation start
variable mdot equal f_mymassflow[3] # current mass flow rate
variable npdot equal f_mymassflow[4] # current particle flow rate

fix pmassout all print 1000 "${time} ${m} ${np} ${mdot} ${npdot}" screen no title "t mass particles massflow particlerate" file massflow.txt

variable time2 equal step*dt # current simulation time (step and dt are built-in variables)
variable m2 equal f_mymassflow2[1] # total mass that has crossed the mesh since simulation start
variable np2 equal f_mymassflow2[2] # total number of particles that have crossed the mesh since simulation start
variable mdot2 equal f_mymassflow2[3] # current mass flow rate
variable npdot2 equal f_mymassflow2[4] # current particle flow rate

fix pmassout all print 1000 "${time2} ${m2} ${np2} ${mdot2} ${npdot2}" screen no title "t mass particles massflow particlerate" file massflow2.txt

###################################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 type x y z vx vy vz fx fy fz radius mass v_Dragz
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

####################################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}