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.

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

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.