[BUG] Inconsistency using fix wall/gran with the granular pairstyle

[kjhanley]

Summary

The granular pairstyle allows for the mixing of coefficients specified via pair_coeff. With the hertz/material normal contact model, a particle with E = 400 in contact with a particle with E = 600, both with the same Poisson’s ratio of 0.2, would have an effective Young’s modulus of 250. That’s fine for particles in contact.

However, doing the same using a particle and a wall (fix wall/gran, granular pairstyle, hertz/material normal contact model) gives an unexpected value of the effective Young’s modulus.

LAMMPS Version and Platform

latest development branch (17 April 2024); reproducible on multiple platforms

Expected Behavior

For particle–wall contacts using pair_style granular, the documentation (https://docs.lammps.org/fix_wall_gran.html) says:

"For granular, fstyle_params are set using the same syntax as for the pair_coeff command of pair_style granular" and "The equation for the force between the wall and particles touching it is the same as the corresponding equation on the pair_style gran/* and pair_style granular doc pages, in the limit of one of the two particles going to infinite radius and mass (flat wall)".

Therefore, the expectation is the same mixing rule would be applied for Young’s modulus as for particle–particle contacts. The documentation indicates that the parameters provided in the fix wall/gran command are applied to the wall, treating it as a particle of infinite radius, and these are mixed appropriately with the equivalent parameters for the contacting particles.

For example, a user would expect the ‘5000’ in the line below to be interpreted as the wall’s Young’s modulus:

fix 1 all wall/gran granular hertz/material 5000 0.2 0.2…

They would also expect the same mixing rule applied for particle–particle contacts to be applied to particle–wall contacts, so that a particle with E = 400 in contact with this wall would have an effective Young’s modulus of ~386.

Actual Behavior

The particle properties are ignored for particle–wall contacts when pair_style granular is used. The simple provided script demonstrates this. In this script, two particles with different Young’s moduli are placed in contact with a wall that is using pair_style granular. Both particles have the same overlap. Exactly the same contact forces are calculated for both. This is not what one would expect. By analogy with the particle–particle contact case, one would expect the contact force on each particle to be different because the particles have different stiffnesses. Evidently the wall is not being treated like a particle of infinite radius as the documentation implies.

The effective Young’s modulus for a particle–wall contact is the value for a wall–wall contact, being calculated only from the values provided in fix wall/gran and not taking the particle properties into consideration. For the example above with a wall Young’s modulus of 5000, the effective Young’s modulus is calculated within LAMMPS as ~2604 for particle–wall contacts (what should be the value for a wall–wall contact).

Knowing this, a user could specify hertz rather than hertz/material. The problem would still arise for jkr or dmt though. A potential workaround may be to have a different fix wall/gran command corresponding to each pair_coeff command for each wall boundary. Each atom type would be assigned to a per-type group and atoms within that group would only interact with a subset of walls. The numbers in each fix wall/gran command could be manipulated to get the correct particle–wall contact stiffness for that atom type. This would be rather awkward though.

Ideally this workaround would be unnecessary if the same mixing rules were implemented as for particle–particle contacts, so that a user could specify a wall Young’s modulus and it would work in the expected way. At least the documentation should be clarified for fix wall/gran — it is currently quite misleading. This may be the more straightforward interim solution.

Steps to Reproduce

The attached script places two particles with different Young’s moduli in contact with a wall that is using pair_style granular and a hertz/material normal contact model. Both particles have the same overlap. Exactly the same contact forces are calculated for both, as shown in the attached dump file.

Further Information, Files, and Links

dump.granular_info.txt
in.granular_coeffs.txt

[dsbolin]

Hi Kevin,
Thanks for the very detailed description.

In short, I agree with everything you pointed out, especially that the solution (short-term, and arguably even long term) is a clarification in the documentation.

Slightly longer answer: the original intent/philosophy was to have material properties (Young's modulus, etc.) not necessarily associated with particle types, but rather associated with pairwise interactions, but use mixing rules based on Hertz theory when needed/unspecified. This is partly for consistency with how all other pair styles work in LAMMPS, and partly for additional flexibility to the user. So, with pair coeff, the intent is to allow the user to manually override the mixing rules for different types, e.g. pair coeff 1 2 would override the mixing of the type 1 and type 2 coefficients, even though that may not be consistent with the effective modulus from Hertz theory. For walls, the idea is similar, where the modulus is not the modulus of the wall, but the effective modulus associated with the interaction of the wall and the particle specified in the wall fix. As you correctly point out, if the effective modulus associated with Hertz contact is desired, the user can simply compute this for a given combination of a wall material modulus and particle type, and use multiple fix wall/gran commands for different particle types. While that can be slightly clunky in a script, it would not incur a noticeable performance penalty (and personally I've never needed more that 2 or 3 particle types in a granular simulation).

Please let me know if you have a specific suggestion about the rewording for the doc page, otherwise I can come up with something to clarify that the wall coefficients are associated with the wall-particle interaction rather than the material of the wall.

[kjhanley]

Hi Dan,

I see the intent — that makes sense. A clarification of the doc page would be the best solution then. I had a go at drafting an addition to the doc page, though it may be a little wordy. Please edit in any way you see fit, of course:

When fstyle granular is specified, the associated fstyle_params are taken as those for a wall–particle interaction. For example, for the hertz/material normal contact model with $E$ = 960 and $\nu$ = 0.2, the effective Young’s modulus for a wall–particle interaction is computed as $\frac{960}{2(1-0.2^2)}$ = 500. A mixing rule is not applied which implicitly takes the particles' properties into consideration.

If the simulation contains particles with different Young's moduli, and the correct effective moduli for wall–particle interactions are sought, the user may define multiple fix wall/gran commands: one (or more) for each atom type. The fix wall/gran command(s) corresponding to each atom type would be specified with the Young's modulus that would yield the desired wall–particle effective Young's modulus.