Particles not interacting with .stl Walls

Submitted by willroc7 on Wed, 08/03/2011 - 19:00

Greetings. I am attempting to run a simple compaction program of particles in a cylinder. I seem to get a successful packing of particles but they do not interact with the topwall as it moves down in the z direction. Using the cohesion and compression examples as templates, I have the following code:

#Contact model example

dimension 3
atom_style granular #diameter, density, angular velocity
atom_modify map array #processor stores lookup table of length N, N atoms. fast but limited by memory
boundary f f f #fixed boundaries in 3 dimensions
newton off #turns Newton's 3rd law off for pairwise/bonded interactions
#should get same answer for on/off
#on yields modest savings in computation at cost of twice communication

communicate single vel yes #communicates velocity information b/w processors

units si

region reg cylinder z 0 0 .005 0 .05 units box #cylindrical simulation space
#r=0.005, h=0.05 (meters)
create_box 1 reg

neighbor 0.005 bin #search for neighbors within 0.005m of particle
neigh_modify delay 0 #build neighbor list ever step

#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.5
fix m4 all property/global coefficientFriction peratomtypepair 1 0.1
fix m5 all property/global characteristicVelocity scalar 2.
fix m6 all property/global cohesionEnergyDensity peratomtypepair 1 500

#New pair style
pair_style gran/hertz/history 1 0 #Hertzian without cohesion
pair_coeff * *

timestep 0.001

fix gravi all gravity 9.81 vector 0.0 0.0 -1.0

#fix zwall all wall/gran/hertz/history 1 0 zplane 0.0 .051 1
fix cylwalls all wall/gran/hertz/history 1 0 zcylinder 0.005 1

#granular walls
fix topwall all mesh/gran/stressanalysis topwall.stl 1 1. 0. 0. 0. 0. 0. 0.
fix bottomwall all mesh/gran/stressanalysis bottomwall.stl 1 1. 0. 0. 0. 0. 0. 0.
fix wall all wall/gran/hertz/history 1 0 mesh/gran 2 topwall bottomwall

fix F_movetop all move/mesh/gran linear 0. 0. -0.025 units box topwall 1
#distributions for insertion
fix pts1 all particletemplate/sphere 1 atom_type 1 density constant 2500 radius constant 0.0005
fix pdd1 all particledistribution/discrete 1. 1 pts1 1.0

#region and insertion
group nve_group region reg
fix ins all pour/dev/packing 1 distributiontemplate pdd1 vol 0.5 1000 region reg

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

dump D_stl all stl 1000 post/Wall-*.stl
dump I_stl all mesh/gran/VTK 1000 post/Stress_file-*.vtk stress

#output settings, include total thermal energy
compute 1 all erotate/sphere
thermo_style custom step atoms ke c_1 vol
thermo 1000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes

#insert the first particles so that dump is not empty
run 1
dump dmp all custom 1000 post/dump.hertztest5 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

The .stl files are:
# ======
#topwall.stl:
solid Mesh
facet normal 0.000000 0.000000 -1.000000
outer loop
vertex -0.00500000 -0.00500000 .051
vertex -0.00500000 .00500000 .051
vertex .00500000 .00500000 .051
endloop
endfacet
facet normal 0.000000 0.000000 -1.000000
outer loop
vertex .00500000 -0.00500000 .051
vertex -0.00500000 -0.00500000 .051
vertex .00500000 .00500000 .051
endloop
endfacet
endsolid Mesh
#=====

and the same for bottom wall except with z=0.0.

I can view the topwall moving in ParaView after converting with pizza.py. Any help would be much appreciated.

Thanks,
Will

ckloss's picture

ckloss | Thu, 08/04/2011 - 16:14

Hi Wil,

that is due to a faulty input script. 2 comments: (a) you are not integrating your particles, (b) the time-step is way too large. If in doubt, use fix chek/timestep gran to check it

Cheers,
Christoph

willroc7 | Thu, 08/04/2011 - 19:35

Thank you very much for the help Christoph. It is now running! I added the following line:

fix integr2 all nve/sphere

Is there a reason the fix nve/sphere command I had in place was not running the integration? I thought "group nve_group region reg" would put all of the particles in the "reg" region into the nve_group which would be integrated by "fix integr nve_group nve/sphere".

I also decreased the time step by two orders of magnitude.

Thanks again,
Will

ckloss's picture

ckloss | Thu, 08/04/2011 - 19:39

>>Is there a reason the fix nve/sphere command I had in place was not running the integration?
The group command is a one-time assignment, so the particles that you inserted were not in the group you were integrating

Cheers,
Christoph