Multisphere model issues

Submitted by mtennant on Tue, 02/12/2019 - 16:38

Hi,

I am currently attempting to model multisphere particles falling onto a bed and a spreader being used to run over the particles. I am having an issue that the particles seem to disappear as they are about to fall onto the bed. I have reduced my timestep and altered the youngs modulus incase it is a force issue between the particles. I do not seem to receive any errors and I was wondering if anyone had any ideas. Thanks, the script is as follows:

log log-000.txt

#Simulation info and manipulated variable
variable run_num string 000
variable surfaceLiquidContentInitial equal 0.00

atom_style granular
atom_modify map array
boundary p f f
newton off
dimension 3

communicate single vel yes

processors * * *

units micro

# Dimensions for simulation boundary
variable regionXhi equal 500
variable regionXlo equal -500
variable regionYhi equal 32000
variable regionYlo equal -40000
variable regionZhi equal 7000
variable regionZlo equal 0

region reg block ${regionXlo} ${regionXhi} ${regionYlo} ${regionYhi} ${regionZlo} ${regionZhi} units box
create_box 1 reg

# Dimensions for particle insertion
variable regionYhi2 equal 24000
variable regionYlo2 equal 0
variable regionZhi2 equal 6000
variable regionZlo2 equal 4000

# Walls for particle insertion
variable inserthi equal ${regionYhi2}+200
variable insertlo equal ${regionYlo2}-200

# Dimensions for volume fraction region
variable regionYhi3 equal 22000
variable regionYlo3 equal -26000
variable regionZhi3 equal 400
variable regionZlo3 equal 0

neighbor 30 bin
neigh_modify every 1 delay 0

variable youngsModulus equal 5e2
variable poissonsRatio equal 0.29 # ref azom
variable coefficientRestitution equal 0.50
variable coefficientFriction equal 0.20
variable coefficientRollingFriction equal 0.005
variable minSeparationDistanceRatio equal 1.01 # Reccomended by Liggghts docs
variable maxSeparationDistanceRatio equal 1.10 # Reccomended by Liggghts docs
variable surfaceTension equal 72.9
variable fluidViscosity equal 890000
variable contactAngle equal 35.00 # Check

variable radius equal 1 # Particle radius in micrometers
variable density equal 7.87 # g/cm^3 ref azom

# Basic material properties
fix m1 all property/global youngsModulus peratomtype ${youngsModulus}
fix m2 all property/global poissonsRatio peratomtype ${poissonsRatio}
fix m3 all property/global coefficientRestitution peratomtypepair 1 ${coefficientRestitution}
fix m4 all property/global coefficientFriction peratomtypepair 1 ${coefficientFriction}
fix m5 all property/global coefficientRollingFriction peratomtypepair 1 ${coefficientRollingFriction}

# Material properties for washino/capillary/viscous cohesion model
fix m6 all property/global minSeparationDistanceRatio scalar ${minSeparationDistanceRatio}
fix m7 all property/global maxSeparationDistanceRatio scalar ${maxSeparationDistanceRatio}
fix m8 all property/global surfaceLiquidContentInitial scalar ${surfaceLiquidContentInitial}
fix m9 all property/global surfaceTension scalar ${surfaceTension}
fix m10 all property/global fluidViscosity scalar ${fluidViscosity}
fix m11 all property/global contactAngle peratomtype ${contactAngle}

# Set initial pair style for insertion (Hertzian without cohesion)
pair_style gran model hertz tangential history rolling_friction cdt
pair_coeff * *

# Include gravity in the z-direction
fix gravi all gravity 9.81e-6 vector 0.0 0.0 -1.0

variable dt equal 0.001

# Set the simulation timestep
timestep ${dt} #Microseconds

# Number of steps for each phase of the simulation
variable time_scale equal 1/${dt}

variable insert_steps equal 5000*${time_scale}
variable check_timestep equal 5000*${time_scale}
variable thermo_dump equal 5000*${time_scale}
variable data_dump equal 5000*${time_scale}
variable insert_p equal 50000*${time_scale}
variable settle_p equal 150000*${time_scale}
variable spread_p equal 1100000*${time_scale}

# Import blade geometry file
fix blade all mesh/surface file blade-narrow.stl type 1 scale 1000 move 0 25500 500
fix bladewall all wall/gran model hertz tangential history rolling_friction cdt mesh n_meshes 1 meshes blade

# Set Wall locations and collision models for boundaries
fix zwalls2 all wall/gran model hertz tangential history rolling_friction cdt primitive type 1 zplane ${regionZlo}
fix ywalls1 all wall/gran model hertz tangential history primitive type 1 yplane ${regionYhi}
fix ywalls2 all wall/gran model hertz tangential history primitive type 1 yplane ${regionYlo}

# Introduce walls to insertion region to contain particles
fix ywalls3 all wall/gran model hertz tangential history primitive type 1 yplane ${inserthi}
fix ywalls4 all wall/gran model hertz tangential history primitive type 1 yplane ${insertlo}

# Region for insertion
region regi block EDGE EDGE ${regionYlo2} ${regionYhi2} ${regionZlo2} ${regionZhi2} units box

group multispheres initialize
fix intMulti multispheres multisphere

# Define a multisphere particle template for each non-spherical particle shape
# that will be used as an input for the fix_particledistribution_discrete command

fix pts1 multispheres particletemplate/multisphere 15485867 atom_type 1 &
density constant ${density} nspheres 12 ntry 1000000 &
spheres file Part.2.lammps scale 1e6 type 1 mass 12.709713683 inertia_tensor 3.5423379326e2 3.7130271494e2 3.9388290335e2 1.9718882783e1 2.6516484147e1 1.2830073431e2 use_density

fix pts2 multispheres particletemplate/multisphere 32452843 atom_type 1 &
density constant ${density} nspheres 12 ntry 1000000 &
spheres file Part.3.lammps scale 1e6 type 2 mass 19.714351546 inertia_tensor 8.8052395874e2 4.6698450869e2 8.5726867291e2 -1.5962785071e2 1.0440694358e1 -1.0678891619e2 use_density

fix pts3 multispheres particletemplate/multisphere 49979687 atom_type 1 &
density constant ${density} nspheres 12 ntry 1000000 &
spheres file Part.4.lammps scale 1e6 type 3 mass 26.878677059 inertia_tensor 1.1805295600e3 9.0126152883e2 1.4545507014e3 1.0582391271e2 3.1134169036e1 1.1961458017e2 use_density

fix pts4 multispheres particletemplate/multisphere 67867967 atom_type 1 &
density constant ${density} nspheres 12 ntry 1000000 &
spheres file Part.6.lammps scale 1e6 type 4 mass 10.273518668 inertia_tensor 2.6522781863e2 2.4705954152e2 2.3031343468e2 5.6700500788e1 -9.4456038033 4.3815860004e1 use_density

fix pts5 multispheres particletemplate/multisphere 86028157 atom_type 1 &
density constant ${density} nspheres 12 ntry 1000000 &
spheres file Part.9.lammps scale 1e6 type 5 mass 13.244806030 inertia_tensor 4.1197447136e2 2.7144104236e2 4.2880065557e2 -9.8056722899e1 -2.5467845640e1 -7.6527593738e1 use_density

# Define a discrete particle distribution that will be used by the insert command
# All templates are added with a weighting factor

fix pdd1 all particledistribution/discrete 37189 5 pts1 0.2 pts2 0.2 pts3 0.2 pts4 0.2 pts5 0.2

fix ins all insert/rate/region seed 57097 distributiontemplate pdd1 &
nparticles 3000 particlerate 0.15 &
insert_every ${insert_steps} overlapcheck yes &
region regi

# ======================================
# Output Settings
# ======================================

# Integrator for multisphere rigid bodies
#fix integr all multisphere

# Check reliability of timestep
fix ts all check/timestep/gran ${check_timestep} 0.25 0.25

# Output settings, include total thermal energy
compute rke all erotate/sphere

thermo_style custom step atoms c_rke vol cpu
thermo ${thermo_dump}
thermo_modify lost ignore norm no

# ======================================
# Particle Insertion
# ======================================

# Insert the first particles so that dump is not empty
run 1

#### Which group should be dumped?
dump dmp1 all custom/vtk ${data_dump} post/moist${run_num}/particlecoords${run_num}_*.vtk id type &
x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
dump dumpstl all mesh/stl ${data_dump} post/moist${run_num}/spreader${run_num}_*.stl

# Insert particles and allow to settle (+- 68k to insert all particles)
# (+- 32k for last pastircles to drop)
run ${insert_p} upto

# Drag used as energy minimisation technique to improve initial pile
fix damping1 all viscous 4000

# Remove insertion walls and allow particles to settle
unfix ywalls3
unfix ywalls4
run ${settle_p} upto # +- 80k for particles to settle
unfix damping1

# Include Drag effects
fix damping2 all viscous 1677.61 # Stoke's Drag on particles

# ======================================
# Blade Motion
# ======================================

# Change pair style to include cohesion (Hertzian with cohesion)
pair_style gran model hertz tangential history cohesion washino/capillary/viscous rolling_friction cdt
pair_coeff * *

# Begin spreader motion
fix movecad1 all move/mesh mesh blade linear 0 -0.06 0.
run ${spread_p} upto

# Computation of final volume fraction
region regvc block EDGE EDGE ${regionYlo3} ${regionYhi3} ${regionZlo3} ${regionZhi3} units box
group groupvc region regvc

dump dmp3 groupvc custom/vtk 1 volume_fraction/moist${run_num}/moist${run_num}.vtk id type &
x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
dump dmp4 groupvc custom 1 volume_fraction/moist${run_num}/moist${run_num}.txt id type &
x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
run 0