DEM exploding with certain input parameters

Hi all. I am currently running simulations of a dry chute avalanche using LIGGGHTS-DEM at various different physical sizes. My test case was using particles with a Gaussian distribution from 0.75-1.95 mm with a timestep of 3.54E-7 seconds, material stiffnesses of 60-200 GPa, and an initial input velocity of 10 m/s during a settlement phase. This settlement phase is necessary for ensuring the right packing of particles in a timely manner. I have run a simulation that is 10 times larger, with particles of distribution 7.5-19.5 mm, a timestep of 1.12E-6 seconds, the same material stiffness, and input velocity of 31.62 m/s. These parameters are scaled such that the slide undergoes the same motion relative to the ramp geometry in the same number of simulation timesteps.

The simulations using both of these sizes have worked fine without any problems for 6 million iterations each. I have tried running a small-scale counterpart with particles 0.075-0.195mm, a timestep of 1.12E-7 seconds, the same material stiffness, and an input velocity of 3.16 m/s. This simulation is now exploding after approximately 300 iterations, with the point of origin being a particle colliding with a channel wall. I have tried using an input velocity of 0.5 m/s, and this resulted in the simulation exploding at 44000 iterations.

I suspect this is probably something to do with the timestep being larger relative to the particle diameter, but since the timestep is scaled with the input velocity (and the particle velocities throughout the actual avalanche) I didn't think this would be a factor. Let me know what your thoughts are. Attached is the relevant part of the input script I've used for LIGGGHTS:

log log.newersupersmallflat

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

communicate single vel yes

units si

region reg block -0.001 0.051 -0.16 0.12 -0.001 0.11 units box
create_box 3 reg

neighbor 0.00027 bin
neigh_modify delay 0

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 #Hertzian without cohesion
pair_coeff * *

variable tstep equal 1.118033989e-7
timestep ${tstep}

fix gravi all gravity 9.81 vector 0.0 0.0 -1.0

fix mesh1 all mesh/surface file meshes/mesh1.stl element_exclusion_list read mesh1fixed.stl type 2 move 0.0 -3.08981 -0.003 scale 0.05 curvature_tolerant yes
fix mesh2 all mesh/surface file meshes/mesh2.stl element_exclusion_list read mesh2fixed.stl type 3 move 0.0 -3.08981 -0.003 scale 0.05 curvature_tolerant yes
fix shutter all mesh/surface file meshes/shuttermidlabtest40.stl type 2 scale 0.1
fix xwalls1 all wall/gran model hertz tangential history rolling_friction epsd primitive type 2 xplane 0.0
fix xwalls2 all wall/gran model hertz tangential history rolling_friction epsd primitive type 2 xplane 0.05
fix ywalls1 all wall/gran model hertz tangential history rolling_friction epsd primitive type 2 yplane 0.083891
fix ywalls2 all wall/gran model hertz tangential history rolling_friction epsd primitive type 2 yplane 0.103338
fix granwalls all wall/gran model hertz tangential history rolling_friction epsd mesh n_meshes 3 meshes mesh1 mesh2 shutter

fix pts1 all particletemplate/sphere 15485863 atom_type 1 density constant 2650 radius gaussian number 0.0000675 0.00001
fix pdd1 all particledistribution/discrete 32452843 1 pts1 1.0 #pts2 0.4 pts3 0.3

group nve_group region reg

region gen1 block 0.0 0.05 0.083891 0.103338 0.086711 0.101711 units box
region del0 block 0.0 0.05 0.083891 0.103338 0.077135 0.086711 units box
#region del1 prism 0.0 0.05 0.055322 0.083891 0.086711 0.094365 0.0 0.0 0.028569 side out
region del2 prism 0.0 0.05 0.083891 0.091926 0.077135 0.086711 0.0 0.0 -0.008035 side out
region del3 prism 0.0 0.05 0.091926 0.103338 0.077135 0.086711 0.0 0.0 0.011412 side out
region wedge intersect 3 del0 del2 del3
region Del1 block 0.0 0.05 0.083891 0.103338 0.086711 0.11 units box
region Del2 prism 0.0 0.05 0.055322 0.083891 0.086711 0.094365 0.0 0.0 0.028569 units box
region wedge2 union 2 Del1 Del2

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

fix integr nve_group nve/sphere
group ke_group region wedge

variable update equal 2000
variable break equal 20000
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

restart 100000 newersupersmallflat.one.restart newersupersmallflat.two.restart
dump dmp1 all mesh/gran/VTK 1000000000 newersupersmallflat-data/mesh1-*.vtk id mesh1
dump dmp2 all mesh/gran/VTK 1000000000 newersupersmallflat-data/mesh2-*.vtk id mesh2
dump dmp3 all custom ${update} newersupersmallflat-data/particles-*.gz id type x y z radius vx vy vz
dump dmp4 all mesh/gran/VTK ${update} newersupersmallflat-data/shutter-*.vtk id shutter