Multisphere particles disappearing unexpectedly after delete_atoms command

Submitted by mattkesseler on Mon, 02/12/2018 - 19:12

Hi all. I am currently adapting an old script of mine, replacing spherical particles being released down a ramp with multisphere particles. This has worked for the most part, but when I try to delete some particles from the simulation when the slide mass has settled and is ready to be released from the shutter, more and more particles keep disappearing between timesteps and some particles become fixed in place. The delete_atoms command successfully deletes the particles outside the region I'm interested in (and any spheres connected to those particles in their respective multisphere). The particles that disappear afterwards seem to disappear without reason or logic from the region I'm interested in. I suspect it may have something to do with the periodic boundary condition that is being used. I've noticed that this issue occurs even if I have fixed sidewalls and the periodic boundary is several particle diameters from those walls. The issue only seems to disappear when the distance between the periodic boundaries is twice that of the distance between the fixed sidewalls. I've attached the script I've been using:

log log.rayleighmidflat_multi

atom_style granular
atom_modify map array
boundary p f f
newton off
#processors 96 1 1

communicate single vel yes

units si

region reg block 0.2 0.3 -1.60 1.20 -0.01 1.10 units box
create_box 3 reg

variable update equal 547
variable break equal 5470

neighbor 0.027 bin
neigh_modify delay 0 every 200 check no

fix m1 all property/global youngsModulus peratomtype 65.e9 210.e9 210.e9
fix m2 all property/global poissonsRatio peratomtype 0.25 0.25 0.25
fix m3 all property/global coefficientRestitution peratomtypepair 3 0.893 0.893 0.893 0.893 0.893 0.893 0.893 0.893 0.893
fix m4 all property/global coefficientFriction peratomtypepair 3 0.5 0.5 0.45 0.5 0.5 0.5 0.45 0.5 0.5
fix m5 all property/global coefficientRollingFriction peratomtypepair 3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3
fix m6 all property/global coefficientRollingViscousDamping peratomtypepair 3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3
fix m7 all property/global coeffRollingStiffness scalar 0.5
fix color all property/atom color scalar yes no no 0

pair_style gran model hertz tangential history rolling_friction epsd
pair_coeff * *

variable tstep equal 1.292699783e-6
timestep ${tstep}

fix gravi all gravity 9.81 vector 0.0 0.0 -1.0

fix mesh1 all mesh/surface file meshes/mesh1_periodic.stl element_exclusion_list read mesh1fixed.stl type 2 move 0.4 -3.08981 -0.003 scale 0.5 curvature_tolerant yes
fix mesh2 all mesh/surface file meshes/mesh2_periodic.stl element_exclusion_list read mesh2fixed.stl type 3 move 0.4 -3.08981 -0.003 scale 0.5 curvature_tolerant yes
fix shutter all mesh/surface file meshes/shutter_periodic.stl type 2 scale 1
fix ywalls1 all wall/gran model hertz tangential history rolling_friction epsd primitive type 2 yplane 0.83891
fix ywalls2 all wall/gran model hertz tangential history rolling_friction epsd primitive type 2 yplane 1.03338
fix granwalls all wall/gran model hertz tangential history rolling_friction epsd mesh n_meshes 3 meshes mesh1 mesh2 shutter

fix pts1 all particletemplate/multisphere 10000169 atom_type 1 density constant 2650 nspheres 4 ntry 100000 spheres 0.001291590 0.0 -0.000913292 0.002742230 -0.001291590 0.0 -0.000913292 0.002742230 0.0 0.001291590 0.000913292 0.002742230 0.0 -0.001291590 0.000913292 0.002742230 type 1

fix pts2 all particletemplate/multisphere 10000223 atom_type 1 density constant 2650 nspheres 4 ntry 100000 spheres 0.001377697 0.0 -0.000974179 0.002925046 -0.001377697 0.0 -0.000974179 0.002925046 0.0 0.001377697 0.000974179 0.002925046 0.0 -0.001377697 0.000974179 0.002925046 type 2

.... etc

fix pdd1 all particledistribution/discrete 32452843 25 pts1 0.001116 pts2 0.002277 pts3 0.004389 pts4 0.007948 pts5 0.013521 pts6 0.021606 pts7 0.032435 pts8 0.045741 pts9 0.060597 pts10 0.075414 pts11 0.088168 pts12 0.096834 pts13 0.099908 pts14 0.096834 pts15 0.088168 pts16 0.075414 pts17 0.060597 pts18 0.045741 pts19 0.032435 pts20 0.021606 pts21 0.013521 pts22 0.007948 pts23 0.004389 pts24 0.002277 pts25 0.001116

group nve_group region reg

region gen1 block 0.2 0.3 0.83891 1.03338 0.86711 1.01711 units box
region del0 block 0.2 0.3 0.83891 1.03338 0.77135 0.86711 units box
region del2 prism 0.2 0.3 0.83891 0.91926 0.77135 0.86711 0.0 0.0 -0.08035 side out
region del3 prism 0.2 0.3 0.91926 1.03338 0.77135 0.86711 0.0 0.0 0.11412 side out
region wedge intersect 3 del0 del2 del3
region wedge2 block 0.2 0.3 0.83891 1.03338 0.86711 1.10 units box

fix ins nve_group insert/pack seed 32452867 distributiontemplate pdd1 insert_every once overlapcheck yes all_in yes vel constant 0.0 0.0 -10.0 orientation random volumefraction_region 0.25 region gen1 ntry_mc 100000

fix integr all multisphere allow_group_and_set yes
group ke_group region wedge

compute k ke_group ke
variable kill equal c_k

thermo ${update}
thermo_style custom spcpu step time atoms c_k ke cpu cpuremain #performance checking
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes

dump dmp1 all mesh/gran/VTK 1000000000 rayleighmidflat_multi-data/mesh1-*.vtk id mesh1
dump dmp2 all mesh/gran/VTK 1000000000 rayleighmidflat_multi-data/mesh2-*.vtk id mesh2
dump dmp3 all custom/vtk ${update} rayleighmidflat_multi-data/particles-*.vtk id id_multisphere type x y z radius vx vy vz f_color
dump dmp4 all mesh/gran/VTK ${update} rayleighmidflat_multi-data/shutter-*.vtk id shutter

label loop
variable a loop 9001
variable b equal ke
run ${break}
variable trigger0 equal step
if "(${trigger0} > 5470)" then "jump SELF break" # && (${kill} < 0.0002)"
uncompute k
group ke_group delete
group ke_group region wedge
set region wedge property/atom color 1
set region wedge2 property/atom color 0
compute k ke_group ke
variable kill equal c_k
next a
jump SELF loop

label break
uncompute k
group ke_group delete
group ke_group region wedge
set region wedge property/atom color 1
set region wedge2 property/atom color 0
compute k ke_group ke
variable kill equal c_k
variable trigger equal ${trigger0}
print "Settlement complete, granular slide triggered."
print "Trigger Timestep = ${trigger}"

delete_atoms region wedge2 mol yes compress yes
velocity all set 0 0 0 units box
unfix ins
unfix ywalls1
unfix ywalls2

variable omegav equal 0
variable omegat equal "step - v_trigger"
variable omegav equal "(-v_omegat * 1) * v_tstep * 0.5 * 9.81 / 0.125"
fix rotate all move/mesh mesh shutter rotate/variable origin 0.25 0.74817 1.00635 axis 1 0 0 omega v_omegav

thermo_style custom spcpu step time atoms ke cpu cpuremain v_omegat v_omegav
thermo ${update}
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes

label loop2
variable c loop 9001
variable d equal ke
run ${break}
if "(${omegat} > 54700) && ($d < 0.00002)" then "jump SELF break2" #
if "${omegat} >= 218800" then "variable omegav equal 0"
next c
jump SELF loop2
label break2
print "Granular slide finished."

Any help with this would be appreciated.

mattkesseler | Tue, 02/13/2018 - 20:39

On closer inspection I've noticed that the velocity command doesn't seem to affect the multisphere particles correctly, even when I use it in a fixed boundary simulation without these issues. Is there a better way of removing all motion from a multisphere particle at a certain tine in the simulation?