$$

pair_style granular command

Syntax

pair_style granular cutoff

• cutoff = global cutoff (optional). See discussion below.

Examples

pair_style granular
pair_coeff * * hooke 1000.0 50.0 tangential linear_nohistory 1.0 0.4 damping mass_velocity

pair_style granular
pair_coeff * * hooke 1000.0 50.0 tangential linear_history 500.0 1.0 0.4 damping mass_velocity

pair_style granular
pair_coeff * * hertz 1000.0 50.0 tangential mindlin 1000.0 1.0 0.4 limit_damping

pair_style granular
pair_coeff * * hertz/material 1e8 0.3 0.3 tangential mindlin_rescale NULL 1.0 0.4 damping tsuji

pair_style granular
pair_coeff 1 * jkr 1000.0 500.0 0.3 10 tangential mindlin 800.0 1.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall
pair_coeff 2 2 hertz 200.0 100.0 tangential linear_history 300.0 1.0 0.1 rolling sds 200.0 100.0 0.1 twisting marshall

pair_style granular
pair_coeff 1 1 dmt 1000.0 50.0 0.3 0.0 tangential mindlin NULL 0.5 0.5 rolling sds 500.0 200.0 0.5 twisting marshall
pair_coeff 2 2 dmt 1000.0 50.0 0.3 10.0 tangential mindlin NULL 0.5 0.1 rolling sds 500.0 200.0 0.1 twisting marshall

pair_style granular
pair_coeff * * hertz 1000.0 50.0 tangential mindlin 1000.0 1.0 0.4 heat area 0.1

pair_style granular
pair_coeff * * mdr 5e6 0.4 1.9e5 2.0 0.5 0.5 tangential linear_history 940.0 0.0 0.7 rolling sds 2.7e5 0.0 0.6 damping none

Description

The granular styles support a variety of options for the normal, tangential, rolling and twisting forces resulting from contact between two granular particles. This expands on the options offered by the pair gran/* pair styles. The total computed forces and torques are the sum of various models selected for the normal, tangential, rolling and twisting modes of motion.

All model choices and parameters are entered in the pair_coeff command, as described below. Unlike e.g. pair gran/hooke, coefficient values are not global, but can be set to different values for different combinations of particle types, as determined by the pair_coeff command. If the contact model choice is the same for two particle types, the mixing for the cross-coefficients can be carried out automatically. This is shown in the last example, where model choices are the same for type 1 - type 1 as for type 2 - type2 interactions, but coefficients are different. In this case, the mixed coefficients for type 1 - type 2 interactions can be determined from mixing rules discussed below. For additional flexibility, coefficients as well as model forms can vary between particle types, as shown in the fourth example: type 1 - type 1 interactions are based on a Johnson-Kendall-Roberts normal contact model and 2-2 interactions are based on a DMT cohesive model (see below). In that example, 1-1 and 2-2 interactions have different model forms, in which case mixing of coefficients cannot be determined, so 1-2 interactions must be explicitly defined via the pair_coeff 1 * command, otherwise an error would result.

---

The first required keyword for the pair_coeff command is the normal contact model. Currently supported options for normal contact models and their required arguments are:

1. hooke : $k_n$, ${\eta }_{n0}$ (or $e$)

2. hertz : $k_n$, ${\eta }_{n0}$ (or $e$)

3. hertz/material : E, ${\eta }_{n0}$ (or $e$), $\nu$

4. dmt : E, ${\eta }_{n0}$ (or $e$), $\nu$, $\gamma$

5. jkr : E, ${\eta }_{n0}$ (or $e$), $\nu$, $\gamma$

6. mdr : $E$, $\nu$, $Y$, $\mathrm{\Delta }\gamma$, ${\psi }_b$, $e$

Here, $k_n$ is spring stiffness (with units that depend on model choice, see below); ${\eta }_{n0}$ is a damping prefactor (or, in its place a coefficient of restitution $e$, depending on the choice of damping mode, see below); E is Young’s modulus in units of force/length^2, i.e. pressure; $\nu$ is Poisson’s ratio and $\gamma$ is a surface energy density, in units of energy/length^2.

For the hooke model, the normal, elastic component of force acting on particle i due to contact with particle j is given by:

$${\mathbf{F}}_{ne,Hooke}=k_n{\delta }_{ij}\mathbf{n}$$

Where ${\delta }_{ij}=R_i+R_j-\Vert {\mathbf{r}}_{ij}\Vert$ is the particle overlap, $R_i,R_j$ are the particle radii, ${\mathbf{r}}_{ij}={\mathbf{r}}_i-{\mathbf{r}}_j$ is the vector separating the two particle centers (note the i-j ordering so that ${\mathbf{F}}_{ne}$ is positive for repulsion), and $\mathbf{n}=\frac{{\mathbf{r}}_{ij}}{\Vert {\mathbf{r}}_{ij}\Vert }$. Therefore, for hooke, the units of the spring constant $k_n$ are force/distance, or equivalently mass/time^2.

For the hertz model, the normal component of force is given by:

$${\mathbf{F}}_{ne,Hertz}=k_n{R}_{eff}^{1/2}{\delta }_{ij}^{3/2}\mathbf{n}$$

Here, $R_{eff}=R=\frac{R_i{R}_j}{R_i+R_j}$ is the effective radius, denoted for simplicity as R from here on. For hertz, the units of the spring constant $k_n$ are force/length^2, or equivalently pressure.

For the hertz/material model, the force is given by:

$${\mathbf{F}}_{ne,Hertz/material}=\frac{4}{3}{E}_{eff}{R}^{1/2}{\delta }_{ij}^{3/2}\mathbf{n}$$

Here, $E_{eff}=E={(\frac{1-{\nu }_i^2}{E_i}+\frac{1-{\nu }_j^2}{E_j})}^{-1}$ is the effective Young’s modulus, with ${\nu }_i,{\nu }_j$ the Poisson ratios of the particles of types i and j. $E_{eff}$ is denoted as E from here on. Note that if the elastic modulus and the shear modulus of the two particles are the same, the hertz/material model is equivalent to the hertz model with $k_n=4/3E$

The dmt model corresponds to the (Derjaguin-Muller-Toporov) cohesive model, where the force is simply Hertz with an additional attractive cohesion term:

$${\mathbf{F}}_{ne,dmt}=(\frac{4}{3}E{R}^{1/2}{\delta }_{ij}^{3/2}-4\pi \gamma R)\mathbf{n}$$

The jkr model is the (Johnson-Kendall-Roberts) model, where the force is computed as:

$${\mathbf{F}}_{ne,jkr}=(\frac{4E{a}^3}{3R}-2\pi a^2\sqrt{\frac{4\gamma E}{\pi a}})\mathbf{n}$$

Here, $a$ is the radius of the contact zone, related to the overlap $\delta$ according to:

$$\delta =a^2/R-2\sqrt{\pi \gamma a/E}$$

LAMMPS internally inverts the equation above to solve for a in terms of $\delta$, then solves for the force in the previous equation. Additionally, note that the JKR model allows for a tensile force beyond contact (i.e. for $\delta <0$), up to a maximum of $3\pi \gamma R$ (also known as the ‘pull-off’ force). Note that this is a hysteretic effect, where particles that are not contacting initially will not experience force until they come into contact $\delta \ge 0$; as they move apart and ($\delta <0$), they experience a tensile force up to $3\pi \gamma R$, at which point they lose contact.

The mdr model is a mechanically-derived contact model designed to capture the contact response between adhesive elastic-plastic particles into large deformation. The theoretical foundations of the mdr model are detailed in the two-part series Zunker and Kamrin Part I and Zunker and Kamrin Part II. Further development and demonstrations of its application to industrially relevant powder compaction processes are presented in Zunker et al..

The model requires the following inputs:

1. Young’s modulus $E>0$ : The Young’s modulus is commonly reported for various powders.

2. Poisson’s ratio $0\le \nu \le 0.5$ : The Poisson’s ratio is commonly reported for various powders.

3. Yield stress $Y\ge 0$ : The yield stress is often known for powders composed of materials such as metals but may be unreported for ductile organic materials, in which case it can be treated as a free parameter.

4. Effective surface energy $\mathrm{\Delta }\gamma \ge 0$ : The effective surface energy for powder compaction applications is most easily determined through its relation to the more commonly reported critical stress intensity factor $K_{Ic}=\sqrt{2\mathrm{\Delta }\gamma E/(1-{\nu }^2)}$.

5. Critical confinement ratio $0\le {\psi }_b\le 1$ : The critical confinement ratio is a tunable parameter that determines when the bulk elastic response is triggered. Lower values of ${\psi }_b$ delay the onset of the bulk elastic response.

6. Coefficient of restitution $0\le e\le 1$ : The coefficient of restitution is a tunable parameter that controls damping in the normal direction.

Note

The values for $E$, $\nu$, $Y$, and $\mathrm{\Delta }\gamma$ (i.e., $K_{Ic}$) should be selected for zero porosity to reflect the intrinsic material property rather than the bulk powder property.

The mdr model produces a nonlinear force-displacement response, therefore the critical timestep $\mathrm{\Delta }t$ depends on the inputs and level of deformation. As a conservative starting point the timestep can be assumed to be dictated by the bulk elastic response such that $\mathrm{\Delta }t=0.35\sqrt{m/k_{\text{bulk}}}$, where $m$ is the mass of the smallest particle and $k_{\text{bulk}}=\kappa R_{\text{min}}$ is an effective stiffness related to the bulk elastic response. Here, $\kappa =E/(3(1-2\nu ))$ is the bulk modulus and $R_{\text{min}}$ is the radius of the smallest particle.

Note

The mdr model requires some specific settings to function properly, please read the following text carefully to ensure all requirements are followed.

The atom_style must be set to sphere 1 to enable dynamic particle radii. The mdr model is designed to respect the incompressibility of plastic deformation and inherently tracks free surface displacements induced by all particle contacts. In practice, this means that all particles begin with an initial radius, however as compaction occurs and plastic deformation is accumulated, a new enlarged apparent radius is defined to ensure that that volume change due to plastic deformation is not lost. This apparent radius is stored as the atom radius meaning it is used for subsequent neighbor list builds and contact detection checks. The advantage of this is that multi-neighbor dependent effects such as formation of secondary contacts caused by radial expansion are captured by the mdr model. Setting atom_style sphere 1 ensures that updates to the particle radii are properly reflected throughout the simulation.

atom_style sphere 1

Newton’s third law must be set to off. This ensures that the neighbor lists are constructed properly for the topological penalty algorithm used to screen for non-physical contacts occurring through obstructing particles, an issue prevalent under large deformation conditions. For more information on this algorithm see Zunker et al..

newton off

The damping model must be set to none. The mdr model already has a built in damping model.

pair_coeff * * mdr 5e6 0.4 1.9e5 2 0.5 0.5 damping none

The definition of multiple mdr models in the pair_style is currently not supported. Similarly, the mdr model cannot be combined with a different normal model in the pair_style. Physically this means that only one homogeneous collection of particles governed by a single mdr model is allowed.

The mdr model currently only supports fix wall/gran/region, not fix wall/gran. If the mdr model is specified for the pair_style any fix wall/gran/region commands must also use the mdr model. Additionally, the following mdr inputs must match between the pair_style and fix wall/gran/region definitions: $E$, $\nu$, $Y$, ${\psi }_b$, and $e$. The exception is $\mathrm{\Delta }\gamma$, which may vary, permitting different adhesive behaviors between particle-particle and particle-wall interactions.

Note

The mdr model has a number of custom property/atom and pair/local definitions that can be called in the input file. The useful properties for visualization and analysis are described below.

In addition to contact forces the mdr model also tracks the following quantities for each particle: elastic volume change, average normal stress components, total surface area involved in contact, and individual contact areas. In the input script, these quantities are initialized by calling run 0 and can then be accessed using subsequent compute commands. The last compute command uses pair/local p13 to calculate the pairwise contact areas for each active contact in the group-ID. Due to the use of an apparent radius in the mdr model, the keyword/arg pair cutoff radius must be specified for pair/local to properly detect existing contacts.

run 0
compute ID group-ID property/atom d_Velas
compute ID group-ID property/atom d_sigmaxx
compute ID group-ID property/atom d_sigmayy
compute ID group-ID property/atom d_sigmazz
compute ID group-ID property/atom d_Acon1
compute ID group-ID pair/local p13 cutoff radius

Note

The mdr model has two example input scripts within the examples/granular directory. The first is a die compaction simulation involving 200 particles named in.tableting.200. The second is a triaxial compaction simulation involving 12 particles named in.triaxial.compaction.12.

---

In addition, the normal force is augmented by a damping term of the following general form:

$${\mathbf{F}}_{n,damp}=-{\eta }_n{\mathbf{v}}_{n,rel}$$

Here, ${\mathbf{v}}_{n,rel}=({\mathbf{v}}_j-{\mathbf{v}}_i)\cdot \mathbf{n}\mathbf{n}$ is the component of relative velocity along $\mathbf{n}$.

The optional damping keyword to the pair_coeff command followed by a keyword determines the model form of the damping factor ${\eta }_n$, and the interpretation of the ${\eta }_{n0}$ or $e$ coefficients specified as part of the normal contact model settings. The damping keyword and corresponding model form selection may be appended anywhere in the pair coeff command. Note that the choice of damping model affects both the normal and tangential damping (and depending on other settings, potentially also the twisting damping). The options for the damping model currently supported are:

1. velocity

2. mass_velocity

3. viscoelastic

4. tsuji

5. coeff_restitution

If the damping keyword is not specified, the viscoelastic model is used by default.

For damping velocity, the normal damping is simply equal to the user-specified damping coefficient in the normal model:

$${\eta }_n={\eta }_{n0}$$

Here, ${\eta }_{n0}$ is the damping coefficient specified for the normal contact model, in units of mass/time.

For damping mass_velocity, the normal damping is given by:

$${\eta }_n={\eta }_{n0}{m}_{eff}$$

Here, ${\eta }_{n0}$ is the damping coefficient specified for the normal contact model, in units of 1/time and $m_{eff}=m_i{m}_j/(m_i+m_j)$ is the effective mass. Use damping mass_velocity to reproduce the damping behavior of pair gran/hooke/*.

The damping viscoelastic model is based on the viscoelastic treatment of (Brilliantov et al), where the normal damping is given by:

$${\eta }_n={\eta }_{n0}a{m}_{eff}$$

Here, a is the contact radius, given by $a=\sqrt{R\delta }$ for all models except jkr, for which it is given implicitly according to $\delta =a^2/R-2\sqrt{\pi \gamma a/E}$. For damping viscoelastic, ${\eta }_{n0}$ is in units of 1/(time*distance).

The tsuji model is based on the work of (Tsuji et al). Here, the damping coefficient specified as part of the normal model is interpreted as a restitution coefficient $e$. The damping constant ${\eta }_n$ is given by:

$${\eta }_n=\alpha (m_{eff}{k}_{nd}{)}^{1/2}$$

where $k_{nd}$ is an effective harmonic stiffness equal to the ratio of the normal force to the overlap. For example, $k_{nd}=4/3Ea$ for a Hertz contact model based on material parameters with $a$ being the contact radius of $\sqrt{\delta R}$. For Hooke, $k_{nd}$ is simply the spring constant or $k_n$. This damping model is not compatible with cohesive normal models such as JKR or DMT. The parameter $\alpha$ is related to the restitution coefficient e according to:

$$\alpha =1.2728-4.2783e+11.087e^2-22.348e^3+27.467e^4-18.022e^5+4.8218e^6$$

The dimensionless coefficient of restitution $e$ specified as part of the normal contact model parameters should be between 0 and 1, but no error check is performed on this.

The coeff_restitution model is useful when a specific normal coefficient of restitution $e$ is required. It operates much like the Tsuji model but, the normal coefficient of restitution $e$ is specified as an input in place of the usual ${\eta }_{n0}$ value in the normal model. Following the approach of (Brilliantov et al), when using the hooke normal model, coeff_restitution then calculates the damping coefficient as:

$${\eta }_n=\sqrt{\frac{4m_{eff}{k}_{nd}}{1+{\left(\frac{\pi }{\log (e)}\right)}^2}},$$

where $k_{nd}$ is the same stiffness defined in the above Tsuji model. For any other normal model, e.g. the hertz and hertz/material models, the damping coefficient is:

$${\eta }_n=-2\sqrt{\frac{5}{6}}\frac{\log (e)}{\sqrt{{\pi }^2+(\log (e){)}^2}}\sqrt{\frac{3}{2}{k}_{nd}{m}_{eff}},$$

Since coeff_restitution accounts for the effective mass, effective radius, and pairwise overlaps (except when used with the hooke normal model) when calculating the damping coefficient, it accurately reproduces the specified coefficient of restitution for both monodisperse and polydisperse particle pairs. This damping model is not compatible with cohesive normal models such as JKR or DMT.

The total normal force is computed as the sum of the elastic and damping components:

$${\mathbf{F}}_n={\mathbf{F}}_{ne}+{\mathbf{F}}_{n,damp}$$

---

The pair_coeff command also requires specification of the tangential contact model. The required keyword tangential is expected, followed by the model choice and associated parameters. Currently supported tangential model choices and their expected parameters are as follows:

1. linear_nohistory : $x_{\gamma ,t}$, ${\mu }_s$

2. linear_history : $k_t$, $x_{\gamma ,t}$, ${\mu }_s$

3. mindlin : $k_t$ or NULL, $x_{\gamma ,t}$, ${\mu }_s$

4. mindlin/force : $k_t$ or NULL, $x_{\gamma ,t}$, ${\mu }_s$

5. mindlin_rescale : $k_t$ or NULL, $x_{\gamma ,t}$, ${\mu }_s$

6. mindlin_rescale/force : $k_t$ or NULL, $x_{\gamma ,t}$, ${\mu }_s$

Here, $x_{\gamma ,t}$ is a dimensionless multiplier for the normal damping ${\eta }_n$ that determines the magnitude of the tangential damping, ${\mu }_t$ is the tangential (or sliding) friction coefficient, and $k_t$ is the tangential stiffness coefficient.

For tangential linear_nohistory, a simple velocity-dependent Coulomb friction criterion is used, which mimics the behavior of the pair gran/hooke style. The tangential force ${\mathbf{F}}_t$ is given by:

$${\mathbf{F}}_t=-min({\mu }_t{F}_{n0},\Vert {\mathbf{F}}_{\mathrm{t},\mathrm{damp}}\Vert )\mathbf{t}$$

The tangential damping force ${\mathbf{F}}_{\mathrm{t},\mathrm{damp}}$ is given by:

$${\mathbf{F}}_{\mathrm{t},\mathrm{damp}}=-{\eta }_t{\mathbf{v}}_{t,rel}$$

The tangential damping prefactor ${\eta }_t$ is calculated by scaling the normal damping ${\eta }_n$ (see above):

$${\eta }_t=-x_{\gamma ,t}{\eta }_n$$

The normal damping prefactor ${\eta }_n$ is determined by the choice of the damping keyword, as discussed above. Thus, the damping keyword also affects the tangential damping. The parameter $x_{\gamma ,t}$ is a scaling coefficient. Several works in the literature use $x_{\gamma ,t}=1$ (Marshall, Tsuji et al, Silbert et al). The relative tangential velocity at the point of contact is given by ${\mathbf{v}}_{t,rel}={\mathbf{v}}_t-(R_i{\mathbf{\Omega }}_i+R_j{\mathbf{\Omega }}_j)\times \mathbf{n}$, where ${\mathbf{v}}_t={\mathbf{v}}_r-{\mathbf{v}}_r\cdot \mathbf{n}\mathbf{n}$, ${\mathbf{v}}_r={\mathbf{v}}_j-{\mathbf{v}}_i$ . The direction of the applied force is $\mathbf{t}={\mathbf{v}}_{\mathbf{t}\mathbf{,}\mathbf{rel}}/\Vert {\mathbf{v}}_{\mathbf{t}\mathbf{,}\mathbf{rel}}\Vert$ .

The normal force value $F_{n0}$ used to compute the critical force depends on the form of the contact model. For non-cohesive models (hertz, hertz/material, hooke), it is given by the magnitude of the normal force:

$$F_{n0}=\Vert {\mathbf{F}}_n\Vert$$

For cohesive models such as jkr and dmt, the critical force is adjusted so that the critical tangential force approaches ${\mu }_t{F}_{pulloff}$, see Marshall, equation 43, and Thornton. For both models, $F_{n0}$ takes the form:

$$F_{n0}=\Vert {\mathbf{F}}_{ne}+2F_{pulloff}\Vert$$

Where $F_{pulloff}=3\pi \gamma R$ for jkr, and $F_{pulloff}=4\pi \gamma R$ for dmt.

The remaining tangential options all use accumulated tangential displacement (i.e. contact history), except for the options mindlin/force and mindlin_rescale/force, that use accumulated tangential force instead, and are discussed further below. The accumulated tangential displacement is discussed in details below in the context of the linear_history option. The same treatment of the accumulated displacement applies to the other options as well.

For tangential linear_history, the tangential force is given by:

$${\mathbf{F}}_t=-min({\mu }_t{F}_{n0},\Vert -k_t\xi +{\mathbf{F}}_{\mathrm{t},\mathrm{damp}}\Vert )\mathbf{t}$$

Here, $\xi$ is the tangential displacement accumulated during the entire duration of the contact:

$$\xi ={\int }_{t0}^t{\mathbf{v}}_{t,rel}(\tau )\mathrm{d}\tau$$

This accumulated tangential displacement must be adjusted to account for changes in the frame of reference of the contacting pair of particles during contact. This occurs due to the overall motion of the contacting particles in a rigid-body-like fashion during the duration of the contact. There are two modes of motion that are relevant: the ‘tumbling’ rotation of the contacting pair, which changes the orientation of the plane in which tangential displacement occurs; and ‘spinning’ rotation of the contacting pair about the vector connecting their centers of mass ($\mathbf{n}$). Corrections due to the former mode of motion are made by rotating the accumulated displacement into the plane that is tangential to the contact vector at each step, or equivalently removing any component of the tangential displacement that lies along $\mathbf{n}$, and rescaling to preserve the magnitude. This follows the discussion in Luding, see equation 17 and relevant discussion in that work:

$$\xi =({\xi }^{\prime }-(\mathbf{n}\cdot {\xi }^{\prime })\mathbf{n})\frac{\Vert {\xi }^{\prime }\Vert }{\Vert {\xi }^{\prime }-(\mathbf{n}\cdot {\xi }^{\prime })\mathbf{n}\Vert }$$

Here, ${\xi }^{\prime }$ is the accumulated displacement prior to the current time step and $\xi$ is the corrected displacement. Corrections to the displacement due to the second mode of motion described above (rotations about $\mathbf{n}$) are not currently implemented, but are expected to be minor for most simulations.

Furthermore, when the tangential force exceeds the critical force, the tangential displacement is re-scaled to match the value for the critical force (see Luding, equation 20 and related discussion):

$$\xi =-\frac{1}{k_t}({\mu }_t{F}_{n0}\mathbf{t}-{\mathbf{F}}_{t,damp})$$

The tangential force is added to the total normal force (elastic plus damping) to produce the total force on the particle. The tangential force also acts at the contact point (defined as the center of the overlap region) to induce a torque on each particle according to:

$${\tau }_i=-(R_i-0.5\delta )\mathbf{n}\times {\mathbf{F}}_t$$

$${\tau }_j=-(R_j-0.5\delta )\mathbf{n}\times {\mathbf{F}}_t$$

For tangential mindlin, the Mindlin no-slip solution is used which differs from the linear_history option by an additional factor of $a$, the radius of the contact region. The tangential force is given by:

$${\mathbf{F}}_t=-min({\mu }_t{F}_{n0},\Vert -k_{t}a\xi +{\mathbf{F}}_{\mathrm{t},\mathrm{damp}}\Vert )\mathbf{t}$$

Here, $a$ is the radius of the contact region, given by $a=\sqrt{R\delta }$ for all normal contact models, except for jkr, where it is given implicitly by $\delta =a^2/R-2\sqrt{\pi \gamma a/E}$, see discussion above. To match the Mindlin solution, one should set $k_t=8G_{eff}$, where $G_{eff}$ is the effective shear modulus given by:

$$G_{eff}={(\frac{2-{\nu }_i}{G_i}+\frac{2-{\nu }_j}{G_j})}^{-1}$$

where $G_i$ is the shear modulus of a particle of type $i$, related to Young’s modulus $E_i$ and Poisson’s ratio ${\nu }_i$ by $G_i=E_i/(2(1+{\nu }_i))$. This can also be achieved by specifying NULL for $k_t$, in which case a normal contact model that specifies material parameters $E_i$ and ${\nu }_i$ is required (e.g. hertz/material, dmt or jkr). In this case, mixing of the shear modulus for different particle types i and j is done according to the formula above.

Note

The radius of the contact region $a$ depends on the normal overlap. As a result, the tangential force for mindlin can change due to a variation in normal overlap, even with no change in tangential displacement.

For tangential mindlin/force, the accumulated elastic tangential force characterizes the contact history, instead of the accumulated tangential displacement. This prevents the dependence of the tangential force on the normal overlap as noted above. The tangential force is given by:

$${\mathbf{F}}_t=-min({\mu }_t{F}_{n0},\Vert {\mathbf{F}}_{te}+{\mathbf{F}}_{\mathrm{t},\mathrm{damp}}\Vert )\mathbf{t}$$

The increment of the elastic component of the tangential force ${\mathbf{F}}_{te}$ is given by:

$$\mathrm{d}{\mathbf{F}}_{te}=-k_{t}a{\mathbf{v}}_{t,rel}\mathrm{d}\tau$$

The changes in frame of reference of the contacting pair of particles during contact are accounted for by the same formula as above, replacing the accumulated tangential displacement $\xi$, by the accumulated tangential elastic force $F_{te}$. When the tangential force exceeds the critical force, the tangential force is directly re-scaled to match the value for the critical force:

$${\mathbf{F}}_{te}=-{\mu }_t{F}_{n0}\mathbf{t}+{\mathbf{F}}_{t,damp}$$

The same rules as those described for mindlin apply regarding the tangential stiffness and mixing of the shear modulus for different particle types.

The mindlin_rescale option uses the same form as mindlin, but the magnitude of the tangential displacement is re-scaled as the contact unloads, i.e. if $a<a_{t_{n-1}}$:

$$\xi ={\xi }_{{\mathbf{t}}_{\mathbf{n}\mathbf{-}\mathbf{1}}}\frac{a}{a_{t_{n-1}}}$$

Here, $t_{n-1}$ indicates the value at the previous time step. This rescaling accounts for the fact that a decrease in the contact area upon unloading leads to the contact being unable to support the previous tangential loading, and spurious energy is created without the rescaling above (Walton ).

Note

For mindlin, a decrease in the tangential force already occurs as the contact unloads, due to the dependence of the tangential force on the normal force described above. By re-scaling $\xi$, mindlin_rescale effectively re-scales the tangential force twice, i.e., proportionally to $a^2$. This peculiar behavior results from use of the accumulated tangential displacement to characterize the contact history. Although mindlin_rescale remains available for historic reasons and backward compatibility purposes, it should be avoided in favor of mindlin_rescale/force.

The mindlin_rescale/force option uses the same form as mindlin/force, but the magnitude of the tangential elastic force is re-scaled as the contact unloads, i.e. if $a<a_{t_{n-1}}$:

$${\mathbf{F}}_{te}={\mathbf{F}}_{te,t_{n-1}}\frac{a}{a_{t_{n-1}}}$$

This approach provides a better approximation of the Mindlin-Deresiewicz laws and is more consistent than mindlin_rescale. See discussions in Thornton et al, 2013, particularly equation 18(b) of that work and associated discussion, and Agnolin and Roux, 2007, particularly Appendix A.

---

The optional rolling keyword enables rolling friction, which resists pure rolling motion of particles. The options currently supported are:

1. none

2. sds : $k_{roll}$, ${\gamma }_{roll}$, ${\mu }_{roll}$

If the rolling keyword is not specified, the model defaults to none.

For rolling sds, rolling friction is computed via a spring-dashpot-slider, using a ‘pseudo-force’ formulation, as detailed by Luding. Unlike the formulation in Marshall, this allows for the required adjustment of rolling displacement due to changes in the frame of reference of the contacting pair. The rolling pseudo-force is computed analogously to the tangential force:

$${\mathbf{F}}_{roll,0}=k_{roll}{\xi }_{roll}-{\gamma }_{roll}{\mathbf{v}}_{roll}$$

Here, ${\mathbf{v}}_{roll}=-R({\mathbf{\Omega }}_i-{\mathbf{\Omega }}_j)\times \mathbf{n}$ is the relative rolling velocity, as given in Wang et al and Luding. This differs from the expressions given by Kuhn and Bagi and used in Marshall; see Wang et al for details. The rolling displacement is given by:

$${\xi }_{roll}={\int }_{t_0}^t{\mathbf{v}}_{roll}(\tau )\mathrm{d}\tau$$

A Coulomb friction criterion truncates the rolling pseudo-force if it exceeds a critical value:

$${\mathbf{F}}_{roll}=min({\mu }_{roll}{F}_{n,0},\Vert {\mathbf{F}}_{roll,0}\Vert )\mathbf{k}$$

Here, $\mathbf{k}={\mathbf{v}}_{roll}/\Vert {\mathbf{v}}_{roll}\Vert$ is the direction of the pseudo-force. As with tangential displacement, the rolling displacement is rescaled when the critical force is exceeded, so that the spring length corresponds the critical force. Additionally, the displacement is adjusted to account for rotations of the frame of reference of the two contacting particles in a manner analogous to the tangential displacement.

The rolling pseudo-force does not contribute to the total force on either particle (hence ‘pseudo’), but acts only to induce an equal and opposite torque on each particle, according to:

$${\tau }_{roll,i}=R\mathbf{n}\times {\mathbf{F}}_{roll}$$

$${\tau }_{roll,j}=-{\tau }_{roll,i}$$

---

The optional twisting keyword enables twisting friction, which resists rotation of two contacting particles about the vector $\mathbf{n}$ that connects their centers. The options currently supported are:

1. none

2. sds : $k_{twist}$, ${\gamma }_{twist}$, ${\mu }_{twist}$

3. marshall

If the twisting keyword is not specified, the model defaults to none.

For both twisting sds and twisting marshall, a history-dependent spring-dashpot-slider is used to compute the twisting torque. Because twisting displacement is a scalar, there is no need to adjust for changes in the frame of reference due to rotations of the particle pair. The formulation in Marshall therefore provides the most straightforward treatment:

$${\tau }_{twist,0}=-k_{twist}{\xi }_{twist}-{\gamma }_{twist}{\mathrm{\Omega }}_{twist}$$

Here ${\xi }_{twist}={\int }_{t_0}^t{\mathrm{\Omega }}_{twist}(\tau )\mathrm{d}\tau$ is the twisting angular displacement, and ${\mathrm{\Omega }}_{twist}=({\mathbf{\Omega }}_i-{\mathbf{\Omega }}_j)\cdot \mathbf{n}$ is the relative twisting angular velocity. The torque is then truncated according to:

$${\tau }_{twist}=min({\mu }_{twist}{F}_{n,0},{\tau }_{twist,0})$$

Similar to the sliding and rolling displacement, the angular displacement is rescaled so that it corresponds to the critical value if the twisting torque exceeds this critical value:

$${\xi }_{twist}=\frac{1}{k_{twist}}({\mu }_{twist}{F}_{n,0}sgn({\mathrm{\Omega }}_{twist})-{\gamma }_{twist}{\mathrm{\Omega }}_{twist})$$

For twisting sds, the coefficients $k_{twist},{\gamma }_{twist}$ and ${\mu }_{twist}$ are simply the user input parameters that follow the twisting sds keywords in the pair_coeff command.

For twisting_marshall, the coefficients are expressed in terms of sliding friction coefficients, as discussed in Marshall (see equations 32 and 33 of that work):

$$k_{twist}=0.5k_t{a}^2$$

$${\eta }_{twist}=0.5{\eta }_t{a}^2$$

$${\mu }_{twist}=\frac{2}{3}a{\mu }_t$$

Finally, the twisting torque on each particle is given by:

$${\tau }_{twist,i}={\tau }_{twist}\mathbf{n}$$

$${\tau }_{twist,j}=-{\tau }_{twist,i}$$

---

If two particles are moving away from each other while in contact, there is a possibility that the particles could experience an effective attractive force due to damping. If the optional limit_damping keyword is used, this option will zero out the normal component of the force if there is an effective attractive force. This keyword cannot be used with the JKR or DMT models.

---

The optional heat keyword enables heat conduction. The options currently supported are:

1. none

2. radius : $k_s$

3. area : $h_s$

If the heat keyword is not specified, the model defaults to none. All heat models calculate an additional pairwise quantity accessible by the single() function (described below) which is the heat conducted between the two particles.

For heat radius, the heat $Q$ conducted between two particles is given by

$$Q=2k_{s}a\mathrm{\Delta }T$$

where $\mathrm{\Delta }T$ is the difference in the two particles’ temperature, $k_s$ is a non-negative numeric value for the conductivity (in units of power/(length*temperature)), and $a$ is the radius of the contact and depends on the normal force model. This is the model proposed by Vargas and McCarthy.

For heat area, the heat $Q$ conducted between two particles is given by

$$Q=h_{s}A\mathrm{\Delta }T$$

where $\mathrm{\Delta }T$ is the difference in the two particles’ temperature, $h_s$ is a non-negative numeric value for the heat transfer coefficient (in units of power/(area*temperature)), and $A=\pi a^2$ is the area of the contact and depends on the normal force model.

Note that the option none must either be used in all or none of of the pair_coeff calls. See fix heat/flow and fix property/atom for more information on this option.

---

The granular pair style can reproduce the behavior of the pair gran/* styles with the appropriate settings (some very minor differences can be expected due to corrections in displacement history frame-of-reference, and the application of the torque at the center of the contact rather than at each particle). The first example above is equivalent to pair gran/hooke 1000.0 NULL 50.0 50.0 0.4 1. The second example is equivalent to pair gran/hooke/history 1000.0 500.0 50.0 50.0 0.4 1. The third example is equivalent to pair gran/hertz/history 1000.0 500.0 50.0 50.0 0.4 1 limit_damping.

---

LAMMPS automatically sets pairwise cutoff values for pair_style granular based on particle radii (and in the case of jkr pull-off distances). In the vast majority of situations, this is adequate. However, a cutoff value can optionally be appended to the pair_style granular command to specify a global cutoff (i.e. a cutoff for all atom types). Additionally, the optional cutoff keyword can be passed to the pair_coeff command, followed by a cutoff value. This will set a pairwise cutoff for the atom types in the pair_coeff command. These options may be useful in some rare cases where the automatic cutoff determination is not sufficient, e.g. if particle diameters are being modified via the fix adapt command. In that case, the global cutoff specified as part of the pair_style granular command is applied to all atom types, unless it is overridden for a given atom type combination by the cutoff value specified in the pair coeff command. If cutoff is only specified in the pair coeff command and no global cutoff is appended to the pair_style granular command, then LAMMPS will use that cutoff for the specified atom type combination, and automatically set pairwise cutoffs for the remaining atom types.

---

Mixing, shift, table, tail correction, restart, rRESPA info

The pair_modify mix, shift, table, and tail options are not relevant for granular pair styles.

Mixing of coefficients is carried out using geometric averaging for most quantities, e.g. if friction coefficient for type 1-type 1 interactions is set to ${\mu }_1$, and friction coefficient for type 2-type 2 interactions is set to ${\mu }_2$, the friction coefficient for type1-type2 interactions is computed as $\sqrt{{\mu }_1{\mu }_2}$ (unless explicitly specified to a different value by a pair_coeff 1 2 … command). The exception to this is elastic modulus, only applicable to hertz/material, dmt and jkr normal contact models. In that case, the effective elastic modulus is computed as:

$$E_{eff,ij}={(\frac{1-{\nu }_i^2}{E_i}+\frac{1-{\nu }_j^2}{E_j})}^{-1}$$

If the i-j coefficients $E_{ij}$ and ${\nu }_{ij}$ are explicitly specified, the effective modulus is computed as:

$$E_{eff,ij}={(\frac{1-{\nu }_{ij}^2}{E_{ij}}+\frac{1-{\nu }_{ij}^2}{E_{ij}})}^{-1}$$

or

$$E_{eff,ij}=\frac{E_{ij}}{2(1-{\nu }_{ij}^2)}$$

These pair styles write their information to binary restart files, so a pair_style command does not need to be specified in an input script that reads a restart file.

These pair styles can only be used via the pair keyword of the run_style respa command. They do not support the inner, middle, outer keywords.

The single() function of these pair styles returns 0.0 for the energy of a pairwise interaction, since energy is not conserved in these dissipative potentials. It also returns only the normal component of the pairwise interaction force. However, the single() function also calculates at least 13 extra pairwise quantities. The first 3 are the components of the tangential force between particles I and J, acting on particle I. The fourth is the magnitude of this tangential force. The next 3 (5-7) are the components of the rolling torque acting on particle I. The next entry (8) is the magnitude of the rolling torque. The next entry (9) is the magnitude of the twisting torque acting about the vector connecting the two particle centers. The next 3 (10-12) are the components of the vector connecting the centers of the two particles (x_I - x_J). If a granular sub-model calculates additional contact information (e.g. the heat sub-models calculate the amount of heat exchanged), these quantities are appended to the end of this list. First, any extra values from the normal sub-model are appended followed by the damping, tangential, rolling, twisting, then heat models. See the descriptions of specific granular sub-models above for information on any extra quantities. If two or more models are defined by pair coefficients, the size of the array is set by the maximum number of extra quantities in a model but the order of quantities is determined by each model’s specific set of sub-models. Any unused quantities are zeroed.

These extra quantities can be accessed by the compute pair/local command, as p1, p2, …, p12.

---

Restrictions

This pair style is part of the GRANULAR package. It is only enabled if LAMMPS was built with that package. See the Build package page for more info.

This pair style requires that atoms store per-particle radius, torque, and angular velocity (omega) as defined by the atom_style sphere.

This pair style requires you to use the comm_modify vel yes command so that velocities are stored by ghost atoms.

This pair style will not restart exactly when using the read_restart command, though it should provide statistically similar results. This is because the forces it computes depend on atom velocities and the atom velocities have been propagated half a timestep between the force computation and when the restart is written, due to using Velocity Verlet time integration. See the read_restart command for more details.

Accumulated values for individual contacts are saved to restart files but are not saved to data files. Therefore, forces may differ significantly when a system is reloaded using the read_data command.

Related commands

pair_coeff pair gran/*

Default

For the pair_coeff settings: damping viscoelastic, rolling none, twisting none.

References

(Brilliantov et al, 1996) Brilliantov, N. V., Spahn, F., Hertzsch, J. M., & Poschel, T. (1996). Model for collisions in granular gases. Physical review E, 53(5), 5382.

(Tsuji et al, 1992) Tsuji, Y., Tanaka, T., & Ishida, T. (1992). Lagrangian numerical simulation of plug flow of cohesionless particles in a horizontal pipe. Powder technology, 71(3), 239-250.

(Johnson et al, 1971) Johnson, K. L., Kendall, K., & Roberts, A. D. (1971). Surface energy and the contact of elastic solids. Proc. R. Soc. Lond. A, 324(1558), 301-313.

(Derjaguin et al, 1975) Derjaguin, B. V., Muller, V. M., & Toporov, Y. P. (1975). Effect of contact deformations on the adhesion of particles. Journal of Colloid and interface science, 53(2), 314-326.

(Zunker and Kamrin, 2024) Zunker, W., & Kamrin, K. (2024). A mechanically-derived contact model for adhesive elastic-perfectly plastic particles, Part I: Utilizing the method of dimensionality reduction. Journal of the Mechanics and Physics of Solids, 183, 105492.

(Zunker and Kamrin, 2024) Zunker, W., & Kamrin, K. (2024). A mechanically-derived contact model for adhesive elastic-perfectly plastic particles, Part II: Contact under high compaction-modeling a bulk elastic response. Journal of the Mechanics and Physics of Solids, 183, 105493.

(Zunker et al, 2025) Zunker, W., Dunatunga, S., Thakur, S., Tang, P., & Kamrin, K. (2025). Experimentally validated DEM for large deformation powder compaction: mechanically-derived contact model and screening of non-physical contacts.

(Luding, 2008) Luding, S. (2008). Cohesive, frictional powders: contact models for tension. Granular matter, 10(4), 235.

(Marshall, 2009) Marshall, J. S. (2009). Discrete-element modeling of particulate aerosol flows. Journal of Computational Physics, 228(5), 1541-1561.

(Silbert, 2001) Silbert, L. E., Ertas, D., Grest, G. S., Halsey, T. C., Levine, D., & Plimpton, S. J. (2001). Granular flow down an inclined plane: Bagnold scaling and rheology. Physical Review E, 64(5), 051302.

(Kuhn and Bagi, 2005) Kuhn, M. R., & Bagi, K. (2004). Contact rolling and deformation in granular media. International journal of solids and structures, 41(21), 5793-5820.

(Wang et al, 2015) Wang, Y., Alonso-Marroquin, F., & Guo, W. W. (2015). Rolling and sliding in 3-D discrete element models. Particuology, 23, 49-55.

(Thornton, 1991) Thornton, C. (1991). Interparticle sliding in the presence of adhesion. J. Phys. D: Appl. Phys. 24 1942

(Mindlin, 1949) Mindlin, R. D. (1949). Compliance of elastic bodies in contact. J. Appl. Mech., ASME 16, 259-268.

(Thornton et al, 2013) Thornton, C., Cummins, S. J., & Cleary, P. W. (2013). An investigation of the comparative behavior of alternative contact force models during inelastic collisions. Powder Technology, 233, 30-46.

(Otis R. Walton) Walton, O.R., Personal Communication

(Mindlin and Deresiewicz, 1953) Mindlin, R.D., & Deresiewicz, H (1953). Elastic Spheres in Contact under Varying Oblique Force. J. Appl. Mech., ASME 20, 327-344.

(Agnolin and Roux 2007) Agnolin, I. & Roux, J-N. (2007). Internal states of model isotropic granular packings. I. Assembling process, geometry, and contact networks. Phys. Rev. E, 76, 061302.

(Vargas and McCarthy 2001) Vargas, W.L. and McCarthy, J.J. (2001). Heat conduction in granular materials. AIChE Journal, 47(5), 1052-1059.