Help jamming a granular specimen and equilibrating it so that its kinetic energy is negligible in 3D with 2 periodic boundary conditions

Hi all,

I have a 3-dimensional granular packing that is enclosed by a rectangular cuboid (l x w x w). 2 of the directions of the cuboid correspond to the width dimension, w, which is periodic and there is a reflective wall bounding the simulation in the l dimension

I’ve added a picture to get a better understanding…

My apologies for it being rudimentary but P denotes periodic in that direction and R denotes the reflective wall on the top and bottom.

I initialize the packing without gravity such that there is a packing fraction close to the critical RCP volume fraction of 0.64 but not exceeding it…i.e pressure is 0.

My goal is to confine this initial state such that there is an initial pressure and the granules are equilibrated such that there is negligible kinetic energy (rotational+translational) when viscosity is removed. I am using pair/gran/hertz/history with tangential effects and use NVE for time integration.

The main approach I used was using fix deform on the periodic dimensions and compressing at a slow strain rate. Then unfix this command and let the granules equilibrate with fix viscous. What I found is that fix viscous does not dissipate translational or rotational kinetic energy even when acting for large times and these energy values remain almost constant although there is a finite pressure.

Pertinent code is as follows

units si

timestep 1e-6

fix 1111 all deform 1 x vel -0.05 y vel -0.05

run 		14000

unfix 1111

fix  12 all viscous 10

run 1000000

Is there another way I can achieve what I want? Any help is much appreciated.

Is it also possible to use fix controller to move the position of the periodic box so that the pressure can be controlled? I tried doing this but kept getting a floating point error.

Thank you!

There is no time integration here. Thus particles will not relax or dissipate kinetic energy.

Please also note that all that fix viscous does is to compute a velocity dependent translational force.

This is by far too vague a description to give meaningful advice.

Sorry Axel, I forgot to include:

fix  		1 all nve/sphere

Ahh, maybe fix viscous/sphere might be useful for rotational damping then? But I don’t understand why the translational kinetic energy won’t dampen.

I tried running a frictionless simulation and here rotational kinetic energy is 0 as expected but here too fix viscous doesn’t dampen kinetic energy to near 0. It should be for this case, right?

About the controller, I was wanting to know if it is possible for me to change the location of the box in the periodic dimensions using a controller while keeping the reflective wall on the top and bottom fixed so that the pressure can be controlled. This is the code I used:

variable       foo internal -0.05 
fix 1111 all deform 1 x vel v_foo y vel v_foo

fix 1111 rw controller 1 1 1 1 1 v_press 100000 foo

I don’t know if this is possible but I thought it might be worth a shot.

Are there any other strategies you might know of to jam a granular system described by my domain?

Unless you provide small and complete test examples that can be easily and quickly run and clearly demonstrate the respective issues at hand, I have no further comments on your inputs or problems. What I do have to say is the following.

The practice of providing only piece-wise inputs, yet expecting detailed and competent discussion of complex physics is something I very much dislike. You need to keep in mind that you are not talking to your adviser or a collaborator that is very familiar with your research domain and - possibly - your specific research project, but with people that have a general interest in MD, use LAMMPS and are willing to help others in order to share their expertise and improve their understanding of LAMMPS and hone their skills at using it.
Thus, in order to give competent answer, people like me need to construct (complete!) input decks, experiment with them, teach ourselves what the issues are and check out what possible solutions can be. That would be - in case of examples like yours - a lot of effort. But with things like “oh, I forgot to add”, you are creating a lot of frustration. How can anybody now trust you that this is the only relevant thing you have left out? Would you spend a few hours of your free(!) time to build some test input in order to be able to discuss some question if it becomes likely that whatever you find out will be invalidated because you started from incomplete assumptions?

In short, you have to make a better effort in helping us to be able to help you.

I am so sorry Dr Axel, I totally understand. My apologies, this has been bugging me and I haven’t been thinking with clarity.

I hope you don’t mind me sharing my input deck where I randomly initialize granules, equilibrate, then deform box, and perform cycles of viscous damping to relax the system.

I’ve attached the read_data data script where i have around 11000 granules.

units si
atom_style  	sphere
boundary 	pp pp fs
newton  	off
comm_modify 	vel yes
read_data	small_3d_63p5.spheres

log 			n_0p63p5_2.lammps

neighbor 	0.02 bin
neigh_modify 	delay 0

#pair_style      gran/hertz/history 8.424908425E+10 NULL 10 10 0.05 1
pair_style      gran/hertz/history 8.424908425E+10 0 0 0 0 0

pair_coeff	* *

timestep 	1e-6
fix  		1 all nve/sphere
#fix  		2 all gravity 9.0 spherical 0.0 -180.0
fix 11 all viscous 5

region 		atom_region block 0 0.1 0 0.1 0 0.8
region 		atom_region2 block 0 0.08 0 0.08 0 0.75
group 		atoms region atom_region
group		atoms2 region atom_region2

compute        peratom atoms stress/atom NULL pair
compute        p atoms reduce sum c_peratom[1] c_peratom[2] c_peratom[3]
compute        p1 atoms reduce sum c_peratom[4] c_peratom[5] c_peratom[6]
variable       press equal -(c_p[1]+c_p[2]+c_p[3])/(3*0.1*0.1*0.8)#((y[896]-1e+5)-(y[895]+1e+5))*((z[893]-1e+5)-(z[894]+1e+5)))
variable       stress_ratio equal c_p[1]/c_p[2]

compute        peratom2 atoms2 stress/atom NULL pair
compute        p2 atoms2 reduce sum c_peratom2[1] c_peratom2[2] c_peratom2[3]
compute        p12 atoms2 reduce sum c_peratom2[4] c_peratom2[5] c_peratom2[6]
variable       press2 equal -(c_p2[1]+c_p2[2]+c_p2[3])/(3*0.08*0.08*0.75)#((y[896]-1e+5)-(y[895]+1e+5))*((z[893]-1e+5)-(z[894]+1e+5)))
variable       stress_ratio2 equal c_p2[1]/c_p2[2]

fix 		mywall1 all wall/reflect zlo EDGE 
fix		mywall2 all wall/reflect zhi EDGE

compute  	rot_kin all erotate/sphere # rotational kinetic energy of the spherical particles
thermo_style 	custom step atoms ke c_rot_kin v_press v_press2
thermo  	1000
thermo_modify 	lost ignore norm no # atom number non-constant
dump myDump all atom 1000 n_0p63p5_2.atom

run  		20000 # final longer run

unfix 11

fix 111 all viscous 0

run 	        4000

fix 1111 all deform 1 x vel -0.05 y vel -0.05

run 		14000

unfix 1111

run 		20000

fix 1112 all viscous 10

run 	400000

unfix	1112

fix 1113 all viscous 0

run	20000

unfix	1113

fix 1114 all viscous 1

run 400000

unfix 1114

fix 1115 all viscous 0

run 100000

write_data      n_0p63p5_2.sim

Apologies once again, and I thank what you have done for the community :slight_smile:

small_3d_63p5.spheres (1.5 MB)

This is the second most irritating thing that people do after what I complained about in my previous response. It is particularly irritating because it is so predictable.

Instead of constructing a specific test case following the guidelines that I just outlined, you just dump your unmodified input on people. That is almost as bad. This is larger than it needs to be, this has parts in it that are not relevant to the problem and you still failed to explain clearly what you are getting and what you expect to be getting instead. So basically, you are still not doing much to help getting an answer.