# On enhanced descent algorithms for solving frictional multicontact problems: application to the discrete element method

Serge Dumont

► To cite this version:

Serge Dumont. On enhanced descent algorithms for solving frictional multicontact problems: application to the discrete element method. International Journal for Numerical Methods in Engineering, 2013, 93 (11), pp.1170-1190. 10.1002/nme.4424. hal-01113707

HAL Id: hal-01113707

https://hal.science/hal-01113707

Submitted on 10 Mar 2023

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On enhanced descent algorithms for solving frictional multicontact problems: application to the discrete element method

S. Dumont*,†

LAMFA, Université de Picardie Jules Verne - CNRS UMR 7352, 33, rue Saint-Leu, France

In this article, various numerical methods to solve multicontact problems within the nonsmooth discrete element method are presented. The techniques considered to solve the frictional unilateral conditions are based both on the bipotential theory introduced by G. de Saxcé and the augmented Lagrangian theory introduced by P. Alart. Following the ideas of Z.-Q. Feng a new Newton method is developed to improve these classical algorithms, and numerical experiments are presented to show that these methods are faster than the previous ones, provide results with a better quality, and are less sensitive to the numerical parameters. Moreover, a stopping criterion that ensures a good mechanical property of the solution is provided.

KEY WORDS: granular materials; contact mechanics; Newton algorithms; bipotential; augmented Lagrangian

## 1. INTRODUCTION

The simulation of granular media concerns a wide range of engineering applications, such as concrete or geomaterials. All these materials can be considered as particles in interaction. In this study, the interaction is the frictional contact.

The numerical simulation of frictional contact problem remains a difficult problem, because the frictional contact law is strongly nonlinear and multivalued. Nevertheless, several approaches exist to solve contact problems. A wide bibliography exists on this subject, and an extensive list of references can be found, for example, in [4] or in [5]. The most popular method is based on the penalty approximation, due to its simplicity. But its important drawback is that the numerical parameters that have to be considered are difficult to choose efficiently. Indeed, if the penalty coefficient is small, then the contact laws are not accurately satisfied, and if the penalty coefficient is larger, then the problem to be solved becomes stiff.

Another important family of approaches is based on the notion of augmented Lagrangian. But because the frictional contact problem is no longer a classical minimization problem, the standard approaches need to be extended. This article deals with two methods of this family, which are the Alart and Curnier [2] method and the De Saxcé and Feng [1] method. These methods have been already adapted to the solving of multicontact problems, in the LMGC90 software for the first method and in the MULTICOR software for the second.

More precisely, the Alart and Curnier method consists in a Newton method to find the saddle point of an augmented Lagrangian method where the convex of constraints depends on the solution.

<footnote>*Correspondence to: S. Dumont, LAMFA, Université de Picardie Jules Verne - CNRS UMR 7352, 33, rue Saint-Leu, 80 000 Amiens, France.</footnote>

<footnote>†E-mail: serge.dumont@u-picardie.fr</footnote>
The method developed by De Saxcé and Feng is based on the theory of implicit standard materials and the writing of a so-called bipotential, which is minimized using a gradient method.

In this article, we will consider the works of Feng et al. [3] to improve the resolution of the multi-body problem, firstly, within the bipotential framework and then within the augmented Lagrangian framework. The fundamental idea of this algorithm is to search the solution of the contact problem no longer as the solution of a minimization problem but as the solution of the Euler equations of the minimization problem. The equation is then solved using a Newton method.

In addition, finding a good stopping criterion remains a difficulty. Very often, algorithms are stopped when the norm of the increment becomes sufficiently small. The main drawback of this technique is that there is no guarantee on the quality of the numerical solutions. Moreover, the solutions obtained by various algorithms are then difficult to compare. In this article, I propose a stopping criterion that can be used with the whole methods presented. The main advantage of this criterion lies in the fact that the mechanical quality of the solution is ensured (penetration, Signorini-Coulomb contact law, and equation of motion), both locally at each contact point and globally on the overall behavior.

The article is organized as follows: in the next part, I present the equations to be solved within the discrete element method context, and the frictional contact law considered. In the third part, I first present two classical methods to numerically solve the full problem, the first one based on the bipotential theory and the second one on the augmented Lagrangian theory. Then, I show how these methods can be enhanced using an appropriate Newton method, and a stopping criterion that works for all these algorithms is provided. The last part of this article is devoted to the numerical experiments to show the main properties of these algorithms.

## 2. PROBLEM SETTING

### 2.1. The equations of motion of a multicontact system

Classically (see, e.g., [6–9]), the motion of a multicontact system is described using a global generalized coordinate $\mathbf{q}$ (for $N_p$ particles, $\mathbf{q} \in \mathbb{R}^{\tilde{d} \times N_p}$, where $\tilde{d} = 6$ for a three-dimensional (3D) problem and $\tilde{d} = 3$ for a two-dimensional (2D) problem). Because of the possible shocks between particles, the equations of motion have to be formulated in terms of differential measure equation:

$$ \mathbb{M}d\dot{\mathbf{q}} + \mathbf{F}^{\text{int}}(t, \mathbf{q}, \dot{\mathbf{q}})dt = \mathbf{F}^{\text{ext}}(t, \mathbf{q}, \dot{\mathbf{q}})dt + dR \quad (1) $$

where

* $\mathbb{M}$ represents the generalized mass matrix;

* $\mathbf{F}^{\text{int}}$ and $\mathbf{F}^{\text{ext}}$ represent the internal and external forces, respectively;

* $d\mathbf{R}$ is a nonnegative real measure, representing the reaction forces and impulses between particles in contact.

For the sake of simplicity and without loss of generality, only the external forces are considered in the following. The internal forces are neglected because the general case can be easily derived through a linearizing procedure.

Then, for the numerics, Equation (1) is integrated on each time interval $[t_k, t_{k+1}]$ and approximated using a $\theta$-method with $\theta \in [\frac{1}{2}, 1]$ for stability reasons (see [10, 11]).

Therefore, the classical approximation of Equation (1) yields

$$ \left\{ \begin{array}{l} \mathbb{M}(\dot{\mathbf{q}}_{k+1} - \dot{\mathbf{q}}_n) = \Delta t(\theta \mathbf{F}_{k+1} + (1-\theta)\mathbf{F}_k) + \mathbf{R}_{k+1} \\ \mathbf{q}_{k+1} = \mathbf{q}_k + \Delta t \theta \dot{\mathbf{q}}_{k+1} + \Delta t (1-\theta) \dot{\mathbf{q}}_k \end{array} \right. \quad (2) $$

We will denote $\dot{\mathbf{q}}_k^{\text{free}} = \dot{\mathbf{q}}_k + \mathbb{M}^{-1}\Delta t(\theta \mathbf{F}_{k+1} + (1-\theta)\mathbf{F}_k)$ the free velocity (velocity when the contact forces vanish). Then, the first equation in (2) becomes

$$ \dot{\mathbf{q}}_{k+1} = \dot{\mathbf{q}}_k^{\text{free}} + \mathbb{M}^{-1}\mathbf{R}_{k+1} \quad (3) $$
To write the contact law, for a contact c between two particles ($1 \le c \le N_c$, where $N_c$ is the total number of contact), we define the local-global mapping

$$ \left\{ \begin{array}{l} \mathbf{u}^c = P^*(\mathbf{q}, c)\dot{\mathbf{q}} \\ \mathbf{R} = P(\mathbf{q}, c)\mathbf{r}^c \end{array} \right. \qquad (4) $$

where $\mathbf{u}^c$ is the local relative velocity between the two bodies in contact and $\mathbf{r}^c$ is the local contact force ($\mathbf{u}^c, \mathbf{r}^c \in \mathbb{R}^d$ where $d$ is the dimension of the problem, and $P^*$ is the transpose of matrix $P$). We also denote $\mathbb{P}(\mathbf{q})$ the total-global mapping, for $\mathbf{u}$ and $\mathbf{r}$ in $\mathbb{R}^{d \times N_c}$ (vectors composed of all relative velocities and contact forces, respectively):

$$ \left\{ \begin{array}{l} \mathbf{u} = \mathbb{P}^*(\mathbf{q})\dot{\mathbf{q}} \\ \mathbf{R} = \mathbb{P}(\mathbf{q})\mathbf{r} \end{array} \right. \qquad (5) $$

In the discretization, a prediction of $\mathbf{q}$ is computed to estimate the mapping $\mathbb{P}(\mathbf{q})$ (see Equations (18) and (19) in the following).

Using Equations (2) and (5), we can write the discretization of the motion of a multicontact system, with frictional contact between particles, as

$$ \left\{ \begin{array}{l} \tilde{\mathbf{u}}_{k+1} = \tilde{\mathbf{u}}_{k}^{\text{free}} + \mathbb{W}\mathbf{r}_{k+1} \\ \text{law}_c (\tilde{\mathbf{u}}_{k+1}^c, \mathbf{r}_{k+1}^c) = \text{true} \quad \forall c \in \{1, 2, \dots, N_c\} \end{array} \right. \qquad (6) $$

where $\mathbb{W} = \mathbb{P}^*\mathbb{M}^{-1}\mathbb{P}$ is the Delassus operator, and $\tilde{\mathbf{u}}_{k}^{\text{free}} = \mathbb{P}^*\dot{\mathbf{q}}_{k}^{\text{free}}$ is the relative free velocity. Notice that a Newton impact law is also considered (see [8] and Equation (20) in the following), which modify $\mathbf{u}_k$ and $\mathbf{u}_k^{\text{free}}$ by $\tilde{\mathbf{u}}_k$ and $\tilde{\mathbf{u}}_k^{\text{free}}$, respectively.

The second equation in (6) is the implicit frictional contact law that is in our case the classical Signorini condition and Coulomb's friction law.

## 2.2. The frictional contact law

In the local coordinates system defined by the local normal vector $\mathbf{n}$ and the tangential vector $\mathbf{t} \perp \mathbf{n}$, any element $\mathbf{u}$ and $\mathbf{r}$ can be uniquely decomposed as $\mathbf{u} = u_n\mathbf{n} + \mathbf{u}_t$ and $\mathbf{r} = r_n\mathbf{n} + \mathbf{r}_t$, respectively. Using these coordinates, we can state the unilateral contact law using the Signorini's conditions (see Figure 1 for a graphical representation):

$$ u_n \ge 0, \quad r_n \ge 0, \quad u_n r_n = 0 \qquad (7) $$

On the other hand, the Coulomb's law of friction can be stated using the algorithmic form (see Figure 2 for a graphical representation):

$$ \left[ \begin{array}{lll} \text{If } r_n = 0 & \text{then } u_n \ge 0 & \text{! No contact} \\ \text{Else if } r_n > 0 \text{ and } \|\mathbf{r}_t\| < \mu r_n & \text{then } \mathbf{u} = 0 & \text{! Sticking} \\ \text{Else } r_n > 0 \text{ and } \|\mathbf{r}_t\| = \mu r_n & \text{then } \exists \lambda \ge 0 \text{ such that } \mathbf{u}_t = \lambda \frac{\mathbf{r}_t}{\|\mathbf{r}_t\|} & \text{! Sliding} \end{array} \right. \qquad (8) $$

<figure><img src="image_2.png" /><figcaption>Figure 1. The Signorini conditions.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 2. The Coulomb’s conditions.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3. The Coulomb’s cone.</figcaption></figure>

For a given friction coefficient $\mu$, let $K_\mu$ be the isotropic Coulomb’s cone, which defines the set of admissible forces (Figure 3):

$$K_{\mu} = \{ \mathbf{r} = r_n \mathbf{n} + \mathbf{r}_t : \|\mathbf{r}_t\| - \mu r_n \le 0 \} \quad (9)$$

The previous law can also be written as follows:

$$ \left[ \begin{array}{lll} \text{If } r_n = 0 & \text{then } u_n \ge 0 & ! \text{ No contact} \\ \text{Else if } \mathbf{r} \in I(K_\mu) & \text{then } \mathbf{u} = 0 & ! \text{ Sticking} \\ \text{Else } r_n > 0 \text{ and } \mathbf{r} \in B(K_\mu) & \text{then } \exists \lambda \ge 0 \text{ such that } \mathbf{u}_t = \lambda \frac{\mathbf{r}_t}{\|\mathbf{r}_t\|} & ! \text{ Sliding} \end{array} \right. \quad (10) $$

where $I(K_\mu)$ and $B(K_\mu)$ are, respectively, the interior and the boundary of the cone $K_\mu$.

### 3. NUMERICAL RESOLUTION OF THE CONTACT/FRICTION PROBLEMS

I will describe in this section the numerical algorithms that will be considered in what follows. Generally, to solve problem (6), the numerical algorithms considered are based on two levels: the global level where the equations of motion are solved and the local level devoted to the resolution of the contact law.

#### 3.1. Resolution of the global problem : the nonlinear Gauss-Seidel (NLGS) method

In this paragraph, I describe the algorithm used at the global level to solve problem (6). Following the ideas of Jean and Moreau [7, 8, 12], I use the NLGS algorithm, which is the most commonly used. It consists in considering successively each contact until the convergence. The numerical criterion used to state the convergence will be studied later in this article.
This method is intrinsically sequential, but it is possible to use a simple multithreading technique that consists in splitting the contact loop onto several threads. This method has been studied in [13] in the case of a local algorithm based on the augmented Lagrangian method.

Notice that it is also possible to consider at this stage more sophisticated methods such as conjugate gradient-type methods (see, e.g., [11]). However, these techniques do not improve significantly the speedup of the convergence. That is why, they are not considered in the following.

## 3.2. The standard bipotential-based method (SBP)

In this paragraph, I provide a first method to solve the contact problem at the local level (contact point between two particles). The method is based on the notion of bipotential, introduced by de Saxcé et al. [1].

Using the bipotential framework, we can easily show (see, e.g., [1, 14–16]) that a couple $(\mathbf{u}, \mathbf{r})$ verifies the Signorini–Coulomb contact rules if and only if

$$b_c(\mathbf{v}, \mathbf{s}) + \mathbf{v} \cdot \mathbf{s} \geq b_c(\mathbf{u}, \mathbf{r}) + \mathbf{u} \cdot \mathbf{r} = 0 \quad \forall \mathbf{v}, \mathbf{s} \quad (11)$$

where $b_c$ is the bipotential

$$b_c(-\mathbf{u}, \mathbf{r}) = \Psi_{\mathbb{R}+}(u_n) + \Psi_{K_{\mu}}(\mathbf{r}) + \mu r_n \|u_t\| \quad (12)$$

and $\Psi_C$ stands for the indicator function of the set $C$: $\Psi_C(x) = 0$ if $x \in C$, $\Psi_C(x) = +\infty$ if $x \notin C$.

Consequently, the contact law can be written in a compact form of an implicit subnormality rule (or a differential inclusion rule):

$$-\mathbf{u} \in \partial_{\mathbf{r}} b_c(-\mathbf{u}, \mathbf{r}) \quad (13)$$

Then, for a contact $c$, at an NLGS iteration $i$, knowing the relative velocity $\tilde{\mathbf{u}}^{c,i}$, the algorithm to compute $\mathbf{r}^{c,i+1}$ from $\mathbf{r}^{c,i}$ is based on the minimization of the bipotential (see, e.g., [14]), using the inequality

$$b_c(-\tilde{\mathbf{u}}^{c,i}, \mathbf{r}) + \tilde{\mathbf{u}}^{c,i} \cdot \mathbf{r} \geq b_c(-\tilde{\mathbf{u}}^{c,i}, \mathbf{r}^{c,i+1}) + \tilde{\mathbf{u}}^{k,i} \cdot \mathbf{r}^{c,i+1} \quad \forall \mathbf{r} \in K_{\mu} \quad (14)$$

or $g(\mathbf{r}) \geq g(\mathbf{r}^{c,i+1})$, $\forall \mathbf{r} \in K_{\mu}$, if we denote

$$g(\mathbf{r}) = \Psi_{\mathbb{R}+}(\tilde{u}_n^{c,i}) + \Psi_{K_{\mu}}(\mathbf{r}) + \mu r_n \|\tilde{u}_t^{c,i}\| + \tilde{u}^{c,i} \cdot \mathbf{r} \quad (15)$$

The minimization of (14) is classically realized using a projected gradient method (Uzawa method) without considering the singular term $\Psi_{\mathbb{R}+}(\tilde{u}_n^{c,i})$. This minimization can also be viewed as the search of the proximal point of the augmented force $\mathbf{r} - \rho\tilde{\mathbf{u}}$, with respect to the function $\mathbf{r} \mapsto \rho b_c(-\tilde{\mathbf{u}}, \mathbf{r})$ (see, e.g., [1, 14, 15]):

$$\mathbf{r} = \text{prox}(\mathbf{r} - \rho \tilde{\mathbf{u}}, \rho b_c(-\tilde{\mathbf{u}}, \mathbf{r}))$$

More precisely, the Uzawa method leads to compute the augmented force $\tau^{c,i+1} = \mathbf{r}^{c,i} - \rho \nabla \tilde{g}(\mathbf{r}^{c,i})$, where $\tilde{g}$ is the differential part of $g$:

$$\nabla \tilde{g}(\mathbf{r}^{c,i}) = \nabla_{\mathbf{r}} (\mu r_n \|\tilde{u}_t^{c,i}\| + \tilde{u}^{c,i} \cdot \mathbf{r}) = \mu \|\tilde{u}_t^{c,i}\| \mathbf{n} + \tilde{u}^{c,i}$$

and to consider the force at the next step as a projection of the augmented force onto the set of admissible force $\mathbf{r}^{c,i+1} = \text{proj}(\tau^{c,i+1}, K_{\mu})$, which provides Equations (21) and (22) in the resolution algorithm of the global problem presented in the following paragraphs. The proj($\tau^{c,i+1}$, $K_{\mu}$) stands for the orthogonal projection over the convex $K_{\mu}$, which can be computed exactly (see [14]).

This algorithm will be referred as the SBP method thereafter.

For the sake of simplicity, we denote hereafter the descent direction

$$\mathbf{D}^{c,i} = \mu \|\tilde{\mathbf{u}}_t^{c,i}\| \mathbf{n} + \tilde{\mathbf{u}}^{c,i}$$
Remark 1

A first improvement of this method could be to compute the optimal step $\rho^{c,i}$. To do so, we have to minimize

$$ \rho \mapsto g(\mathbf{r}^{c,i} - \rho \mathbf{D}^{c,i}) \qquad (16) $$

or, more precisely,

$$ \begin{aligned} \rho &\mapsto \Psi_{\mathbb{R}+} (\tilde{u}_n^{c,i}) + \Psi_{K_\mu} (\mathbf{r}^{c,i} - \rho \mathbf{D}^{c,i}) + \mu (r_n^{c,i} - \rho \mathbf{D}^{c,i} \cdot \mathbf{n}) \|\tilde{\mathbf{u}}_t^{c,i}\| + \tilde{\mathbf{u}}^{c,i} \cdot (\mathbf{r}^{c,i} - \rho \mathbf{D}^{c,i}) \\ &= \Psi_{\mathbb{R}+} (\tilde{u}_n^{c,i}) + \Psi_{K_\mu} (\mathbf{r}^{c,i} - \rho \mathbf{D}^{c,i}) - \rho \mathbf{D}^{c,i} \cdot (\mu \|\tilde{\mathbf{u}}_t^{c,i}\| \mathbf{n} + \tilde{\mathbf{u}}^{c,i}) + Cte \\ &= \Psi_{\mathbb{R}+} (\tilde{u}_n^{c,i}) + \Psi_{K_\mu} (\mathbf{r}^{c,i} - \rho \mathbf{D}^{c,i}) - \rho \|\mathbf{D}^{c,i}\|^2 + Cte \end{aligned} \qquad (17) $$

We can observe that this method does not permit to choose an optimal parameter $\rho$ because $g$, as a function of $\rho$, is linear, excepted in the case where it exists $\rho > 0$ such that $\mathbf{r}^{c,i} - \rho \mathbf{D}^{c,i} \notin K_\mu$. A solution could be to modify the function $g$, for example, by replacing $\tilde{\mathbf{u}}^{c,i}$ by a prediction of $\tilde{\mathbf{u}}^{c,i+1}$ using the equations of the dynamics. Unfortunately, this method does not provide satisfying numerical results.

Then, the SBP algorithm can be written (see, e.g., [16]):

* Loop on the step time $k$

- Prediction of a position (for the computation of the local-global mapping):

$$ \mathbf{q}_{k+\frac{1}{2}} = \mathbf{q}_k + \frac{\Delta t}{2} \dot{\mathbf{q}}_k \qquad (18) $$

- Initialization of the motion: $\dot{\mathbf{q}}_{k+1}^0 = \dot{\mathbf{q}}_k^{\text{free}}$ (initialization of the contact forces with $\mathbf{R} = 0$).

- Loop on $i \ge 0$ (NLGS), until convergence

* Loop on the contacts $c$:

- Computation of the local-global mapping

$$ \dot{\mathbf{u}}^{-} = P^{*} (\mathbf{q}_{k+\frac{1}{2}}, c) \dot{\mathbf{q}}_{k} : \qquad \dot{\mathbf{u}}^{c,+i} = P^{t} (\mathbf{q}_{k+\frac{1}{2}}, c) \dot{\mathbf{q}}_{k+1}^{i} \qquad (19) $$

- Newton shock law

$$ \tilde{u}_n^{c,i} = \frac{u_n^{c,+i} + e_n \bar{u}_n^{-}}{1 + e_n} ; \qquad \tilde{u}_t^{c,i} = \frac{\mathbf{u}_t^{c,+i} + e_n \bar{\mathbf{u}}_t^{-}}{1 + e_t} \qquad (20) $$

- Prediction of the reaction:

$$ \tau^{c,i+1} = \mathbf{r}^{c,i} - \rho \left[ \tilde{\mathbf{u}}_t^{c,i} + \left( \tilde{u}_n^{c,i} + \mu \|\tilde{\mathbf{u}}_t^{c,i}\| \right) \mathbf{n} \right] \qquad (21) $$

- Correction of the reaction:

$$ \mathbf{r}^{c,i+1} = \text{proj}(\tau^{c,i+1}, K_\mu) \qquad (22) $$

- Actualization of the generalized displacement:

$$ \dot{\mathbf{q}}_{k+1}^{i+1} = \dot{\mathbf{q}}_{k}^{\text{free}} + \mathbb{M}^{-1} \left( \sum_{\alpha \le c} P(\mathbf{q}_{k+\frac{1}{2}}, \alpha) \mathbf{r}^{\alpha,i+1} + \sum_{\alpha > c} P(\mathbf{q}_{k+\frac{1}{2}}, \alpha) \mathbf{r}^{\alpha,i} \right) \qquad (23) $$

* End of the loop on contacts $c$.

- End of the loop on $i$ of NLGS. When the convergence is reached, actualization of the velocity $\dot{\mathbf{q}}_{k+1} = \dot{\mathbf{q}}_{k+1}^{i+1}$

- Actualization of the generalized displacements: $\mathbf{q}_{k+1} = \mathbf{q}_{k+\frac{1}{2}} + \frac{\Delta t}{2} \dot{\mathbf{q}}_{k+1}$

* End of the loop on the step time $k$.
### Remark 2

Notice that only one iteration of the Uzawa algorithm at the local level is considered. Various previous studies (see, e.g., [17]) show that there is no significant improvement of the method if several iterations of the Uzawa algorithm are considered at this stage.

### 3.3. Newton method and enhanced bipotential method (EBP)

I introduce in this section a Newton method to speed up the convergence in the computation of the solution. This method has been already used, especially in the case of the augmented Lagrangian method developed by Alart et al. [2], and the ideas presented in this article follow those of Feng et al. [3] and have been adapted to the problem of the discrete element method. The main idea of this technique is to find the solution of the optimization problem, not as a minimum but as a zero of a function, using the Euler equation (or the saddle point) of the problem. Then, a standard Newton method can be developed to solve this Euler equation.

The technique is first described in the case of the bipotential framework and will be adapted to the augmented Lagrangian method farther.

We recall that the local problem to be solved, for each contact c, can be written as

$$ egin{cases} 
\tilde{\mathbf{u}}_{k+1}^c = \tilde{\mathbf{u}}_k^{c, \text{free}} + \sum_{\alpha=1}^{N_c} W_{c\alpha} \mathbf{r}^{\alpha} \\ 
\mathbf{r}^c = \text{proj}(\tau^c, K_\mu) 
\end{cases} 
\quad \forall c = 1, \dots, N_c \qquad (24) $$

where $\tau^c = \mathbf{r}^c - \rho (\mu \|\tilde{\mathbf{u}}_t^c\|\mathbf{n} + \tilde{\mathbf{u}})$ is the augmented reaction (see (21)), and $W_{c\alpha} = P^*(\mathbf{q}_{k+\frac{1}{2},c}) \mathbb{M}^{-1} P(\mathbf{q}_{k+\frac{1}{2},\alpha})$ is the local Delassus operator.

This problem can be written equivalently

$$ egin{cases} 
\tilde{\mathbf{u}}_{k+1}^c - \tilde{\mathbf{u}}_k^{c, \text{free}} - \sum_{\alpha=1}^{N_c} W_{c\alpha} \mathbf{r}^{\alpha} = 0 \\ 
\mathbf{r}^c - \text{proj}(\tau^c, K_\mu) = 0 
\end{cases} 
\quad \forall c = 1, \dots, N_c \qquad (25) $$

Reminding now that we want to use a Newton algorithm to solve these equations inside the NLGS loop on the variable $i$, we define now, for each contact $c = 1, \dots, N_c$, the function

$$ f_c^i(\chi) = \begin{pmatrix} \tilde{\mathbf{u}}^{c,i} - \tilde{\mathbf{u}}_k^{c, \text{free}} - \sum_{\alpha=1}^{N_c} W_{c\alpha} \mathbf{r}^{\alpha,i} \\ \mathbf{Z}^{c,i} \end{pmatrix} $$

where

* the vector $\mathbf{Z}^c$ is the error on the prediction of the reaction

$$ \mathbf{Z}^{c,i}(\mathbf{r}^{c,i}, \tilde{\mathbf{u}}^{c,i}) = \mathbf{r}^{c,i} - \text{proj}(\tau^{c,i}, K_\mu) \qquad (26) $$

* $\chi_c = (\mathbf{r}^{c,i}, \tilde{\mathbf{u}}^{c,i})^t$

* $\chi = (\chi_1, \chi_2, \dots, \chi_{N_c})^t$

### Remark 3

The first equality in the relation $f(\chi) = 0$ is the equation of motion for the bodies in contact, and the second relation is the frictional Coulomb law between the bodies in contact, written within the bipotential framework.
Then, we have to write a Newton algorithm to solve the problem $f(\chi) = 0$. This algorithm can be written, for a contact c, by substituting Equations (21) and (22) in algorithm (SBP) by the following:

* Initialization:

$$ \chi_c^0 = (\mathbf{r}^0 = \mathbf{r}^{c,i}, \mathbf{v}^0 = \tilde{\mathbf{u}}^{c,i})^t, \quad \ell = 0 $$

* Loop on $\ell$, until convergence:

$$ - \tau_{\ell}^{c} = \mathbf{r}^{\ell} - \rho (\mu \|\mathbf{v}_{t}^{\ell} \|\mathbf{n} + \mathbf{v}^{\ell}) $$

- Resolution:

$$ \left[ \frac{\partial f_c}{\partial \chi^c} (\chi^\ell) \right] \Delta \chi_c = -f_c(\chi^\ell) \qquad (27) $$

- Actualization: $\chi_c^{\ell+1} = \chi_c^{\ell} + \Delta\chi_c$

* End of the loop on $\ell$ until convergence, $\tilde{\mathbf{u}}^{c,i+1} = \mathbf{v}^{\ell}$ and $\mathbf{r}^{c,i+1} = \mathbf{r}^{\ell}$.

### Remark 4

This algorithm needs more than one iteration at each NLGS iteration to be efficient. As a consequence and compared with the Uzawa algorithm, the solution in the Newton algorithm is controlled by both the local (iteration $\ell$) and global convergence criteria (iteration $i$, see [3, 17]).

The local convergence criterion for the Newton algorithm is defined by

$$ \epsilon_{\text{Newt}}^{c}(\chi_{\ell}) = \|\mathbf{v}^{\ell} - \mathbf{u}_{k}^{c, \text{free}} - W\mathbf{r}^{\ell}\| + \|\mathbf{r}^{\ell} - \text{proj}(\mathbf{r}^{\ell}, K_{\mu})\| \qquad (28) $$

This criterion measures $f_c(\chi_\ell)$ that has to be sufficiently small.

The matrix $[\partial f_c / \partial \chi_c(\chi)]$ represents the tangential matrix of the local equations for the contact c. This matrix is of dimension $6 \times 6$ for a 3D problem and $4 \times 4$ for a 2D problem. For a 3D problem, the general form of this matrix is the following:

$$ \left[ \frac{\partial f_c}{\partial \chi_c}(\chi) \right] = \begin{bmatrix} -W & Id_{3 \times 3} \\ A_c & B_c \end{bmatrix} \qquad (29) $$

where

$$ A_c = \left[ \frac{\partial Z_c}{\partial r_n} \middle| \frac{\partial Z_c}{\partial r_{t_1}} \middle| \frac{\partial Z_c}{\partial r_{t_2}} \right] \quad B_c = \left[ \frac{\partial Z_c}{\partial v_n} \middle| \frac{\partial Z_c}{\partial v_{t_1}} \middle| \frac{\partial Z_c}{\partial v_{t_2}} \right] \qquad (30) $$

The matrices $A_c$ and $B_c$ take different forms according to the contact status:

* First case: sliding contact.

In that case, we have

$$ \mu \|\tau_t\| \geq -\tau_n \qquad \|\tau_t\| \geq \mu \tau_n $$

then

$$ \text{Proj}(\tau, K_{\mu}) = \tau - \left( \frac{\|\tau_t\| - \mu \tau_n}{1 + \mu^2} \right) \left( \frac{\tau_t}{\|\tau_t\|} - \mu \mathbf{n} \right) $$

and

$$ \mathbf{Z}_c = \rho (\mu \|\mathbf{v}_t^k\|\mathbf{n} + \mathbf{v}^k) + \left( \frac{\|\tau_t\| - \mu \tau_n}{1 + \mu^2} \right) \left( \frac{\tau_t}{\|\tau_t\|} - \mu \mathbf{n} \right) $$
The computation of the derivatives of $\mathbf{Z}_c$ provides the matrices $A_c$ and $B_c$:

$$
\begin{aligned}
-\frac{\partial \mathbf{Z}_c}{\partial r_n} &= -\frac{\mu}{1+\mu^2} \left( \frac{\tau_t}{\|\tau_t\|} - \mu \mathbf{n} \right) \\
-\frac{\partial \mathbf{Z}_c}{\partial r_{t_1}} &= \frac{\tau_{t_1}}{(1+\mu^2)\|\tau_t\|} \left( \frac{\tau_t}{\|\tau_t\|} - \mu \mathbf{n} \right) + \frac{\|\tau_t\| - \mu \tau_n}{1+\mu^2} \left( \frac{\mathbf{t}_1}{\|\tau_t\|} - \frac{\tau_{t_1}}{\|\tau_t\|^3} \tau_t \right) \\
-\frac{\partial \mathbf{Z}_c}{\partial r_{t_2}} &= \frac{\tau_{t_2}}{(1+\mu^2)\|\tau_t\|} \left( \frac{\tau_t}{\|\tau_t\|} - \mu \mathbf{n} \right) + \frac{\|\tau_t\| - \mu \tau_n}{1+\mu^2} \left( \frac{\mathbf{t}_2}{\|\tau_t\|} - \frac{\tau_{t_2}}{\|\tau_t\|^3} \tau_t \right) \\
-\frac{\partial \mathbf{Z}_c}{\partial v_n} &= \rho \mathbf{n} + \frac{\rho\mu}{1+\mu^2} \left( \frac{\tau_t}{\|\tau_t\|} - \mu \mathbf{n} \right) \\
-\frac{\partial \mathbf{Z}_c}{\partial v_{t_1}} &= \rho \left( \mathbf{t}_1 + \mu \frac{v_{t_1}}{\|\mathbf{v}_t\|} \mathbf{n} \right) - \frac{\rho}{1+\mu^2} \left( \begin{pmatrix} \frac{\tau_{t_1}}{\|\tau_t\|} - \frac{\mu^2 v_{t_1}}{\|\mathbf{v}_t\|} \\ (\|\tau_t\| - \mu \tau_n) \left( \frac{\mathbf{t}_1}{\|\tau_t\|} - \frac{\tau_{t_1}}{\|\tau_t\|^3} \tau_t \right) \end{pmatrix} \right) \\
-\frac{\partial \mathbf{Z}_c}{\partial v_{t_2}} &= \rho \left( \mathbf{t}_2 + \mu \frac{v_{t_2}}{\|\mathbf{v}_t\|} \mathbf{n} \right) - \frac{\rho}{1+\mu^2} \left( \begin{pmatrix} \frac{\tau_{t_2}}{\|\tau_t\|} - \frac{\mu^2 v_{t_2}}{\|\mathbf{v}_t\|} \\ (\|\tau_t\| - \mu \tau_n) \left( \frac{\mathbf{t}_2}{\|\tau_t\|} - \frac{\tau_{t_2}}{\|\tau_t\|^3} \tau_t \right) \end{pmatrix} \right)
\end{aligned}
$$

For a 2D problem, these computations yield the following:

$$
\begin{aligned}
-\frac{\partial \mathbf{Z}_c}{\partial r_n} &= \frac{\mu}{1+\mu^2} (\mu \mathbf{n} - \theta_r \mathbf{t}) \\
-\frac{\partial \mathbf{Z}_c}{\partial r_t} &= \frac{1}{(1+\mu^2)} (-\mu \theta_r \mathbf{n} + \mathbf{t}) \\
-\frac{\partial \mathbf{Z}_c}{\partial v_n} &= \frac{\rho}{1+\mu^2} (\mathbf{n} + \mu \theta_r \mathbf{t}) \\
-\frac{\partial \mathbf{Z}_c}{\partial v_t} &= \frac{\rho \mu}{1+\mu^2} ((\theta_v + \theta_r)\mathbf{n} + \mu(1-\theta_r\theta_v)\mathbf{t})
\end{aligned}
$$

where $\theta_v = \text{sign}(\mathbf{v}_t)$ and $\theta_r = \text{sign}(\tau_t)$.

• Second case: sticking contact.

In that case, we have

$$ \mu \|\tau_t\| \ge -\tau_n \qquad \|\tau_t\| < \mu \tau_n $$

then

$$ \mathbf{Z}_c = \rho (\mu \|\mathbf{v}_t^k\| \mathbf{n} + \mathbf{v}^k) $$

and the computation of the derivatives of $\mathbf{Z}_c$ reads

$$
\begin{aligned}
-A_c &= 0_{3 \times 3} \\
-\frac{\partial \mathbf{Z}_c}{\partial v_n} &= \rho \mathbf{n} \\
-\frac{\partial \mathbf{Z}_c}{\partial v_{t_1}} &= \rho \mu \frac{v_{t_1}}{\|\mathbf{v}_t\|} \mathbf{n} + \rho \mathbf{t}_1 \\
-\frac{\partial \mathbf{Z}_c}{\partial v_{t_2}} &= \rho \mu \frac{v_{t_2}}{\|\mathbf{v}_t\|} \mathbf{n} + \rho \mathbf{t}_2
\end{aligned}
$$
For a 2D problem, these computations lead to the following:

$$
\begin{align*}
-A_c &= 0_{2 \times 2} \\
-\frac{\partial \mathbf{Z}_c}{\partial v_n} &= \rho \mathbf{n} \\
-\frac{\partial \mathbf{Z}_c}{\partial v_t} &= \rho \mu \theta_v \mathbf{n} + \rho \mathbf{t}
\end{align*}
$$

• Third case: no contact.

In that case, the matrices $A_c = \text{Id}_{3 \times 3}$ and $B_c$ vanish, and $\chi_c^{\ell+1} = \begin{cases} 0 \\ \mathbf{v}^k \end{cases}$

3.4. Resolution of the linear system

Generally, the drawback of a Newton method is the computational cost of the linear system to be solved at each iteration. Here, the particular form of the tangent matrix allows the use of a condensation technique. More precisely, the linear system to be solved can be written as

$$
\begin{bmatrix} -W & Id_{3 \times 3} \\ A_c & B_c \end{bmatrix} \begin{pmatrix} \delta \mathbf{r} \\ \delta \mathbf{v} \end{pmatrix} = \begin{pmatrix} -\mathbf{f} \\ -\mathbf{g} \end{pmatrix} \qquad (31)
$$

The first equation yields $\delta \mathbf{v} = -\mathbf{f} + W \delta \mathbf{r}$, and introducing this equality in the second equation leads to solve the linear system

$$
(A_c + B_c W) \delta \mathbf{r} = -\mathbf{g} + B_c \mathbf{f} \quad (32)
$$

This property halves the size of the linear system to be solved.

Remark 5

A drawback of the bipotential framework is that, due to its specificity, it is rather difficult to con-
sider fully coupled problems, where the contact law and another phenomena, such as electricity or
thermal effects, are strongly coupled. The other method presented in this article has a better prop-
erty from this point of view because it is based on a more standard mathematical background in the
theory of optimization (see, e.g., [18]).

3.5. Standard augmented Lagrangian (SAL) and enhanced augmented Lagrangian (EAL) method

In [2], Alart *et al.* proposed another method to solve the frictional contact problem. This method has been also used with various improvement (e.g., parallelization and conjugate gradient method) to solve multicontact problems (see, e.g., [11, 13, 19–22]). Even if the coupled frictional contact problem is not an optimization problem anymore, it is always possible to formally formulate a ‘quasi’-optimization problem, for which the constraint set depends on the normal components of the solution as a parameter. The solution is then searched as a saddle point of a quasi-augmented Lagrangian of the problem.

More precisely, the global problem on all unknowns that has to be solved at each time step
(in place of Equation (24)) has the following form:

$$
\left\{
\begin{array}{l}
\mathbf{u} = \mathbf{u}^{\text{free}} + W\mathbf{r} \\
\mathbf{r} \ge 0, \mathbf{u} \ge 0, \mathbf{r} \cdot \mathbf{u} = 0
\end{array}
\right.
\qquad (33)
$$

To solve this problem, for a given **r** ∈ ℝ³×N_c, one can define the Cartesian product of infinite
half-cylinder with section equal to the ball B(0, μr_c) of radius μr_c by:

$$
C(\mu \mathbf{r}) = \prod_{c=1}^{N_c} \mathbb{R}^{+} \times \mathcal{B}(0, \mu r_c)
$$
and then, the granular-type frictional contact problem is given by

$$ \mathbf{r} \in \operatorname{argmin}_{\mathbf{r} \in C(\mu\mathbf{r})} \frac{1}{2} \mathbf{r} \cdot \mathbf{W}\mathbf{r} + \mathbf{u}^{\text{free}} \cdot \mathbf{r} = \operatorname{argmin}_{\mathbf{r} \in C(\mu\mathbf{r})} J(\mathbf{r}) \quad (34) $$

and the projected gradient method to minimize this problem reads (for each iteration $i$ of the NLGS algorithm) as follows:

$$ \mathbf{r}^{i+1} = \operatorname{proj}(\mathbf{r}^i - \rho(\mathbf{u}^{\text{free}} + \mathbf{W}\mathbf{r}^i), C(\mu\mathbf{r}^{i+1})) \quad (35) $$

or $\mathbf{r}^{i+1} = \operatorname{proj}(\tau^{i+1}, C(\mu\mathbf{r}^{i+1}))$, with $\tau^{i+1} = \mathbf{r}^i - \rho\mathbf{u}^i$, $\mathbf{u}^i = \mathbf{u}^{\text{free}} + \mathbf{W}\mathbf{r}^i$. This algorithm will be referred to hereinafter as the SAL method.

Notice that this method is very close to the SBP method. More precisely, for a contact $c$, only the descent direction $\tilde{\mathbf{u}}^{c,i} + \mu \|\tilde{\mathbf{u}}^{c,i}\| \mathbf{n}$ in (21) is replaced by $\tilde{\mathbf{u}}^{c,i}$, and the projection $\mathbf{r}^{c,i+1} = \operatorname{proj}(\tau^{c,i+1}, K_\mu)$ in (22) is replaced by

$$ \left\{ \begin{array}{l} r_n^{c,i+1} = \max(0, \tau_n^{c,i+1}) \\ \\ \mathbf{r}_t^{c,i+1} = \frac{\tau_t^{c,i+1}}{\|\tau_t^{c,i+1}\|} \min(\mu r_n^{c,i+1}, \|\tau_t^{c,i+1}\|) \end{array} \right. $$

### Remark 6

It is also possible to see the algorithm developed from the bipotential formalism as a slight modification of the algorithm defined previously. Indeed, it is only necessary to change the set $C(\mathbf{r})$ by

$\mathcal{K} = \prod_{c=1}^{N_c} K_\mu$, and to change the descent direction $\tilde{\mathbf{u}}^{c,i}$ by $\tilde{\mathbf{u}}^{c,i} + \mu \|\tilde{\mathbf{u}}^{c,i}\| \mathbf{n}$, which remains a descent direction for the SAL method, because

$$ \nabla J(\mathbf{r}^{c,i+1}) \cdot \mathbf{D}^{c,i} = -\|\tilde{\mathbf{u}}^{c,i}\|^2 - \tilde{\mathbf{u}}^{c,i} \cdot (\mu \|\tilde{\mathbf{u}}^{c,i}\|) \mathbf{n} = -\|\tilde{\mathbf{u}}^{c,i}\|^2 - \mu u_n^{c,i} \|\tilde{\mathbf{u}}^{c,i}\| $$

which is negative because $\mu \ge 0$.

Then, acting by analogy, we can develop a Newton method to find the minimum of $J$ by seeking the solution as a zero of the function $\tilde{f}(\chi)$ where, for a contact $c$,

$$ \tilde{f}_c(\chi) = \left( \begin{array}{l} \tilde{\mathbf{u}}_{k+1}^c - \tilde{\mathbf{u}}_k^{c,\text{free}} - \sum_{\alpha=1}^{N_c} W_{c\alpha} \mathbf{r}^\alpha \\ \\ \frac{\tilde{\mathbf{Z}}^c}{\tilde{\mathbf{Z}}^c} \end{array} \right) $$

the vector $\tilde{\mathbf{Z}}^c$ is the error on the prediction of the reaction

$$ \tilde{\mathbf{Z}}^c(\mathbf{r}^c, \tilde{\mathbf{u}}_{k+1}^c) = \mathbf{r}^c - \operatorname{proj}(\tau_{k+1}^c, C_c(\mu\mathbf{r}^c)) \quad (36) $$

and the set $C_c(\mu\mathbf{r}^c)$ is the set of admissible forces $C_c(\mu\mathbf{r}^c) = \mathbb{R}^+ \times \mathcal{B}(0, r_c)$. This method will be referred as the EAL method hereafter.

Then, in what follows, we have three cases in the computation of the tangent matrix $[\partial \tilde{f}/\partial \chi^c(\chi^l)]$:

* First case: sliding contact ($\tau_n > 0$, $\tau_t \ge \mu r_n$)

We have $\operatorname{proj}(\tau^c, C_c(\mu\mathbf{r}^c)) = \tau_n \mathbf{n} + \frac{\tau_t}{\|\tau_t\|} \mu r_n$ and $\tilde{\mathbf{Z}}_c = \rho v_n \mathbf{n} - \frac{\tau_t}{\|\tau_t\|} \mu r_n + \mathbf{r}_t$.
The computation of the derivatives of $\tilde{Z}_c$ provides the matrices $A_c$ and $B_c$:

$$
\begin{align*}
-\frac{\partial \tilde{Z}_c}{\partial r_n} &= -\mu \frac{\tau_t}{\|\tau_t\|} \\
-\frac{\partial \tilde{Z}_c}{\partial r_{t_1}} &= \mathbf{t}_1 - \mu r_n \left( \frac{\mathbf{t}_1}{\|\tau_t\|} - \frac{\tau_{t_1}}{\|\tau_t\|^3} \tau_t \right) \\
-\frac{\partial \tilde{Z}_c}{\partial r_{t_2}} &= \mathbf{t}_2 - \mu r_n \left( \frac{\mathbf{t}_2}{\|\tau_t\|} - \frac{\tau_{t_2}}{\|\tau_t\|^3} \tau_t \right) \\
-\frac{\partial \tilde{Z}_c}{\partial v_n} &= \rho \mathbf{n} \\
-\frac{\partial \tilde{Z}_c}{\partial v_{t_1}} &= -\rho \mu r_n \left( \frac{\mathbf{t}_1}{\|\tau_t\|} - \frac{\tau_{t_1}}{\|\tau_t\|^3} \tau_t \right) \\
-\frac{\partial \tilde{Z}_c}{\partial v_{t_2}} &= -\rho \mu r_n \left( \frac{\mathbf{t}_2}{\|\tau_t\|} - \frac{\tau_{t_2}}{\|\tau_t\|^3} \tau_t \right)
\end{align*}
$$

For a 2D problem, these computations lead to

$$
A_c = \begin{pmatrix} 0 & 0 \\ -\mu\theta_r & 1 \end{pmatrix} \quad B_c = \rho \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix}
$$

* Second case: sticking contact ($\tau_n > 0$, $\tau_t < \mu r_n$)
  $\mathrm{proj}(\tau^c, C_c(\mu\mathbf{r}^c)) = \tau^c$, and the computation of the derivatives of $Z_c$ reads

$$
-A_c = 0_{3 \times 3}
$$

$$
-B_c = \rho \mathrm{Id}_{3 \times 3}
$$

* Third case: no contact ($\tau_n \le 0$)
  $\mathrm{proj}(\tau^c, C_c(\mu\mathbf{r}^c)) = 0$, then the matrices $A_c = \mathrm{Id}_{3 \times 3}$ and $B_c$ vanishes, and $\chi_c^{\ell+1} = \begin{Bmatrix} 0 \\ \mathbf{v}^k \end{Bmatrix}$

3.6. The global stopping (convergence) criterion

I present in this paragraph the convergence criterion on the global NLGS iterations. This criterion,
developed from that proposed in [23], has been extended in the case of the Newton and bipoten-
tial (EBP) method, where some terms are naturally vanishing in the original Uzawa and bipotential
(SBP) method. This criterion $\epsilon_{glob}$ has been written in such a way that if the solution verifies that
$\epsilon_{glob}$ is sufficiently small, then this solution provides good properties of the equation of motion and
Signorini–Coulomb contact law. Consequently, this criterion remains valid for methods developed
using an augmented Lagrangian (e.g., SAL and EAL methods).

This criterion can be stated as:

$$
\epsilon_{\text{glob}} = \frac{1}{N_c} \sum_{c=1}^{N_c} \left[ \epsilon_{\text{motion}}^{c} + \epsilon_{\text{proj}}^{c} + \epsilon_{b_c} + \epsilon_{\text{pen}}^{c} \right] \quad (37)
$$

where

* $\epsilon_{\text{motion}}^c = \|\tilde{\mathbf{u}}^c - \tilde{\mathbf{u}}_m^c\|$ where $\tilde{\mathbf{u}}_m^c = \tilde{\mathbf{u}}^{c,i} + \sum_{\alpha=1}^{N_c} W_{c\alpha} \mathbf{r}^{\alpha}$, so $\epsilon_{\text{motion}}$ measures the error on the equation of motion (see Equation (24), this term vanishes for the SBP and SAL methods);

* $\epsilon_{\text{proj}}^c = \|\mathbf{r}^c - \text{proj}(\mathbf{r}^c, K_\mu)\|$ is the error for the projection onto the Coulomb cone (vanishing for the SBP method);

* $\epsilon_{b_c} = |\tilde{\mathbf{u}}^c \cdot \mathbf{r}^c + \mu r_n^c ||\tilde{\mathbf{u}}_t^c||$ is the absolute value of the bi-potential that has to vanish if and only if the couple $(\tilde{\mathbf{u}}^c, \mathbf{r}^c)$ verifies the Signorini Coulomb contact law (see formula 11);

* $\epsilon_{\text{pen}}^c = -\min(0, \tilde{u}_n^c)$ is the value of the penetration.
### Remark 7

One can notice that it is absolutely necessary to verify in the criterion that there is no penetration, because nothing in the presented algorithm ensures that this condition is satisfied at the end of the loop. Moreover, if this condition is not satisfied, the other part of the bipotential can be negative or can vanish, even if the couple $(\tilde{u}, \tilde{r})$ is not a solution.

## 4. NUMERICAL RESULTS

I present in this section three numerical examples with an increasing complexity.

In these computations, the descent parameter $\rho$ is taken in such a way that the result is optimal in terms of time computing. Denoting $\bar{\rho} = \frac{m_i m_j}{m_i+m_j} \frac{1}{\Delta t}$, for the SBP and the SAL methods, we have chosen $\rho = 0.6\bar{\rho}$, whereas for the EBP and the EAL methods, it is better the take $\rho = \bar{\rho}$. We recall that it has been shown that, for the bipotential method (see, e.g., [3]) and the augmented Lagrangian method (see, e.g., [11]), the parameter $\rho$ has to verify $\rho < 2\bar{\rho}$ in order to ensure the convergence. Generally, for these two methods, the convergence is very sensitive to this parameter. We will show in the last paragraph of this study that for the EBP method, the parameter $\rho$ can be taken in a large range around the value $\bar{\rho}$ without changing dramatically the convergence of the method.

At each iteration of the NLGS algorithm, the Newton algorithm is stopped either if the convergence is obtained ($\epsilon_{\text{Newt}}^{c} \le 10^{-5}$) or if the number of iterations of the Newton algorithm reaches 100 when there is no convergence.

### 4.1. Ball sliding on a plane

In this first example, we consider a ball placed on a table with an initial horizontal velocity equal to $1.5 \text{ m} \cdot \text{s}^{-1}$. The radius of the ball is equal to $5 \cdot 10^{-3} \text{ m}$, and the friction coefficient between the table and the ball is equal to $\mu = 0.7$. The time step of discretization is equal to $10^{-4} \text{ s}$. In this experiment, the ball first slides on the table, and then, the ball rolls without sliding (Figure 4). The global stopping criterion is equal to $\epsilon_{\text{glob}} = 10^{-10}$.

We can observe from these numerical results that the error coming from the projection is very small for the four methods. The SBP method and the SAL method give very close results, both in terms of quality (Figures 5 and 6) and in terms of time computing (Table I). Nevertheless, we can notice that the time computing is smaller with the SAL method, because there is less computations at each iteration (e.g., there is no term such as $\|\tilde{u}_r\|$ and projection is easier to compute). The EBP method provides better results, both in terms of quality (Figure 7) and in terms of time computing (6.5% better). The EAL method converges after the first NLGS iteration for every time step, and consequently, this is the faster method on this example (7% faster than the SBP method).

### 4.2. Sedimentation of four balls in a box

In this second experiment, we consider the sedimentation of four balls of radii ranging from $4 \cdot 10^{-4}$ to $5 \cdot 10^{-4} \text{ m}$ (Figure 8). For the computations, the time step of discretization is equal to $\Delta t = 10^{-4} \text{ s}$, and the NLGS loop is stopped either if the global stopping criterion on the NLGS method is equal to $\epsilon_{\text{glob}} = 10^{-10}$ or after 5000 iterations if there is no convergence (this case never occurs in this experiment). The friction coefficient between the balls and between the balls and the walls is equal to $\mu = 0.3$.

<figure><img src="image_5.png" /><figcaption>Figure 4. Example 1—A ball is launched with an initial horizontal velocity (left). First, the ball slides. Then, the ball rolls without slipping (right).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5. Example 1—Convergence for the standard bipotential-based (SBP) method, fifth iteration. The last two curves collapse. Thus, the total error is essentially due to the penetrations.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6. Example 1—Convergence for the standard augmented Lagrangian (SAL) method, fifth iteration. The last two curves collapse. Thus, the total error is essentially due to the penetrations.</figcaption></figure>

<table><caption>Table I. Comparison of the results obtained by the four methods on the first example (after the 2000th time step).</caption><thead><tr><th>Method</th><th>Number of NLGS iterations (last time step)</th><th>Error ε<sub>glob</sub> (last time step)</th><th>Total CPU time (s)</th></tr></thead><tbody><tr><td>SBP</td><td>18</td><td>0.384 · 10<sup>-10</sup></td><td>9.44</td></tr><tr><td>SAL</td><td>18</td><td>0.384 · 10<sup>-10</sup></td><td>9.28</td></tr><tr><td>EBP</td><td>1</td><td>0</td><td>8.83</td></tr><tr><td>EAL</td><td>1</td><td>0.175 · 10<sup>-13</sup></td><td>8.78</td></tr></tbody></table>

NLGS, nonlinear Gauss-Seidel; SBP, standard bipotential-based method; SAL, standard augmented Lagrangian; EBP, enhanced bipotential method; EAL, enhanced augmented Lagrangian.
<figure><img src="image_8.png" /><figcaption>Figure 7. Example 1—Convergence for the Newton and enhanced bipotential (EBP) method, fifth iteration. The curve corresponding to the penetrations does not appear because the penetrations are lower than $10^{-14}$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8. Example 2—Sedimentation of four balls under the gravity effect.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9. Example 2—Convergence of the nonlinear Gauss-Seidel iterations for the standard bipotential-based (SBP) method (1000th time step). The last two curves overlap, showing that the global error is governed by the error of penetration.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10. Example 2—Convergence of the nonlinear Gauss-Seidel iterations for the standard augmented Lagrangian (SAL) method (1000th time step). The last two curves overlap, showing that the global error is governed by the error of penetration.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 11. Example 2—Convergence of the nonlinear Gauss-Seidel iterations for the Newton and enhanced bipotential (EBP) method (1000th time step). The last two curves overlap, which shows that the global error is governed by the error on the equations of motion.</figcaption></figure>

<table><caption>Table II. Comparison of the results obtained by the four methods on the second example (after the 1000th time step).</caption><thead><tr><th>Method</th><th>Number of NLGS iterations (last time step)</th><th>Error ε<sub>glob</sub> (last time step)</th><th>Maximal penetration (last time step)</th><th>Total CPU time (s)</th></tr></thead><tbody><tr><td>SBP</td><td>305</td><td>0.949 · 10<sup>-12</sup></td><td>0.310 · 10<sup>-11</sup></td><td>2.92</td></tr><tr><td>SAL</td><td>301</td><td>0.980 · 10<sup>-12</sup></td><td>0.340 · 10<sup>-11</sup></td><td>2.87</td></tr><tr><td>EBP</td><td>161</td><td>0.635 · 10<sup>-12</sup></td><td>0.641 · 10<sup>-12</sup></td><td>2.59</td></tr><tr><td>EAL</td><td>158</td><td>0.973 · 10<sup>-12</sup></td><td>0.208 · 10<sup>-19</sup></td><td>2.43</td></tr></tbody></table>

NLGS, nonlinear Gauss-Seidel; SBP, standard bipotential-based method; SAL, standard augmented Lagrangian; EBP, enhanced bipotential method; EAL, enhanced augmented Lagrangian.
Like in the previous simulations, the SBP and SAL methods provide very similar results (Figures 9 and 10). For these two methods, we can notice that here the global error is essentially due to the penetrations. The SAL method is 2% faster than the SBP method (Table II).

Results obtained by the EBP method are better (Figure 11), and here, the overall error is governed by the error on the equations of motion. The EBP method is 11% faster than the SBP method, and the penetration is five times smaller. In this example, the EAL is the faster method (16.8% faster than the SBP method), and the penetration is very small (Figure 12).

### 4.3. Sedimentation of 500 balls

In this example, we consider the sedimentation of 500 balls (Figure 13) of radii ranging from $2.5 \cdot 10^{-4}$ to $5 \cdot 10^{-4}$ m, the time step of discretization is equal to $\Delta t = 5 \cdot 10^{-5}$ s, and the NLGS loop is stopped if the global estimator (37) verifies $\varepsilon_{glob} \le 10^{-12}$ or after 5000 iterations if there is no convergence. The friction coefficient between the balls, and between the balls and the walls is equal to $\mu = 0.3$.

The results in Table III are obtained after 1000 time steps ($t = 5 \cdot 10^{-2}$ s).

<figure><img src="image_13.png" /><figcaption>Figure 12. Example 2—Convergence of the nonlinear Gauss-Seidel iterations for the Newton and enhanced augmented Lagrangian (EAL) method (1000th time step). The last two curves overlap, and the other one does not appear on the figure because the corresponding errors are lower than $10^{-16}$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 13. Example 3—Zoom-in images of balls falling under the gravity effect. Initial configuration on the left, final configuration on the right.</figcaption></figure>
<table><caption>Table III. Comparison of the results obtained by the four methods on the third example (after the 1000th iteration, $N_{max} = 5000$ iterations).</caption><thead><tr><th>Method</th><th>Number of NLGS iterations (last time step)</th><th>Error ε<sub>glob</sub> (last time step)</th><th>Maximal penetration (last time step)</th><th>Total CPU time (s)</th></tr></thead><tbody><tr><td>SBP</td><td>5000</td><td>0.119 · 10<sup>-6</sup></td><td>0.213 · 10<sup>-5</sup></td><td>1092.95</td></tr><tr><td>SAL</td><td>5000</td><td>0.135 · 10<sup>-6</sup></td><td>0.533 · 10<sup>-5</sup></td><td>973.31</td></tr><tr><td>EBP</td><td>5000</td><td>0.156 · 10<sup>-6</sup></td><td>0.286 · 10<sup>-6</sup></td><td>854.31</td></tr><tr><td>EAL</td><td>5000</td><td>0.101 · 10<sup>-6</sup></td><td>0.390 · 10<sup>-17</sup></td><td>916.65</td></tr></tbody></table>

NLGS, nonlinear Gauss-Seidel; SBP, standard bipotential-based method; SAL, standard augmented Lagrangian; EBP, enhanced bipotential method; EAL, enhanced augmented Lagrangian.

<table><caption>Table IV. Comparison of the results obtained for various values of $ρ = α\bar{ρ}$ on the third example (after the 500th iteration, $N_{max} = 5000$ iterations)</caption><thead><tr><th>α</th><th>Number of NLGS iterations (last time step)</th><th>Maximal penetration (last time step)</th><th>Total CPU time (s)</th></tr></thead><tbody><tr><td>5</td><td>652</td><td>0.110 · 10<sup>-6</sup></td><td>65.08</td></tr><tr><td>2</td><td>414</td><td>0.177 · 10<sup>-6</sup></td><td>55.86</td></tr><tr><td>1</td><td>750</td><td>0.149 · 10<sup>-6</sup></td><td>53.95</td></tr><tr><td>1/2</td><td>812</td><td>0.634 · 10<sup>-6</sup></td><td>78.43</td></tr><tr><td>1/5</td><td>667</td><td>0.219 · 10<sup>-5</sup></td><td>176.14</td></tr></tbody></table>

NLGS, nonlinear Gauss-Seidel.

<figure><img src="image_15.png" /><figcaption>Figure 14. Example 3—Convergence of the nonlinear Gauss-Seidel iterations for the standard bipotential-based (SBP) method (1000th time step). The last two curves collapse.</figcaption></figure>

In this example, the difference between methods SBP and SAL on the one hand, and the methods EBP and EAL on the other hand is larger (Table III). We can notice that the SAL method is 10.95% faster than the SBP method, and the EBP is 21.83% faster than the SBP method. Here, the EAL method is no longer the faster one, but the penetration is very small (Figure 17). Again, for the first two methods the global error is essentially due to the penetrations (Figures 14 and 15), whereas for the last two methods, the error is essentially due to the failure to follow precisely the equations of motion (Figures 16 and 17).
### 4.4. Discussion on the descent parameter $\rho$

We consider again the third example solved by the Newton and bipotential method ($\epsilon_{tot} = 10^{-8}$, maximal number of iterations of Newton method equal to 100, $\epsilon_{Newt} = 10^{-5}$, 500th time step). Here, we take $\bar{\rho} = \frac{m_i m_j}{m_i + m_j} \frac{1}{\Delta t}$, and we consider $\rho = \alpha \bar{\rho}$, for various values of $\alpha$.

These results show one of the main advantage of the EBP method. Indeed, one can notice that in Table IV, the CPU time and the quality of the solution are very similar if $\alpha$ is equal to 1 or 2. Even if $\alpha$ is equal to 5, the convergence is not badly damaged. One remains that for this value of $\alpha$, the SBP and the SAL methods are no longer convergent. If the parameter $\alpha$ is small, the method converges, but the convergence rate is very small. One can notice that the EAL method is much more sensitive to the parameter $\alpha$, essentially in the convergence of the Newton method. This phenomena can be explained by the fact that the EAL method differs from the SAL method only if the augmented reaction $\tau^{i+1}$ is outside the admissible set $C_c(\mu\mathbf{r}^{i+1})$.

<figure><img src="image_16.png" /><figcaption>Figure 15. Example 3—Convergence of the nonlinear Gauss–Seidel iterations for the standard augmented Lagrangian (SAL) method (1000th time step). The last two curves collapse.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 16. Example 3—Convergence of the nonlinear Gauss–Seidel iterations for the Newton and enhanced bipotential method (EBP) (1000th time step). The last two curves collapse.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure 17. Example 3—Convergence of the nonlinear Gauss-Seidel iterations for the Newton and enhanced augmented Lagrangian (EAL) method (1000th time step). The last two curves collapse. The curve corresponding to the penetrations does not appear because the penetrations are lower than $10^{-14}$.</figcaption></figure>

## 5. CONCLUSION

The results presented show that, using an appropriate Newton method, it is possible to improve the computational time to more than 20% compared with the standard methods. Moreover, one principal drawback of that type of methods, that is, the dependence of the results on the parameter $\rho$, does not exist anymore.

On the other hand, the stopping criterion developed appears to be efficient. It permits to have a unified criterion to compare the methods. Moreover, it ensures the mechanical quality of the solutions.

In the future, this method will be extended to the case of a contact law with adhesion. The difficulty here is the writing of appropriate local equations to model the adhesion. Even in the case of a frictional contact law, one can notice that the local equations that have to be solved come from a choice, and the efficiency of the method strongly depends on this choice.

## ACKNOWLEDGEMENTS

This work was supported by the French CNRS grant ‘PEPS Opticontact’ and by the project ‘TRIBAL’ supported by the Picardy region and the European funds FEDER.

## REFERENCES

1. De Saxcé G, Feng Z-Q. New inequality and functional for contact with friction: the implicit standard material approach. *Mechanics of Structures and Machines* 1991; **19**:301–325.

2. Alart P, Curnier A. A mixed formulation for frictional contact problems prone to Newton like solution method. *Computer Methods in Applied Mechanics and Engineering* 1991; **92**(3):353–375.

3. Feng Z-Q, Joli P, Cros J-M, Magnain B. The bi-potential method applied to the modeling of dynamics problems with friction. *Computational Mechanics* 2005; **36**(5):375–383.

4. Klarbring A. Mathematical programming in contact problems. In *Computational Methods in Contact Mechanics*, Aliabdali MH, Brebbia CA (eds). Computational Mechanics. Publications: Southampton, 1993; 233-263.

5. Wriggers P. Finite element algorithms for contact problems. *Archives of Computational Methods in Engineering* 1995; **2**:1–49.

6. Jean M, Moreau JJ. Unilaterality and dry friction in the dynamics of rigid bodies collection. In *Contact Mechanics International Symposium*, Curnier A (ed.). Presses Polytechniques et Universitaires Romanes: Lausanne, 1992; 31–48.
7. Jean M. The non smooth contact dynamics method. *Computer Methods Applied Mathematics Engineering* 1999; **177**:235-257.

8. Moreau JJ. Unilateral contact and dry friction in finite freedom dynamics. In *Non Smooth Mechanics and Applications, CISM Courses and Lectures*, Vol. 302, Moreau JJ, Panagiotopoulos P-D (eds). Springer-Verlag, Wien: New York, 1988; 1-82.

9. Moreau JJ. Some numerical methods in multibody dynamics: application to granular materials. *European Journal of Mechanics - A/Solids* 1994; **13**(4-suppl.):93-114.

10. Moreau JJ. Numerical aspect of sweeping process. *Computer Methods Applied Mathematics Engineering* 1999; **177**:329-349.

11. Renouf M, Alart P. Conjugate gradient type algorithms for frictional multicontact problems: applications to granular materials. *Computer Methods in Applied Mechanics and Engineering* 2004; **194**(18-20):2019-2041.

12. Jourdan F, Alart P, Jean M. A Gauss-Seidel like algorithm to solve frictional contact problems. *Computer Methods in Applied Mechanics and Engineering* 1998; **155**:31-47.

13. Renouf M, Dubois F, Alart P. A parallel version of the non smooth contact dynamics algorithm applied to the simulation of granular media. *Journal of Computational and Applied Mathematics* 2004; **168**:375-38.

14. Fortin J. Simulation numérique de la dynamique des systèmes multicorps appliquée aux milieux granulaires. *PhD Thesis (in french)*, University of Lille 1, 2000.

15. Fortin J, Millet O, de Saxcé G. Numerical simulation of granular materials by an improved discrete element method. *International Journal for Numerical Methods in Engineering* 2005; **62**(5):639-663.

16. Sanni I. Modélisation et simulation bi et tri-dimensionnelles de la dynamique unilatérale des systèmes multi-corps de grandes tailles: application aux milieux granulaires. *PhD Thesis (in french)*, University of Picardy, 2006.

17. Joli P, Feng Z-Q. Uzawa and Newton algorithms to solve frictional contact problems within the bi-potential framework. *International Journal for Numerical Methods in Engineering* 2008; **73**(3):317-330.

18. Ramaniraka NA. Thermomécanique des contacts entre deux solides déformables. *PhD thesis no. 1651 (in french)*, Ecole Polytechnique Fédérale de Lausanne, 1991.

19. Renouf M, Alart P. Gradient type algorithms for 2d/3d frictionless/frictional multicontact problems. In *ECCOMAS 2004*, Majava K, Neittaanmaki P, Rossi T, Pironneau O, Owen R, Mikkola M (eds): Jyvaskyla, July 2004; 24-28.

20. Renouf M, Acary V, Dumont G. Comparison of algorithms for collisions, contact and friction in view of real-time applications. In *Multibody Dynamics 2005 proceedings, International Conference on Advances in Computational Multibody*, Goicolea JM, Cuadrado J, García Orden JC (eds), Madrid, June 2005; 21-24.

21. Renouf M, Bonamy D, Dubois F, Alart P. Numerical simulation of two-dimensional steady granular flows in rotating drum: on surface flow rheology. *Physics of Fluids* 2005; **17**:103-303.

22. Rockafellar RT. *Convex Analysis*, Vol. 28, Princeton Math. Series. Princeton Univ. Press: Princeton, 1970.

23. Fortin J, Hjiaj M, de Saxcé G. An improved discrete element method based on a variational formulation of the contact law. *Computers and Geotechnics* 2002; **29**:609-640.