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
Not to my understanding of the physics of DEM as a non-expert. Why should the granular particles not move into z direction? They are not point particles but extended particles with a radius a surface and friction. So if you shake them along in x, they will start rotating and if two particles touch each other the rotation could lead to one lifting the other in z direction from the friction.
@jtclemm may need to correct my rationalization and have a different explanation.
But in any case, without seeing a usable input example input deck, it will be difficult to verify your claims and give an accurate assessment.
I agree, an input script is needed. There’s too little information to comment on. E.g., are grains monodisperse? Are they poured into the container so they have slight variations in their z position? Am I correct to assume gravity is along z?
I know grains should move upward naturally. I mean it should not move in z direction in LAMMPS since no net z direction force exists in the calculation.
To my knowledge, LAMMPS regards grains as points. The radius is used to calculate the force when two grains interact each other by using pair_style granular. Thus If no net displacement in z direction at t=0, we should not see grains move to the air in the simulation. Although the rotation is involved in the calculation, it should not affect the displacement in z direction.
Sorry about that. I pasted my script in the original paste.
Grains are monodisperse. The initial state is prepared by another C code. So they are just suddenly show up at equilibrium positions. Such that no net z-displacement between two grains at initial. The gravity is along the z axis.
Thanks for posting an input script. If your horizontal wall has friction and particles are perpendicularly falling in the z direction, will this not produce frictional forces along z producing variations in the particles’ z coordinates?
I agree with @jtclemm. What you said in your initial post here:
is not a general rule of physics. It is an explanation of specific circumstances, namely when an object is observed to be at rest on a flat plane. Normal forces do not automatically equal gravitational forces, or else jumping (and rockets, airplanes and basketball) would be physically impossible.
I still do not see anything inconsistent from your description/script. Like @srtee said, the wall and gravitational forces are not guaranteed to cancel. The wall is inherently soft so gravity can pull particles downward creating some overlap. Are you trying to initialize the system with a z coordinate such that gravity exactly cancels the wall’s force? As the particles move laterally, are there not also frictional forces from the floor which induce rotation creating more frictional forces between particles?
If I’m misunderstanding, do you have a small data file to go with the input file that illustrates exactly what your concern is?
Thanks for replying. I will try my best to explain your questions.
I have another script which just makes granules rest on the floor, and the z-positions of them stay the same. Since the coefficient of restitution is 1, I believe the gravity and the hertz force are balanced at initial.
The rotation is not included. There is no parameter of rotation inside the script.
I still think the difference of z-positions should come from the floating-point-error. However, I can’t prove it.
Now I just replace the mathematical wall by the real wall. So I think I can just move on.
If you are settling particles on the floor, I think you’ll find their vertical velocities decay but are unlikely to ever exactly reach zero. This small residual velocity would be enough to create additional tangential damping forces (even without friction) along the walls creating variation the vertical position of particles. However, you also do have rotation in your script: atom style sphere, fix nve/sphere, and a granular pair style are all used for extended finite-sized particles (see 8.5.1. Finite-size spherical and aspherical particles — LAMMPS documentation).
I"m not sure what you mean by replacing the mathematical wall with a real wall, but if you don’t want any vertical motion (which is physically expected in your setup), you could always run a 2D simulation.
