units si # SI atom_style sphere # spheric particles boundary p f f # peiodic boundary condition in x-direction newton off communicate single vel yes hard_particles yes # Computational domain region reg block -0.01 0.01 -0.005 0.005 -0.001 0.2 units box create_box 3 reg neighbor 0.01 bin neigh_modify delay 0 # Particle properties fix m1 all property/global youngsModulus peratomtype 9.7e10 9.7e10 9.7e10 # Young's modulus fix m2 all property/global poissonsRatio peratomtype 0.078 0.078 0.078 # Poisson's ratio fix m3 all property/global coefficientRestitution peratomtypepair 3 0.051 0.051 0.051 0.051 0.051 0.051 0.051 0.051 0.051 # restitution coefficient fix m4 all property/global coefficientFriction peratomtypepair 3 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 # friction coefficient pair_style gran model hertz tangential history # Hertz's contacting model pair_coeff * * # Particle distribution fix pts1 all particletemplate/sphere 10007 atom_type 1 density constant 2650 radius constant 0.0005 fix pts2 all particletemplate/sphere 10009 atom_type 2 density constant 2650 radius constant 0.00075 fix pts3 all particletemplate/sphere 10037 atom_type 3 density constant 2650 radius constant 0.001 fix pdd1 all particledistribution/discrete 10039 3 pts1 0.25 pts2 0.5 pts3 0.25 # Region of insertion fix ins_mesh all mesh/surface/planar file particle.stl type 1 scale 0.001 fix y_wall_1 all wall/gran model hertz tangential no_history primitive type 1 yplane -0.005 # Frictionless wall, perpendicular to y-axis fix y_wall_2 all wall/gran model hertz tangential no_history primitive type 2 yplane 0.005 # Frictionless wall, perpendicular to y-axis fix cad1 all mesh/surface/stress file bottom.stl type 1 scale 0.001 stress on fix geometry all wall/gran model hertz tangential history mesh n_meshes 1 meshes cad1 fix integrator all nve/sphere # Gravitational acceleration fix grav all gravity 9.81 vector 0.0 0.0 -1.0 timestep 1.0e-7 # Preserve the state of caluclation and particles thermo_style custom step atoms ke cpu thermo 100000 thermo_modify lost ignore norm no #---------------------------------------------------------------- #---------------------------------------------------------------- # step : time-steps # atoms : number of particles which have been already inserted # ke : sum of kinematic energy of all particles # cpu : elapsed CPU time (s) #---------------------------------------------------------------- #---------------------------------------------------------------- # Insertion of particles fix ins1 all insert/stream seed 10061 distributiontemplate pdd1 nparticles 1700 particlerate 10000 overlapcheck yes vel constant 0. 0. -2. all_in yes insertion_face ins_mesh extrude_length 0.1 # Output the state of particles dump dmp all custom 100000 post_w/dump/dump*.shearcell id type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius #----------------------------------------------------------------------------------------------------------------- #----------------------------------------------------------------------------------------------------------------- # id : id of the particle # type : an integer defined as 'atom_type' before # x, y, z : coordinate of the particle (m) # ix, iy, iz : how many times the particle stepped across periodic boundaries in positive direction # vx, vy, vz : velocity of the particle (m/s) # fx, fy, fz : force acting on the particle (N) # omegax, omegay, omegaz: angular velocity of the particle (/s) # radius : radius of the particle (m) # # For more details, please see online documentation for LIGGGHTS(R)-PUBLIC. # https://www.cfdem.com/media/DEM/docu/Manual.html #----------------------------------------------------------------------------------------------------------------- #----------------------------------------------------------------------------------------------------------------- # Run (Particle insertion) run 3000000 # Change restitution coefficient from 0.051 to 0.5 unfix m3 fix m3 all property/global coefficientRestitution peratomtypepair 3 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 # Run (Change of restitution coefficient) run 3000000 # Prepare a plate on which vertical force acts (stress plate) fix cad2 all mesh/surface/stress/servo file stress.stl type 1 scale 0.001 stress on com 0 0 0.03 axis 0. 0. 1. ctrlPV force target_val -5000 vel_max 3.0 unfix geometry fix geometry all wall/gran model hertz tangential history mesh n_meshes 2 meshes cad1 cad2 # Preserve the state of caluclation and particles thermo_style custom step atoms ke cpu f_cad2[1] f_cad2[2] f_cad2[3] f_cad2[4] f_cad2[5] f_cad2[6] thermo 1000000 #--------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------- # f_cad2[1] f_cad2[2] f_cad2[3]: total force acting on the staress plate (N) # f_cad2[4] f_cad2[5] f_cad2[6]: total torque by particles on the stress plate (N m) #--------------------------------------------------------------------------------------- #--------------------------------------------------------------------------------------- # Output time-step and the state of plates variable n equal step # time-steps variable x equal f_cad2[1] # x-th component of force acting on cad2 variable y equal f_cad2[3] # z-th component of force acting on cad2 variable z equal f_cad2[9] # z-th corrdinate of cad2 variable u equal f_cad1[1] # x-th component of force acting on cad1 variable v equal f_cad1[3] # z-th component of force acting on cad1 fix thermo all print 100000 "$n $x $y $z $u $v" file thermo_w # Output these values on the file 'thermo_w' every 100000 time-steps # Output state of the stress plate (cad2) dump dmpstl all mesh/stl 100000 post_w/lid/lid*.stl cad2 # Make restart-file every 100000 time-steps write_restart restart/restart*.shearcell restart 100000 restart/restart*.shearcell # Run (Pre-compaction) run 6000000