Recently, I studied a system with a granular monolayer. The container is vibrated in the x direction by using the “fix wall/gran/region” command.

Although the gravity exists in the system, it should be balanced with the normal force given by the botton wall.

The simulation shows that the granules will leave the floor if the amplitude and the frequency satisfies (roughly) Aω²/g>1. This is purely straight forward since we can imagine what will happen if we squeeze a ball by another two balls on the floor.

However, the command “pair_style granular” should not provide any extra force in z direction in this case.

My opinion is the floating point error provide some fluctuation in z-position such that cause the excitation of granules.

Is my opinion right? Does LAMMPS remain trustworthy in this case?

Following is my script:

```
variable freq index 3 4 5 6 7 8 9 10 11 12
units cgs
atom_style sphere
boundary f f f
newton off
comm_modify vel yes
timestep 0.000001
# Global parameters
variable PI equal 3.14159265358979
# Input the initial condition and some setting
read_data ./C_code/init.data
fix 1 all nve/sphere
thermo_style custom step atoms ke
thermo 1000
# Force
## Parameters for forces
variable r1 equal 0.106
variable r2 equal 0.106
variable reff equal ${r1}*${r2}/(${r1}+${r2})
neighbor ${r2} bin
neigh_modify every 1 delay 0
variable rho1 equal 0.2
variable rho2 equal 8.8
variable m1 equal 4/3*${PI}*${r1}*${r1}*${r1}*${rho1}
variable m2 equal 4/3*${PI}*${r2}*${r2}*${r2}*${rho2}
variable meff equal ${m1}*${m2}/(${m1}+${m2})
variable mu1 equal 0.2 # Poisson's ratio
variable mu2 equal 0.2
variable E1 equal 1e11 # Young's modulus
variable E2 equal 1e11
variable Eeff11 equal 1/((1-${mu1}*${mu1})/${E1}+(1-${mu1}*${mu1})/${E1}) # Effective E for the ij pair
variable Eeff22 equal 1/((1-${mu2}*${mu2})/${E2}+(1-${mu2}*${mu2})/${E2})
variable Eeff12 equal 1/((1-${mu1}*${mu1})/${E1}+(1-${mu2}*${mu2})/${E2})
variable Geff11_8 equal 8/(2*(2-${mu1})*(1+${mu1})/${E1}+2*(2-${mu1})*(1+${mu1})/${E1}) # Effective G times 8 for ij
variable Geff22_8 equal 8/(2*(2-${mu2})*(1+${mu2})/${E2}+2*(2-${mu2})*(1+${mu2})/${E2})
variable Geff12_8 equal 8/(2*(2-${mu1})*(1+${mu1})/${E1}+2*(2-${mu2})*(1+${mu2})/${E2})
variable CoeRes equal 1
variable xmu11 equal 0.9 # Glass-glass
variable xmu22 equal 1.1 # Iron-iron
variable xmu12 equal 0.6 # Glass-iron
## Setting
# pair_coeff * * hertz/material 1e8 0.3 0.3 tangential mindlin 1e5 1.0 0.4 damping tsuji
pair_style granular
pair_coeff * * hertz/material ${Eeff11} ${CoeRes} ${mu1} tangential mindlin ${Geff11_8} 2.0 ${xmu11} damping tsuji
fix grav all gravity 980 vector 0 0 -1
variable amp equal 1.0
variable T equal 1/${freq}
variable shake_x equal ${amp}*sin(2*${PI}*${freq}*elaplong*dt)
region container block -9 9 -2 2 -1 1
region shake_box block -9 9 -2 2 -1 1 side in move v_shake_x NULL NULL
fix 4 all wall/gran/region granular hertz/material ${Eeff11} ${CoeRes} ${mu1} tangential mindlin ${Geff11_8} 2.0 ${xmu11} damping tsuji region shake_box
dump 1 all custom 100000 ./dump_1.0_${freq}/*.dump id type mass diameter x y z vx vy vz
run 70000000 # 70 seconds
undump 1
# -------------------------
clear
next freq
jump ./Hvibro_collide_10.lmp
```