# Particle packing by insertion by gravity # DST packing, with respect to Ni et al. 2000, Leighton Buzzard sand, Shearbox 60 mm x 60 mm x 20 mm # Constant velocity can be applied to move the plane and box # Constant settings atom_style granular atom_modify map array boundary f f f # will change to f f f boundary to save time newton off #echo both communicate single vel yes units si processors 2 2 1 read_restart liggghts.packing_gravity1 neighbor 0.00064 bin neigh_modify delay 0 #Material properties required for new pair styles fix m1 all property/global youngsModulus peratomtype 2.6e8 2.6e8 fix m2 all property/global poissonsRatio peratomtype 0.3 0.3 fix m3 all property/global coefficientRestitution peratomtypepair 2 0.2 0.2 0.2 0.2 fix m4 all property/global coefficientFriction peratomtypepair 2 0.5 0.36 0.36 0.5 fix m5 all property/global coefficientRollingFriction peratomtypepair 2 0.2 0.36 0.36 0.2 #fix m6 all property/global characteristicVelocity scalar 2 fix m7 all property/global cohesionEnergyDensity peratomtypepair 2 0 0 0 0 #New pair style pair_style gran model hertz tangential history cohesion sjkr rolling_friction epsd2 pair_coeff * * timestep 1e-6 # 10% tcrit ## Variables for moving upper_plane meshes and apply servo stresses variable sigma_1 equal 100000 # Normal stress starting with 100kPa variable sample_width equal 0.03 # m variable sample_height equal 0.01 # m variable f_top equal ${sigma_1}*${sample_width}*${sample_width} # equal to 45 N variable nstep1 equal 10000 variable vel equal 10 ## Insert CAD models fix bot_wall all mesh/surface/stress file meshes/plane1.stl type 2 scale 0.01 stress on fix top_wall all mesh/surface/stress/servo file meshes/plane7.stl type 2 scale 0.01 com 0.0 0.0 0.0095 ctrlPV force axis 0.0 0.0 1.0 target_val -${f_top} vel_max ${vel} fix botx_plane1 all mesh/surface/stress file meshes/plane2.stl type 2 scale 0.01 stress on fix botx_plane2 all mesh/surface/stress file meshes/plane3.stl type 2 scale 0.01 stress on fix boty_plane1 all mesh/surface/stress file meshes/plane4.stl type 2 scale 0.01 stress on fix boty_plane2 all mesh/surface/stress file meshes/plane5.stl type 2 scale 0.01 stress on fix left_wing all mesh/surface file meshes/plane6.stl type 2 scale 0.01 fix right_wing all mesh/surface file meshes/plane12.stl type 2 scale 0.01 fix topx_plane1 all mesh/surface/stress file meshes/plane8.stl type 2 scale 0.01 stress on fix topx_plane2 all mesh/surface/stress file meshes/plane9.stl type 2 scale 0.01 stress on fix topy_plane1 all mesh/surface/stress file meshes/plane10.stl type 2 scale 0.01 stress on fix topy_plane2 all mesh/surface/stress file meshes/plane11.stl type 2 scale 0.01 stress on fix wall all wall/gran model hertz tangential history mesh n_meshes 12 meshes bot_wall top_wall botx_plane1 botx_plane2 boty_plane1 boty_plane2 left_wing right_wing topx_plane1 topx_plane2 topy_plane1 topy_plane2 # Output forces on walls fix force1 all ave/time ${nstep1} 1 ${nstep1} f_top_wall[1] f_top_wall[2] f_top_wall[3] f_top_wall[7] f_top_wall[8] f_top_wall[9] file graph/force_twall.dat fix force2 all ave/time ${nstep1} 1 ${nstep1} f_bot_wall[1] f_bot_wall[2] f_bot_wall[3] f_bot_wall[7] f_bot_wall[8] f_bot_wall[9] file graph/force_bwall.dat fix force3 all ave/time ${nstep1} 1 ${nstep1} f_botx_plane1[1] f_botx_plane1[2] f_botx_plane1[3] f_botx_plane1[7] f_botx_plane1[8] f_botx_plane1[9] file graph/force_botx_plane1.dat fix force4 all ave/time ${nstep1} 1 ${nstep1} f_botx_plane2[1] f_botx_plane2[2] f_botx_plane2[3] f_botx_plane2[7] f_botx_plane2[8] f_botx_plane2[9] file graph/force_botx_plane2.dat fix force5 all ave/time ${nstep1} 1 ${nstep1} f_boty_plane1[1] f_boty_plane1[2] f_boty_plane1[3] f_boty_plane1[7] f_boty_plane1[8] f_boty_plane1[9] file graph/force_boty_plane1.dat fix force6 all ave/time ${nstep1} 1 ${nstep1} f_boty_plane2[1] f_boty_plane2[2] f_boty_plane2[3] f_boty_plane2[7] f_boty_plane2[8] f_boty_plane2[9] file graph/force_boty_plane2.dat #apply nve integration to all particles that are inserted as single particles fix integr all nve/sphere fix gravi all gravity 9.81 vector 0.0 0.0 -1.0 #output settings, include total thermal energy compute 1 all erotate/sphere thermo_style custom step atoms ke c_1 vol thermo 10000 thermo_modify lost ignore norm no # Compute contact forces for the measurement region region measurement_reg block -0.01 0.01 -0.01 0.01 0.001 0.007 units box group measured_group region measurement_reg compute c1 all pair/gran/local pos id force delta # Atom stresses compute peratom all stress/atom compute p1 all reduce sum c_peratom[1] c_peratom[2] c_peratom[3] variable press equal -(c_p1[1]+c_p1[2]+c_p1[3])/(3*vol) #insert the first particles run 1 dump dmp all custom/vtk 10000 post1_shearing_12walls/Particle/packing_*.vtk id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius #dump dmp_upperbox all mesh/stl 10000 post1_shearing/Mesh_Box/Upperbox*.stl upperbox dump dmp_lowerbox all mesh/stl 10000 post1_shearing/Mesh_Box/Lowerbox*.stl botx_plane1 botx_plane2 boty_plane1 boty_plane2 left_wing bot_wall dump dmp_stress all mesh/gran/VTK 10000 post1_shearing_12walls/Mesh_servo/mesh_top_*.vtk stress top_wall ## Dump contact information dump dmp_contactforce all local 10000 post1_shearing_12walls/Contact_forces/Contact*.vtk c_c1[1] c_c1[2] c_c1[3] c_c1[4] c_c1[5] c_c1[6] c_c1[7] c_c1[8] c_c1[9] c_c1[10] c_c1[11] c_c1[12] c_c1[13] #dump dmp_contactforce all local/gran/vtk 10000 post1_shearing/Contact_forces/Contact*.vtk c1 run 200000 fix movemesh1 all move/mesh mesh botx_plane1 linear 0.002 0 0 fix movemesh2 all move/mesh mesh botx_plane2 linear 0.002 0 0 run 800000 write_restart liggghts.packing1_loading