Unexpected bouncing in closed packing

Submitted by rahulsoni on Fri, 10/25/2013 - 15:12

I am attempting to run drum mixing simulation of packed granular bed. The trial outputs in two cases of simple cubic and HCP packing are attached with the link below:

Even after trying several combination of properties I am getting similar undesired bouncing of particles. At the prima facie it appears to be caused by the force between the drum wall and particles nearby it.

Anyone please help how can I get rid of that.

rahulsoni | Fri, 10/25/2013 - 16:19

# Moving mesh example

atom_style granular
boundary f f f
newton off

communicate single vel yes
units si

#region reg block -.1 0.6 -.1 0.2 -0.3 0.25 units box
#region reg block -1 1 -.5 .5 -.5 .5 units box
region reg block -0.08 0.08 -0.15 0.15 -0.15 0.15 units box
create_box 2 reg

neighbor 0.005 bin
neigh_modify delay 0

#Material properties
fix m1 all property/global youngsModulus peratomtype 3.e8 3.e8
fix m2 all property/global poissonsRatio peratomtype 0.5 0.5
fix m3 all property/global coefficientRestitution peratomtypepair 2 &
0.7 0.7 &
0.7 0.7
fix m4 all property/global coefficientFriction peratomtypepair 2 &
0.5 0.5 &
0.5 0.5
fix m5 all property/global characteristicVelocity scalar 0.5.

#pair style
pair_style gran/hooke/history #Hooke without cohesion
pair_coeff * *

timestep 0.00005

fix 1 all nve/sphere

#box walls
fix boxwalls_x1 all wall/gran/hooke/history primitive type 1 xplane -0.075
fix boxwalls_x2 all wall/gran/hooke/history primitive type 1 xplane +0.075
fix boxwalls_y1 all wall/gran/hooke/history primitive type 1 yplane -0.145
fix boxwalls_y2 all wall/gran/hooke/history primitive type 1 yplane +0.145
fix boxwalls_z1 all wall/gran/hooke/history primitive type 1 zplane -0.145
fix boxwalls_z2 all wall/gran/hooke/history primitive type 1 zplane +0.145

#import mesh from cad:
group cad_group_1 region reg
fix cad1 cad_group_1 mesh/surface file Drum_No_lifter.stl type 1 move 0. 0. 0. &
scale 0.001 rotate axis 0. 1. 0. angle 0. rotate axis 0. 0. 1. angle 0.

#import splitter cad:
#group cad_group_2 region reg
#fix cad2 cad_group_2 mesh/surface file Drum_Splitter.stl type 1 move 0. 0. 0. &
# scale 0.001 rotate axis 0. 1. 0. angle 0. rotate axis 0. 0. 1. angle 0.

#region and insertion (taken from mesh_tet code)

#region mesh_1 mesh/tet file test_1.vtk scale 1. move 0. 0. 0. rotate 0. 0. 0. units box
#group nve_group_1 region mesh_1
#create_atoms 1 region mesh_1
#set group nve_group_1 density 500 diameter $d
#region mesh_2 mesh/tet file test_2.vtk scale 1. move 0. 0. 1.2. rotate 0. 0. 0. units box
#group nve_group_2 region mesh_2
#create_atoms 2 region mesh_2
#set group nve_group_2 density 500 diameter $d

#particle insertion (old code)
#region bc_1 block -0.05 0.05 -0.10 -0.005 0.075 0.095 units box
#group nve_group_1 region reg
#fix ins_1 nve_group_1 pour/legacy 32000 1 1 vol 0.5 5000 diam 0.005 0.005 dens 2700 2700 vel 0. 0. 0. 0. -0.0. region bc_1

#particle insertion
#region bc_2 block -0.05 0.05 0.005 0.10 0.075 0.095 units box
#group nve_group_2 region reg
#fix ins_2 nve_group_2 pour/legacy 32000 2 1 vol 0.5 5000 diam 0.005 0.005 dens 2700 2700 vel 0. 0. 0. 0. -0.0. region bc_2

variable d equal 0.005
lattice sc $d
#reg_red mesh/tet file Packing_Drum_Red.vtk scale 0.001. move 0. 0. 0. rotate 0. 0. 0. units box
#particle insertion_red
region reg_red mesh/tet file Packing_Drum_Red_Volume.vtk scale 0.001. move 0. 0. 0. rotate 0. 0. 0. units box
#group nve_group_1 region reg_red
create_atoms 1 region reg_red
group nve_group_1 region reg_red
#set group nve_group_1 density 2700 diameter $d
set region reg_red density 2700 diameter $d

#particle insertion_gree
region reg_green mesh/tet file Packing_Drum_Green_Volume.vtk scale 0.001. move 0. 0. 0. rotate 0. 0. 0. units box
#group nve_group_2 region reg_green
create_atoms 2 region reg_green
group nve_group_2 region reg_green
#set group nve_group_2 density 2700 diameter $d
set region reg_green density 2700 diameter $d

#thermo settings
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

#make a dump of particles and the stl file
dump dmp_1 nve_group_1 custom 200 post/dump*.red id type type x y z ix iy iz vx vy vz fx fy fz omegay omegaz radius
dump dmp_2 nve_group_2 custom 200 post/dump*.green id type type x y z ix iy iz vx vy vz fx fy fz omegay omegaz radius
dump dumpstl cad_group_1 mesh/stl 200 post/dump*.stl cad1

#restart command

restart 500 tmp.restart

#use the imported mesh (drum) as granular wall
fix drum_mixer all wall/gran/hooke/history mesh n_meshes 1 meshes cad1
run 1 upto

#use the imported mesh (drum_splitter) as granular wall
#fix splitter all wall/gran/hooke/history mesh n_meshes 1 meshes cad2
#run 2 upto

#run for particle insertion
#run 100000 upto
run 5000 upto
#stopping particle insertion
#unfix ins_1
#run 1
#unfix ins_2
#run 1

fix 2 all gravity 9.81 vector 0.0 0.0 -1

#run to let particle settle
run 2000

#removing drum_splitter
#unfix splitter
run 1
#unfix cad2
#run 15002 upto

#now lift the particles up
fix rotate_drum all move/mesh mesh cad1 rotate origin 0. 0. 0. axis 1. 0. 0. period 14.0.
#run 2800000
run 20000

--
Regards

Rahul Kumar Soni
Scientist, CSIR - IMMT, India

richti83's picture

richti83 | Fri, 10/25/2013 - 17:59


fix m1 all property/global youngsModulus peratomtype 3.e8 3.e8

and

timestep 0.00005

seems not to be a reasonable combination for 5mm particles ...
pls read some literature about critical timestep (Rayleigh time) and try fix check/timestep/gran

I'm not an associate of DCS GmbH and not a core developer of LIGGGHTS®
but I have 7 years user experience using LIGGGHTS®, ParaView and LINUX

rahulsoni | Fri, 10/25/2013 - 18:55

Okay Richter, a nice piece of information. Should 5e-07 would be a reasonable value for the time-step.

One more thing, as you can see that I have imported a mesh file representing the drum. I have assigned the type 1 in fix command to import stl. Therefore, youngmodulus, coeff of friction, coeff of restitution and poission's ratio are assigned to stl file. But I also want to assign density for this, how can I do this.

Also, please see the mixing for randomly distributed particles in the drum in the link below:

Really unexpected results are seen on the later half of the simulation. Please tell if timestep is the reason for same or anything else. (Material properties are as above)

--
Regards

Rahul Kumar Soni
Scientist, CSIR - IMMT, India

richti83's picture

richti83 | Fri, 10/25/2013 - 19:27

As far as I know the "density" for walls is infinitely and the contact force depends only from particle stiffness -maybe christoph can correct me when I'm wrong.
I would suggest 1e-6 as timestep or lower for Y>1e6, but its only a feeling. We use a commercial DEM software to determine a good timestep for a given particle/Y/density combination. Use fix check/timestep to get a feeling (you will get a lot of warnings if timestep to big, than make it smaller).
About your video I've never seen such a behaviour it can be timestep related problem or the fact that you make two dumpfiles (maybe they arent synchron) try dump id all (all particles are in group all by default !)
best,
Christian.

I'm not an associate of DCS GmbH and not a core developer of LIGGGHTS®
but I have 7 years user experience using LIGGGHTS®, ParaView and LINUX

cstoltz | Fri, 10/25/2013 - 21:12

To echo Christian's comments, timestep is the big problem with the way the simulation is running. A decent rule of thumb for stability is to set the timestep equal to ~20% of the Rayleigh timestep, which for your simulation, would give dt ~ 1.7e-5 sec. You can sometimes get away with larger values depending on the simulation, and sometimes need smaller; just depends on the simulation.

I don't think that the timestep could be responsible for the odd video. That looks like two image sequences that aren't matched up for some reason. Any reason to create two separate particle dumps, rather than simply dump all to one file and color by type?

Regards,
Chris

rahulsoni | Sat, 10/26/2013 - 06:07

Actually, I wanted the two side particles to be colored differently and therefore, I assigned nve_group_1 and nve_group_2 and created two set of dump files. I am really unaware how to color the two type particles differently in paraview and therefore I did it. When I import set dump files in paraview then it shows as collection. Therefore, if two sets results in two collection and I can color them differently. But it would be great to know if this can be done with type ID itself.

--
Regards

Rahul Kumar Soni
Scientist, CSIR - IMMT, India

cstoltz | Sat, 10/26/2013 - 14:31

In Paraview, go to the Display tab in the Object Inspector, and in the section on Color, select Color by = Type.

You may also want to consider dumping the particles straight to .vtu format (dump atom/vtk) as these can be directly loaded into Paraview and includes type, id, radius, mass, force, and velocity as attributes.

Regards,
Chris

rahulsoni | Sun, 10/27/2013 - 06:53

Thank you and sorry Chris. Meanwhile I got another way to set color in paraview. Glyph>> Set type in scalars and coloring. Anyway the other way told by you is useful.

Thank you Richter for always helping in to troubleshoot the issue.

--
Regards

Rahul Kumar Soni
Scientist, CSIR - IMMT, India

cstoltz | Sun, 10/27/2013 - 12:36

One other thing you might want to do is to use the Point Sprite plugin instead of using the Glyphs. It allows for much faster rendering when you have a large number of particles.

Regards,
Chris

rahulsoni | Sun, 10/27/2013 - 14:19

Chris and Richter, I have reduced the timestep but still the simulation is facing similar bouncing. What I am suspecting is that it is happening because when the particles are created in packing and simulation starts then since they were already near to each other the forces get calculated and implemented. In the dump file I can see the non-zero values for forces, velocity and angular velocity even before the gravity was switched on. And it because this bouncing happens.

I think I can turn all above attribute values to zero and then switch gravity on followed by starting the simulation. I can set force values to zero by either
fix no_force all setforce 0.0 0.0 0.0
or
fix no all freeze

But I am not getting any option to set zero values for linear and angular velocities.

--
Regards

Rahul Kumar Soni
Scientist, CSIR - IMMT, India

rahulsoni | Sun, 10/27/2013 - 15:40

Okay, I got how to do set velocity and forces to zero. See the below part and correct me if I am wrong:

set group all vx 0.0
set group all vy 0.0
set group all vz 0.0
set group all omegax 0.0
set group all omegay 0.0
set group all omegaz 0.0
fix no all freeze

--
Regards

Rahul Kumar Soni
Scientist, CSIR - IMMT, India

rahulsoni | Mon, 10/28/2013 - 17:30

Chris and Richter, I have tried several combinations and options but could not get what I want. See below the intial pahse of animation which clearly shows the upward bouncing:

Dump files clearly show that forces are zero after giving fix to set all forces zero. However, right after unfixing the same dump files witness forces of huge magnitude. And I think it is because forces are calculated for neighbor list. I just want a low energy system; even if settling happens particles should move downwards.

PS: I gone through the dump file of one other simulation which shows that forces, velocity and angular velocity values are zero when the particles are settled and are in no motion. Inspired by this I turn these attribute values to zero. Unfortunately, it could not work.

Desperately need your help.

--
Regards

Rahul Kumar Soni
Scientist, CSIR - IMMT, India

richti83's picture

richti83 | Mon, 10/28/2013 - 19:38

I dont know what you are trying and what is the final goal.

BUT I think you can not generate particles with large initial overlap !
I would generate partilcle-packaging as you did, but use only half the diamter so that there is no overlap at step 0.
Do not enable gravity BUT enable integration (nve/sphere)
Than grow the particles (fix adapt) in 100 steps with 1000 runsteps between every growstep (see example packing for details)
than enable gravity.

fix freeze will only turn off particle-force interaction as long as it is fixed, thats why your bed explodes when unfixing it.

btw: my mill looks like this at the moment:
http://www.youtube.com/watch?v=A37B_5L-5Ao
Why do you use only one particle fraction ?

best,
Christian

I'm not an associate of DCS GmbH and not a core developer of LIGGGHTS®
but I have 7 years user experience using LIGGGHTS®, ParaView and LINUX

cstoltz | Mon, 10/28/2013 - 19:53

Can't say for sure as I don't have your CAD geometry, but from looking at your video, it appears to me that at the outset of the simulation, there is a rapid radial contraction of the particle bed. I'm guessing that this is due to your insertion region placing a particle too close to the wall of the geometry. Even if you've already zeroed all forces (which shouldn't be required anyway), as soon as you turn everything on, the high degree of overlap between some particles and the wall causes an enormous force to drive those particles radially inward and then transmit that energy through the bed.

Suggest ensuring the insertion domain is more than a particle radius smaller in all dimensions than the actual geometry.

Also, be aware that placing particles in such a structured lattice can lead to some unintended consequences. It may take quite a while to 'melt' the lattice structure during the simulation.

Regards,
Chris

rahulsoni | Wed, 10/30/2013 - 08:02

Thanks Chris and Richter, your suggestions have worked. Especially of using the fix adapt command. Interestingly, one part of the problem was related to particles being too close to wall, as apprehended by Chris.

Thanks again to both of you.

--
Regards

Rahul Kumar Soni
Scientist, CSIR - IMMT, India