Having problem during running a cohesive force model

Submitted by mmkamyabi on Tue, 11/07/2017 - 14:25

Hello everybody,
I'm trying to run an example of cohesive force model of washino/capillary/viscous. When I run this script (as below) the memory usage of the computer is increased illogically so that the run is dead. I decreased the amount of particles but it doesn't solve the problem. May you please help me to know what's the reason and what I can do to fix it?
Thank you

the input script is:

#Contact model example

atom_style granular
atom_modify map array

boundary m m m
newton off

communicate single vel yes

units si

region reg block -0.05 0.05 -0.05 0.05 0. 0.15 units box
create_box 1 reg

neighbor 0.002 bin
neigh_modify delay 0

#Material properties required for new pair styles

fix m1 all property/global youngsModulus peratomtype 5.e6
fix m2 all property/global poissonsRatio peratomtype 0.45
fix m3 all property/global coefficientRestitution peratomtypepair 1 0.9
fix m4 all property/global coefficientFriction peratomtypepair 1 0.05
fix m5 all property/global characteristicVelocity scalar 2.
fix m6 all property/global cohesionEnergyDensity peratomtypepair 1 300000

fix m7 all property/global minSeparationDistanceRatio scalar 1.01
#(value=value for the minimum separation distance, recommended as 1.01)
fix m8 all property/global maxSeparationDistanceRatio scalar 1.1
#(value=value for the maximum separation distance, recommended as 1.1)
fix m9 all property/global surfaceLiquidContentInitial scalar 0.1
#(value=value for the initial surface liquid volume in % of the solid volume)
fix m10 all property/global surfaceTension scalar 0.073
#(value=value for the surface tension of liquid (SI units N/m))
fix m11 all property/global fluidViscosity scalar 1000
#(value=value for the fluidViscosity (SI units Pas))
fix m12 all property/global contactAngle peratomtype 25
#(value_i=value for contact angle of atom type i and fluid)
fix m13 all property/global maxLiquidContent peratomtype 0.2
#(value_i=value for maximum liquid content of an atom of type i in % of particle volume)

#New pair style
pair_style gran model hertz tangential history cohesion washino/capillary/viscous #Hertzian with cohesion
pair_coeff * *

timestep 0.00001

fix gravi all gravity 9.81 vector 0.0 0.0 -1.0

fix zwalls1 all wall/gran model hertz tangential history primitive type 1 zplane 0.0
fix zwalls2 all wall/gran model hertz tangential history primitive type 1 zplane 0.15
fix cylwalls all wall/gran model hertz tangential history primitive type 1 zcylinder 0.05 0. 0.

#region for insertion
region bc cylinder z 0.01 0.01 0.025 0.05 0.0603 units box
group nve_group region reg

#distributions for insertion
fix pts1 all particletemplate/sphere 15485863 atom_type 1 density constant 2500 radius constant 0.0015
fix pdd1 all particledistribution/discrete 15485867 1 pts1 1.0

#particle insertion
fix ins nve_group insert/pack seed 32452843 distributiontemplate pdd1 &
maxattempt 100 insert_every 3000 overlapcheck yes all_in yes vel constant 0.0 0.0 -0.2 &
region bc particles_in_region 250 ntry_mc 10000

#apply nve integration to all particles that are inserted as single particles
fix integr nve_group nve/sphere

#output settings, include total thermal energy
compute rke all erotate/sphere
thermo_style custom step atoms ke c_rke vol
thermo 1000
thermo_modify lost ignore norm no

#insert the first particles so that dump is not empty
run 1
dump dmp all custom/vtk 800 post/cohesion_*.vtk id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius

#insert particles
run 10000
unfix ins

#run
run 50000 upto

aaigner's picture

aaigner | Wed, 11/08/2017 - 00:00

Hello!

Your simulation is exploding and since you are using 'boundary m m m', the simulation domain increases to contain all particles.
This is very important: Always use 'boundary f f f' if nothing else is necessary. In this case you would have lost the particles and LIGGGHTS would not stop with 'too many neighbor bins'

Sry, I have currently not the time to search for the reason of the really fast motion. Maybe it is just the timestep width.

Best wishes
Andreas