[lammps-users] non-zero shear stress (z-components) in a 2D (x-y) problem

Hi everyone,

I am performing a biaxial (x-y directions) compression of frictional granules in a rectangular box using rigid walls in LAMMPs. I have used fix enforce2d to enforce 2-dimensionality and my simulation procedure is as follows:

  1. compress granules from a loose state by moving the top and right walls inward till it forms a dense sample (beyond jamming transition)
  2. make walls stationary and let the sample equilibrate.

During step 1 when I compute stress, I have shear stress only in xy but I see shear stress in xz and yz after the onset of step 2. Why does this happen? I’ve attached a code snippet of my input script which emphasizes step 1 and 2.

###Step 1#####
variable zero equal 0
fix 1 granules nve/sphere
fix 2 granules viscous 1
fix 3 stationary move linear 0 0 0

fix 4 granules langevin 2 0 0.001 100000 omega yes

fix 5 topwall move linear 0 -0.055 0
fix 6 rightwall move linear -0.03056 0 0

run 100000

######## Step 2#########
unfix 2
unfix 5
unfix 6

fix 8 topwall move variable v_zero v_zero NULL v_zero v_zero NULL
fix 9 rightwall move variable v_zero v_zero NULL v_zero v_zero NULL
velocity walls zero angular
velocity walls zero linear

fix 10 granules viscous 10
run 2000000

I compute stresses using the following command:

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])/(20.040.4*4)

I don’t understand why step 2 will generate stresses in z-direction.

Any help will be much appreciated.


Are you trying to make stress-free equilibrium structure following these two steps?

I am wondering why you see xy shear stress at step 1?

If you do NPT at ambient conditions (say 1atm and 300K), all stress components should be zero once equilibrium is reached.


In fixes 8 and 9, your “NULL” keywords will not keep the z-coordinates of those particles fixed at zero – it will integrate the z-coordinates as if under “fix nve”. Try replacing those with actual "0"s and see what happens.




I am trying to jam the granules so I want to stabilize the granules with a non-zero stress. I do that by confining the granules so that it tips beyond the jamming transition point.

Yes, you are correct, the shear stress is very small when compared with the other magnitudes like nearly three orders of magnitude less than the normal stress components. I think that is because I calculate stress about a boundary that doesn’t encompass all the granules (just excludes the ones near the wall) so I think it might be an artifact as a result of that. Without friction I know that my two normal stress components will be equal but with friction, they will be unequal. However for the latter, can there be a slight build up of shear stress due to friction?

Right now I have been using NVE because my system is athermal… Is it common to use NPT for granular systems as I thought it is mainly for thermal systems which I don’t think is my case?

I know this is not LAMMPs related but is there any literature on what I am trying to achieve?


Hi Dr Shern,

Since I am doing a 2d simulation, I cannot specify 0 for those components. I get the following error:
Fix move cannot define z or vz variable for 2d problem (…/fix_move.cpp:192)



Thanks for your suggestion. I will try that :slight_smile:

I know this is not LAMMPs related but is there any literature on what I am trying to achieve?

I don’t know.

Anyway, I would prefer NPT to develop a stress-free system at the right mass density.

Once you have the correct initial structure, you could use NVE.


See this paper. Molecules are randomly packed at low density. The model is then simulated with NPT at high temp and pressure to get rid of voids/empty spaces followed by NPT at ambient conditions to have the right denisty.