[lammps-users] A problem with friction? - Hertzian contact

Dear Lammps users,

I have been running a granular simulation with Hertz contact on the
May 21 build with updates through Aug 4, 2008. The problem is that I
have a random jam packed group of atoms (packing density 63.99%), but
when subjected to various loads, such as an indenter, it acts almost
like a fluid. There appears to be almost no resistance to an indenter
at all since the atoms can very easily rearrange themselves. I even
implemented cohesion, but the problem remains. I'm starting to get
suspicious about the friction coefficient.

To test friction, I wrote a simple input file to press one particle
down on another one, which rests on a frictional wall. The friction
coefficient for the atoms and the wall is 0.9 (normally I use 0.5 but
I was testing several cases). I applied a normal compressive force,
F, to the top particle so that it presses down on the bottom particle.
I ran the simulation for enough timesteps to ensure equilibrium (no
oscillations). I also checked that the displacement of both particles
was correct according to Hertz theory (as expected, it was correct).
Then I tried adding a tangential force to either the top or bottom
particle while the compressive force was still in effect. From the
2001 publication by Silbert et al. in Phys Rev E referenced in the
pair_style gran_hertzian page in the manual, I would expect that the
particle "slips" with a tangential force exceeding 0.9*F (since mu in
the case was 0.9). However, even with only a tangential force of
0.05*F, the particle still slips. So in this simple 2-particle
simulation, there seems to be no tangential force small enough so that
the static yield criterion is not met. In other words the particle
always slips.

If anyone can tell me where I've gone wrong, I'd be much appreciative.
The 2-particle input file is attached below.

Sincerely,
Priscilla Fonseca
PhD candidate
Northwestern University

Input file:
# Apply external compressive force to top sphere, then test friction

atom_style granular
units cgs
boundary f f f
newton off
neighbor 1.5E-7 bin

region regionID block -5.5E-7 10.5E-7 -5.5E-7 10.5E-7 0
15E-7 units box
create_box 2 regionID
create_atoms 1 single 2.5E-7 2.5E-7 2.5E-7 units box
create_atoms 2 single 2.5E-7 2.5E-7 7.5E-7 units box

group bottom type 1
group top type 2

set group all diameter 0.0000005
set group bottom density 2.3
set group top density 2.3

pair_style gran/hertzian 160588353.237 0.5 0.9 1
timestep 1E-16

fix 1 all wall/gran zplane 0 NULL 10 0.9

fix 2 all nve/sphere
fix 3 top addforce 0 0 -0.0001605873
fix 4 all viscous 0.00000006
fix 5 all coord/original
compute c1 all displace/atom 5

thermo_style custom step atoms ke press temp pzz
thermo 10000

dump id1 all atom 10000 presscsh1.dump
dump_modify id1 scale no
dump id2 all custom 10000 presscsh2.dump tag z vz fz

run 2000000
fix 6 bottom addforce 0.000008029365 0 0
#Also, I have tried: fix 6 top addforce 0.000008029365 0 0
run 1000000

Leo or Gary - any comments on this post someone
made to the LAMMPS mailing list?

Thanks,
Steve

Your script specifies cgs units. In principle that should
be OK, but everyone I know of uses LJ units for granular
simulations. So I suggest you try to get your model
working in LJ units, using published parameter values from
the literature. Once it is OK, you can convert to CGS units.
The process for doing so can often be a bit tricky.

Steve