Models for tangential friction in granular interactions
Tangential models apply a force \(\mathbf{F}_t = \mathbf{F}_{i,t} = -\mathbf{F}_{j,t}\) on each particle \(i,j\) in a direction tangential to the surfaces of the particles, and induce torques \(\tau_i,\tau_j\). The tangential force is added to the total normal force to produce the total force on each particle.
The general form of the tangential contact force for all models is:
where \(\mathbf{F}_{te}\) and \(\mathbf{F}_\mathrm{t,damp}\) are the elastic and damping components of the tangential force, respectively, and \(\mathbf{t}\) is the direction of the tangential force given by:
Critical force
The normal force value \(F_{n0}\) used to compute the critical force at which slip occurs 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:
For cohesive models such as jkr and dmt , the critical force is adjusted so that the critical tangential force approaches \(\mu_tF_\text{pulloff}\), see Marshall, equation 43, and Thornton. For these models, the normal force for purposes of computing the critical force is given by:
Where \(F_\text{pulloff} = 3\pi \gamma R\) for jkr, \(F_\text{pulloff} = 4\pi \gamma R\) for dmt.
For cohesive models such as epa/linear and epa/nonlinear, the critical force is incremented by the value of the cohesive force, see Luding2008, Morrisey2013 <Morrisey2013> .. math:
F_{n0} = \|\mathbf{F}_{ne} + F_\text{cohesive}\|
Where \(F_\text{cohesive} = f0 + k_c\delta\) for epa/linear, \(F_\text{cohesive} = f0 + k_c\delta^n\) for epa/nonlinear.
Tangential damping
The tangential damping force is the same for all models and is given by:
The tangential damping prefactor \(\eta_t\) is calculated by scaling the normal damping \(\eta_n\) (see above):
The normal damping prefactor \(\eta_n\) is determined by the choice of the damping keyword, see Models for normal damping in granular interactions. 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\boldsymbol{\Omega}_i + R_j\boldsymbol{\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 elastic tangential force \(\mathbf{F}_{te}\) depends on the tangential model and is detailed for each tangential model.
Tangential displacement
History-dependent tangential models (linear_history*, mindlin*) rely on accumulated tangential displacement (i.e. contact history), except for the options mindlin/force and mindlin_rescale/force, which use accumulated tangential force instead, and are discussed further below. The tangential displacement accumulated during the entire duration of the contact is denoted \(\mathbf{\xi}\), and is given by:
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:
Here, \(\mathbf{\xi'}\) is the accumulated displacement prior to the current time step and \(\mathbf{\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 a critical force \(F_{t,crit}\) (see critical force), the magnitude of the tangential force is truncated to the critical force, i.e., \(F_t = F_{t,crit}\). Additionally, the tangential displacement is re-scaled to match the value for the critical force (see Luding, equation 20 and related discussion):
linear_nohistory tangential model
Parameters: \(x_{\gamma,t}\), \(\mu_t\)
Example:
pair_style granular
pair_coeff * * hooke 1000.0 50.0 tangential linear_nohistory 1.0 0.4 damping mass_velocity
For the linear_nohistory tangential model, 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:
linear_history tangential model
Parameters: \(k_t\), \(x_{\gamma,t}\), \(\mu_t\)
Example:
pair_style granular
pair_coeff * * hooke 1000.0 50.0 tangential linear_history 500.0 1.0 0.4 damping mass_velocity
For tangential linear_history, the tangential force is given by:
mindlin tangential model
Parameters: \(k_t\) or NULL, \(x_{\gamma,t}\), \(\mu_t\)
Example:
pair_style granular
pair_coeff * * hertz/material 1.0e8 50.0 0.3 tangential mindlin NULL 1.0 0.4 damping viscoelastic
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:
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 jkr. To match the Mindlin solution, one should set \(k_t = 8G_{eff}\), where \(G_{eff}\) is the effective shear modulus given by:
where \(G\) is the shear modulus, related to Young’s modulus \(E\) and Poisson’s ratio \(\nu\) by \(G = E/(2(1+\nu))\). This can also be achieved by specifying NULL for \(k_t\), in which case a normal contact model that specifies material parameters \(E\) and \(\nu\) 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.
mindlin/force tangential model
Parameters: \(k_t\) or NULL, \(x_{\gamma,t}\), \(\mu_t\)
Example:
pair_style granular
pair_coeff * * hertz/material 1.0e8 50.0 0.3 tangential mindlin/force NULL 1.0 0.4 damping viscoelastic
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 for the mindlin model. The tangential force is given by:
The increment of the elastic component of the tangential force \(\mathbf{F}_{te}\) is given by:
The changes in frame of reference of the contacting pair of particles during contact are accounted for by the same formula as tangential displacement, 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:
The same rules as those described for mindlin apply regarding the tangential stiffness and mixing of the shear modulus for different particle types.
mindlin/rescale tangential model
Parameters: \(k_t\) or NULL, \(x_{\gamma,t}\), \(\mu_t\)
Example:
pair_style granular
pair_coeff * * hertz/material 1.0e8 50.0 0.3 tangential mindlin_rescale NULL 1.0 0.4 damping viscoelastic
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}}\):
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.
mindlin_rescale/force tangential model
Parameters: \(k_t\) or NULL, \(x_{\gamma,t}\), \(\mu_t\)
Example:
pair_style granular
pair_coeff * * hertz/material 1.0e8 50.0 0.3 tangential mindlin_rescale/force NULL 1.0 0.4 damping viscoelastic
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}}\):
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.
linear_history/static tangential model
Parameters: \(k_t\), \(x_{\gamma,t}\), \(\mu_{t,\mathrm{static}}\), \(\mu_{t,\mathrm{dynamic}}\)
Example:
pair_style granular
pair_coeff * * hooke 1000.0 50.0 tangential linear_history/static 500.0 1.0 0.5 0.3 damping mass_velocity
The linear_history/static option mirrors the linear_history model, but allows for different static and dynamic friction coefficients, \(\mu_{t,static}\) and \(\mu_{t,dynamic}\).
The static and dynamic critical friction forces are defined as:
and
The state of a contact is designated as static or dynamic, and tracked throughout a simulation. The initial state of a contact is initialized to static at the start of a simulation (or read from a restart file, similar to other history entries).
While the contact is designated as static, the critical force for purposes of determining slip and truncating the tangential displacement is given by \(F_\text{t,crit} = F_\text{t,crit,static}\). If the tangential force exceeds the static critical force while the contact is static, i.e. \(F_t > F_\text{t,crit,static}\), the contact state is changed to dynamic.
Conversely, while the contact is designated as dynamic, the critical force is given by \(F_\text{t,crit} = F_\text{t,crit,dynamic}\). If the tangential force drops below the dynamic critical force while the contact is dynamic, i.e. \(F_t < F_\text{t,crit,dynamic}\), the contact state is changed back to static.
In most physical scenarios, \(\mu_\text{t,dynamic}\le\mu_\text{t,static}\). If the user specifies \(\mu_\text{t,dynamic}\ge\mu_\text{t,static}\), the static friction coefficient will have no effect, and the model recovers the behavior of a linear_history model with \(\mu = \mu_\text{t,dynamic}\). A warning is issued to this effect.
mindlin/static tangential model
Parameters: \(k_t\) or NULL, \(x_{\gamma,t}\),
\(\mu_{t,\text{static}}\), \(\mu_{t,\text{dynamic}}\)
Example:
pair_style granular
pair_coeff * * hertz/material 1.0e8 50.0 0.3 tangential mindlin/static NULL 1.0 0.5 0.3 damping viscoelastic
The mindlin/static option mirrors the mindlin model, but allows different static and dynamic friction coefficients, \(\mu_\text{t,static}\) and \(\mu_\text{t,dynamic}\). The determination of the critical force and designation of a contact as static or dynamic are the same as described for linear_history/static.
References
(Otis R. Walton) Walton, O.R., Personal Communication
(Mindlin, 1949) Mindlin, R. D. (1949). Compliance of elastic bodies in contact. J. Appl. Mech., ASME 16, 259-268.
(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.
(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.