Discrepancy with particle forces and pair forces

LAMMPS LAMMPS Mailing List Mirror
1

Dear LAMMPS users,

I am colliding two granular particles (LAMMPS downloaded April 2013) in normal (along X axis) and oblique (just off X axis - see script below) configurations and outputting:

  1. The particle FX FY FZ components via dump custom fx fy fz etc

  2. The pair FX FY FZ components from a compute pair/local fx fy fz etc command

With no other forces acting i.e. no gravity etc, the only force acting on the particles is the contact force, and therefore I would expect 1) and 2) to give the same values. I have investigated gran/hooke and gran/hooke/history pair styles and I find that:

a) In the case of damping = 0 and friction coeff = 0, the two outputs are identical

b) In the case of damping = 0 and friction coeff > 0, the two outputs are identical for the normal collision, but for the oblique collision I find a small discrepancy (approx 1%) in the FX components between output methods, and a large discrepancy (approx 50%) in the FY components.

c) In the case of damping > 0 and friction coeff = 0, there is a small discrepancy (<1%) in FX for the normal case, and a small discrepancy (<1%) in FX and FY for the oblique case.

d) In the case of damping > 0 and friction coeff > 0, there is a small discrepancy (<1%) in FX for the normal case, and for the oblique collision I find a small discrepancy (approx 1%) in the FX components between output methods, and a large discrepancy (approx 50%) in the FY components.

Please see an input script below that demonstrates the oblique case.

Has anyone experienced this or have an explanation for why the FX FY components are not always the same when output as pairs or as atom properties from dump custom?

Many thanks,

Chris

atom_style sphere

boundary f f f

newton off

communicate single vel yes

region reg block -10 10 -10 10 -10 10 units box

create_box 2 reg

create_atoms 1 single -0.5 0.1 0.0 units box

create_atoms 2 single 0.5 0.0 0.0 units box

group left type 1

group right type 2

neighbor 0.2 bin

neigh_modify delay 0

pair_style gran/hooke/history 200000 NULL 100 NULL 1 1

pair_coeff * *

timestep 3.6317e-5

velocity left set 5.0 0. 0. units box

velocity right set -5.0 0. 0. units box

compute pairs all pair/local fx fy fz

fix 1 all nve/sphere

compute 1 all erotate/sphere

thermo_style custom step atoms ke c_1 vol

thermo 1

thermo_modify lost ignore norm no

compute_modify thermo_temp dynamic yes

dump pairdump all local 1 dump.pairs c_pairs[1] c_pairs[2] c_pairs[3]

dump id all custom 1 dump.two_elastic id type fx fy fz

run 100

2

Jeremy might be able to comment on this.

Steve

3

Chris,

Did you happen to compare these values to the analytic solutions? In any case, I believe that pair/local includes tangential forces, while the particle dump only returns normal forces, see http://lammps.sandia.gov/doc/pair_gran.html in the Mixing, shift, table, tail correction, restart, rRESPA info entry. This makes sense why you see 50% difference, when your coefficient of friction is set to 1.

4

Hi Eric,

Thanks for the helpful response

  1. I think the dump custom fx fy fz part is calling atom->f[i][0], atom->f[i][1] etc and should therefore be the total force on each particle, including normal and tangential?

  2. As for pair/local, this is calling the “single” part of the gran pair style, so, referencing http://lammps.sandia.gov/doc/pair_gran.html, I think a command of pair/local fx fy fz will only output the normal components, and I must include p1 p2 p3 p4 to get the tangential part - perhaps this is the source of the discrepancy - I will have a look

Many thanks,

Chris

5

Eric and others,

Thanks for your help so far -

I am now colliding two granular particles (LAMMPS downloaded April 2013) in normal (along X axis) and oblique (just off X axis - see script below) configurations and outputting:

  1. The particle FX FY FZ components via dump custom fx fy fz etc

  2. The pair force components from a compute pair/local fx fy fz p1 p2 p3 command

Including p1 p2 p3 in pair/local has solved the discrepancy associated with the tangential force but not those associated with damping.

See my original post - cases a) and b) are now solved, thanks to including p1 p2 p3 in the pair/local output. However, with non-zero damping, I still see the discrepancies described before in c), d).

Any further comment would be appreciated.

Chris

6

Dear Chris,

Thanks for pointing this out.

Would it be possible for you to download and run with the newest version so as to make sure that the difference is reproducible there?

Also is the difference still approximately 1%? Or is it much smaller? Sorry I cannot run LAMMPS at the moment to check myself.

The difference probably comes about due to the fact that one value is calculated by the compute function and the other by the single. So there is always a small danger that the code in single and compute does not match 100%.

As for your original point maybe the manual for compute_pair_local should be changed slightly to highlight this?

Steve I suggest this, but please feel free to modify:

“The output force is the normal force acting between the pair of atoms, which is positive for a repulsive force and negative for an attractive force. The outputs fx, fy, and fz are the xyz components of the normal force on atom I. Note that for most pair styles the only force that acts between atoms is the normal force, so that normal force is equal to the total force. Exceptions to this are some granular pairs, where in order to calculate total forces between atoms the shear forces should be added too (also available through arguments p1 p2 p3).”

Thanks,

George

http://users.otenet.gr/~marketop/

7

Hi George, thanks for your reply.

I have just run the same case (input script below) on the June 17 version of LAMMPS and find the same discrepancy.

I would expect e.g. the FX component output from compute stress/atom to be equal to the FX+p1 component output from pair/local.

In the case below:

For damping = 0, there is a discrepancy of approx 0.0001%
For damping = 100, there is a discrepancy of approx 0.2%
For damping = 1000, there is a discrepancy of approx 2%
For damping = 10000, there is a discrepancy of approx 20%

Yes, I agree that the difference between these calculations is that pair/local is using the “single” function - I was hoping someone might shed some light on exactly where the difference arises.

Thanks again for your comments and I look forward to your reply,

Chris

atom_style sphere

boundary f f f

newton off

communicate single vel yes

region reg block -10 10 -10 10 -10 10 units box

create_box 2 reg

create_atoms 1 single 0. 0.1 0 units box

create_atoms 2 single 0.5 0.0 0.0 units box

group left type 1

group right type 2

neighbor 0.2 bin

neigh_modify delay 0

pair_style gran/hooke/history 20000 NULL 1000 NULL 1 1

pair_coeff * *

timestep 3.6317e-5

velocity left set 5.0 0. 0. units box

#velocity right set -5.0 0. 0. units box

compute pairs all pair/local fx fy fz p1 p2 p3 p4

fix 1 all nve/sphere

compute 1 all erotate/sphere

thermo_style custom step atoms ke c_1 vol

thermo 1

thermo_modify lost ignore norm no

compute_modify thermo_temp dynamic yes

dump id all custom 1 dump.two_elastic id type fx fy fz

dump pairdump all local 1 dump.pairs c_pairs[1] c_pairs[2] c_pairs[3] c_pairs[4] c_pairs[5] c_pairs[6] c_pairs[7]

run 100

8

I didn’t understand your first email. I thought you
were comparing forces on the 2 particles.

George’s point is correct:

The difference probably comes about due to the fact that one value is calculated by the compute function and the other by the
single. So there is always a small danger that the code in single and compute does not match 100%.

This is not a “danger” or a bug. It is due to the fact that the pair style
is velocity dependent. The velocities are updated at the end
of the timestep (Verlet), so the computation of forces in the middle
of the timestep (compute function) is different than at the end
of the timestep (single function invoked for output for the dump).

Steve

9

Steve,

Does this go for the pressure compute as well?

10

Not sure what you mean. It is true
that the pressure for a granular model
(e.g. printed with thermo output) will have
a virial contribution from forces in the middle
of the timetstep (using velocities at the mid pt),
and a kinetic contribution from the velocities
at the end of the timestep. Thus it is slightly
inconsitent. Emphasis on slightly.

Steve

11

You got at the question I was asking. Thanks again.