Rolling Resistance Linear Model

The rolling resistance linear model is based on the linear model, to which a rolling-resistance mechanism is added. It is a linear-based model that can be installed at both ball-ball and ball-facet contacts, and is referred to in commands and FISH by the name rrlinear.

Introduction

The effect of rolling resistance at contacts between particles, and associated energy dissipation, may be of major importance for many granular applications in both dense, quasistatic and dynamic regimes. In real granular systems, these mechanisms may have different micro-mechanical origins, such as adhesion of the contact area, or the steric effect due to surface roughness or non-sphericity about the contact point. The rolling resistance contact model provided in PFC is a simple model, based on the linear model, that incorporates a torque acting on the contacting pieces to counteract rolling motion. This model is based on the review paper [Ai2011] and on the work presented in [Wensrich2012]. The user is encouraged to consult these papers and references therein for a thorough introduction to the importance of rolling-resistance mechanisms in granular systems along with their implementation in DEM models.

Behavior Summary

The behavior of the rolling resistance linear contact model is similar to the linear model, except that the internal moment is incremented linearly with the accumulated relative rotation of the contacting pieces at the contact point. This accumulation occurs up to a maximum limiting value equal to the product of the current normal force with a rolling friction coefficient, \(\mu_r\), and an effective contact radius \(\bar{R}\). The constant that relates the increment of internal moment to the increment of relative rotation at the contact point is denoted as the rolling stiffness, \(k_r\). This value is related to the shear stiffness \(k_s\) and to the effective radius of the contact following [Ishiwata1998]. As discussed in [Wensrich2012], this choice sets the nominal rotational natural frequency due to the rolling stiffness equal to the nominal rotational natural frequency due to the shear stiffness, leading to a well-behaved and well-damped rolling-resistance mechanism without the need for any additional damping parameters. The rolling resistance linear model is therefore similar to the one described in [Wensrich2012], except that it is based on a linear model instead of a Hertz model. Note that this model only resists relative bending at the contact point. More complex models, that include resistance to bending and twisting have been proposed (see for instance [Jiang2015]).

Activity-Deletion Criteria

Similar to the linear model, a contact with the rolling resistance linear model is active if and only if the surface gap is less than or equal to zero. The force-displacement law is skipped for inactive contacts. The surface gap is shown in this figure of the linear formulation. When the reference gap is zero, the notional surfaces coincide with the piece surfaces.

Force-Displacement Law

The force-displacement law for the rolling resistance linear model updates the contact force and moment as:

(1)\[\mathbf{F_{c}} =\mathbf{F^{l}} +\mathbf{F^{d}} ,\quad \mathbf{M_{c}} = \mathbf{M^r}\]

where \(\mathbf{F^{l}}\) is the linear force, \(\mathbf{F^{d}}\) is the dashpot force, and \(\mathbf{M^r}\) is the rolling resistance moment. The linear and dashpot forces are updated as in the linear model, while the rolling resistance moment is updated with the following steps.

First, the rolling resistance moment is incremented as:

(2)\[\mathbf{M^r} :=\mathbf{M^r} - k_r {\kern 1pt} \Delta \pmb{θ}_\mathbf{b}\]

where \(\Delta \pmb{θ}_\mathbf{b}\) is the relative bend-rotation increment of this equation of the “Contact Resolution” section.

Note that the rolling resistance moment has no twisting component. It can be expressed in the contact plane coordinate system as:

(3)\[\mathbf{M^{r}} =M_{bs}^{r} {\kern 1pt} \hat{\mathbf{s}}_\mathbf{c} +M_{bt}^{r} {\kern 1pt} \hat{\mathbf{t}}_\mathbf{c} \; \; \left({\rm 2D\; model:\; }M_{bt}^{r} \equiv 0\right).\]

The rolling resistance stiffness \(k_r\) is defined as:

(4)\[k_r = k_s {\kern 1pt} \bar{R}^{2}\]

with \(\bar{R}\) the contact effective radius, defined as:

(5)\[\frac{1}{\bar{R}} = \frac{1}{R^{(1)}} + \frac{1}{R^{(2)}}.\]

In the above equation \(R^{(1)}\) and \(R^{(2)}\) are the radii of end (1) and end (2) of the contact respectively (with \(R^{(2)} = \infty\) for ball-facet contacts).

The magnitude of the updated rolling resistance moment is then checked against a threshold limit:

(6)\[\begin{split}\mathbf{M^r} = \left\{ \begin{array}{rl} \mathbf{M^r} , & \| \mathbf{M^r} \| \le M^* \\ M^* \bigl( \mathbf{M^r} / \| \mathbf{M^r} \| \bigr), & \mbox{otherwise.} \end{array} \right.\end{split}\]

The limiting torque is defined as:

(7)\[ M^{*} = \mu_r \bar{R} F_n^l.\]

The rolling slip state is also updated as:

(8)\[\begin{split}s_r = \left\{ \begin{array}{rl} \mbox{true}, & \| \mathbf{M^r} \| = M^* \\ \mbox{false}, & \mbox{otherwise.} \end{array} \right.\end{split}\]

As discussed in [Ai2011], the rolling friction coefficient \(\mu_r\) corresponds to the tangent of the maximum angle of a slope on which the rolling resistance torque counterbalances the torque produced by gravity acting on the body (as illustrated in the verification example Rolling Resistance Linear Contact Model: Single Ball on a Flat Surface).

Energy Partitions

In addition to the energy partitions of the linear model, the rolling resistance linear model provides the following two energy partitions:

  • rolling strain energy, \(E_{k_r}\), stored in the linear spring; and
  • rolling slip energy, \(E_{\mu_r}\), defined as the total energy dissipated by rolling slip.

If energy tracking is activated, these energy partitions are updated as follows.

  1. Update the rolling strain energy:

    (9)\[E_{k_r} = \frac{1}{2} \frac{\left\| \mathbf{M^{r}} \right\| ^{2} }{k_{r}}.\]

  2. Update the rolling slip energy:

    (10)\[\begin{split}\begin{array}{l} {E_{\mu_r} := E_{\mu_r} - {\tfrac{1}{2}} \biggl(\left(\mathbf{M^{r}} \right)_{o} +\mathbf{M^{r}} \biggr)\cdot \Delta \pmb{θ}_\mathbf{b}^{\mu_r} } \\ {{\rm with}\qquad \Delta \pmb{θ}_\mathbf{b}^{\mu_r} = \Delta \pmb{θ}_\mathbf{b} -\Delta \pmb{θ}_\mathbf{b}^{k} =\Delta \pmb{θ}_\mathbf{b} -\left(\frac{\mathbf{M^{r}} -\left(\mathbf{M^{r}} \right)_{o} }{k_{r} } \right)} \end{array}\end{split}\]

    where \(\left(\mathbf{M^{r}} \right)_{o}\) is the rolling resistance moment at the beginning of the timestep, and the adjusted relative bend-rotation increment of Equation (2) has been decomposed into an elastic \(\Delta \pmb{θ}_\mathbf{b}^{k}\) and a slip \(\Delta \pmb{θ}_\mathbf{b}^{\mu_r}\) component. The dot product of the slip component and the rolling resistance moment occurring during the timestep gives the increment of rolling slip energy.

Additional information, including the keywords by which these partitions are referred to in commands and FISH, is provided in the table below.

Table 1: Rolling Resistance Linear Model Energy Partitions
Keyword Symbol Description Range Accumulated
Linear Group:
energy‑strain \(E_{k}\) strain energy \([0.0,+\infty)\) NO
energy‑slip \(E_{\mu}\) total energy dissipated by slip \((-\infty,0.0]\) YES
Dashpot Group:
energy‑dashpot \(E_{\beta}\) total energy dissipated by dashpots \((-\infty,0.0]\) YES
Rolling-Resistance Group:
energy‑rrstrain \(E_{k_r}\) Rolling strain energy \([0.0,+\infty)\) NO
energy‑rrslip \(E_{\mu_r}\) Total energy dissipated by rolling slip \((-\infty,0.0]\) YES

Properties

The properties defined by the rolling resistance linear contact model are listed in the table below as a concise reference; see the Contact Properties section for a description of the information in the table columns. The mapping from the surface inheritable properties to the contact model properties is also discussed below.

Table 2: Rolling Resistance Linear Model Properties
Keyword Symbol Description Type Range Default Modifiable Inheritable
rrlinear Model name
Linear Group:
kn \(k_n\) Normal stiffness [force/length] FLT \([0.0,+\infty)\) 0.0 YES YES
ks \(k_s\) Shear stiffness [force/length] FLT \([0.0,+\infty)\) 0.0 YES YES
fric \(\mu\) Friction coefficient [-] FLT \([0.0,+\infty)\) 0.0 YES YES
rgap \(g_r\) Reference gap [length] FLT \(\mathbb{R}\) 0.0 YES NO
lin_mode \(M_l\) Normal-force update mode [-] INT \(\{0,1\}\) 0 YES NO
    \(\;\;\;\;\;\;\begin{cases} \mbox{0: update is absolute} \\ \mbox{1: update is incremental} \end{cases}\)          
emod \(E^*\) Effective modulus [force/area] FLT \([0.0,+\infty)\) 0.0 NO N/A
kratio \(\kappa^*\) Normal-to-shear stiffness ratio [-] FLT \([0.0,+\infty)^*\) 0.0 * NO N/A
    \(\kappa^* \equiv \frac{k_n}{k_s}\)          
lin_slip \(s\) Slip state [-] BOOL {false,true} false NO N/A
    \(\;\;\;\;\;\;\begin{cases} \mbox{true: slipping} \\ \mbox{false: not slipping} \end{cases}\)          
lin_force \(\mathbf{F^l}\) Linear force (contact plane coord. system) VEC \(\mathbb{R}^3\) \(\mathbf{0}\) YES NO
    \(\left( -F_n^l,F_{ss}^l,F_{st}^l \right) \quad \left(\mbox{2D model: } F_{ss}^l \equiv 0 \right)\)          
user_area \(A\) Constant area [length*length] FLT \((0.0,+\infty)\) 0.0 YES NO
Dashpot Group:
dp_nratio \(\beta_n\) Normal critical damping ratio [-] FLT \([0.0,1.0]\) 0.0 YES NO
dp_sratio \(\beta_s\) Shear critical damping ratio [-] FLT \([0.0,1.0]\) 0.0 YES NO
dp_mode \(M_d\) Dashpot mode [-] INT {0,1,2,3} 0 YES NO
    \(\;\;\;\;\;\;\begin{cases} \mbox{0: full normal & full shear} \\ \mbox{1: no-tension normal & full shear} \\ \mbox{2: full normal & slip-cut shear} \\ \mbox{3: no-tension normal & slip-cut shear} \end{cases}\)          
dp_force \(\mathbf{F^d}\) Dashpot force (contact plane coord. system) VEC \(\mathbb{R}^3\) \(\mathbf{0}\) NO NO
    \(\left( -F_n^d,F_{ss}^d,F_{st}^d \right) \quad \left(\mbox{2D model: } F_{ss}^d \equiv 0 \right)\)          
Rolling-Resistance Group:
rr_kr \(k_r\) Rolling resistance stiffness [force*length] FLT \([0.0,+\infty)\) 0.0 NO NO
rr_fric \(\mu_r\) Rolling friction coefficient [-] FLT \([0.0,+\infty)\) 0.0 YES YES
rr_slip \(s_r\) Rolling slip state [-] BOOL {false,true} false NO N/A
rr_moment \(\mathbf{M^r}\) Rolling resistance moment (contact plane coord. system) VEC \(\mathbb{R}^3\) \(\mathbf{0}\) YES NO
    \(\left( 0,M_{bs}^r,M_{bt}^r \right) \quad \left(\mbox{2D model: } M_{bt}^r \equiv 0 \right)\)          
* By convention, \(\kappa^* = 0.0\) if either the normal or the shear stiffness is 0.

Note

Out-of-balance forces acting on bodies are accumulated during force-displacement calculations. Modifying the forces stored in the contact models will not alter them instantaneously. Therefore, any change to \(\mathbf{F^l}\) may only be effective during the next force-displacement calculation. When \(M_l = 0\), the normal component of the linear force is automatically overridden during the next force-displacement calculation.

Surface Property Inheritance

The linear stiffnesses, \(k_n\) and \(k_s\), and the friction coefficient, \(\mu\), may be inherited from the contacting pieces. See this section from the linear formulation for details. Additionally, the rolling resistance stiffness \(k_r\) is updated if the shear stiffness \(k_s\) is modified.

The rolling friction coefficient is inherited using the minimum of the values set for the pieces:

(11)\[\mu_r = \min \left( \mu_r^{(1)},\mu_r^{(2)} \right).\]

Methods

The methods of the rolling resistance linear model are identical to those of the linear model.

Table 3: Rolling Resistance Linear Model Methods
Method Arguments Symbol Type Range Default Description
area Set user_area to the area
deformability Set deformability
  emod \(E^*\) FLT \([0.0,+\infty)\) N/A Effective modulus
  kratio \(\kappa^*\) FLT \([0.0,+\infty)^*\) N/A Normal-to-shear stiffness ratio
* By convention, setting \(\kappa^* = 0.0\) sets \(k_s = 0.0\) but does not alter \(k_n\).

Area

Set the user_area property via the current contact area. This operation means that the contact area stays constant and is fixed independent of changes to the piece sizes/geometries. In order for the stiffnesses to be recomputed accounting for this area, one should subsequently call the deformabilty method.

Deformability

See this section from the linear formulation for details on these methods.

Callback Events

The callback events defined by the rolling resistance linear model are identical to those defined by the linear model.

Table 4: Rolling Resistance Linear Model Callback Events
Event Array Slot Value Type Range Description
contact_activated Contact has become active
  1 C_PNT N/A Contact pointer
slip_change Slip state has changed
  1 C_PNT N/A Contact pointer
  2 INT {0,1} Slip change mode
        \(\;\;\;\;\;\;\begin{cases} \mbox{0: slip has initiated} \\ \mbox{1: slip has ended} \end{cases}\)

Usage and Verification Examples

The verification example Rolling Resistance Linear Contact Model: Single Ball on a Flat Surface exercises the model for a system consisting of one ball sliding and rolling on a flat surface. The example Rolling Resistance Linear Contact Model: Repose Angle studies the influence of the friction and rolling friction coefficients on the repose angle of a granular heap.

Model Summary

An alphabetical list of the rolling resistance linear model methods is given here. An alphabetical list of the rolling resistance linear model properties is given here.

References

[Ai2011](1, 2) Ai, J., J. F. Chen, J. M. Rotter, and J. Y. Ooi, “Assessment of Rolling Resistance Models in Discrete Element Simulations,” Powder Technology, 206, 269-282, 2011.
[Wensrich2012](1, 2, 3) Wensrich, C. M., and A. Katterfeld, “Rolling Friction as a Technique for Modelling Particle Shape in DEM,” Powder Technology, 217, 409-417, 2012.
[Ishiwata1998]Iwashita, K., and M. Oda, “Rolling Resistance at Contacts in Simulation of Shear Band Development by DEM,” Journal of Engineering Mechanics ASCE, 124, 285-292, 1998.
[Jiang2015]Jiang, M., Z. Shen, and J. Wang, “A Novel Three-Dimensional Contact Model for Granulates Incorporating Rolling and Twisting Resistances,” Computers and Geotechnics, 65, 147-163, 2015.