<figure><img src="image_1.png" /></figure>

HAL
open science

# Coulomb friction with rolling resistance as a cone complementarity problem

Vincent Acary, Franck Bourrier

► To cite this version:

Vincent Acary, Franck Bourrier. Coulomb friction with rolling resistance as a cone complementarity problem. European Journal of Mechanics - A/Solids, 2021, 85, pp.1-12. 10.1016/j.euromechsol.2020.104046. hal-02382213v2

HAL Id: hal-02382213

https://inria.hal.science/hal-02382213v2

Submitted on 11 Jun 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Coulomb friction with rolling resistance as a cone complementarity problem

Vincent Acary<sup>a</sup>, Franck Bourrier<sup>b</sup>

<sup>a</sup>Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP*, LJK, 38000 Grenoble, France

*Institute of Engineering Univ. Grenoble Alpes

<sup>b</sup>Univ. Grenoble Alpes, INRAE, ETNA, 38000 Grenoble, France

## Abstract

Coulomb friction model with unilateral contact is a basic, but reliable, model to represent the resistance to sliding between solid bodies. It is nowadays well-known that this model can be formulated as a second-order cone complementarity problem, or equivalently, as a variational inequality. In this article, the Coulomb friction model is enriched to take into account the resistance to rolling, also known as rolling friction. Introducing the rolling friction cone, an extended Coulomb's cone, and its dual, a formulation of the Coulomb friction with rolling resistance as a cone complementarity problem is shown to be equivalent to the standard formulation of the Coulomb friction with rolling resistance. Based on this complementarity formulation, the maximum dissipation principle and the bi-potential function are derived. Several iterative numerical methods based on projected fixed point iterations for variational inequalities and block-splitting techniques are given. The efficiency of these methods strongly relies on the computation of the projection onto the rolling friction cone. In this article, an original closed-form formulae for the projection on the rolling friction cone is derived. The abilities of the model and the numerical methods are illustrated on the examples of a single sphere sliding and rolling on a plane, and of the evolution of spheres piles under gravity.

## 1. Introduction and rolling friction model

Rolling friction has a long history; see for instance [14] and references therein. Like the standard Coulomb friction, the threshold effect of rolling friction is the main physical phenomenon that is required to be modeled for a basic friction law as described in [34, Chapter 1, §20]. Although it may be considered as a naive model from a tribological point of view, the simplest rolling friction model postulates the existence of an upper bound of the norm of the rolling resistive moment $m_R$ that is proportional to the normal pressure $r_N$:

$$ \|m_R\| \le \mu_r r_N \qquad (1) $$

<footnote>Email addresses: vincent.acary@inria.fr (Vincent Acary), franck.bourrier@inriae.fr (Franck Bourrier)</footnote>
where $\mu_r$ is the rolling friction coefficient. Such a simple model has been used successfully in analytical nonlinear dynamics [30, 29] and in [37, 23, 11, 10, 17] for numerical computations of granular materials. The standard Coulomb friction law is a set-valued force law that generates a resistive force to sliding, i.e, opposite to the sliding velocity. The standard rolling friction law, considered in this article, is also a set-valued force law that generates a resistive moment to rolling, i.e, opposite to the rolling velocity denoted by $\omega_R$. The same type of law can also be formulated for the pivoting friction, which will not be treated in this article although the proposed approach could also be applied.

Our goal is to formulate the rolling friction model as a set-valued force law in a nonsmooth setting following the seminal work of Moreau [31]. For a good account on set-valued force laws approach and nonsmooth mechanics, we refer to [33, 13, 7]. This approach enables to compute solutions with efficient numerical methods that come mainly from mathematical programming techniques, especially complementarity problems, normal cone inclusions, and variational inequalities [12]. Our paper is motivated by the numerical evaluation of the classical rolling friction in the isotropic case rather than a novel friction law based on microscopic tribological considerations as it can be found in [20, 21, 25, 24].

As we mentioned above, the standard rolling friction model is the simplest macroscopic model of resistance to rolling. For a discussion on more complex models of rolling and pivoting friction, we refer to [28, 26, 27, 17]. In [28, 26, 27], the Coulomb-Contensou model for pivoting resistance is formulated in terms of normal cones inclusions. This model considers the coupling between the sliding friction and the pivoting friction. A more accurate model of rolling friction should also consider a coupling between sliding and rolling friction based on tribological considerations. In this article, we consider a simple decoupled isotropic rolling friction model. The authors in [26, 27] introduced also the concept of contour friction as a special type of resistance to rolling. The main difference lies in the definition of the rolling friction $\omega_R$ for an inclined disk rolling on a plane. Indeed, the resistance to rolling requires the definition of the rolling velocity $\omega_R$ between two contacting bodies. In general configurations, this is not a trivial question to define $\omega_R$ from the kinematics of the bodies. Some definitions may yield unrealistic behaviors and the study of the elastic rolling contact problem [21] should be considered. In this article, the development of the formulation as cone complementarity problem and the associated numerical technique are not dependent of the definition of the rolling velocity, so that the contour friction can also be taken into account. In [16, 36], the standard decoupled rolling friction is addressed with the help of second-order cone complementarity problem, but the friction flow rule is associated. In other words, the relative velocities at contact are chosen such that they are orthogonal to the friction cone yielding a dilatance in case of sliding and no dissipation. The problem becomes a convex relaxation of the initial model, but does not represent standard Coulomb friction without dilatancy.
The article is organized as follows. After the presentation of the model in Section 2, a formulation of the Coulomb friction with unilateral contact and rolling resistance with no-dilatancy as a cone complementarity problem is given in Section 3. Following the work on the bi-potential method [9, 35], the modified velocity at contact is introduced to avoid dilatancy. The formulation is proved to be equivalent to disjuntive (or modal) description of the rolling friction. In Section 4, optimization formulations are exhibited that lead to the principle of maximum dissipation for a given normal reaction and to the definition of a bipotential for the rolling friction. Numerical schemes based on projected fixed iterations and block-splitting techniques are presented in Section 5 to solve the problem with multiple contact points. Especially, the explicit formulae for the projection onto the rolling friction cone, that is crucial for the efficiency of the numerical methods, is derived in Section 5.3. The article is concluded by two simple numerical applications that emphasize the capacities of the rolling friction model : the motion of a sphere rolling on a plane and the creation of a pile of spheres using a pluviation process.

**Notation:** For a set $X \subset \mathbb{R}^n$, $\bar{X}, \hat{X}$ and $\partial X$ will respectively denote the closure, the interior and the boundary of $X$. For a convex set $C \subset \mathbb{R}^n$, the normal cone to $C$ at $x \in \mathbb{R}^n$ is defined by $N_C(x) := \{y | y^\top(z-x) \le 0, \forall z \in C\}$. The indicator function $i_C$ is defined by $i_C(x) = 0$, if $x \in C$ and $i_C(x) = +\infty$, if $x \notin C$. For a cone $K$, its dual cone $K^*$ is defined by $K^* := \{y | y^\top p \ge 0, \forall p \in K\}$ and its polar by $K^o = -K^*$. For a convex function $f : \mathbb{R}^n \to \mathbb{R}$, the subdifferential of $f$ is defined by $\partial f = \{\gamma \in \mathbb{R}^n | f(y) - f(x) \ge \gamma^T(y - x) \text{ for all } y \in \mathbb{R}^n\}$.

## 2. Rolling friction model

### 2.1. Standard Coulomb friction model with unilateral contact.

Coulomb friction model with unilateral contact is a basic, but reliable model to represent the mechanical behavior of interfaces between solid bodies. Let us recall its formulation. In Figure 1, we represent two bodies, labeled by $A$ and $B$, potentially in contact at points $C_A$ and $C_B$ in a three-dimensional configuration space<fnref n="1" />. Let assume that we can uniquely define an orthonormal contact frame at $C_A$ denoted by $(C_A, \mathbf{N}, \mathbf{T}_1, \mathbf{T}_2)$, where $\mathbf{N} \in \mathbb{R}^3$ defines an outward unit normal vector to $A$ at point $C_A$ and $\mathbf{T}_1 \in \mathbb{R}^3, \mathbf{T}_2 \in \mathbb{R}^3$ are unit tangent vectors. The gap between $C_A$ and $C_B$ reads as

$$g_N := (C_B - C_A)^\top \mathbf{N}. \qquad (2)$$

The reaction force exerted by $A$ on $B$ is denoted by $r \in \mathbb{R}^3$. It can be decomposed in the contact frame as

$$r := r_N \mathbf{N} + r_{T_1} \mathbf{T}_1 + r_{T_2} \mathbf{T}_2, \quad \text{with } r_N \in \mathbb{R} \text{ and } r_T := [r_{T_1}, r_{T_2}]^\top \in \mathbb{R}^2. \qquad (3)$$

<footnote n="1">¹The two-dimensional case can be treated as well by reducing the dimension of the problem</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 1: Local frame at contact.</figcaption></figure>

The Signorini condition expresses the unilateral contact as a complementarity condition between the gap and the normal reaction force as

$$0 \le g_N \perp r_N \ge 0. \qquad (4)$$

Since we want to deal with friction, the relative velocity at contact $u \in \mathbb{R}^3$ is used as natural way to formulate friction. It is also decomposed as

$$u := u_N N + u_{T_1} T_1 + u_{T_2} T_2 \quad \text{with } u_N \in \mathbb{R} \text{ and } u_T = [u_{T_1}, u_{T_2}]^\top \in \mathbb{R}^2. \qquad (5)$$

At the velocity level, the Signorini condition is written

$$\left\{ \begin{array}{ll} 0 \le u_N \perp r_N \ge 0 & \text{if } g_N \le 0 \\ r_N = 0 & \text{otherwise.} \end{array} \right. \qquad (6)$$

The Moreau's viability Lemma [32] ensures that (6) implies (4) if $g_N \ge 0$ holds in the initial configuration. In the sequel, we will focus our attention on closed contact ($g_N \le 0$) since the other case is trivial. Equivalently, the complementarity condition can be written as a normal cone inclusion, or a variational inequality:

$$-u_N \in N_{\mathbb{R}_+}(r_N) \iff u_N(s - r_N) \ge 0, \forall s \in \mathbb{R}_+. \qquad (7)$$

If needed, an impact law must be added to complete the interface model. This is particularly the case when we deal with rigid solid bodies. The Newton impact law might be chosen as a basic impact model. Considering the post-impact velocity $u_N^+$, and the pre-impact velocity $u_N^-$, the Newton impact law can be written as

$$u_N^+ = -e_N u_N^- \qquad (8)$$
where $e_n \ge 0$ is the coefficient of restitution. As in [32], the Newton impact law can be included in the complementary condition at the velocity level as

$$0 \le \bar{u}_N \perp r_N \ge 0, \qquad (9)$$

with $\bar{u}_N := u_N^+ + e_n u_N^-$. In order to simplify the presentation, we will consider in the sequel that $e_n = 0$ and $\bar{u}_N = u_N^+ = u_N^-$.

In the tangent plane at contact, Coulomb friction is based on the postulate of an admissible set for the reaction force, called the Coulomb's cone:

$$K = \{r \in \mathbb{R}^3 \mid \|r_T\| \le \mu r_N\}, \qquad (10)$$

where $\mu$ is the coefficient of friction.

In the sticking case, Coulomb's law states

$$u_T = 0, \quad r \in K, \qquad (11)$$

and for the sliding case

$$u_T \neq 0, \quad \|r_T\| = \mu r_N \quad \text{and} \quad \|r_T\| u_T = -\|u_T\| r_T. \qquad (12)$$

Introducing the disk $D(c) := \{x \in \mathbb{R}^2 \mid \|x\| \le c\}$ of radius $c \ge 0$, the Coulomb model can be written as an inclusion

$$-u_T \in N_{D(\mu r_N)}(r_T). \qquad (13)$$

Using definition of the normal cone, an equivalent variational inequality can be written

$$u_T^\top(z - r_T) \ge 0, \quad \forall z \in D(\mu r_N). \qquad (14)$$

The latter formulation is often related to Moreau's maximum dissipation principle of the frictional behavior:

$$r_T \in \underset{\|z\| \le \mu r_N}{\operatorname{argmax}} -z^\top u_T. \qquad (15)$$

which is completely equivalent to the inclusion (13). This principle postulates that, for a given relative tangent velocity $u_T$, the reaction $r_T$ maximizes the dissipation $-r_T u_T$ over all the admissible tangent reactions. It is noteworthy to remark that this principle does not claim that $r$ is chosen such that is maximizes the dissipation $r^\top u$. This is related to the non-associated character of Coulomb friction.
## 2.2. Rolling friction model.

In order to formulate the rolling friction at contact, we introduce the relative angular velocity $\omega_R \in \mathbb{R}^2$ and the rolling friction reaction moment $m_R \in \mathbb{R}^2$ at contact. Following [34], a rolling friction model can be defined with two following modes

$$ \left\{ \begin{array}{ll} \omega_R = 0, \|m_R\| \le \mu_r r_N & \text{(no-rolling)} \\ \omega_R \ne 0, \|m_R\| = \mu_r r_N, \text{ and } \|m_R\|\omega_R = -\|\omega_R\|m_R & \text{(rolling)}. \end{array} \right. \qquad (16) $$

where $\mu_r$ is the rolling friction coefficient which has the dimension of a length.

As for the Coulomb friction, the model (16) can be expressed as a normal cone inclusion into the disk of rolling friction as

$$ -\omega_R \in N_D(\mu_r r_N)(m_T). \qquad (17) $$

Again, a maximum dissipation principle can be written as

$$ m_R \in \underset{\|z\| \le \mu_r r_N}{\operatorname{argmax}} -z^\top \omega_R. \qquad (18) $$

## 2.3. Compact formulation

In order to obtain a compact formulation of the model, we denote the local variables at contact by

$$ p := \begin{bmatrix} r_N \\ r_T \\ m_R \end{bmatrix} = \begin{bmatrix} r \\ m_R \end{bmatrix} \quad \text{and} \quad y := \begin{bmatrix} u_N \\ u_T \\ \omega_R \end{bmatrix} = \begin{bmatrix} u \\ \omega_R \end{bmatrix}. \qquad (19) $$

The rolling friction cone, as the cone of admissible reaction forces and torques, is denoted by

$$ K_r = \{p \in \mathbb{R}^5 \mid \|r_T\| \le \mu r_N, \|m_R\| \le \mu_r r_N\} \subset \mathbb{R}^5. \qquad (20) $$

which is a closed convex cone. The Coulomb friction model with unilateral contact and rolling resistance is formulated in a disjunctive form for $g_N \le 0$ as

$$ \left\{ \begin{array}{ll} p = 0, u_N \ge 0, & \text{(take-off)} \\ p \in K_r, u = 0, \omega_R = 0, & \text{(sticking and no-rolling)} \\ p \in \partial K_r, u_N = 0, \|r_T\| = \mu r_N, \|m_R\| < \mu_r r_N, & \|r_T\|u_T = -\|u_T\|r_T, \omega_R = 0, \text{(sliding and no-rolling)} \\ p \in \partial K_r, u_N = 0, \|r_T\| < \mu r_N, \|m_R\| = \mu_r r_N, & \|m_R\|\omega_R = -\|\omega_R\|m_R, u_T = 0, \text{(sticking and rolling)} \\ p \in \partial K_r, u_N = 0, \|r_T\| = \mu r_N, \|m_R\| = \mu_r r_N, & \|r_T\|u_T = -\|u_T\|r_T, \\ \|m_R\|\omega_R = -\|\omega_R\|m_R. & \text{(sliding and rolling)} \end{array} \right. \qquad (21) $$

In the disjunctive form (21), the model is hardly tractable from the computational point of view. In the sequel, we propose a formulation of the model as a single normal cone inclusion to the rolling friction cone $K_r$, and we prove its equivalence with the disjunctive formulation of the model.
### 3. Formulation as a cone complementarity problem

One of the objective of the article is to formulate the rolling friction as a unique inclusion to the normal cone to $K_r$. Using the normal cone inclusions (7), (13) and (17), the rolling friction model is given by the following set of inclusions,

$$ \begin{cases} -u_N \in N_{\mathbb{R}_+}(r_N) \\ -u_T \in N_{D(\mu_r r_N)}(r_T) \\ -\omega_R \in N_{D(\mu_r r_N)}(m_T). \end{cases} \quad (22) $$

that appears to be a quasi-variational inequality. Unfortunately, it is not possible to deduce from (22) a complementarity problem, since the disks $D(\mu_r r_N)$ and $D(\mu_r r_N)$ are not cones. A single normal cone inclusion to $K_r$ enables to write a complementarity problem, and then, to take benefit from the theory and numerical methods for complementarity problems [12].

In the case of standard Coulomb friction and following the seminal work of De Saxcé and Feng [35], both normal cone inclusions (7) and (13) can be joined into one normal cone inclusion as

$$ -\hat{u} \in N_K(r), \quad (23) $$

where the modified relative velocity is defined by

$$ \hat{u} := u + [\mu \|u_T\|, 0, 0]^\top. \quad (24) $$

The modified relative velocity takes into account the non-associated character of the friction rule with no dilatancy. Since $K$ is the second-order cone, the inclusion can be expressed as a second-order complementarity condition as

$$ K^* \ni \hat{u} \perp r \in K \quad (25) $$

where $K^*$ is the dual cone of $K$. For more details on this reformulation, we refer to [5].

As we remark in (25), the first step is to characterize the dual cone of $K_r$.

**Lemma 3.1** ($K_r^*$ dual cone of $K_r$). The dual cone of $K_r$ is given by

$$ K_r^* = \{y = \begin{bmatrix} u_N \\ u_T \\ \omega_R \end{bmatrix} \mid \mu \|u_T\| + \mu_r \|\omega_R\| \le u_N\}. \quad (26) $$

**Proof**: Let us introduce the following convex cone

$$ L = \{y = [u_N, u_T^\top, \omega_R^\top]^\top \mid \mu \|u_T\| + \mu_r \|\omega_R\| \le u_N\}. \quad (27) $$
The goal is to prove that $K_r^* = L$. Let us choose $y \in L$ and $p \in K_r$. We obtain

$$y^\top p = u_N r_N + u_T^\top r_T + \omega_R^\top m_R \ge u_N r_N - \|u_T\| \|r_T\| - \|\omega_R\| \|m_R\|, \quad (28)$$

by the reverse Cauchy-Schwarz inequality that ensures $u_T^\top r_T \ge -\|u_T\| \|r_T\|$ and $\omega_R^\top m_R \ge -\|\omega_R\| \|m_R\|$. Since $y \in L$ and $p \in K_r$, we obtain

$$y^\top p \ge u_N r_N - \mu r_N \|u_T\| - \mu_r r_N \|\omega_R\| \ge r_N (u_N - \mu \|u_T\| - \mu_r \|\omega_R\|) \ge 0, \quad (29)$$

and then $r_N \ge 0$ implies $L \subset K_r^*$. Let us choose $y \notin L$ such that $\mu \|u_T\| + \mu_r \|\omega_R\| > u_N$. Let us choose a special element on the boundary of $K_r$ but in the opposite direction to $y$, that is, $p$ such that

$$\|r_T\| = \mu r_N, \quad \|m_R\| = \mu_r r_N, \quad \|r_T\| \|u_T\| = -\|u_T\| \|r_T\|, \quad \|m_R\| \|\omega_R\| = -\|\omega_R\| \|m_R\|. \quad (30)$$

Let us now compute the scalar product

$$\begin{aligned} y^\top p &= u_N r_N + u_T^\top r_T + \omega_R^\top m_R = u_N r_N - \|r_T\| \frac{u_T^\top}{\|u_T\|} u_T - \|m_R\| \frac{\omega_T^\top}{\|\omega_R\|} \omega_R \\ &= u_N r_N - \|r_T\| \|u_T\| - \|m_R\| \|\omega_R\| \\ &= u_N r_N - \mu r_N \|u_T\| - \mu_r r_N \|\omega_R\| = r_N (u_N - \mu \|u_T\| - \mu_r \|\omega_R\|) < 0 \end{aligned} \quad (31)$$

so $y \notin K_r^*$. $\square$

Inspired by the definition of the dual cone in Lemma 3.1, the following proposition defines the rolling friction model as a complementarity problem over the cone $K_r$ and proves the equivalence between the formulations (21) and (22).

**Proposition 3.1.** Let us define the modified velocity as

$$\hat{y} := y + [\mu \|u_T\| + \mu_r \|\omega_R\|, 0, 0, 0, 0]^\top. \quad (32)$$

Let us consider the following cone complementarity problem.

$$K_r^* \ni \hat{y} \perp p \in K_r \iff -\hat{y} \in N_{K_r}(p). \quad (33)$$

Then, the models in (21), (22) and (33) are equivalent.

**Proof:** Let us introduce the notation $\hat{u}_N := u_N + \mu \|u_T\| + \mu_r \|\omega_R\|$. If $\hat{y} \in K_r^*$, then $\hat{u}_N \ge \mu \|u_T\| + \mu_r \|\omega_R\|$ or equivalently $u_N \ge 0$, so we have

$$\hat{y} \in K_r^* \iff u_N \ge 0. \quad (34)$$

Let us show that (21) implies (22). From (21), we have either $r_N = 0$ and $u_N \ge 0$ or $r_N \ge 0$ and $u_N = 0$ so the complementarity between $u_N$ and $r_N$ in (22) is satisfied. If $r_N = 0$, the model in (21)
implies that $r_T = 0$, $m_R = 0$ and the velocities are free to evolve in $\mathbb{R}^2$. In (22), the normal cone to the disk is given by $N_{D(0)} = \mathbb{R}^2$ for $r_N = 0$. For $r_N > 0$, let us recall the definition of $N_{D(c)}$ with $c > 0$:

$$ N_{D(c)}(x) = \begin{cases} 0 & \text{if } x < c, \\ \lambda \frac{x}{\|x\|}, \lambda \ge 0 & \text{if } x = c. \end{cases} \qquad (35) $$

In the sticking and no-rolling case ($r \in K, u = 0$), it is trivial since $u = 0$ belongs to the normal cone $N_{D(\mu r_N)}$. For the sticking cases, we have $\|r_T\| < \mu r_N$ and $u_T = 0$ so $-u_T \in N_{D(\mu r_N)}(r_T)$ is satisfied. For the sliding cases, we have $\|r_T\| = \mu r_N$ and $\|r_T\|u_T = -\|u_T\|\|r_T\|$ so $-u_T \in N_{D(\mu r_N)}(r_T)$ is satisfied. The same argumentation applies to show that $-\omega_R \in N_{D(\mu r_N)}(m_T)$.

Let us show that (22) implies (33). Since $-u_N \in N_{\mathbb{R}_+(r_N)}$ holds, we have $u_N \ge 0$ and then $\hat{y} \in K_r^*$. If (22) is satisfied, then $p \in K_r$. It remains to check the orthogonality in the different cases.

Starting from the constraints $p \in K_r$, let us enumerate the different cases:

a) $r_N = 0$. In that case, the model (22) implies $r_T = 0$ and $m_R = 0$ that is $p = 0$. The orthogonality condition is then trivially satisfied for all $\hat{y} \in K_r^*$. We conclude that (33) is satisfied.

b) $r_N > 0$, $\|r_T\| < \mu r_N$, $\|m_R\| < \mu r_N$. In that case, it is obvious that $p \in K_r^*$. The relations in (22) imply that $u_N = 0$, $u_T = 0$ and $\omega_R = 0$. So we have $\hat{y} = y = 0$, the orthogonality is also trivially satisfied.

c) $r_N > 0$, $\|r_T\| = \mu r_N$, $\|m_R\| < \mu r_N$. From (22), we conclude that $u_N = 0$ and $\omega_R = 0$. Since $\|r_T\| = \mu r_N$, from the definition of the normal cone in (35), we have $-u_T = \lambda \frac{r_T}{\|r_T\|}$, $\lambda \ge 0$. Taking the norm of the last expression yields $\lambda = \|u_T\|$, and then $\|r_T\|u_T = -\|u_T\|\|r_T\|$ and $u_T^\top r_T = -\|u_T\|\|r_T\|$. The orthogonality condition is

$$ \hat{y}^\top p = \hat{u}_N r_N + u_T r_T = u_N r_N + \mu \|u_T\| \|r_N - \|u_T\| \|r_T\|\|u_T\| = u_N r_N + (\mu r_N - \|r_T\|) \|u_T\| = 0, \quad (36) $$

and we conclude that (33) is satisfied.

d) $r_N > 0$, $\|r_T\| < \mu r_N$, $\|m_R\| = \mu r_N$. From (22), we conclude that $u_N = 0$ and $u_T = 0$. Since $\|m_R\| = \mu r_N$, we have $-\omega_R = \lambda \frac{m_R}{\|m_R\|}$, $\lambda \ge 0$. Similarly, as the previous case, we have $\lambda = \|\omega_R\|$ and then $\omega_R^\top m_R = -\|\omega_R\|\|m_R\|$. The orthogonality condition is

$$ \hat{y}^\top p = \hat{u}_N r_N + \omega_R m_R = u_N r_N + \mu_r \|\omega_R\| \|r_N - \|\omega_R\| \|m_R\|\|u_N\| = u_N r_N + (\mu_r r_N - \|m_R\|) \|\omega_R\| = 0, \quad (37) $$

and we conclude that (33) is satisfied.

e) $r_N > 0$, $\|r_T\| = \mu r_N$, $\|m_R\| = \mu r_N$. From (22), we have $u_N = 0$, $u_T^\top r_T = -\|u_T\|\|r_T\|$ and
$\omega_R^\top m_R = -\|\omega_R\|\|m_R\|$. The orthogonality condition is

$$ \hat{y}^\top p = \hat{u}_N r_N + u_T r_T + \omega_R m_R = u_N r_N + \mu \|u_T\| r_N - \|u_T\| \|r_T\| + \mu_r \|\omega_R\| r_N - \|\omega_R\| \|m_R\| = 0, \quad (38) $$

we conclude that (33) is satisfied.

Let us finally show that (33) implies (21). From (33), we have $r \in K_r$ and $u_N \ge 0$. Let us again enumerate the cases.

a) $p = 0$. Since we have $u_N \ge 0$ so it implies the take-off condition in (21).

b) $p \in \mathring{K}_r$. In that case, the complementarity implies $\hat{y} = 0$, and then $u_T = \omega_R = 0$ and finally $u_N = 0$ so it implies the sticking and no-rolling case.

c) $p \in \partial K_r$, $\|r_T\| = \mu r_N$, $\|m_R\| < \mu_r r_N$. Let us compute the normal cone to $K_r$ at $p$:

$$ N_{K_r}(p) = \left\{ y = \lambda_1 \begin{bmatrix} -\mu \\ \frac{r_T}{\|r_T\|} \\ 0 \end{bmatrix}, \lambda_1 \ge 0 \right\}. \quad (39) $$

From $-\hat{y} \in N_{K_r}(p)$, we conclude that

$$ u_T = \lambda_1 \frac{r_T}{\|r_T\|}, \lambda_1 \ge 0 \text{ and } \omega_R = 0. \quad (40) $$

Taking the norm of the first expression in (40) yields $\lambda_1 = \|u_T\|$, and then $\|r_T\|u_T = -\|u_T\|r_T$.

d) $p \in \partial K_r$, $\|r_T\| < \mu r_N$, $\|m_R\| = \mu_r r_N$. This case is similar to the previous one.

e) $p \in \partial K_r$, $\|r_T\| = \mu r_N$, $\|m_R\| = \mu_r r_N$. Let us compute the normal cone to $K_r$ at $p$:

$$ N_{K_r}(p) = \left\{ y = \lambda_1 \begin{bmatrix} -\mu \\ \frac{r_T}{\|r_T\|} \\ 0 \end{bmatrix} + \lambda_2 \begin{bmatrix} -\mu_r \\ 0 \\ \frac{m_R}{\|m_R\|} \end{bmatrix}, \lambda_1 \ge 0, \lambda_2 \ge 0 \right\}. \quad (41) $$

From $-\hat{y} \in N_K(p)$, we conclude

$$ u_T = \lambda_1 \frac{r_T}{\|r_T\|}, \lambda_1 \ge 0 \text{ and } \omega_R = \lambda_2 \frac{m_R}{\|m_R\|}, \lambda_2 \ge 0. \quad (42) $$

Taking the norm of the expressions in (42) yields $\lambda_1 = \|u_T\|$ and $\lambda_2 = \|\omega_R\|$, and then $\|r_T\|u_T = -\|u_T\|r_T$ and $\|m_R\|\omega_R = -\|\omega_R\|m_R$. $\square$

## 4. Maximum dissipation principle, optimization formulation and the bipotential method

Proposition 3.1 shows that the model of rolling friction is equivalent to a cone complementarity problem. In this section, we are interested in showing that the maximum dissipation principle can be deduced easily from the complementarity condition. In the formulation of the maximum
dissipation principle for standard Coulomb friction, only the tangent variables $u_T, r_T$ are involved. With a single inclusion into the normal cone to $K_r$, it is also possible to write a maximum dissipation principle over the whole variable $y := [u_N, u_T^\top, \omega_R^\top]^\top$ and $p := [r_N, r_T^\top, m_R^\top]^\top$.

To state dissipation principles, the power of the contact reaction $p$ with the relative velocity $y$ is computed as $y^\top p$. In the sequel, we show that this power is equal to the dissipated power at contact and is maximal over the choice of $y$ and $p$. This principle is closely related to the notion of bipotential function that is introduced and computed at the end of the section.

If the problem (33) has a solution, then we have $\hat{y}^\top p = 0$ by the complementarity condition. Coming back to the original variable $y$ and $p$, we get

$$ -y^\top p = (\mu \|u_T\| + \mu_r \|\omega_R\|) r_N, \quad (43) $$

which corresponds to the energy dissipation at contact by the friction forces and moments.

By maximizing the inner product of $\hat{y}$ and $p$ in (33), the following optimization problem can be written

$$ \begin{array}{ll} \max_{\hat{y},p} & -\hat{y}^\top p \\ \text{s.t.} & \hat{y} \in K_r^* \\ & p \in K_r \end{array} \quad (44) $$

The problem is feasible since $K_r \neq \emptyset$. Since $\hat{y} \in K_r^*$ and $p \in K_r$, we have $-\hat{y}^\top p \le 0$. The maximizer exists since the cost function is bounded from above. If the problem (33) has a solution, then it is also a solution of the optimization problem (44) with $-\hat{y}^\top p = 0$.

Let us now have a look at the rolling friction law in original variables $y$ and $p$. Using (34) and substituting $\hat{y}$ in the cost function, we get an equivalent optimization problem as

$$ \begin{array}{ll} \max_{y,p} & -(\mu \|u_T\| + \mu_r \|\omega_R\|)r_N - y^\top p \\ \text{s.t.} & u_N \ge 0, \\ & p \in K_r. \end{array} \quad (45) $$

For a given value of $y$ such that $u_N = 0$ and a given value of $r_N \ge 0$, the maximization problem (45) reduces to

$$ \begin{array}{ll} \max_{r_T, m_R} & -r_T u_T - m_R \omega_R \\ \text{s.t.} & \|r_T\| \le \mu r_N \\ & \|m_R\| \le \mu_r r_N \end{array} \quad (46) $$

The problem (46) corresponds to the maximum dissipation principle for the rolling friction problem.

**Remark 1.** The assumption of an associated law of friction yields the following model of rolling friction

$$ K^* \ni y \perp p \in K \iff -y \in N_{K_r}(p). \quad (47) $$
This model has several drawbacks from the modeling point of view that can be listed :

* The first drawback is the dilatancy of the model when the contact slides and/or rolls ($\|u_T\| > 0, \|\omega_R\| > 0$). Indeed, we get from the definition of the normal cone in (39) that

$$u_N = \mu \|u_t\| + \mu_r \|\omega_R\| > 0 \qquad (48)$$

In other words, the normal relative velocity is positive with sliding and/or rolling.

* The second drawback is the fact that the model does not dissipate energy at contact. Indeed, the orthogonality condition in (47), $y^\top p = 0$ implies that the dissipation vanishes.

Let us now have a look at how to define the bipotential for the rolling friction model. Firstly, let us recall now the definition of a bipotential in $\mathbb{R}^n$ with the standard scalar product.

**Définition 4.1** (Bipotential mapping and extremal points [8]). Let $b : \mathbb{R}^n \times \mathbb{R}^n \to \mathbb{R} \cup +\infty$ such that $v \mapsto b(v, f)$ and $f \mapsto b(v, f)$ are convex. This function is called a bipotential if

$$b(v', f') \ge v'^\top f', \quad v' \in \mathbb{R}^n, f' \in \mathbb{R}^n \qquad (49)$$

A couple $(v, f)$ is said to be extremal if it satisfies

$$b(v, f) = v^\top f \qquad (50)$$

For an extremal couple $(v, f)$, we verify easily that

$$b(v', f) - b(v, f) \ge f'^\top (v' - v), \quad \forall v' \qquad (51)$$

$$b(v, f') - b(v, f) \ge v'^\top (f' - f), \quad \forall f' \qquad (51)$$

Using the definition of the subdifferential, the extremal couple satisfies

$$v \in \partial_f b(v, f), \quad \text{and} \quad f \in \partial_v b(v, f) \qquad (52)$$

The relations (51) and (52) can also be related to optimization problems. An extremal point is also a solution of

$$(v, f) = \underset{p}{\operatorname{argmin}} \quad b(v, f) - v^\top f \qquad (53)$$

since the cost function $b(v, f) - v^\top f$ is positive from (49) and vanishes at the extremal points from (50). Taking the optimality conditions of (53) we get (52).

Coming back to the rolling friction problem and introducing $i_{\mathbb{R}_+}$ and $i_{K_r}$, respectively the indicator functions of $\mathbb{R}_+$ and $K_r$, the optimization problem (45) can be written as

$$\min_{y,p} (\mu \|u_T\| + \mu_r \|\omega_R\|) r_N + y^\top p + i_{\mathbb{R}_+}(u_N) + i_{K_r}(p). \qquad (54)$$

By inspecting (53) and (54), we propose to define the bipotential for the rolling friction by taking $v = -y$ and $f = p$ and $b(-y, p) := (\mu \|u_T\| + \mu_r \|\omega_R\|)r_N + i_{\mathbb{R}_+}(u_N) + i_{K_r}(p)$.
**Proposition 4.2.** *The mapping defined by*

$$b(-y, p) := (\mu\|u_T\| + \mu_r\|\omega_R\|)r_N + i_{\mathbb{R}_+}(u_N) + i_{K_r}(p) \quad (55)$$

is a bi-potential. Furthermore, an extremal point of b satisfies the rolling friction problem (33).

**Proof:** The mapping $p \mapsto b(-y, p)$ is the sum of a convex function as an indicator function of a convex set and a linear function. The mapping $-y \mapsto b(-y, p)$ is the sum of a convex function $(\mu\|u_T\| + \mu_r\|\omega_R\|)r_N$ for $r_N > 0$ and the indicator function of $\mathbb{R}_+$. It remains to prove that

$$b(-y', q') \geq -y'q', \quad \forall y', p' \quad (56)$$

If $y' \notin K_r^*$ i.e. $u_N < 0$ or $p' \notin K_r$, it is trivial. Let us assume that $y' \in K_r^*$ and $p' \in K_r$. The bipotential is then reduced to

$$b(-y, p) = (\mu\|u_T'\| + \mu_r\|\omega_T'\|)r_N \quad (57)$$

Since $p' \in K$, we obtain

$$(\mu\|u_T'\| + \mu_r\|\omega_T'\|)r_N' \geq \|u_T'\|\|r_T'\| + \|\omega_T'\|\|m_T'\| \geq \|u_T'\|\|-r_T'\| + \|\omega_T'\|\|-m_T'\| \geq -u_T'^{\top}r_T' - \omega_T'^{\top}m_T' \quad (58)$$

Since $y' \in K_r^*$ and $p' \in K_r$, we have also $-u_n' r_n' \geq 0$. The inequality (56) is then proved and b is a bipotential.

Let us consider a solution of (33). Then the equation (43) is satisfied, and we have $b(-y, p) = -y'p$ so it is an extremal point. The first inclusion in (52) yields

$$-y \in N_{K_r}(p) + [\mu\|u_T\| + \mu_r\|\omega_T\|, 0, 0, 0, 0]^T \iff -\hat{y} \in N_{K_r}(p). \quad (59)$$

□

**5. Projected based numerical methods for solving rolling friction problems**

Projected based numerical methods for solving standard Coulomb friction problems are based on a reformulation of the complementarity condition as a projection operator. Let us first recall the following result.

**Proposition 5.1 ([15]).** Let K a closed convex cone. Then $p = \text{proj}_K(\bar{p})$ if and only if

$$K^* \ni p - \bar{p} \perp p \in K \iff -(p - \bar{p}) \in N_K(p). \quad (60)$$

From the definition of the complementarity problem in Proposition 3.1 and the characterization of the projection in Proposition 5.1, we get the following equivalence

$$K_r^* \ni \hat{y} \perp p \in K_r \iff -\hat{y} \in N_{K_r}(p) \iff p = \text{proj}_{K_r}(p - \rho\hat{y}), \rho > 0. \quad (61)$$
The last term is an equality formulation of the complementarity problem based on the natural map [12]. This equation is prone to the development of numerical methods based on fixed-point iterations or semi-smooth Newton techniques [2].

Let us assume that the local relative velocity $y$ is related to the reaction torque $p$ by a linear relation of the form

$$y = Wp + q. \tag{62}$$

This relation can be obtained from the linearization of quasi-static balance equations or the time-discretization and a linearization of equations of motion and a reduction to local variables (see [3, 4])

Altogether, we obtain a nonlinear cone complementarity to solve:

$$\left\{
\begin{array}{l}
y = Wp + q, \\
\hat{y} = y + g(y) \\
K_r^* \ni \hat{y} \perp p \in K_r
\end{array}
\right.
\tag{63}$$

where $g(y) = [\mu \|u_T\| + \mu_r \|\omega_R\|, 0, 0, 0]^T$. This problem is usually a nonsmooth and nonlinear complementarity problem for which it is rather difficult to state a very general result of existence and uniqueness. Nevertheless, existence results should follow the line in [5].

For a finite set of $n_c$ contact points and their associated local frames, it is possible to write exactly the same problem as in (63) where $y \in \mathbb{R}^{5n_c}$, $p \in \mathbb{R}^{5n_c}$ collects the contact variables for each contact, $K_r^* \subset \mathbb{R}^{5n_c}$ and $K_r \subset \mathbb{R}^{5n_c}$ are the Cartesian product of local contact cones and $W \in \mathbb{R}^{5n_c \times 5n_c}$, $q \in \mathbb{R}^{5n_c}$, $\mu \in \mathbb{R}^{n_c}$ are the problem data.

### 5.1. Projected fixed point techniques

Projected fixed point techniques are standard methods to solve non-linear variational inequalities, and especially, complementarity problems. The method is based on the following iterations

$$p_{k+1} \leftarrow \operatorname{proj}_{K_r}(p_k - \rho_k(Wp_k + q + g(Wp_k + q))), \tag{64}$$

with $\rho_k > 0$. The fact that the projection is given in closed form in the sequel, is crucial for the implementation of the method that will not be related to an optimization at each iteration. The method can be improved by using the so-called extra-gradient algorithm:

$$\begin{align}
\tilde{p}_k &\leftarrow \operatorname{proj}_{K_r}(p_k - \rho(Wp_k + q + g(Wp_k + q))) \tag{65} \\
p_{k+1} &\leftarrow \operatorname{proj}_{K_r}(p_k - \rho_k(W\tilde{p}_k + q + g(W\tilde{p}_k + q))), \nonumber
\end{align}$$

with $\rho_k > 0$. The choice of $\rho_k$ is also important for the convergence and the efficiency of the method. Self-adaptive rules, that try to guess the Lipschitz constant of $p \mapsto Wp + q + g(Wp + q)$
are available in the literature. For more details in the context of contact friction problems, we refer to [2].

## 5.2. Block-splitting techniques

As in the standard Gauss-Seidel technique for frictional contact problems [19, 22], a way to improve the efficiency of the projected fixed-point methods is to perform a block splitting of the matrix $W$ contact by contact. The linear relation (62) for a contact $\alpha \in [1, n_c]$ reads as

$$y^{\alpha} = W^{\alpha\alpha} p^{\alpha} + q^{\alpha} + \sum_{\beta \in [1, n_c], \beta \neq \alpha} W^{\alpha\beta} p^{\beta}. \qquad (66)$$

With provisional values for $p^{\beta}, \beta \in [1, n_c], \beta \neq \alpha$, a one-contact problem

$$\left\{
\begin{array}{l}
y^{\alpha} = W^{\alpha\alpha} p^{\alpha} + \bar{q}^{\alpha}, \\
\hat{y}^{\alpha} = y^{\alpha} + g(y^{\alpha}) \\
K_{r}^{\alpha,*} \ni \hat{y}^{\alpha} \perp p^{\alpha} \in K_{r}^{\alpha}
\end{array}
\right.
\qquad (67)$$

is formulated with $\bar{q}^{\alpha} = q^{\alpha} + \sum_{\beta \in [1, n_c], \beta \neq \alpha} W^{\alpha\beta} p^{\beta}$. This problem can be solved with fixed-point solvers described in Section 5.1. The Gauss-Seidel technique is based on iterating over the contact points.

**Remark 2.** Since the problem (63) is formulated as a standard cone complementarity problem, or more generally as a variational inequality, this opens the door to use a lot of other mathematical programming techniques to solve it. To cite a few, we can apply second order techniques as semi-smooth Newton methods or interior point methods. In semi-smooth Newton methods, the definition of the projection is of great help if the complementarity function is based on the natural map [2, 12]. Other first order techniques can also be applied such as accelerated projected gradient method or Alternating Direction Method of Multipliers (ADMM) [6]. In [26], the proposed numerical method for the rolling contact problem is a fixed point method with projection based on a Quasi-Variational Inequality (QVI) formulation that considers separated projections for $r_N, r_T$ and $m_R$ onto the positive orthant, and the associated disks. In the QVI formulation, the set defining the variational inequality depends on the solution. This renders the problem far more difficult, and the number of available algorithms quite restricted. As an illustration, the method used in [26] has several $\rho$ parameters as in (64) but there is no standard rules to size it. As the authors mentioned, the convergence of the algorithm depends strongly on the choice of these parameters.

Finally, we want to stress that the problem (63) is formulated on cones, and not on generic convex sets. This means that the complementarity condition can be written, and that the computation of the dual cone is simple. On the contrary, in quasi-variational inequalities, the evaluation of the support function of a convex set is not straightforward.
### 5.3. Closed-form formulae for the projection onto the rolling friction cone

As it can be observed in the previous sections, a key element of the projected based numerical methods is the computation of the projection on $K_r$. This computation must be efficient and a closed-form formulae for the projection is not straightforward. In this section, we derive it in the following proposition using Convex Analysis tools.

**Proposition 5.2.** The projection onto $K_r$ of a vector $\bar{p} = [\bar{r}_N, \bar{r}_T^\top, \bar{m}_R^\top]^\top$ denoted by $p = \text{proj}_{K_r}(\bar{p})$ is given by

1. if $\bar{p} \in K_r$, then $p = \text{proj}_{K_r}(\bar{p}) = \bar{p}$,

2. if $-\bar{p} \in K_r^*$ then $p = \text{proj}_{K_r}(\bar{p}) = 0$,

3. if $\mu r_N - \|\bar{r}_T\| < 0$ and $\mu r_N - \|\bar{m}_R\| < 0$ with $r_N = \frac{\bar{r}_N + \mu \|\bar{r}_T\| + \mu_r \|\bar{m}_R\|}{1 + \mu^2 + \mu_r^2}$, then

$$ p = \text{proj}_{K_r}(\bar{p}) = \begin{bmatrix} r_N = \frac{\bar{r}_N + \mu \|\bar{r}_T\| + \mu_r \|\bar{m}_R\|}{1 + \mu^2 + \mu_r^2} \\ r_T = \mu r_N \frac{\bar{r}_T}{\|\bar{r}_T\|} \\ m_R = \mu_r r_N \frac{\bar{m}_R}{\|\bar{m}_R\|} \end{bmatrix}. \quad (68) $$

4. if $\mu r_N - \|\bar{r}_T\| < 0$ and $\mu r_N - \|\bar{m}_R\| > 0$ with $r_N = \frac{\bar{r}_N + \mu \|\bar{r}_T\|}{1 + \mu^2}$, then

$$ p = \text{proj}_{K_r}(\bar{p}) = \begin{bmatrix} r_N = \frac{\bar{r}_N + \mu \|\bar{r}_T\|}{1 + \mu^2} \\ r_T = \mu r_N \frac{\bar{r}_T}{\|\bar{r}_T\|} \\ m_R = \bar{m}_R \end{bmatrix}. \quad (69) $$

5. if $\mu r_N - \|\bar{m}_R\| < 0$ and $\mu r_N - \|\bar{r}_T\| > 0$ with $r_N = \frac{\bar{r}_N + \mu_r \|\bar{m}_R\|}{1 + \mu_r^2}$, then

$$ p = \text{proj}_{K_r}(\bar{p}) = \begin{bmatrix} r_N = \frac{\bar{r}_N + \mu_r \|\bar{m}_R\|}{1 + \mu_r^2} \\ r_T = \bar{r}_T \\ m_R = \mu_r r_N \frac{\bar{m}_R}{\|\bar{m}_R\|} \end{bmatrix}. \quad (70) $$

**Proof:** For the proof, we extensively use the characterization of the projection by complementarity given in Proposition 5.1. Let us distinguish the various cases:

1) $\bar{p} \in K_r$, then the projection $p = \bar{p}$ satisfies the condition (60) trivially.

2) $\bar{p} \in K_r^o$, then the projection $p = 0$ satisfies the condition (60) trivially.

Let us now consider the case where $\bar{p} \notin K_r$ and $\bar{p} \notin K_r^o$. We know that the projection will lie the boundary of $K_r$ denoted by $\partial K_r$ with $r_N > 0$. The case $r_N = 0$ is excluded by contradiction : if $r_n = 0$, then $p = 0$, and by (28), $\bar{p} \in K_r^*$.

We have three options depending on the faces of the boundary on which we project:

* $r_N > 0$, $\|r_T\| = \mu r_N$ and $\|m_R\| < \mu_r r_N$,
• $r_N > 0$, $\|m_R\| = \mu_r r_N$ and $\|r_T\| < \mu r_N$,

• $r_N > 0$, $\|r_T\| = \mu r_N$ and $\|m_R\| = \mu r_N$.

3) $r_N > 0$, $\mu r_N - \|\bar{r}_T\| < 0$ and $\mu r_N - \|\bar{m}_R\| < 0$. Let us assume for a while that $\|r_T\| = \mu r_N$ and $\|m_R\| = \mu r_N$. From the normal cone definition, the inclusion $-(p - \bar{p}) \in N_{K_r}(p)$ implies

$$ \begin{aligned} (r_T - \bar{r}_T) &= \lambda_1 \frac{\bar{r}_T}{\|r_T\|}, \lambda_1 \ge 0, \\ (m_R - \bar{m}_R) &= \lambda_2 \frac{\bar{m}_R}{\|m_R\|}, \lambda_2 \ge 0. \end{aligned} \tag{71} $$

From the first expression in (71), we get

$$ \bar{r}_T = \left(1 + \frac{\lambda_1}{\|r_T\|}\right) r_T \tag{72} $$

Since $\lambda_1 \ge 0$, the vectors are colinear and in the same direction. Since $\|r_T\| = \mu r_N$, we obtain $r_T = \mu r_N \frac{\bar{r}_T}{\|\bar{r}_T\|}$. The same manipulations can be made for the second expression in (71) leading to

$$ r_T = \mu r_N \frac{\bar{r}_T}{\|\bar{r}_T\|} \quad \text{and} \quad m_R = \mu_r r_N \frac{\bar{m}_R}{\|\bar{m}_R\|}. \tag{73} $$

The orthogonality constraint in the complementarity condition (60) yields:

$$ \begin{aligned} 0 = p^\top (p - \bar{p}) &= r_N(r_N - \bar{r}_N) + r_T^\top (r_T - \bar{r}_T) + m_R^\top (m_R - \bar{m}_R) \\ &= r_N^2 - r_N \bar{r}_N + \|r_T\|^2 - r_T^\top \bar{r}_T + \|m_R\|^2 - m_R^\top \bar{m}_R \\ &= r_N^2(1 + \mu^2 + \mu_r^2) - r_N \bar{r}_N - r_T^\top \bar{r}_T - m_R^\top \bar{m}_R. \end{aligned} \tag{74} $$

Substituting (73) in (74), we obtain

$$ \begin{aligned} 0 = p^\top (p - \bar{p}) &= r_N^2(1 + \mu^2 + \mu_r^2) - r_N \bar{r}_N - \mu r_N \|\bar{r}_T\| - \mu_r r_N \|\bar{m}_R\| \\ &= r_N(r_N((1 + \mu^2 + \mu_r^2)) - (\bar{r}_N + \mu \|\bar{r}_T\| + \mu_r \|\bar{m}_R\|)), \end{aligned} \tag{75} $$

and then, for $r_N > 0$,

$$ r_N = \frac{\bar{r}_N + \mu \|\bar{r}_T\| + \mu_r \|\bar{m}_R\|}{1 + \mu^2 + \mu_r^2}. \tag{76} $$

Let us check now that $p - \bar{p} \in K_r^*$. Since $0 = p^\top (p - \bar{p})$, we must have $p - \bar{p} \in \partial K_r^*$, that is

$$ r_N - \bar{r}_N = \mu \|r_T - \bar{r}_T\| + \mu_r \|m_R - \bar{m}_R\|. \tag{77} $$

On one side, we have

$$ \begin{aligned} r_N - \bar{r}_N &= -(\mu^2 + \mu_r^2)r_N + \mu \|\bar{r}_T\| + \mu_r \|\bar{m}_R\| \\ &= \mu(\|\bar{r}_T\| - \mu r_N) + \mu_r (\|\bar{m}_R\| - \mu_r r_N). \end{aligned} \tag{78} $$

On the other side, we get using (73)

$$ \|r_T - \bar{r}_T\| = \|\mu r_N \frac{\bar{r}_T}{\|\bar{r}_T\|} - \bar{r}_T\| = \left| \frac{\mu r_N}{\|\bar{r}_T\|} - 1 \right| \|\bar{r}_T\| = |\mu r_N - \|\bar{r}_T\||, \tag{79} $$
and

$$ \|m_R - \bar{m}_R\| = \|\mu_r r_N \frac{\bar{m}_R}{\|\bar{m}_R\|} - \bar{m}_R\| = \left| \frac{\mu_r r_N}{\|\bar{m}_R\|} - 1 \right| \|\bar{m}_R\| = |\mu_r r_N - \|\bar{m}_R\||, \quad (80) $$

so we have

$$ \mu \|r_T - \bar{r}_T\| + \mu_r \|m_R - \bar{m}_R\| = \mu |\mu r_N - \|\bar{r}_T\| | + \mu_r |\mu_r r_N - \|\bar{m}_R\| |. \quad (81) $$

If the assumption

$$ \mu r_N - \|\bar{r}_T\| < 0 \text{ and } \mu_r r_N - \|\bar{m}_R\| < 0, \quad (82) $$

is satisfied, we deduce from (81) that the equality (77) holds. This assumption means that the value $\|\bar{m}_R\|$ and $\|\bar{r}_T\|$ must be outside the disks of radius $\mu r_N$ and $\mu_r r_N$, where $r_N$ is the forecast value given by (76).

4) $r_N > 0$, $\|r_T\| = \mu r_N$ and $\|m_R\| < \mu_r r_N$. Let us assume for a while that $\|r_T\| = \mu r_N > 0$ and $m_R = \bar{m}_R$. From $-(p - \bar{p}) \in N_{K_r}(p)$, we conclude as before

$$ r_T = \mu r_N \frac{\bar{r}_T}{\|\bar{r}_T\|} \quad \text{and} \quad m_R = \bar{m}_R. \quad (83) $$

Let us remark that in that case the condition $\|m_R\| < \mu_r r_N$ is equivalent to $\|\bar{m}_R\| < \mu_r r_N$. The orthogonality $p^\top(p - \bar{p}) = 0$ with (83) yields

$$ r_N = \frac{\bar{r}_N + \mu \|\bar{r}_T\|}{1 + \mu^2}. \quad (84) $$

Let us check that $p - \bar{p} \in \partial K_r$, that is

$$ r_N - \bar{r}_N = \mu \|r_T - \bar{r}_T\| + \mu_r \|m_R - \bar{m}_R\|. \quad (85) $$

Let us derive the formulae as we did before. On one side, we have

$$ r_N - \bar{r}_N = \mu (\|\bar{r}_T\| - \mu r_N). \quad (86) $$

On the other side we get

$$ \mu \|r_T - \bar{r}_T\| + \mu_r \|m_R - \bar{m}_R\| = \mu \|r_T - \bar{r}_T\| = \mu |\mu r_N - \|\bar{r}_T\|| \quad (87) $$

If the assumption

$$ \mu r_N - \|\bar{r}_T\| < 0 \quad (88) $$

is satisfied, then the equality (60) holds.

5) $r_N > 0$, $\|r_T\| < \mu r_N$ and $\|m_R\| = \mu_r r_N$. Doing as in the previous case, we get the definition of the projection. □
# 6. Applications of the rolling friction model

In this section, the abilities of the rolling friction model to dissipate energy by friction and/or rolling friction are illustrated on simple numerical examples. The simulations have been realized with the Moreau–Jean time–stepping scheme [18, 3] based on a $\theta$-method for the smooth terms and the projected Gauss-Seidel presented in Section 5.2 for the discrete frictional contact problems. The numerical methods are implemented in the Siconos software v4.3.0 [1] and are freely available together with the simulation scripts.

## 6.1. A sphere rolling on a plane

The problem of a sphere of radius $R$ rolling on a plane with initial tangential and rotational velocities is chosen to illustrate the predictive capacity of the rolling friction model compared with the classical friction one. The sphere is initially located on an horizontal planar surface with a tangential velocity $U_X(t = 0)$ and a rotational velocity $\Omega_Y(t = 0)$. The simulation has been made with $\theta = 1/2$ for the parameter of the Moreau–Jean scheme, a time-step of $10^{-4}$s and a relative tolerance of the Gauss–Seidel algorithm of $10^{-12}$.

Under these conditions, the sphere does not stop to roll with a classical friction model. In the case of initial sliding at the contact point $u_T(0) > 0$, the norm of the tangential reaction $\|r_T\|$ is equal to $\mu r_n$ and is in the opposite direction to $u_T$, which induces a decrease of the sliding velocity till the slip-free rolling condition $u_T(t) = U_X(t) - R\Omega_Y(t) = 0$ is reached. Once this condition is reached, $r_T$ jumps to 0 and the sphere velocities $U_X$ and $\Omega_Y$ are constant, which prevents the sphere to stop.

The rolling friction model is first investigated for initial sphere velocities fulfilling the slip-free rolling conditions ($U_X(0) - R\Omega_Y(0) = 0$). In this case, the maximum propagation distance $X_{max}$ reached by the sphere depends on both the friction ($\mu$) and of the rolling friction ($\mu_r$) coefficients (see Figure 2). For $\mu_r \to 0$, as expected, the stopping distance goes to infinity. Similarly, large propagation distances are observed for $\mu \to 0$. Indeed, although the resistive moment $m_T$ progressively reduces the sphere rotational velocity, the small values of $r_T$ do not allow reducing $U_X$ which finally entails block sliding without rotational velocity. The propagation distance is strongly influenced by the values of $\mu$ and $\mu_r$ for $\mu \in [0, 0.15]$ and $\mu_r \in [0, 0.01]$. In this domain, $X_{max}$ strongly decreases if $\mu$ or $\mu_r$ increase. For larger values of both $\mu$ and $\mu_r$, the propagation distance is much less sensitive to these coefficients.

We choose to illustrate the details of the rolling friction model behavior through three representative configurations. The first configuration (Figure 3a) corresponds to initial slip-free rolling conditions with large values of both $\mu$ and $\mu_r$. The resistive moment $m_T$ entails a decrease in both $U_X$ and $\Omega_Y$ and the slip-free condition is conserved (Figure 3a) till the sphere stops. The second configuration is very similar to the first one, except that slippage initially occurs at the contact
<figure><img src="image_4.png" /><figcaption>Figure 2: Evolution of the maximum propagation distance $X_{max}$ reached by the sphere depending on the friction ($\mu$) and rolling friction ($\mu_r$) coefficients for initial slip-free rolling conditions. The sphere radius and density were set at 0.5m and 2500kg/m$^3$, respectively, while the initial velocities were $U_X(t=0) = 2.5$m/s and $\Omega_Y(t=0) = 5$rad/s</figcaption></figure>

point (Figure 3b), which results in velocities changes before reaching the slip-free condition. Afterwards, the decrease in both $U_X$ and $\Omega_Y$ entails sphere stopping. At the transition from slippage to slip-free rolling condition, the tangential force jumps from $\|r_T\| = \mu r_n$ to $\|r_T\| < \mu r_n$. Let us note that $r_T$ does no vanish as in the case of standard friction without rolling resistance. The third example (Figure 3c) illustrates the case where the rotational velocity $\Omega_Y$ vanishes before the sphere stopping because $\mu$ is small compared to $\mu_r$. In this case, the vanishing of $\Omega_Y$ induces a jump of the resitive moment from $\|m_T\| = \mu_r r_n$ to $\|m_T\| < \mu_r r_n$.

## 6.2. Spheres piles

In this section, the interest of the rolling friction model for the simulation of a granular material made of spheres is illustrated. To this end, we consider a small sample of 250 spheres of diameter $\phi = 0.02m$, dropped at a constant rate from a constant height of 10$\phi$ on a square fixed rigid plane of 60$\phi$ width. To avoid the constitution of artificial columns in the sample, the initial horizontal coordinates of the spheres are randomly sampled in the range $[-\phi/2; \phi/2]$. The density of the spheres is 1300kg.m$^{-3}$, the gravity is equal to 9.81m.s$^{-2}$ and the coefficient of restitution is equal to 0. The simulation campaign was carried out with three values of the coefficient of friction $\mu \in \{0.1, 0.3, 1.0\}$ and with values of the coefficient of rolling friction ranging in $[10^{-5}, 1]$. The simulation has been made with $\theta = 1$ for the parameter of the Moreau-Jean scheme, a time-step of $10^{-3}$s over a time interval of 20s, and a relative tolerance of the Gauss-Seidel algorithm of $10^{-4}$.

In Figure 4, the final state of the pile when the kinetic energy vanishes is depicted for different values of the coefficients of friction. We can observe that the height and the width of the pile
<figure><img src="image_5.png" /><figcaption>Figure 3: Evolution of the sphere position X, tangential velocity $U_X$, and rotational velocity $\Omega_Y$ as well as tangential force $r_T$ and resitive moment $m_T$. The sphere radius and density were set at 0.5m and 2500kg/m$^3$, respectively. Three typical initial conditions are presented: a) slip-free rolling conditions ($U_X(t=0) = 2.5\,\text{m/s}$, $\Omega_Y(t=0) = 5\,\text{rad/s}$, $\mu=0.2$, $\mu_r=0.04$), b) sliding ($U_X(t=0) = 2.5\,\text{m/s}$, $\Omega_Y(t=0) = 2.5\,\text{rad/s}$, $\mu=0.2$, $\mu_r=0.04$), c) slip-free rolling conditions for small $\mu$ values ($U_X(t=0) = 2.5\,\text{m/s}$, $\Omega_Y(t=0) = 5\,\text{rad/s}$, $\mu=0.05$, $\mu_r=0.04$)</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 4: Final state of spheres piles when the kinetic energy vanishes.</figcaption></figure>

depend strongly on both coulomb's friction and rolling friction. In Figure 5, some measurements of the height of the pile and its mean radius are reported. The mean radius of the pile is computed by selecting the spheres that are in contact with the horizontal surface, which constitute the bottom of the pile. The centroid of the assembly of these spheres is computed in the horizontal plane and the mean radius is the average distance to the centroid of these spheres. In Figure 5(a), the height of the pile is plotted with respect to the coefficient of rolling friction $\mu_r$ for the three values of the coefficient of friction $\mu$. As expected, the height of the pile is generally increasing with the coefficient of rolling friction since the spheres are subjected to the rolling resistance torque. Three regimes can be observed for all the values of the coefficient of friction $\mu$. In a first regime, for low coefficients of rolling friction, the ratio between the height and the diameter is equal to one, meaning that there is no superposition of spheres. In a second regime, after a critical value of the rolling friction coefficient, the height of the pile quickly increased up to a third regime where the height of the pile is stabilized around a value that does not depend on the coefficient of rolling friction. In this third regime, the height of the pile is higher when the friction coefficient increases. This can be explained by the fact that, for too small friction coefficients, the spheres in contact with the plane slide under the action of the weight of the pile. If the coefficient of friction increases, the slope angle of the pile is higher as it is illustrated on Figures 4(a)-(c)-(e). Let us note that the graph of the height with respect to the rolling friction coefficient is not strictly monotone. This is mainly due to the random character of the pluviation process that may generate avalanches of the spheres (as it can be seen in Figure 4(d)). The mean radius of the piles plotted in Figure 5(a) corroborates the observation on the height of the pile. We observe that the mean radius of the bottom of the pile quickly decreases in the second regime and is stable in the third regime.
<figure><img src="image_7.png" /><figcaption>Figure 5: Sphere piles simulation. Normalized height and mean radius versus the coefficient of rolling friction for three values of $\mu \in \{0.3, 0.5, 0.7\}$.</figcaption></figure>

## 7. Conclusion

In this article, the model of Coulomb's friction with unilateral contact is enriched with a resistance to rolling that is proportional to the normal pressure. With the goal of designing a tractable model from the computational point of view, the model is formulated as a cone complementarity problem over a single rolling friction cone. To this aim a modified relative velocity at contact has been introduced following the seminal work in [9]. The formulation is proved to be equivalent to the disjunctive form of the model and the model with three separated inclusions. The principle of maximal dissipation and the bipotential function have been delineated for the rolling friction, that allows us to understand the dissipation generated by the model. Two projected based numerical methods are given that are based on the projection onto the rolling friction cone. We give a closed-form formulae of the projection that is crucial for the efficiency of these numerical methods. The abilities of the model and the numerical techniques are shown on two examples of a sphere rolling on a plane and the study of a granular sample made of spherical particles.

## References

[1] V. Acary, O. Bonnefon, M. Brémond, O. Huber, F. Pérignon, and S. Sinclair. An introduction to Siconos. Technical Report RT-0340, INRIA, November 2019. http://siconos.gforge.inria.fr, http://github.com/siconos/siconos.
[2] V. Acary, M. Brémond, and O. Huber. *Advanced Topics in Nonsmooth Dynamics*, chapter On solving frictional contact problems: formulations and comparisons of numerical methods. Acary, V. and Brüls. O. and Leine, R. (eds). Springer Verlag, 2018.

[3] V. Acary and B. Brogliato. *Numerical methods for nonsmooth dynamical systems. Applications in mechanics and electronics*. Lecture Notes in Applied and Computational Mechanics 35. Berlin: Springer. xxi, 525 p. , 2008.

[4] V. Acary and F. Cadoux. *Recent Advances in Contact Mechanics, Stavroulakis, Georgios E. (Ed.)*, volume 56 of *Lecture Notes in Applied and Computational Mechanics*, chapter Applications of an existence result for the Coulomb friction problem. Springer Verlag, 2013.

[5] V. Acary, F. Cadoux, C. Lemaréchal, and J. Malick. A formulation of the linear discrete coulomb friction problem via convex optimization. *ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik*, 91(2):155–175, 2011.

[6] S. Boyd, N. Parikh, S. Chu, B. Peleato, and J. Eckstein. Distributed optimization and statistical learning via the alternating direction method of multipliers. *Found. Trends Mach. Learn.*, 3(1):1–122, January 2011.

[7] B. Brogliato. *Nonsmooth Mechanics: Models, Dynamics and Control*. Communications and Control Engineering. Springer-Verlag, London, 3rd edition, 1999.

[8] G. De Saxcé. A generalization of Fenchel's inequality and its applications to the constitutive laws. *Comptes Rendus de l'académie des Sciences. Paris*, 314, II:125-129, 1992.

[9] G. De Saxcé and Z. Q. Feng. New inequality and functional for contact with friction: The implicit standard material approach. *Mechanics of Structures and Machines*, 19(3):301–325, 1991.

[10] N. Estrada, E. Azéma, F. Radjai, and A. Taboada. Identification of rolling resistance as a shape parameter in sheared granular media. *Phys. Rev. E*, 84:011306, Jul 2011.

[11] N. Estrada, A. Taboada, and F. Radjaï. Shear strength and force transmission in granular media with rolling resistance. *Phys. Rev. E*, 78:021301, Aug 2008.

[12] F. Facchinei and J. S. Pang. *Finite-dimensional Variational Inequalities and Complementarity Problems*, volume I & II of Springer Series in Operations Research. Springer Verlag NY. Inc., 2003.

[13] Ch. Glocker. *Set-Valued Force Laws: Dynamics of Non-Smooth systems*, volume 1 of *Lecture Notes in Applied Mechanics*. Springer Verlag, 2001.
[14] M. D. Hersey. Rolling friction, i—historical introduction. *Journal of Lubrication Technology*, 91(2):260-263, 1969.

[15] J.B. Hiriart-Urruty and C. Lemaréchal. *Fundamentals of Convex Analysis*. Springer Verlag, 2001.

[16] J. Huang, M. Vicente da Silva, and K. Krabbenhoft. Three-dimensional granular contact dynamics with rolling resistance. *Computers and Geotechnics*, 49:289–298, 2013.

[17] J. Irazábal, F. Salazar, and E. Oñate. Numerical modelling of granular materials with spherical discrete particles and the bounded rolling friction model. application to railway ballast. *Computers and Geotechnics*, 85:220–229, 2017.

[18] M. Jean and J.J. Moreau. Dynamics in the presence of unilateral contacts and dry friction: a numerical approach. In G. Del Pietro and F. Maceri, editors, *Unilateral problems in structural analysis. II*, pages 151–196. CISM 304, Springer Verlag, 1987.

[19] M. Jean and J.J. Moreau. Dynamics of elastic or rigid bodies with frictional contact and numerical methods. In R. Blanc P. Suquet M. Raous, editor, *Publications du LMA*, pages 9–29, 1991.

[20] K.L. Johnson. Rolling resistance of a rigid cylinder on an elastic-plastic surface. *International Journal of Mechanical Sciences*, 14(2):145–148, 1972.

[21] K.L. Johnson. *Contact mechanics*. Cambridge University Press, 1985.

[22] F. Jourdan, P. Alart, and M. Jean. A Gauss Seidel like algorithm to solve frictional contact problems. *Computer Methods in Applied Mechanics and Engineering*, 155(1):31–47, 1998.

[23] D. Kadau, G. Bartels, L. Brendel, and D.E. Wolf. Pore stabilization in cohesive granular systems. *Phase Transitions*, 76(4-5):315–331, 2003.

[24] J. J. Kalker. The computation of three-dimensional rolling contact with dry friction. *International Journal for Numerical Methods in Engineering*, 14(9):1293–1307, 1979.

[25] J.J. Kalker. A fast algorithm for the simplified theory of rolling contact. *Vehicle System Dynamics*, 11(1):1–13, 1982.

[26] C. Le Saux, R. I. Leine, and Ch. Glocker. Dynamics of a rolling disk in the presence of dry friction. *Journal of Nonlinear Science*, 15(1):27–61, Feb 2005.

[27] R.I. Leine. Experimental and theoretical investigation of the energy dissipation of a rolling disk during its final stage of motion. *Arch Appl Mech*, 2009.
[28] R.I. Leine and Ch. Glocker. A set-valued force law for spatial coulomb-contensou friction. *European Journal of Mechanics - A/Solids*, 22(2):193–216, 2003.

[29] D. Ma, C. Liu, Z. Zhao, and Zhang. H. Rolling friction and energy dissipation in a spinning disc. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 470(2169):20140191, 2014.

[30] D. B. Marghitu and E. D. Stoenescu. Rigid body impact with moment of rolling friction. *Nonlinear Dynamics*, 50(3):597–608, Nov 2007.

[31] J.J. Moreau. On unilateral constraints, friction and plasticity. In G. Capriz and G. Stam-pacchia, editors, *New Variational Techniques in Mathematical Physics, CIME II ciclo 1973*, pages 175-322. Edizioni Cremonese, 1974.

[32] J.J. Moreau. Unilateral contact and dry friction in finite freedom dynamics. In J.J. Moreau and Panagiotopoulos P.D., editors, *Nonsmooth Mechanics and Applications*, pages 1-82. CISM 302, Spinger Verlag, 1988.

[33] P.D Panagiotopoulos. *Inequality problems in Mechanics and applications*. Birkhäuser., 1985.

[34] J. Pérès. *Mécanique Générale*. Masson, 1953.

[35] G. De Saxcé and Z.-Q. Feng. The bipotential method: A constructive approach to design the complete contact law with friction and improved numerical algorithms. *Mathematical and Computer Modelling*, 28(4):225–245, 1998. Recent Advances in Contact Mechanics.

[36] A. Tasora and M. Anitescu. A complementarity-based rolling friction model for rigid contacts. *Meccanica*, 48(7):1643–1659, Sep 2013.

[37] Y.C. Zhou, B.D. Wright, R.Y. Yang, B.H. Xu, and A.B. Yu. Rolling friction in the dynamic simulation of sandpile formation. *Physica A: Statistical Mechanics and its Applications*, 269(2):536–553, 1999.