# Unified primal-dual active set method for dynamic frictional contact problems

Stéphane Abide, Mikaël Barboteu, Soufiane Cherkaoui, Serge Dumont

► To cite this version:

Stéphane Abide, Mikaël Barboteu, Soufiane Cherkaoui, Serge Dumont. Unified primal-dual active set method for dynamic frictional contact problems. Fixed Point Theory and Algorithms for Sciences and Engineering, 2022, 2022 (1), pp.19. 10.1186/s13663-022-00729-4. hal-04053755

HAL Id: hal-04053755

https://hal.science/hal-04053755

Submitted on 22 May 2024

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>
<figure><img src="image_1.png" /></figure>

# Unified primal-dual active set method for dynamic frictional contact problems

Stéphane Abide¹* ID, Mikaël Barboteu¹, Soufiane Cherkaoui¹ and Serge Dumont² ID

stephane.abide@univ-perp.fr

¹Laboratoire de Mathématiques et Physique, Université de Perpignan Via Domitia, 52 Avenue Paul Alduy, 66860 Perpignan, France

Full list of author information is available at the end of the article

**Abstract**

In this paper, we propose a semi-smooth Newton method and a primal-dual active set strategy to solve dynamical contact problems with friction. The conditions of contact with Coulomb's friction can be formulated in the form of a fixed point problem related to a quasi-optimization one thanks to the semi-smooth Newton method. This method is based on the use of the primal-dual active set (PDAS) strategy. The main idea here is to find the correct subset $\mathcal{A}$ of nodes that are in contact (active) opposed to those which are not in contact (inactive). For each case, the nonlinear boundary condition is replaced by a suitable linear one. Numerical experiments on both hyper-elastic problems and rigid granular materials are presented to show the efficiency of the proposed method.

**MSC:** 70F40; 70-08; 70E55; 35Q70

**Keywords:** Granular media; Elasticity; Unilateral constraint; Friction; Rigid body; Deformable body; Discrete element method; Nonsmooth contact dynamics; Semi-smooth Newton method; Primal-dual active set; Numerical simulations

## 1 Introduction

Nowadays, understanding dynamic contact problems and their behaviour, both de-formable and rigid, is a major challenge for industries in several sectors. The scientific community is particularly interested in studying these problems because of the complex nature of the interactions governing their dynamics and their simulations, which involve a large number of degrees of freedom. Some derivations for a model of contact dynamics and their specific numerical methods can illustrated by these works [1–4].

According to the theory of sweeping process of Moreau [5–9], the dynamic contact problems can be written as a measure differential inclusion due to multi-valued mappings between contact reactions/impulses and the contact distances/velocities that introduce a nonsmoothness contact law. In this way, convenient formulations for nonsmooth dynamic problems can be derived. This appears to be well adapted for problems involving rigid bodies and elastodynamic problems.

The *nonsmooth contact dynamics* (NSCD) approach is an alternative simulation strategy that works for both rigid body interactions and deformable problems [3, 4, 8–11]. According to Moreau's sweeping process, it leads to implicit numerical schemes based

<footnote>*Correspondence:</footnote>
<figure><img src="image_2.png" /></figure>

Check for
updates

# Unified primal-dual active set method for dynamic frictional contact problems

Stéphane Abide<fnref n="1" />* ID, Mikaël Barboteu<fnref n="1" />, Soufiane Cherkaoui<fnref n="1" /> and Serge Dumont² ID

**Abstract**

In this paper, we propose a semi-smooth Newton method and a primal-dual active set strategy to solve dynamical contact problems with friction. The conditions of contact with Coulomb's friction can be formulated in the form of a fixed point problem related to a quasi-optimization one thanks to the semi-smooth Newton method. This method is based on the use of the primal-dual active set (PDAS) strategy. The main idea here is to find the correct subset $\mathcal{A}$ of nodes that are in contact (active) opposed to those which are not in contact (inactive). For each case, the nonlinear boundary condition is replaced by a suitable linear one. Numerical experiments on both hyper-elastic problems and rigid granular materials are presented to show the efficiency of the proposed method.

**MSC:** 70F40; 70-08; 70E55; 35Q70

**Keywords:** Granular media; Elasticity; Unilateral constraint; Friction; Rigid body; Deformable body; Discrete element method; Nonsmooth contact dynamics; Semi-smooth Newton method; Primal-dual active set; Numerical simulations

## 1 Introduction

Nowadays, understanding dynamic contact problems and their behaviour, both de-formable and rigid, is a major challenge for industries in several sectors. The scientific community is particularly interested in studying these problems because of the complex nature of the interactions governing their dynamics and their simulations, which involve a large number of degrees of freedom. Some derivations for a model of contact dynamics and their specific numerical methods can illustrated by these works [1–4].

According to the theory of sweeping process of Moreau [5–9], the dynamic contact problems can be written as a measure differential inclusion due to multi-valued mappings between contact reactions/impulses and the contact distances/velocities that introduce a nonsmoothness contact law. In this way, convenient formulations for nonsmooth dynamic problems can be derived. This appears to be well adapted for problems involving rigid bodies and elastodynamic problems.

The *nonsmooth contact dynamics* (NSCD) approach is an alternative simulation strategy that works for both rigid body interactions and deformable problems [3, 4, 8–11]. According to Moreau's sweeping process, it leads to implicit numerical schemes based

<footnote>*Correspondence:
stephane.abide@univ-perp.fr</footnote>

<footnote n="1">¹Laboratoire de Mathématiques et Physique, Université de Perpignan Via Domitia, 52 Avenue Paul Alduy, 66860 Perpignan, France</footnote>

<footnote>Full list of author information is available at the end of the article</footnote>
on a nonregular formulation of the contact laws. There are many numerical methods in the literature for solving these problems. In the case of rigid body problems, one can refer to the well-known nonlinear Gauss-Seidel method (NLGS) developed by M. Jean and J. J. Moreau [8–10, 12]. Moreover, there are several sophisticated approaches based on the bi-potential theory and the augmented Lagrangian theory, see [13–17] for more details. Concerning the deformable bodies, several methods have also been successfully tested, namely, the penalty method (cf. [18, 19]), the quasi-augmented Lagrangian [20], the bi-potential method (cf. [13, 17]), the conjugate gradient method (cf. [1, 2]), the Uzawa method (cf. [21, 22]), and the Nitsche finite element method (cf. [23–25]). Recently proposed methods of semi-smooth Newton and primal-dual active (PDAS) type appear to be some of the most relevant methods for solving friction contact problems [26–34]. They are based on the following principle: the conditions of contact and friction are reformulated in terms of nonlinear complementarity functions whose solution is provided by the semi-smooth Newton method. In practice, the conditions of contact with Coulomb’s friction can be formulated in the form of a problem of fixed point related to quasi-optimisation problem [26, 27, 34]. In fact, the boundary conditions are directly imposed thanks to a semi-smooth Newton method, and therefore their implementation could be achieved without much effort. On the basis of these prerequisites, the objective of this work is to provide a generalisation of the semi-smooth Newton method-PDAS approach for hyper-elastic contact problems [34, 35] and for those in rigid multi-body dynamics [33, 36]. Particular care is paid to the development of these algorithms for solving contact and friction laws in the nonregular framework. Several numerical experiments are reported for verification and validation purposes, but also to assess the efficiency and performance of PDAS methods compared to other numerical methods [13, 17, 20].

The paper is organised as follows. Section 2 presents the sweeping process in the context of deformable and rigid bodies. Section 3 is devoted to the semi-smooth Newton approach for solving the frictional contact conditions. In Sect. 4, we present the numerical treatment of frictional contact conditions by using two primal-dual active set methods. This section is divided into two parts. First, an exact PDAS strategy for unilateral contact with Coulomb’s law of friction is considered. Next, an iterative PDAS with a fixed point method related to the friction bound is introduced. In the last section, the results of some simulations are presented in order to illustrate the numerical PDAS.

## 2 Sweeping process for frictional contact problems

According to Moreau’s sweeping process, the dynamic frictional contact problem can be written as a measure differential inclusion due to multi-valued mappings between contact reactions and the contact velocities that introduce a nonsmoothness in the contact law. This formulation allows us to write nonsmooth dynamic problems such as dynamics contact problems for rigid bodies and elastodynamic problems.

### 2.1 Dynamics of deformable bodies

Let us consider a deformable body occupying a domain $\Omega \subset \mathbb{R}^d$ ($d = 1, 2, 3$). The boundary $\Gamma$ of $\Omega$ is assumed to be Lipschitz continuous and consists in three disjoint measurable parts $\Gamma_1, \Gamma_2$ and $\Gamma_3$ such that $\text{meas}(\Gamma_1) > 0$. The notations $\mathbf{x} = (x_i)$ and $\mathbf{n} = (n_i)$ stand for a point in $\Omega \cup \Gamma$ and the inner unit normal at $\Gamma$, respectively. The dynamic of a deformable
on a nonregular formulation of the contact laws. There are many numerical methods in the literature for solving these problems. In the case of rigid body problems, one can refer to the well-known nonlinear Gauss-Seidel method (NLGS) developed by M. Jean and J. J. Moreau [8–10, 12]. Moreover, there are several sophisticated approaches based on the bi-potential theory and the augmented Lagrangian theory, see [13–17] for more details. Concerning the deformable bodies, several methods have also been successfully tested, namely, the penalty method (cf. [18, 19]), the quasi-augmented Lagrangian [20], the bi-potential method (cf. [13, 17]), the conjugate gradient method (cf. [1, 2]), the Uzawa method (cf. [21, 22]), and the Nitsche finite element method (cf. [23–25]). Recently proposed methods of semi-smooth Newton and primal-dual active (PDAS) type appear to be some of the most relevant methods for solving friction contact problems [26–34]. They are based on the following principle: the conditions of contact and friction are reformulated in terms of nonlinear complementarity functions whose solution is provided by the semi-smooth Newton method. In practice, the conditions of contact with Coulomb’s friction can be formulated in the form of a problem of fixed point related to quasi-optimisation problem [26, 27, 34]. In fact, the boundary conditions are directly imposed thanks to a semi-smooth Newton method, and therefore their implementation could be achieved without much effort. On the basis of these prerequisites, the objective of this work is to provide a generalisation of the semi-smooth Newton method-PDAS approach for hyper-elastic contact problems [34, 35] and for those in rigid multi-body dynamics [33, 36]. Particular care is paid to the development of these algorithms for solving contact and friction laws in the nonregular framework. Several numerical experiments are reported for verification and validation purposes, but also to assess the efficiency and performance of PDAS methods compared to other numerical methods [13, 17, 20].

The paper is organised as follows. Section 2 presents the sweeping process in the context of deformable and rigid bodies. Section 3 is devoted to the semi-smooth Newton approach for solving the frictional contact conditions. In Sect. 4, we present the numerical treatment of frictional contact conditions by using two primal-dual active set methods. This section is divided into two parts. First, an exact PDAS strategy for unilateral contact with Coulomb’s law of friction is considered. Next, an iterative PDAS with a fixed point method related to the friction bound is introduced. In the last section, the results of some simulations are presented in order to illustrate the numerical PDAS.

## 2 Sweeping process for frictional contact problems

According to Moreau’s sweeping process, the dynamic frictional contact problem can be written as a measure differential inclusion due to multi-valued mappings between contact reactions and the contact velocities that introduce a nonsmoothness in the contact law. This formulation allows us to write nonsmooth dynamic problems such as dynamics contact problems for rigid bodies and elastodynamic problems.

### 2.1 Dynamics of deformable bodies

Let us consider a deformable body occupying a domain $\Omega \subset \mathbb{R}^d$ ($d = 1, 2, 3$). The boundary $\Gamma$ of $\Omega$ is assumed to be Lipschitz continuous and consists in three disjoint measurable parts $\Gamma_1, \Gamma_2$ and $\Gamma_3$ such that $\text{meas}(\Gamma_1) > 0$. The notations $\mathbf{x} = (x_i)$ and $\mathbf{n} = (n_i)$ stand for a point in $\Omega \cup \Gamma$ and the inner unit normal at $\Gamma$, respectively. The dynamic of a deformable
<figure><img src="image_2.png" /><figcaption>Figure 1 Physical setting of a deformable body against a foundation</figcaption></figure>

body is a time-dependent problem with $t > 0$ denoting the time variable. The primary variable is the displacement field $\mathbf{u}$ from which one can define the normal and the tangential components of $\mathbf{u}$ on $\Gamma$ by $u_n = \mathbf{u} \cdot \mathbf{n}$ and $u_t = \mathbf{u} - u_n \mathbf{n}$.

In the following, we consider the dynamic frictional contact of a deformable body with a perfectly rigid obstacle, the so-called foundation. Figure 1 shows the configuration and the considered notations (see Fig. 1). The displacements are driven by the body force and the boundary conditions. The displacement and the surface traction are prescribed on the boundary $\Gamma_1$ and $\Gamma_2$, respectively. The contact condition is assumed to occur on $\Gamma_3$. The dynamic of the body is investigated on the time interval $[0, T]$ with $T > 0$. For the sake of clarity, the dependence of the functions on $\mathbf{x}$ and $t$ is not formally expressed, and the derivatives with respect to the time read $\dot{\mathbf{x}}$.

### 2.1.1 Equations of motion

According to the fundamental principle of dynamics of deformable bodies, after a spatial discretization using a standard finite element method, the equations of motion formulated in terms of differential measures can be written as follows:

$$ \mathbb{M} d\dot{\mathbf{U}} + \mathbf{A}(\mathbf{U}) dt = \mathbf{F}^{\text{ext}}(t) dt + d\mathbf{\Lambda}, \quad (2.1) $$

where

- $\mathbf{U} \in \mathbb{R}^{d \times N}$ is the generalised displacement vector (where $d = 3$ for a 3D problem and $d = 2$ for a 2D problem and $N$ is the total number of degree of freedom);

- $\mathbb{M}$ represents the generalised mass matrix;

- $\mathbf{A}$ represents the internal force vector;

- $\mathbf{F}^{\text{ext}}$ represent external forces;

- $d\mathbf{\Lambda}$ is a nonnegative real measure, representing the reaction forces between the body in contact.

### 2.1.2 Frictional contact laws

Here we expressed the contact law by using the local contact variables with the vertex $\alpha \in [1, N_\alpha]$ (where $N_\alpha$ is the total number of potential contacts between the body to the foundation). For each $\alpha \in [1, N_\alpha]$, we have to consider the local relative displacements $u_n^\alpha$ and $u_t^\alpha$ between the body to the foundation involved in the contact $\alpha$, as a part of the
<figure><img src="image_3.png" /><figcaption>**Figure 1** Physical setting of a deformable body against a foundation</figcaption></figure>

body is a time-dependent problem with $t > 0$ denoting the time variable. The primary variable is the displacement field $\mathbf{u}$ from which one can define the normal and the tangential components of $\mathbf{u}$ on $\Gamma$ by $u_n = \mathbf{u} \cdot \mathbf{n}$ and $u_t = \mathbf{u} - u_n\mathbf{n}$.

In the following, we consider the dynamic frictional contact of a deformable body with a perfectly rigid obstacle, the so-called foundation. Figure 1 shows the configuration and the considered notations (see Fig. 1). The displacements are driven by the body force and the boundary conditions. The displacement and the surface traction are prescribed on the boundary $\Gamma_1$ and $\Gamma_2$, respectively. The contact condition is assumed to occur on $\Gamma_3$. The dynamic of the body is investigated on the time interval $[0, T]$ with $T > 0$. For the sake of clarity, the dependence of the functions on $\mathbf{x}$ and $t$ is not formally expressed, and the derivatives with respect to the time read $(\cdot)$.

### 2.1.1 Equations of motion

According to the fundamental principle of dynamics of deformable bodies, after a spatial discretization using a standard finite element method, the equations of motion formulated in terms of differential measures can be written as follows:

$$ \mathbb{M} d\dot{\mathbf{U}} + \mathbf{A}(\mathbf{U}) dt = \mathbf{F}^{\text{ext}}(t) dt + d\mathbf{\Lambda}, \qquad (2.1) $$

where

- $\mathbf{U} \in \mathbb{R}^{d \times N}$ is the generalised displacement vector (where $d = 3$ for a 3D problem and $d = 2$ for a 2D problem and $N$ is the total number of degree of freedom);

- $\mathbb{M}$ represents the generalised mass matrix;

- $\mathbf{A}$ represents the internal force vector;

- $\mathbf{F}^{\text{ext}}$ represent external forces;

- $d\mathbf{\Lambda}$ is a nonnegative real measure, representing the reaction forces between the body in contact.

### 2.1.2 Frictional contact laws

Here we expressed the contact law by using the local contact variables with the vertex $\alpha \in [1, N_\alpha]$ (where $N_\alpha$ is the total number of potential contacts between the body to the foundation). For each $\alpha \in [1, N_\alpha]$, we have to consider the local relative displacements $u_n^\alpha$ and $u_t^\alpha$ between the body to the foundation involved in the contact $\alpha$, as a part of the
generalised vector $\mathbf{u} \in \mathbb{R}^{d \times N}$. Likewise, the local contact reactions $\lambda^{\alpha} \in \mathbb{R}^d$ are represented by a vector $\Lambda \in \mathbb{R}^{d \times N}$. We recall that $\lambda^{\alpha}$ can be decomposed into the sum of a normal component $\lambda_n^{\alpha}$ and a tangential component $\lambda_t^{\alpha}$ as follows: $\lambda^{\alpha} = \lambda_n^{\alpha} \mathbf{n} + \lambda_t^{\alpha} \mathbf{t}$. From there, the idea would be to project the global vector from $\Lambda$ and $\mathbf{u}$ to the local vector $\lambda^{\alpha}$ and $\mathbf{u}^{\alpha}$. In other words, $\lambda^{\alpha}$ and $\mathbf{u}^{\alpha}$ denote the coefficient vector of $\Lambda$ and $\mathbf{u}$ associated with the vertex $\alpha$.

Consider for a while that the normal reaction exerted between the deformable body against a foundation at the contact $\alpha$ of the boundary $\Gamma_3$, and with these considerations, the unilateral conditions on the contact boundary point are

$$u_n^{\alpha} \geq g, \lambda_n^{\alpha} \geq 0, \quad (u_n^{\alpha} - g)\lambda_n^{\alpha} = 0, \qquad (2.2)$$

where the gap $g$ measures the distance between a point on $\Gamma_3$ and its projection onto the rigid obstacle.

Now for the friction, we considered Coulomb's law as follows:

$$\left\{ \begin{array}{l} \|\lambda_t^\alpha\| \le \mu |\lambda_n^\alpha|, \\ \|\lambda_t^\alpha\| < \mu |\lambda_n^\alpha| \Rightarrow \dot{\mathbf{u}}_t^\alpha = 0, \\ \|\lambda_t^\alpha\| = \mu |\lambda_n^\alpha| \Rightarrow \exists \beta \ge 0 : \lambda_t^\alpha = \beta \dot{\mathbf{u}}_t^\alpha, \end{array} \right. \qquad (2.3)$$

where $\mu$ is the friction coefficient.

Note that conditions (2.2) are equivalent to the following subdifferential inclusion (cf. [6, 9]):

$$\lambda_n^\alpha \in \partial I_{\mathbb{R}^-} (u_n^\alpha - g), \qquad (2.4)$$

where $\partial$ represents the subdifferential operator in the sense of the convex analysis and $I_A$ denotes the indicator function of the set $A \subset \mathbb{R}$. A similar consideration for the frictional stress leads to

$$\lambda_t^\alpha \in \partial I_{C(\lambda_n^\alpha)}^*(\dot{\mathbf{u}}_t^\alpha), \qquad (2.5)$$

which is equivalent to (2.3). $I^*$ represents the Fenchel conjugate of the indicator function.

Finally, the objective frictional contact constitutive law can be written jointly in the sub-differential inclusion form as

$$\lambda^{\alpha} \in \partial I_{\mathbb{R}^{-}}(u_{n}^{\alpha} - g) \mathbf{n} + \partial I_{C(\lambda_{n}^{\alpha})}^{*}(\dot{u}_{t}^{\alpha}). \qquad (2.6)$$

As a result, the generalised frictional contact forces $d\Lambda$ can also take this form. The so-called Moreau's sweeping process [5, 9, 37] is a special kind of differential inclusion with a maximal monotone operator, which appears to be a very nice formulation for the unilateral dynamics with friction.

### 2.1.3 Nonlinear elastodynamic contact problems

First, we recall some preliminary material concerning the time discretization step. Let $N$ be an integer, let $\Delta\tau = \frac{T}{N}$ be the time step and define

$$t^n = n \Delta t, \quad 0 \le n \le N.$$
generalised vector $u \in \mathbb{R}^{d \times N}$. Likewise, the local contact reactions $\lambda^\alpha \in \mathbb{R}^d$ are represented by a vector $\Lambda \in \mathbb{R}^{d \times N}$. We recall that $\lambda^\alpha$ can be decomposed into the sum of a normal component $\lambda_n^\alpha$ and a tangential component $\lambda_t^\alpha$ as follows: $\lambda^\alpha = \lambda_n^\alpha \mathbf{n} + \lambda_t^\alpha \mathbf{t}$. From there, the idea would be to project the global vector from $\Lambda$ and $u$ to the local vector $\lambda^\alpha$ and $u^\alpha$. In other words, $\lambda^\alpha$ and $u^\alpha$ denote the coefficient vector of $\Lambda$ and $u$ associated with the vertex $\alpha$.

Consider for a while that the normal reaction exerted between the deformable body against a foundation at the contact $\alpha$ of the boundary $\Gamma_3$, and with these considerations, the unilateral conditions on the contact boundary point are

$$u_n^\alpha \geq g, \lambda_n^\alpha \geq 0, \quad (u_n^\alpha - g)\lambda_n^\alpha = 0, \qquad (2.2)$$

where the gap $g$ measures the distance between a point on $\Gamma_3$ and its projection onto the rigid obstacle.

Now for the friction, we considered Coulomb's law as follows:

$$\left\{ \begin{array}{l} \|\lambda_t^\alpha\| \le \mu |\lambda_n^\alpha|, \\ \|\lambda_t^\alpha\| < \mu |\lambda_n^\alpha| \Rightarrow \dot{u}_t^\alpha = 0, \\ \|\lambda_t^\alpha\| = \mu |\lambda_n^\alpha| \Rightarrow \exists \beta \ge 0 : \lambda_t^\alpha = \beta \dot{u}_t^\alpha, \end{array} \right. \qquad (2.3)$$

where $\mu$ is the friction coefficient.

Note that conditions (2.2) are equivalent to the following subdifferential inclusion (cf. [6, 9]):

$$\lambda_n^\alpha \in \partial I_{\mathbb{R}^{-}}(u_n^\alpha - g), \qquad (2.4)$$

where $\partial$ represents the subdifferential operator in the sense of the convex analysis and $I_A$ denotes the indicator function of the set $A \subset \mathbb{R}$. A similar consideration for the frictional stress leads to

$$\lambda_t^\alpha \in \partial I_{C(\lambda_n^\alpha)}^*(\dot{u}_t^\alpha), \qquad (2.5)$$

which is equivalent to (2.3). $I^*$ represents the Fenchel conjugate of the indicator function.

Finally, the objective frictional contact constitutive law can be written jointly in the sub-differential inclusion form as

$$\lambda^\alpha \in \partial I_{\mathbb{R}^{-}}(u_n^\alpha - g) \mathbf{n} + \partial I_{C(\lambda_n^\alpha)}^*(\dot{u}_t^\alpha). \qquad (2.6)$$

As a result, the generalised frictional contact forces $d\Lambda$ can also take this form. The so-called Moreau's sweeping process [5, 9, 37] is a special kind of differential inclusion with a maximal monotone operator, which appears to be a very nice formulation for the unilateral dynamics with friction.

### 2.1.3 Nonlinear elastodynamic contact problems

First, we recall some preliminary material concerning the time discretization step. Let $N$ be an integer, let $\Delta\tau = \frac{T}{N}$ be the time step and define

$$t^n = n \Delta t, \quad 0 \le n \le N.$$
Finally, for a sequence {$w^n$}$_{n=1}^N$, we denote the midpoint divided differences by

$$ \dot{w}^{n-\frac{1}{2}} = (w^n - w^{n-1})/\Delta t = \frac{1}{2}(\dot{w}^n + \dot{w}^{n-1}), \qquad (2.7) $$

and, equivalently, we have $\dot{w}_n = -\dot{w}^{n-1} + \frac{2}{\Delta t}(w^n - w^{n-1})$. Hereinafter, we use the notation $\Box^{n-\frac{1}{2}} = \frac{1}{2}(\Box^n + \Box^{n-1})$, where $\Box^n$ represents the approximation of $\Box(t^n)$. Note that the time integration scheme we use is based on the implicit second order midpoint rule given in (2.7).

**Problem P<sub>S</sub>.** Find a global displacement vector {$u^{n,\Delta t}$}$_{n=0}^N$, a global normal stress vector {$\Lambda_n^{n,\Delta t}$}$_{n=0}^N$ and a global tangential stress vector {$\Lambda_t^{n,\Delta t}$}$_{n=0}^N$ such that, for all $n = 1$,

$$ \mathbb{M}\ddot{u}^{n-\frac{1}{2},\Delta t} + \mathbb{A}(u^{n-\frac{1}{2},\Delta t}) - \mathbb{F}^{\text{ext},n-\frac{1}{2},\Delta t} - \Lambda^{n-\frac{1}{2},\Delta t} = \mathbf{0}, \qquad (2.8) $$

$$ \lambda_n^{\alpha, n-\frac{1}{2}, \Delta t} \in \partial I_{\mathbb{R}^{-}}(u_n^{\alpha, n-\frac{1}{2}, \Delta t} - g), \qquad (2.9) $$

$$ \lambda_t^{\alpha, n-\frac{1}{2}, \Delta t} \in \partial I_{C(\lambda_n^{\alpha, n-\frac{1}{2}, \Delta t})}^*(\dot{u}_t^{\alpha, n-\frac{1}{2}, \Delta t}), \qquad (2.10) $$

$$ u^{0, \Delta t} = u^0, \quad \dot{u}^{0, \Delta t} = \dot{u}^0, $$

where $\ddot{u}^{n-\frac{1}{2},\Delta t} = \frac{\dot{u}^{n,\Delta t} - \dot{u}^{n-1,\Delta t}}{\Delta t}$ is the midpoint time approximation of the acceleration $\ddot{u}$ at the time $t^{n-\frac{1}{2}}$.

To solve the global Problem P<sub>S</sub>, we can use the so-called Newton-Raphson method. For more details about this method, see the books [1, 2].

## 2.2 Dynamics of rigid bodies

### 2.2.1 Equations of motion

To describe the motion of a multi-contact system between rigid bodies, we use the following writing conventions. Assuming that a particle P among $N_p$ particles is described by the position of its centre of gravity and its rotation, we will denote by $q$ the generalised coordinate describing its position in space ($q \in \mathbb{R}^{d \times N_p}$, where $\bar{d} = 6$ for a 3D problem and $\bar{d} = 3$ for a 2D problem). As a consequence of the possible shocks between particles, we introduce the generalised velocity denoted by $\dot{q}$ as a function of bounded variations and its associated differential $d\dot{q}$. According to the fundamental principle of rigid dynamics, the equations of motion formulated in terms of differential measures can be written as follows:

$$ \mathbb{M}d\dot{q} = \mathbf{F}(t, q, \dot{q})dt + d\mathbf{P}, \qquad (2.11) $$

where

- $\mathbb{M}$ represents the generalised mass matrix;

- $\mathbf{F}$ represents external forces;

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

### 2.2.2 Frictional contact laws

As the contact law is expressed using the contact variables, we have to express the equations of motion using these same variables. Thus, a local-global mapping is defined from
Finally, for a sequence {$w^n$}$_{n=1}^N$, we denote the midpoint divided differences by

$$ \dot{w}^{n-\frac{1}{2}} = (w^n - w^{n-1})/\Delta t = \frac{1}{2}(\dot{w}^n + \dot{w}^{n-1}), \qquad (2.7) $$

and, equivalently, we have $\dot{w}_n = -\dot{w}^{n-1} + \frac{2}{\Delta t}(w^n - w^{n-1})$. Hereinafter, we use the notation $\Box^{n-\frac{1}{2}} = \frac{1}{2}(\Box^n + \Box^{n-1})$, where $\Box^n$ represents the approximation of $\Box(t^n)$. Note that the time integration scheme we use is based on the implicit second order midpoint rule given in (2.7).

**Problem P<sub>S</sub>.** Find a global displacement vector {$u^{n,\Delta t}$}$_{n=0}^N$, a global normal stress vector {$\Lambda_n^{n,\Delta t}$}$_{n=0}^N$ and a global tangential stress vector {$\Lambda_t^{n,\Delta t}$}$_{n=0}^N$ such that, for all $n = 1$,

$$ M\ddot{u}^{n-\frac{1}{2},\Delta t} + A(u^{n-\frac{1}{2},\Delta t}) - F^{\text{ext},n-\frac{1}{2},\Delta t} - \Lambda^{n-\frac{1}{2},\Delta t} = 0, \qquad (2.8) $$

$$ \lambda_n^{\alpha,n-\frac{1}{2},\Delta t} \in \partial I_{\mathbb{R}^{-}}(u_n^{\alpha,n-\frac{1}{2},\Delta t} - g), \qquad (2.9) $$

$$ \lambda_t^{\alpha,n-\frac{1}{2},\Delta t} \in \partial I_{C(\lambda_n^{\alpha,n-\frac{1}{2},\Delta t})}^*(\dot{u}_t^{\alpha,n-\frac{1}{2},\Delta t}), \qquad (2.10) $$

$$ u^{0,\Delta t} = u^0, \qquad \dot{u}^{0,\Delta t} = \dot{u}^0, $$

where $\ddot{u}^{n-\frac{1}{2},\Delta t} = \frac{\dot{u}^{n,\Delta t} - \dot{u}^{n-1,\Delta t}}{\Delta t}$ is the midpoint time approximation of the acceleration $\ddot{u}$ at the time $t^{n-\frac{1}{2}}$.

To solve the global Problem P<sub>S</sub>, we can use the so-called Newton-Raphson method. For more details about this method, see the books [1, 2].

## 2.2 Dynamics of rigid bodies

### 2.2.1 Equations of motion

To describe the motion of a multi-contact system between rigid bodies, we use the following writing conventions. Assuming that a particle P among $N_p$ particles is described by the position of its centre of gravity and its rotation, we will denote by $q$ the generalised coordinate describing its position in space ($q \in \mathbb{R}^{d \times N_p}$, where $\bar{d} = 6$ for a 3D problem and $\bar{d} = 3$ for a 2D problem). As a consequence of the possible shocks between particles, we introduce the generalised velocity denoted by $\dot{q}$ as a function of bounded variations and its associated differential $d\dot{q}$. According to the fundamental principle of rigid dynamics, the equations of motion formulated in terms of differential measures can be written as follows:

$$ M d\dot{q} = F(t, q, \dot{q}) dt + dP, \qquad (2.11) $$

where

- $M$ represents the generalised mass matrix;

- $F$ represents external forces;

- $dP$ is a nonnegative real measure, representing the reaction forces and impulses between particles in contact.

### 2.2.2 Frictional contact laws

As the contact law is expressed using the contact variables, we have to express the equations of motion using these same variables. Thus, a local-global mapping is defined from
the global frame to the local one for each contact node $\alpha \in [1, N_{\alpha}]$ (where $N_{\alpha}$ is the total number of potential contacts between two rigid particles) to write the corresponding local contact law. Local relative velocities $v_n^{\alpha}$ and $v_t^{\alpha}$ between the two rigid bodies involved in the contact $\alpha$ can be generalised as a vector $\dot{q} = v \in \mathbb{R}^{d \times N_c}$ (where $d = 3$ for a 3D problem and $d = 2$ for a 2D problem). Likewise, the local contact impulses $p^{\alpha}$ are represented by a vector $P \in \mathbb{R}^{d \times N_c}$. We recall that $p^{\alpha}$ can be decomposed into the sum of a normal component $p_n^{\alpha}$ and a tangential component $p_t^{\alpha}$ as follows: $p^{\alpha} = p_n^{\alpha} n + p_t^{\alpha}$. From there, the idea would be to transform the equations of motion seen in (2.11) from $P$ and $\dot{q}$ to $p^{\alpha}$ and $v^{\alpha}$.

The gap $D^{\alpha} = D(q_i, q_j)$ represents the signed distance from particle $i$ to particle $j$.

The contact and particles velocities ($v^{\alpha}$ and $\dot{q}$) being linearly related, the local-global mapping to compute $v$ to $\dot{q}$ at a contact node $\alpha$ is

$$v^{\alpha} = H^{*}(q, \alpha)\dot{q}, \qquad (2.12)$$

where $H^{*}(q, \alpha)$ is a local-global mapping matrix $dN_c \times \bar{d}N_p$ carrying the information on the geometry of the contact network. Identically, we have for the impulses

$$P = H(q, \alpha)p^{\alpha}, \qquad (2.13)$$

where $H(q, \alpha)$ is the transpose of matrix $H^{*}(q, \alpha)$.

The inelastic shock law can be written as follows:

$$\text{if } D^{\alpha} > 0 \text{ then } p_n^{\alpha} = 0, \qquad (2.14)$$

$$\begin{cases} \text{if } D^{\alpha} \le 0 \text{ then} \\ \left\{ \begin{array}{l} v_n^{\alpha} \ge 0, \\ p_n^{\alpha} \ge 0, \\ v_n^{\alpha} p_n^{\alpha} = 0. \end{array} \right. \end{cases} \qquad (2.15)$$

Coulomb's law of friction is written as follows:

$$p_t^{\alpha} \neq 0, \qquad (2.16)$$

$$\left\{ \begin{array}{l} p_n^{\alpha} = 0 \Rightarrow v_n^{\alpha} \ge 0, \\ p_n^{\alpha} > 0 \text{ and } \|p_t^{\alpha}\| < \mu p_n^{\alpha} \Rightarrow v_t^{\alpha} = 0, \\ p_n^{\alpha} > 0 \text{ and } \|p_t^{\alpha}\| = \mu p_n^{\alpha} \Rightarrow \exists \beta \ge 0, v_t^{\alpha} = \beta \frac{p_t^{\alpha}}{\|p_t^{\alpha}\|}, \end{array} \right. \qquad (2.17)$$

where $\mu$ is the friction coefficient.

The frictional contact laws (2.15) and (2.17) can be also formulated in subdifferential inclusions as

$$p^{\alpha} \in \partial I_{\mathbb{R}^{-}}(v_n^{\alpha}) n + \partial I_{C(p_n^{\alpha})}(v_t^{\alpha}). \qquad (2.18)$$

### 2.2.3 Nonlinear Gauss-Seidel (NLGS) method for nonsmooth contact dynamics (NSCD)
This section is devoted to the description of the method used to solve the rigid body dynamic contact problems [14]. Following the ideas of Jean and Moreau (see for example
the global frame to the local one for each contact node $\alpha \in [1, N_\alpha]$ (where $N_\alpha$ is the total number of potential contacts between two rigid particles) to write the corresponding local contact law. Local relative velocities $v_n^\alpha$ and $v_t^\alpha$ between the two rigid bodies involved in the contact $\alpha$ can be generalised as a vector $\dot{\mathbf{q}} = \mathbf{v} \in \mathbb{R}^{d \times N_c}$ (where $d = 3$ for a 3D problem and $d = 2$ for a 2D problem). Likewise, the local contact impulses $\mathbf{p}^\alpha$ are represented by a vector $\mathbf{P} \in \mathbb{R}^{d \times N_c}$. We recall that $\mathbf{p}^\alpha$ can be decomposed into the sum of a normal component $p_n^\alpha$ and a tangential component $p_t^\alpha$ as follows: $\mathbf{p}^\alpha = p_n^\alpha \mathbf{n} + \mathbf{p}_t^\alpha$. From there, the idea would be to transform the equations of motion seen in (2.11) from $\mathbf{P}$ and $\dot{\mathbf{q}}$ to $\mathbf{p}^\alpha$ and $\mathbf{v}^\alpha$.

The gap $D^\alpha = D(q_i, q_j)$ represents the signed distance from particle $i$ to particle $j$.

The contact and particles velocities ($\mathbf{v}^\alpha$ and $\dot{\mathbf{q}}$) being linearly related, the local-global mapping to compute $\mathbf{v}$ to $\dot{\mathbf{q}}$ at a contact node $\alpha$ is

$$ \mathbf{v}^\alpha = H^*(\mathbf{q}, \alpha)\dot{\mathbf{q}}, \qquad (2.12) $$

where $H^*(\mathbf{q}, \alpha)$ is a local-global mapping matrix $dN_c \times \bar{d}N_p$ carrying the information on the geometry of the contact network. Identically, we have for the impulses

$$ \mathbf{P} = H(\mathbf{q}, \alpha)\mathbf{p}^\alpha, \qquad (2.13) $$

where $H(\mathbf{q}, \alpha)$ is the transpose of matrix $H^*(\mathbf{q}, \alpha)$.

The inelastic shock law can be written as follows:

$$ \text{if } D^\alpha > 0 \text{ then } p_n^\alpha = 0, \qquad (2.14) $$

$$ \text{if } D^\alpha \le 0 \text{ then } \left\{ \begin{array}{l} v_n^\alpha \ge 0, \\ p_n^\alpha \ge 0, \\ v_n^\alpha p_n^\alpha = 0. \end{array} \right. \qquad (2.15) $$

Coulomb's law of friction is written as follows:

$$ \mathbf{p}_t^\alpha \neq \mathbf{0}, \qquad (2.16) $$

$$ \left\{ \begin{array}{l} p_n^\alpha = 0 \Rightarrow v_n^\alpha \ge 0, \\ p_n^\alpha > 0 \text{ and } \|\mathbf{p}_t^\alpha\| < \mu p_n^\alpha \Rightarrow \mathbf{v}_t^\alpha = 0, \\ p_n^\alpha > 0 \text{ and } |\mathbf{p}_t^\alpha| = \mu p_n^\alpha \Rightarrow \exists \beta \ge 0, \mathbf{v}_t^\alpha = \beta \frac{\mathbf{p}_t^\alpha}{\|\mathbf{p}_t^\alpha\|}, \end{array} \right. \qquad (2.17) $$

where $\mu$ is the friction coefficient.

The frictional contact laws (2.15) and (2.17) can be also formulated in subdifferential inclusions as

$$ \mathbf{p}^\alpha \in \partial I_{\mathbb{R}^-}(v_n^\alpha) \mathbf{n} + \partial I_{C(p_n^\alpha)}(v_t^\alpha). \qquad (2.18) $$

### 2.2.3 Nonlinear Gauss-Seidel (NLGS) method for nonsmooth contact dynamics (NSCD)
This section is devoted to the description of the method used to solve the rigid body dynamic contact problems [14]. Following the ideas of Jean and Moreau (see for example
[6, 8, 38]), we use the NLGS method, which consists in considering successively each contact until the convergence.

Considering that the time interval of interest is [0, T], we discretize for numerical purposes the previous equation (2.11). For that, we introduce uniform time instants $t^n$ defined by $t^{n+1} = t^n + \Delta t$ for $n = 0, \dots, N_T - 1$, where

* $\Delta t = T/N_T$ is the time step;

* $N_T$ is the total number of time steps.

Equation (2.11) is then integrated on each time interval $[t^n, t^{n+1}]$ and approximated using a $\theta$-scheme, with $\theta \in [\frac{1}{2}, 1]$ for stability reasons (see [9, 39]). Therefore, the classical approximation of equation (2.11) yields

$$ egin{cases} 
\mathbf{M}(\dot{\mathbf{q}}^{n+1} - \dot{\mathbf{q}}^n) = \Delta t(\theta \mathbf{F}^{n+1} + (1-\theta)\mathbf{F}_n) + \mathbf{P}^{n+1}, \\ 
\mathbf{q}^{n+1} = \mathbf{q}^n + \Delta t \theta \dot{\mathbf{q}}^{n+1} + \Delta t(1-\theta)\dot{\mathbf{q}}^n, 
\end{cases} 
\quad (2.19) $$

where

* $\mathbf{P}^{n+1}$ represents the value of the total impulsion over the time step $[t^n, t^{n+1}]$;

* $\mathbf{F}^n$ (resp. $\mathbf{F}^{n+1}$) is the external force computed at time $t^n$ (resp. $t^{n+1}$).

Denoting by $\dot{\mathbf{q}}^{n,\text{free}} = \dot{\mathbf{q}}^n + \mathbf{M}^{-1}\Delta t(\theta\mathbf{F}^{n+1} + (1-\theta)\mathbf{F}^n)$ the velocity when the contact impulses vanish, the first equation of (2.19) becomes

$$ \dot{\mathbf{q}}^{n+1} = \dot{\mathbf{q}}^{n,\text{free}} + \mathbf{M}^{-1}\mathbf{P}^{n+1}. \quad (2.20) $$

By combining (2.15), (2.12), (2.13), and (2.19), the discretization of the multi-contact system's motion with contact between rigid bodies can be written for each contact $\alpha$:

$$ egin{cases} 
\tilde{\mathbf{v}}^{\alpha,n+1} = \tilde{\mathbf{v}}^{\alpha,n,\text{free}} + \mathbf{W}\mathbf{p}^{\alpha,n+1}, \\ 
\text{contact\_law}(\tilde{\mathbf{v}}_{n}^{\alpha,n+1}, \mathbf{p}_{n}^{\alpha,n+1}) = \text{.true.} \quad \forall \alpha \in [1, \dots, N_{\alpha}], \\ 
\text{friction\_law}(\tilde{\mathbf{v}}_{n}^{\alpha,n+1}, \mathbf{p}_{n}^{\alpha,n+1}) = \text{.true.} \quad \forall \alpha \in [1, \dots, N_{\alpha}]. 
\end{cases} 
\quad (2.21) $$

The matrix $\mathbf{W}$ is called the Delassus matrix, which represents the mass matrix reduced to local variables.

The time-stepping method combined with the NLGS algorithm takes the following form:

* Loop on the step time $n$

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

$$ \mathbf{q}^{n+\frac{1}{2}} = \mathbf{q}^n + \frac{\Delta t}{2} \dot{\mathbf{q}}^n; \quad (2.22) $$

- Initialization of the motion: $\dot{\mathbf{q}}^{0,n+1} = \dot{\mathbf{q}}^{n,\text{free}}$ (initialization of the contact impulses with $\mathbf{P} = 0$).

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

* Loop on the contacts $\alpha$:

+ Computation of the local-global mapping

$$ \mathbf{v}^{-} = H^{*}(\mathbf{q}^{n+\frac{1}{2}}, \alpha)\dot{\mathbf{q}}^{n}; \quad (2.23) $$
[6, 8, 38]), we use the NLGS method, which consists in considering successively each contact until the convergence.

Considering that the time interval of interest is [0, T], we discretize for numerical purposes the previous equation (2.11). For that, we introduce uniform time instants $t^n$ defined by $t^{n+1} = t^n + \Delta t$ for $n = 0, \dots, N_T - 1$, where

* $\Delta t = T/N_T$ is the time step;

* $N_T$ is the total number of time steps.

Equation (2.11) is then integrated on each time interval $[t^n, t^{n+1}]$ and approximated using a $\theta$-scheme, with $\theta \in [\frac{1}{2}, 1]$ for stability reasons (see [9, 39]). Therefore, the classical approximation of equation (2.11) yields

$$ 
\begin{cases}
\mathbf{M}(\dot{\mathbf{q}}^{n+1} - \dot{\mathbf{q}}^n) = \Delta t(\theta \mathbf{F}^{n+1} + (1-\theta)\mathbf{F}_n) + \mathbf{P}^{n+1}, \\
\mathbf{q}^{n+1} = \mathbf{q}^n + \Delta t \theta \dot{\mathbf{q}}^{n+1} + \Delta t(1-\theta)\dot{\mathbf{q}}^n,
\end{cases}
\quad (2.19) $$

where

* $\mathbf{P}^{n+1}$ represents the value of the total impulsion over the time step $[t^n, t^{n+1}]$;

* $\mathbf{F}^n$ (resp. $\mathbf{F}^{n+1}$) is the external force computed at time $t^n$ (resp. $t^{n+1}$).

Denoting by $\dot{\mathbf{q}}^{n,\text{free}} = \dot{\mathbf{q}}^n + \mathbb{M}^{-1}\Delta t(\theta\mathbf{F}^{n+1} + (1-\theta)\mathbf{F}^n)$ the velocity when the contact impulses vanish, the first equation of (2.19) becomes

$$ \dot{\mathbf{q}}^{n+1} = \dot{\mathbf{q}}^{n,\text{free}} + \mathbb{M}^{-1}\mathbf{P}^{n+1}. \quad (2.20) $$

By combining (2.15), (2.12), (2.13), and (2.19), the discretization of the multi-contact system's motion with contact between rigid bodies can be written for each contact $\alpha$:

$$ 
\begin{cases}
\tilde{\mathbf{v}}^{\alpha,n+1} = \tilde{\mathbf{v}}^{\alpha,n,\text{free}} + \mathbb{W}\mathbf{p}^{\alpha,n+1}, \\
\text{contact\_law}(\tilde{\mathbf{v}}_{n}^{\alpha,n+1}, \mathbf{p}_{n}^{\alpha,n+1}) = .\text{true.} \quad \forall \alpha \in [1, \dots, N_{\alpha}], \\
\text{friction\_law}(\tilde{\mathbf{v}}_{n}^{\alpha,n+1}, \mathbf{p}_{n}^{\alpha,n+1}) = .\text{true.} \quad \forall \alpha \in [1, \dots, N_{\alpha}].
\end{cases}
\quad (2.21) $$

The matrix $\mathbb{W}$ is called the Delassus matrix, which represents the mass matrix reduced to local variables.

The time-stepping method combined with the NLGS algorithm takes the following form:

* Loop on the step time $n$

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

$$ \mathbf{q}^{n+\frac{1}{2}} = \mathbf{q}^n + \frac{\Delta t}{2} \dot{\mathbf{q}}^n; \quad (2.22) $$

- Initialization of the motion: $\dot{\mathbf{q}}^{0,n+1} = \dot{\mathbf{q}}^{n,\text{free}}$ (initialization of the contact impulses with $\mathbf{P} = 0$).

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

* Loop on the contacts $\alpha$:

+ Computation of the local-global mapping

$$ \mathbf{v}^{-} = H^{*}(\mathbf{q}^{n+\frac{1}{2}}, \alpha)\dot{\mathbf{q}}^{n}; \quad (2.23) $$
$$ \mathbf{v}^{\alpha,j+} = H^*(\mathbf{q}^{n+\frac{1}{2}}, \alpha) \dot{\mathbf{q}}^{j,n+1}. \qquad (2.24) $$

+ Newton shock law (using formal Moreau velocity)

$$ \tilde{\mathbf{v}}_{n}^{\alpha,j+1} = \frac{\mathbf{v}_{n}^{\alpha,j+} + e_{n} \mathbf{v}_{n}^{-}}{1 + e_{n}}; \qquad (2.25) $$

$$ \tilde{\mathbf{v}}_{t}^{\alpha,j+1} = \frac{\mathbf{v}_{t}^{\alpha,j+} + e_{t} \mathbf{v}_{t}^{-}}{1 + e_{t}}. \qquad (2.26) $$

+ Computation of the frictional contact law:

$$ \left\{ \begin{array}{l} \text{contact\_law}(\tilde{\mathbf{v}}^{\alpha,j+1}, \mathbf{p}^{\alpha,j+1}) = \text{.true.} \\ \text{friction\_law}(\tilde{\mathbf{v}}^{\alpha,j+1}, \mathbf{p}^{\alpha,j+1}) = \text{.true.} \end{array} \right. \qquad (2.27) $$

+ Actualisation of the generalised displacement:

$$ \dot{\mathbf{q}}^{j+1,n+1} = \dot{\mathbf{q}}^{j,n+1} + \mathbb{M}^{-1}P(\mathbf{q}^{n+\frac{1}{2}}, \alpha)\mathbf{p}^{\alpha,j+1}. $$

* End of the loop on contacts $\alpha$.

- End of the loop on $j$ of NLGS. When the convergence is reached, actualisation of the velocity: $\dot{\mathbf{q}}^{n+1} = \dot{\mathbf{q}}^{n+1,j+1}$.

- Actualisation of the generalised displacements: $\mathbf{q}^{n+1} = \mathbf{q}^{n+\frac{1}{2}} + \frac{\Delta t}{2} \dot{\mathbf{q}}^{n+1}$.

* End of the loop on the step time $k$.

### 3 Semi-smooth Newton approach for solving the frictional contact conditions

The semi-smooth Newton-PDAS methods appear to be one of the most relevant methods for solving friction contact problems (cf. [27, 29, 30]). These methods are based on the following principle: the conditions of contact and friction are reformulated in terms of nonlinear complementarity functions whose solution is provided by the semi-smooth iterative method of Newton [26, 27]. To this end, we need the generalised derivative of complementary functions for contact and friction. In practice, the conditions of contact with Coulomb's friction can be formulated in terms of a fixed point problem related to a quasi-optimisation one. From a purely algorithmic point of view, the main goal of these methods is to separate the nodes potentially in contact into two subsets (active and inactive) and to find the correct subset of all the nodes actually in active contact (subset $\mathcal{A}$), as opposed to those that are inactive (subset $\mathcal{I}$). In practice, the semi-smooth Newton-PDAS methods do not require the use of Lagrange multipliers. In fact, the boundary conditions on the subsets $\mathcal{A}$ and $\mathcal{I}$ are directly imposed thanks to a semi-smooth Newton method, and consequently, their implementation could be achieved without much effort.

In order to generalise the frictional contact conditions for deformable and rigid problems, we will make a change of variables. Let us take the following sthenic and kinematic variables:

$$ \text{For deformable problems: } \phi_n = \lambda_n^{\alpha}, \phi_t = \lambda_t^{\alpha}, w_n = u_n^{\alpha}, w_t = \dot{u}_t^{\alpha}, \qquad (3.1) $$

$$ \text{For rigid problems: } \phi_n = p_n^{\alpha}, \phi_t = \mathbf{p}_t^{\alpha}, w_n = v_n^{\alpha}, w_t = \mathbf{v}_t^{\alpha}. \qquad (3.2) $$
$$ \mathbf{v}^{\alpha,j^+} = H^*(\mathbf{q}^{n+\frac{1}{2}}, \alpha) \dot{\mathbf{q}}^{j,n+1}. \qquad (2.24) $$

+ Newton shock law (using formal Moreau velocity)

$$ \tilde{\mathbf{v}}_n^{\alpha,j+1} = \frac{\mathbf{v}_n^{\alpha,j^+} + e_n \mathbf{v}_n^-}{1 + e_n}; \qquad (2.25) $$

$$ \tilde{\mathbf{v}}_t^{\alpha,j+1} = \frac{\mathbf{v}_t^{\alpha,j^+} + e_t \mathbf{v}_t^-}{1 + e_t}. \qquad (2.26) $$

+ Computation of the frictional contact law:

$$ \left\{ \begin{array}{l} \text{contact\_law}(\tilde{\mathbf{v}}^{\alpha,j+1}, \mathbf{p}^{\alpha,j+1}) = \text{.true.} \\ \text{friction\_law}(\tilde{\mathbf{v}}^{\alpha,j+1}, \mathbf{p}^{\alpha,j+1}) = \text{.true.} \end{array} \right. \qquad (2.27) $$

+ Actualisation of the generalised displacement:

$$ \dot{\mathbf{q}}^{j+1,n+1} = \dot{\mathbf{q}}^{j,n+1} + \mathbb{M}^{-1}P(\mathbf{q}^{n+\frac{1}{2}}, \alpha)\mathbf{p}^{\alpha,j+1}. $$

* End of the loop on contacts α.

- End of the loop on j of NLGS. When the convergence is reached, actualisation of the velocity: $\dot{\mathbf{q}}^{n+1} = \dot{\mathbf{q}}^{n+1,j+1}$.

- Actualisation of the generalised displacements: $\mathbf{q}^{n+1} = \mathbf{q}^{n+\frac{1}{2}} + \frac{\Delta t}{2} \dot{\mathbf{q}}^{n+1}$.

* End of the loop on the step time k.

### 3 Semi-smooth Newton approach for solving the frictional contact conditions

The semi-smooth Newton-PDAS methods appear to be one of the most relevant methods for solving friction contact problems (cf. [27, 29, 30]). These methods are based on the following principle: the conditions of contact and friction are reformulated in terms of nonlinear complementarity functions whose solution is provided by the semi-smooth iterative method of Newton [26, 27]. To this end, we need the generalised derivative of complementary functions for contact and friction. In practice, the conditions of contact with Coulomb's friction can be formulated in terms of a fixed point problem related to a quasi-optimisation one. From a purely algorithmic point of view, the main goal of these methods is to separate the nodes potentially in contact into two subsets (active and inactive) and to find the correct subset of all the nodes actually in active contact (subset A), as opposed to those that are inactive (subset I). In practice, the semi-smooth Newton-PDAS methods do not require the use of Lagrange multipliers. In fact, the boundary conditions on the subsets A and I are directly imposed thanks to a semi-smooth Newton method, and consequently, their implementation could be achieved without much effort.

In order to generalise the frictional contact conditions for deformable and rigid problems, we will make a change of variables. Let us take the following sthenic and kinematic variables:

$$ \text{For deformable problems: } \phi_n = \lambda_n^\alpha, \phi_t = \lambda_t^\alpha, w_n = u_n^\alpha, \mathbf{w}_t = \dot{\mathbf{u}}_t^\alpha, \qquad (3.1) $$

$$ \text{For rigid problems: } \phi_n = p_n^\alpha, \phi_t = \mathbf{p}_t^\alpha, w_n = v_n, \mathbf{w}_t = \mathbf{v}_t^\alpha. \qquad (3.2) $$
## 3.1 Complementary function for frictional contact

Signorini's conditions (2.2) or (2.15) can be formulated from the following nonlinear complementary function $C_n^\phi(w_n, \phi_n)$:

$$C_n^\phi(w_n, \phi_n) = \phi_n - [\phi_n - \gamma_n w_n]_+ \quad \forall \alpha \in S, \qquad (3.3)$$

where $S$ is the set of all potential contact nodes and $\gamma_n > 0$ is the normal active set parameter. Let us state a first result.

**Proposition 3.1** The unilateral contact conditions (2.2) or (2.15) for each potential contact are equivalent to $C_n^\phi(w_n, \phi_n) = 0$, where $\phi_n$ is the normal reaction (or impulse) force.

For the proof, see [35] for deformable problems and see [34] for rigid problems.

Frictional Coulomb's conditions (2.3) or (2.17) can be formulated from the following nonlinear complementary function $C_t^\phi(w_n, \phi_t, \phi_n, \phi_t)$:

$$C_t^\phi(w_n, \mathbf{w}_t, \phi_n, \mathbf{\phi}_t) = \max(\mu\phi_n, \|\mathbf{\phi}_t - \gamma_t \mathbf{w}_t\|)\mathbf{\phi}_t - \mu\phi_n(\mathbf{\phi}_t - \gamma_t \mathbf{w}_t) \quad \forall \alpha \in S, \qquad (3.4)$$

where $\gamma_t > 0$ is the tangential active set parameter. Then, let us present this result.

**Proposition 3.2** The frictional contact conditions (2.3) or (2.17) for each contact are equivalent to $C_t^\phi(w_n, \mathbf{w}_t, \phi_n, \mathbf{\phi}_t) = 0$, where $\mathbf{\phi}_t$ is the tangential reaction (or impulse) force between two particles in contact.

For the proof, see [34] for deformable problems and see [35] for rigid problems.

## 3.2 Generalised derivative of complementary functions

Now, we provide the generalised derivative of the complementary functions in the gap, stick and slip cases.

* Gap case: $\phi_n - \gamma_n w_n \le 0$

According to the complementary functions $C_n^\phi(w_n, \phi_n) = \phi_n$ and $C_t^\phi(w_n, \mathbf{w}_t, \phi_n, \mathbf{\phi}_t) = \|\mathbf{\phi}_t - \gamma_t \mathbf{w}_t\|\mathbf{\phi}_t$, we have the following derivatives:

$$d_{w_n} C_n^\phi = 0, \qquad (3.5)$$

$$d_{\phi_n} C_n^\phi = d\phi_n, \qquad (3.6)$$

$$d_{w_n} C_t^\phi = 0, \qquad (3.7)$$

$$d_{\mathbf{w}_t} C_t^\phi = -\gamma_t \mathbf{\phi}_t \frac{(\mathbf{\phi}_t - \gamma_t \mathbf{w}_t)^T}{\|\mathbf{\phi}_t - \gamma_t \mathbf{w}_t\|} d\mathbf{w}_t = 0, \qquad (3.8)$$

$$d_{\phi_n} C_t^\phi = 0, \qquad (3.9)$$

$$d_{\phi_t} C_t^\phi = \left( \mathbf{\phi}_t \frac{(\mathbf{\phi}_t - \gamma_t \mathbf{w}_t)^T}{\|\mathbf{\phi}_t - \gamma_t \mathbf{w}_t\|} + \|\mathbf{\phi}_t - \gamma_t \mathbf{w}_t\| I_2 \right) d\mathbf{\phi}_t. \qquad (3.10)$$

* Stick case: $\mu\phi_n \ge \|\mathbf{\phi}_t - \gamma_t \mathbf{w}_t\| > 0$

Given the complementary functions $C_n^\phi(w_n, \phi_n) = \gamma_n w_n$ and $C_t^\phi(w_n, \mathbf{w}_t, \phi_n, \mathbf{\phi}_t) = \mu\gamma_t \phi_n \mathbf{w}_t$, we have

$$d_{w_n} C_n^\phi = \gamma_n d w_n, \qquad (3.11)$$
## 3.1 Complementary function for frictional contact

Signorini's conditions (2.2) or (2.15) can be formulated from the following nonlinear complementary function $C_n^\phi(w_n, \phi_n)$:

$$C_n^\phi(w_n, \phi_n) = \phi_n - [\phi_n - \gamma_n w_n]_+ \quad \forall \alpha \in S, \quad (3.3)$$

where $S$ is the set of all potential contact nodes and $\gamma_n > 0$ is the normal active set parameter. Let us state a first result.

**Proposition 3.1** The unilateral contact conditions (2.2) or (2.15) for each potential contact are equivalent to $C_n^\phi(w_n, \phi_n) = 0$, where $\phi_n$ is the normal reaction (or impulse) force.

For the proof, see [35] for deformable problems and see [34] for rigid problems.

Frictional Coulomb's conditions (2.3) or (2.17) can be formulated from the following nonlinear complementary function $C_t^\phi(w_n, \phi_t, \phi_n, \phi_t)$:

$$C_t^\phi(w_n, \mathbf{w}_t, \phi_n, \phi_t) = \max(\mu\phi_n, \|\phi_t - \gamma_t \mathbf{w}_t\|)\phi_t - \mu\phi_n(\phi_t - \gamma_t \mathbf{w}_t) \quad \forall \alpha \in S, \quad (3.4)$$

where $\gamma_t > 0$ is the tangential active set parameter. Then, let us present this result.

**Proposition 3.2** The frictional contact conditions (2.3) or (2.17) for each contact are equivalent to $C_t^\phi(w_n, \mathbf{w}_t, \phi_n, \phi_t) = 0$, where $\phi_t$ is the tangential reaction (or impulse) force between two particles in contact.

For the proof, see [34] for deformable problems and see [35] for rigid problems.

## 3.2 Generalised derivative of complementary functions

Now, we provide the generalised derivative of the complementary functions in the gap, stick and slip cases.

• Gap case: $\phi_n - \gamma_n w_n \le 0$

According to the complementary functions $C_n^\phi(w_n, \phi_n) = \phi_n$ and $C_t^\phi(w_n, \mathbf{w}_t, \phi_n, \phi_t) = \|\phi_t - \gamma_t \mathbf{w}_t\|\phi_t$, we have the following derivatives:

$$d_{w_n} C_n^\phi = 0, \quad (3.5)$$

$$d_{\phi_n} C_n^\phi = d\phi_n, \quad (3.6)$$

$$d_{w_n} C_t^\phi = 0, \quad (3.7)$$

$$d_{\mathbf{w}_t} C_t^\phi = -\gamma_t \phi_t \frac{(\phi_t - \gamma_t \mathbf{w}_t)^T}{\|\phi_t - \gamma_t \mathbf{w}_t\|} d\mathbf{w}_t = 0, \quad (3.8)$$

$$d_{\phi_n} C_t^\phi = 0, \quad (3.9)$$

$$d_{\phi_t} C_t^\phi = \left( \phi_t \frac{(\phi_t - \gamma_t \mathbf{w}_t)^T}{\|\phi_t - \gamma_t \mathbf{w}_t\|} + \|\phi_t - \gamma_t \mathbf{w}_t\| I_2 \right) d\phi_t. \quad (3.10)$$

• Stick case: $\mu\phi_n \ge \|\phi_t - \gamma_t \mathbf{w}_t\| > 0$

Given the complementary functions $C_n^\phi(w_n, \phi_n) = \gamma_n w_n$ and $C_t^\phi(w_n, \mathbf{w}_t, \phi_n, \phi_t) = \mu\gamma_t\phi_n\mathbf{w}_t$, we have

$$d_{w_n} C_n^\phi = \gamma_n d w_n, \quad (3.11)$$
$$d_{\phi_n} C_n^{\phi} = 0, \qquad (3.12)$$

$$d_{w_n} C_t^{\phi} = 0, \qquad (3.13)$$

$$d_{w_t} C_t^{\phi} = \mu \gamma_t \phi_n d w_t, \qquad (3.14)$$

$$d_{\phi_n} C_t^{\phi} = \mu \gamma_t w_t d\phi_n, \qquad (3.15)$$

$$d_{\phi_t} C_t^{\phi} = 0. \qquad (3.16)$$

• Slip case: $||\phi_t - \gamma_t w_t|| > \mu \phi_n > 0$

From $C_n^{\phi}(w_n, \phi_n) = \gamma_n w_n$ and $C_t^{\phi}(w_n, w_t, \phi_n, \phi_t) = ||\phi_t - \gamma_t w_t||\phi_t - \mu \phi_n(\phi_t - \gamma_t w_t)$, it comes

$$d_{w_n} C_n^{\phi} = \gamma_n dw_n, \qquad (3.17)$$

$$d_{\phi_n} C_n^{\phi} = 0, \qquad (3.18)$$

$$d_{w_n} C_t^{\phi} = 0, \qquad (3.19)$$

$$d_{w_t} C_t^{\phi} = \left( -\gamma_t \phi_t \frac{(\phi_t - \gamma_t w_t)^T}{\|\phi_t - \gamma_t w_t\|} + \mu \gamma_t \phi_n I_2 \right) dw_t, \qquad (3.20)$$

$$d_{\phi_n} C_t^{\phi} = -\mu (\phi_t - \gamma_t w_t) d\phi_n, \qquad (3.21)$$

$$d_{\phi_t} C_t^{\phi} = \left( \phi_t \frac{(\phi_t - \gamma_t w_t)^T}{\|\phi_t - \gamma_t w_t\|} + \|\phi_t - \gamma_t w_t\| I_2 - \mu \phi_n I_2 \right) d\phi_t. \qquad (3.22)$$

By combining (3.5)–(3.22) with $D_{C_n^{\phi}}$ and $D_{C_t^{\phi}}$ the generalised derivatives of $C_n^{\phi}$ and $C_t^{\phi}$, respectively, we obtain

$$D_{C_n^{\phi}}(w_n, \phi_n)(\delta w_n, \delta \phi_n) = \gamma_n(1_{\text{Stick}} + 1_{\text{Slip}})\delta w_n + 1_{\text{Gap}}\delta \phi_n, \qquad (3.23)$$

$$D_{C_t^{\phi}}(w_n, w_t, \phi_n, \phi_t)(\delta w_n, \delta w_t, \delta \phi_n, \delta \phi_t) \\ = 1_{\text{Gap}} \|\phi_t - \gamma_t w_t\| \delta \phi_t + 1_{\text{Stick}} (\mu \gamma_t \phi_n \delta w_t + \mu \gamma_t w_t \delta \phi_n) \\ + 1_{\text{Slip}} \left( \left( -\gamma_t \phi_t \frac{(\phi_t - \gamma_t w_t)^T}{\|\phi_t - \gamma_t w_t\|} + \mu \gamma_t \phi_n I_2 \right) \delta w_t - \mu (\phi_t - \gamma_t w_t) \delta \phi_n \right) \\ + \left( \phi_t \frac{(\phi_t - \gamma_t w_t)^T}{\|\phi_t - \gamma_t w_t\|} + \|\phi_t - \gamma_t w_t\| I_2 - \mu \phi_n I_2 \right) \delta \phi_t, \qquad (3.24)$$

where

$$1_{\text{Gap}} = 1, \quad 1_{\text{Stick}} = 0, \quad 1_{\text{Slip}} = 0 \quad \text{if } \phi_n - \gamma_n w_n \le 0,$$

$$1_{\text{Gap}} = 0, \quad 1_{\text{Stick}} = 1, \quad 1_{\text{Slip}} = 0 \quad \text{if } \mu \phi_n \ge \|\phi_t - \gamma_t w_t\| > 0,$$

$$1_{\text{Gap}} = 0, \quad 1_{\text{Stick}} = 0, \quad 1_{\text{Slip}} = 1 \quad \text{if } \|\phi_t - \gamma_t w_t\| > \mu \phi_n > 0.$$

### 3.3 Fixed point conditions from Newton's semi-smooth approach

Using now the semi-smooth Newton formalism (indexed by the subscript k) at the current fixed point $(w_n^k, w_t^k, \phi_n^k, \phi_t^k)$ of the complementary functions $C_n^{\phi}$ and $C_t^{\phi}$, one can derive the new iterate $(w_n^{k+1}, w_t^{k+1}, \phi_n^{k+1}, \phi_t^{k+1})$

$$D_{C_n^{\phi}}(w_n^k, \phi_n^k)(\delta w_n^{k+1}, \delta \phi_n^{k+1}) = -C_n^{\phi}(w_n^k, \phi_n^k), \qquad (3.25)$$
$$d_{\phi_n} C_n^{\phi} = 0, \qquad (3.12)$$

$$d_{w_n} C_n^{\phi} = 0, \qquad (3.13)$$

$$d_{w_t} C_t^{\phi} = \mu \gamma_t \phi_n d w_t, \qquad (3.14)$$

$$d_{\phi_n} C_n^{\phi} = \mu \gamma_t w_t d\phi_n, \qquad (3.15)$$

$$d_{\phi_t} C_t^{\phi} = 0. \qquad (3.16)$$

• Slip case: $||\phi_t - \gamma_t w_t|| > \mu \phi_n > 0$

From $C_n^{\phi}(w_n, \phi_n) = \gamma_n w_n$ and $C_t^{\phi}(w_n, w_t, \phi_n, \phi_t) = ||\phi_t - \gamma_t w_t||\phi_t - \mu \phi_n(\phi_t - \gamma_t w_t)$, it comes

$$d_{w_n} C_n^{\phi} = \gamma_n d w_n, \qquad (3.17)$$

$$d_{\phi_n} C_n^{\phi} = 0, \qquad (3.18)$$

$$d_{w_n} C_n^{\phi} = 0, \qquad (3.19)$$

$$d_{w_t} C_t^{\phi} = \left( -\gamma_t \phi_t \frac{(\phi_t - \gamma_t w_t)^T}{\|\phi_t - \gamma_t w_t\|} + \mu \gamma_t \phi_n \mathbf{I}_2 \right) d w_t, \qquad (3.20)$$

$$d_{\phi_n} C_n^{\phi} = -\mu(\phi_t - \gamma_t w_t) d\phi_n, \qquad (3.21)$$

$$d_{\phi_t} C_t^{\phi} = \left( \phi_t \frac{(\phi_t - \gamma_t w_t)^T}{\|\phi_t - \gamma_t w_t\|} + \|\phi_t - \gamma_t w_t\| \mathbf{I}_2 - \mu \phi_n \mathbf{I}_2 \right) d\phi_t. \qquad (3.22)$$

By combining (3.5)–(3.22) with $D_{C_n^{\phi}}$ and $D_{C_t^{\phi}}$ the generalised derivatives of $C_n^{\phi}$ and $C_t^{\phi}$, respectively, we obtain

$$D_{C_n^{\phi}}(w_n, \phi_n)(\delta w_n, \delta \phi_n) = \gamma_n(1_{\text{Stick}} + 1_{\text{Slip}})\delta w_n + 1_{\text{Gap}}\delta \phi_n, \qquad (3.23)$$

$$D_{C_t^{\phi}}(w_n, w_t, \phi_n, \phi_t)(\delta w_n, \delta w_t, \delta \phi_n, \delta \phi_t) \\ = 1_{\text{Gap}} \|\phi_t - \gamma_t w_t\| \delta \phi_t + 1_{\text{Stick}} (\mu \gamma_t \phi_n \delta w_t + \mu \gamma_t w_t \delta \phi_n) \\ + 1_{\text{Slip}} \left( \left( -\gamma_t \phi_t \frac{(\phi_t - \gamma_t w_t)^T}{\|\phi_t - \gamma_t w_t\|} + \mu \gamma_t \phi_n \mathbf{I}_2 \right) \delta w_t - \mu (\phi_t - \gamma_t w_t) \delta \phi_n \right) \\ + \left( \phi_t \frac{(\phi_t - \gamma_t w_t)^T}{\|\phi_t - \gamma_t w_t\|} + \|\phi_t - \gamma_t w_t\| \mathbf{I}_2 - \mu \phi_n \mathbf{I}_2 \right) \delta \phi_t, \qquad (3.24)$$

where

$$1_{\text{Gap}} = 1, \quad 1_{\text{Stick}} = 0, \quad 1_{\text{Slip}} = 0 \quad \text{if } \phi_n - \gamma_n w_n \le 0,$$

$$1_{\text{Gap}} = 0, \quad 1_{\text{Stick}} = 1, \quad 1_{\text{Slip}} = 0 \quad \text{if } \mu \phi_n \ge \|\phi_t - \gamma_t w_t\| > 0,$$

$$1_{\text{Gap}} = 0, \quad 1_{\text{Stick}} = 0, \quad 1_{\text{Slip}} = 1 \quad \text{if } \|\phi_t - \gamma_t w_t\| > \mu \phi_n > 0.$$

### 3.3 Fixed point conditions from Newton’s semi-smooth approach

Using now the semi-smooth Newton formalism (indexed by the subscript $k$) at the current fixed point $(w_n^k, w_t^k, \phi_n^k, \phi_t^k)$ of the complementary functions $C_n^{\phi}$ and $C_t^{\phi}$, one can derive the new iterate $(w_n^{k+1}, w_t^{k+1}, \phi_n^{k+1}, \phi_t^{k+1})$

$$D_{C_n^{\phi}}(w_n^k, \phi_n^k)(\delta w_n^{k+1}, \delta \phi_n^{k+1}) = -C_n^{\phi}(w_n^k, \phi_n^k), \qquad (3.25)$$
$$ \mathcal{D}_{C_t^\phi} (w_n^k, \mathbf{w}_t^k, \phi_n^k, \phi_t^k) (\delta w_n^{k+1}, \delta \mathbf{w}_t^{k+1}, \delta \phi_n^{k+1}, \delta \phi_t^{k+1}) \\ = -C_t^\phi (w_n^k, \mathbf{w}_t^k, \phi_n^k, \phi_t^k), \qquad (3.26) $$

$$ (w_n^{k+1}, \mathbf{w}_t^{k+1}, \phi_n^{k+1}, \phi_t^{k+1}) = (w_n^k, \mathbf{w}_t^k, \phi_n^k, \phi_t^k) + (\delta w_n^{k+1}, \delta \mathbf{w}_t^{k+1}, \delta \phi_n^{k+1}, \delta \phi_t^{k+1}). \qquad (3.27) $$

• Gap case: $1_{\text{Gap}} = 1, 1_{\text{Stick}} = 0, 1_{\text{Slip}} = 0$

From equations (3.25) and (3.26) we have

$$ \phi_n^{k+1} - \phi_n^k = -\phi_n^k, \qquad (3.28) $$

$$ \| \boldsymbol{\phi}_t^k - \gamma_t \mathbf{w}_t^k \| (\boldsymbol{\phi}_t^{k+1} - \boldsymbol{\phi}_t^k) = - \| \boldsymbol{\phi}_t^k - \gamma_t \mathbf{w}_t^k \| \boldsymbol{\phi}_t^k. \qquad (3.29) $$

Next, the gap conditions of the semi-smooth Newton formalism are as follows:

$$ \phi_n^{k+1} = 0, \qquad (3.30) $$

$$ \boldsymbol{\phi}_t^{k+1} = \mathbf{0}, \qquad (3.31) $$

since $\|\boldsymbol{\phi}_t^k - \gamma_t \mathbf{w}_t^k\| > 0$.

• Stick case: $1_{\text{Gap}} = 0, 1_{\text{Stick}} = 1, 1_{\text{Slip}} = 0$

From equations (3.25) and (3.26) we have

$$ \gamma_n (w_n^{k+1} - w_n^k) = -\gamma_n w_n^k, \qquad (3.32) $$

$$ \mu \gamma_n \phi_n^k (\mathbf{w}_t^{k+1} - \mathbf{w}_t^k) + \mu \gamma_t \mathbf{w}_t^k (\phi_n^{k+1} - \phi_n^k) = -\mu \gamma_t \phi_n^k \mathbf{w}_t^k. \qquad (3.33) $$

Next,

$$ w_n^{k+1} = 0, \qquad (3.34) $$

$$ \mathbf{w}_t^{k+1} - \mathbf{w}_t^k = -\mathbf{w}_t^k \frac{\phi_n^{k+1}}{\phi_n^k}. \qquad (3.35) $$

For dynamics of rigid bodies and for a given contact node $\alpha$, the fundamental principle of dynamics can be written as follows:

$$ \mathbf{w}^{\alpha} = \mathbf{w}^{\alpha, \text{free}} + \mathcal{W}^{\alpha \alpha} \boldsymbol{\phi}^{\alpha} + \sum_{\beta \neq \alpha} \mathcal{W}^{\beta \alpha} \boldsymbol{\phi}^{\beta}. \qquad (3.36) $$

Finally, the stick conditions for rigid bodies of the semi-smooth Newton formalism are as follows:

$$ w_n^{k+1} = 0, \qquad (3.37) $$

$$ \boldsymbol{\phi}_t^{k+1} + \frac{\phi_n^{k+1} \mathbf{w}_t^k}{\mathcal{W}_{tt}^{\alpha\alpha}} = \boldsymbol{\phi}_t^k. \qquad (3.38) $$

• Slip case: $1_{\text{Gap}} = 0, 1_{\text{Stick}} = 0, 1_{\text{Slip}} = 1$
$$ \begin{aligned} & \mathcal{D}_{C_t^\phi} (w_n^k, \mathbf{w}_t^k, \phi_n^k, \boldsymbol{\phi}_t^k) (\delta w_n^{k+1}, \delta \mathbf{w}_t^{k+1}, \delta \phi_n^{k+1}, \delta \boldsymbol{\phi}_t^{k+1}) \\ & = -C_t^\phi (w_n^k, \mathbf{w}_t^k, \phi_n^k, \boldsymbol{\phi}_t^k), \end{aligned} \quad (3.26) $$

$$ \begin{aligned} & (w_n^{k+1}, \mathbf{w}_t^{k+1}, \phi_n^{k+1}, \boldsymbol{\phi}_t^{k+1}) \\ & = (w_n^k, \mathbf{w}_t^k, \phi_n^k, \boldsymbol{\phi}_t^k) + (\delta w_n^{k+1}, \delta \mathbf{w}_t^{k+1}, \delta \phi_n^{k+1}, \delta \boldsymbol{\phi}_t^{k+1}). \end{aligned} \quad (3.27) $$

• Gap case: $1_{\text{Gap}} = 1, 1_{\text{Stick}} = 0, 1_{\text{Slip}} = 0$

From equations (3.25) and (3.26) we have

$$ \phi_n^{k+1} - \phi_n^k = -\phi_n^k, \quad (3.28) $$

$$ \| \boldsymbol{\phi}_t^k - \gamma_t \mathbf{w}_t^k \| (\boldsymbol{\phi}_t^{k+1} - \boldsymbol{\phi}_t^k) = - \| \boldsymbol{\phi}_t^k - \gamma_t \mathbf{w}_t^k \| \boldsymbol{\phi}_t^k. \quad (3.29) $$

Next, the gap conditions of the semi-smooth Newton formalism are as follows:

$$ \phi_n^{k+1} = 0, \quad (3.30) $$

$$ \boldsymbol{\phi}_t^{k+1} = \mathbf{0}, \quad (3.31) $$

since $\|\boldsymbol{\phi}_t^k - \gamma_t \mathbf{w}_t^k\| > 0$.

• Stick case: $1_{\text{Gap}} = 0, 1_{\text{Stick}} = 1, 1_{\text{Slip}} = 0$

From equations (3.25) and (3.26) we have

$$ \gamma_n (w_n^{k+1} - w_n^k) = -\gamma_n w_n^k, \quad (3.32) $$

$$ \mu \gamma_n \phi_n^k (\mathbf{w}_t^{k+1} - \mathbf{w}_t^k) + \mu \gamma_t \mathbf{w}_t^k (\phi_n^{k+1} - \phi_n^k) = -\mu \gamma_t \phi_n^k \mathbf{w}_t^k. \quad (3.33) $$

Next,

$$ w_n^{k+1} = 0, \quad (3.34) $$

$$ \mathbf{w}_t^{k+1} - \mathbf{w}_t^k = -\mathbf{w}_t^k \frac{\phi_n^{k+1}}{\phi_n^k}. \quad (3.35) $$

For dynamics of rigid bodies and for a given contact node $\alpha$, the fundamental principle of dynamics can be written as follows:

$$ \mathbf{w}^{\alpha} = \mathbf{w}^{\alpha, \text{free}} + \mathcal{W}^{\alpha \alpha} \boldsymbol{\phi}^{\alpha} + \sum_{\beta \neq \alpha} \mathcal{W}^{\beta \alpha} \boldsymbol{\phi}^{\beta}. \quad (3.36) $$

Finally, the stick conditions for rigid bodies of the semi-smooth Newton formalism are as follows:

$$ w_n^{k+1} = 0, \quad (3.37) $$

$$ \boldsymbol{\phi}_t^{k+1} + \frac{\phi_n^{k+1} \mathbf{w}_t^k}{\mathcal{W}_{tt}^{\alpha\alpha}} = \boldsymbol{\phi}_t^k. \quad (3.38) $$

• Slip case: $1_{\text{Gap}} = 0, 1_{\text{Stick}} = 0, 1_{\text{Slip}} = 1$
For $C_n^\phi$, we get once again

$$w_n^{k+1} = 0. \qquad (3.39)$$

For $C_n^\phi$, we obtain

$$ \begin{aligned} & \left(-\gamma_t \boldsymbol{\phi}_t^k \frac{(\boldsymbol{\phi}_t^k - \gamma_t \boldsymbol{w}_t^k)^T}{\|\boldsymbol{\phi}_t^k - \gamma_t \boldsymbol{w}_t^k\|} + \mu \gamma_t \boldsymbol{\phi}_n^k \mathbf{I}_2\right) (\boldsymbol{w}_t^{k+1} - \boldsymbol{w}_t^k) \\ & - \mu (\boldsymbol{\phi}_t^k - \gamma_t \boldsymbol{w}_t^k) (\boldsymbol{\phi}_n^{k+1} - \boldsymbol{\phi}_n^k) \\ & + \left(\boldsymbol{\phi}_t^k \frac{(\boldsymbol{\phi}_t^k - \gamma_t \boldsymbol{w}_t^k)^T}{\|\boldsymbol{\phi}_t^k - \gamma_t \boldsymbol{w}_t^k\|} + \|\boldsymbol{\phi}_t^k - \gamma_t \boldsymbol{w}_t^k\| \mathbf{I}_2 - \mu \boldsymbol{\phi}_n^k \mathbf{I}_2\right) (\boldsymbol{\phi}_t^{k+1} - \boldsymbol{\phi}_t^k) \\ & = - \|\boldsymbol{\phi}_t^k - \gamma_t \boldsymbol{w}_t^k\| \boldsymbol{\phi}_t^k + \mu \boldsymbol{\phi}_n^k (\boldsymbol{\phi}_t^k - \gamma_t \boldsymbol{w}_t^k). \end{aligned} \qquad (3.40) $$

Let us introduce

$$ \mathbf{F}^{(k)} = \frac{\boldsymbol{\phi}_t^k (\boldsymbol{\phi}_t^k - \gamma_t \boldsymbol{w}_t^k)^T}{\|\boldsymbol{\phi}_t^k - \gamma_t \boldsymbol{w}_t^k\|}, \qquad E^{(k)} = \frac{1}{\|\boldsymbol{\phi}_t^k - \gamma_t \boldsymbol{w}_t^k\|}. \qquad (3.41) $$

Therefore, after an elementary computation

$$ - \gamma_t E^{(k)} (\mathbf{F}^{(k)} - \mu \boldsymbol{\phi}_n^k \mathbf{I}_2) (\boldsymbol{w}_t^{k+1} - \boldsymbol{w}_t^k) - \mu E^{(k)} (\boldsymbol{\phi}_t^k - \gamma_t \boldsymbol{w}_t^k) \boldsymbol{\phi}_n^{k+1} \\ + (E^{(k)} (\mathbf{F}^{(k)} - \mu \boldsymbol{\phi}_n^k \mathbf{I}_2) + \mathbf{I}_2) \boldsymbol{\phi}_t^{k+1} - E^{(k)} (\mathbf{F}^{(k)} - \mu \boldsymbol{\phi}_n^k \mathbf{I}_2) \boldsymbol{\phi}_t^k = 0. $$

Now, let us introduce the following operators:

$$ \mathbf{M}_{\alpha}^{*(k)} = E^{(k)} (\mathbf{F}^{(k)} - \mu \boldsymbol{\phi}_n^k \mathbf{I}_2), $$

$$ \mathbf{h}_{\alpha}^{(k)} = E^{(k)} \mathbf{F}^{(k)} (\boldsymbol{\phi}_t^k - \gamma_t \boldsymbol{w}_t^k) = \boldsymbol{\phi}_t^k. $$

$$ \mathbf{L}_{\alpha}^{*(k)} = -\gamma_t (\mathbf{I}_2 + \mathbf{M}_{\alpha}^{*(k)})^{-1} \mathbf{M}_{\alpha}^{*(k)}, $$

$$ \mathbf{r}_{\alpha}^{*(k)} = (\mathbf{I}_2 + \mathbf{M}_{\alpha}^{*(k)})^{-1} \mathbf{h}_{\alpha}^{(k)}, $$

$$ \mathbf{v}_{\alpha}^{(k)} = \mu (\mathbf{I}_2 + \mathbf{M}_{\alpha}^{*(k)})^{-1} E^{(k)} (\boldsymbol{\phi}_t^k - \gamma_t \boldsymbol{w}_t^k). $$

And, at last after computation, we have

$$ \mathbf{L}_{\alpha}^{*(k)} \boldsymbol{w}_t^{k+1} - \mathbf{v}_{\alpha}^{(k)} \boldsymbol{\phi}_n^{k+1} + \boldsymbol{\phi}_t^{k+1} = \mathbf{r}_{\alpha}^{*(k)} - \mathbf{v}_{\alpha}^{(k)} \boldsymbol{\phi}_n^k. $$

In this specific problem, and for a two-dimensional case, one can obtain a simplified equiv- alent version of the algorithm. Let $\mathcal{D}_{C_t^\phi}^{\text{slip}}$ be the generalised derivative of $C_t^\phi$ in the slip case

$$ \mathcal{D}_{C_t^\phi}^{\text{slip}} (\boldsymbol{w}_t^\alpha, \boldsymbol{\phi}_t) (\delta \boldsymbol{w}_t^\alpha, \delta \boldsymbol{\phi}_t) = \boldsymbol{\phi}_t \frac{(\boldsymbol{\phi}_t - \gamma_t \boldsymbol{w}_t^\alpha)^T}{\|\boldsymbol{\phi}_t - \gamma_t \boldsymbol{w}_t^\alpha\|} (\delta \boldsymbol{\phi}_t - \delta \boldsymbol{w}_t^\alpha) \\ - \mu \boldsymbol{\phi}_n (\delta \boldsymbol{\phi}_t - \gamma_t \boldsymbol{w}_t^\alpha) + \|\boldsymbol{\phi}_t - \gamma_t \boldsymbol{w}_t^\alpha\| \delta \boldsymbol{\phi}_t. \qquad (3.42) $$
For $C_n^\phi$, we get once again

$$w_n^{k+1} = 0. \qquad (3.39)$$

For $C_n^\phi$, we obtain

$$ \begin{aligned} & \left(-\gamma_t \boldsymbol{\phi}_t^k \frac{(\boldsymbol{\phi}_t^k - \gamma_t \boldsymbol{w}_t^k)^T}{\|\boldsymbol{\phi}_t^k - \gamma_t \boldsymbol{w}_t^k\|} + \mu \gamma_t \boldsymbol{\phi}_n^k \mathbf{I}_2\right) (\boldsymbol{w}_t^{k+1} - \boldsymbol{w}_t^k) \\ & \quad - \mu (\boldsymbol{\phi}_t^k - \gamma_t \boldsymbol{w}_t^k) (\boldsymbol{\phi}_n^{k+1} - \boldsymbol{\phi}_n^k) \\ & \quad + \left(\boldsymbol{\phi}_t^k \frac{(\boldsymbol{\phi}_t^k - \gamma_t \boldsymbol{w}_t^k)^T}{\|\boldsymbol{\phi}_t^k - \gamma_t \boldsymbol{w}_t^k\|} + \|\boldsymbol{\phi}_t^k - \gamma_t \boldsymbol{w}_t^k\| \mathbf{I}_2 - \mu \boldsymbol{\phi}_n^k \mathbf{I}_2\right) (\boldsymbol{\phi}_t^{k+1} - \boldsymbol{\phi}_t^k) \\ & = -\|\boldsymbol{\phi}_t^k - \gamma_t \boldsymbol{w}_t^k\| \boldsymbol{\phi}_t^k + \mu \boldsymbol{\phi}_n^k (\boldsymbol{\phi}_t^k - \gamma_t \boldsymbol{w}_t^k). \end{aligned} \qquad (3.40) $$

Let us introduce

$$ \mathbf{F}^{(k)} = \frac{\boldsymbol{\phi}_t^k (\boldsymbol{\phi}_t^k - \gamma_t \boldsymbol{w}_t^k)^T}{\|\boldsymbol{\phi}_t^k - \gamma_t \boldsymbol{w}_t^k\|}, \qquad E^{(k)} = \frac{1}{\|\boldsymbol{\phi}_t^k - \gamma_t \boldsymbol{w}_t^k\|}. \qquad (3.41) $$

Therefore, after an elementary computation

$$ \begin{aligned} & -\gamma_t E^{(k)} (\mathbf{F}^{(k)} - \mu \boldsymbol{\phi}_n^k \mathbf{I}_2) (\boldsymbol{w}_t^{k+1} - \boldsymbol{w}_t^k) - \mu E^{(k)} (\boldsymbol{\phi}_t^k - \gamma_t \boldsymbol{w}_t^k) \boldsymbol{\phi}_n^{k+1} \\ & \quad + (E^{(k)} (\mathbf{F}^{(k)} - \mu \boldsymbol{\phi}_n^k \mathbf{I}_2) + \mathbf{I}_2) \boldsymbol{\phi}_t^{k+1} - E^{(k)} (\mathbf{F}^{(k)} - \mu \boldsymbol{\phi}_n^k \mathbf{I}_2) \boldsymbol{\phi}_t^k = 0. \end{aligned} $$

Now, let us introduce the following operators:

$$ \mathbf{M}_{\alpha}^{*(k)} = E^{(k)} (\mathbf{F}^{(k)} - \mu \boldsymbol{\phi}_n^k \mathbf{I}_2), $$

$$ \mathbf{h}_{\alpha}^{(k)} = E^{(k)} \mathbf{F}^{(k)} (\boldsymbol{\phi}_t^k - \gamma_t \boldsymbol{w}_t^k) = \boldsymbol{\phi}_t^k. $$

$$ \mathbf{L}_{\alpha}^{*(k)} = -\gamma_t (\mathbf{I}_2 + \mathbf{M}_{\alpha}^{*(k)})^{-1} \mathbf{M}_{\alpha}^{*(k)}, $$

$$ \mathbf{r}_{\alpha}^{*(k)} = (\mathbf{I}_2 + \mathbf{M}_{\alpha}^{*(k)})^{-1} \mathbf{h}_{\alpha}^{(k)}, $$

$$ \mathbf{v}_{\alpha}^{(k)} = \mu (\mathbf{I}_2 + \mathbf{M}_{\alpha}^{*(k)})^{-1} E^{(k)} (\boldsymbol{\phi}_t^k - \gamma_t \boldsymbol{w}_t^k). $$

And, at last after computation, we have

$$ \mathbf{L}_{\alpha}^{*(k)} \boldsymbol{w}_t^{k+1} - \mathbf{v}_{\alpha}^{(k)} \boldsymbol{\phi}_n^{k+1} + \boldsymbol{\phi}_t^{k+1} = \mathbf{r}_{\alpha}^{*(k)} - \mathbf{v}_{\alpha}^{(k)} \boldsymbol{\phi}_n^k. $$

In this specific problem, and for a two-dimensional case, one can obtain a simplified equiv- alent version of the algorithm. Let $\mathcal{D}_{C_t^\phi}^{\text{slip}}$ be the generalised derivative of $C_t^\phi$ in the slip case

$$ \begin{aligned} \mathcal{D}_{C_t^\phi}^{\text{slip}} (\boldsymbol{w}_t^\alpha, \boldsymbol{\phi}_t) (\delta \boldsymbol{w}_t^\alpha, \delta \boldsymbol{\phi}_t) = & \boldsymbol{\phi}_t \frac{(\boldsymbol{\phi}_t - \gamma_t \boldsymbol{w}_t^\alpha)^T}{\|\boldsymbol{\phi}_t - \gamma_t \boldsymbol{w}_t^\alpha\|} (\delta \boldsymbol{\phi}_t - \delta \boldsymbol{w}_t^\alpha) \\ & - \mu \boldsymbol{\phi}_n (\delta \boldsymbol{\phi}_t - \gamma_t \boldsymbol{w}_t^\alpha) + \|\boldsymbol{\phi}_t - \gamma_t \boldsymbol{w}_t^\alpha\| \delta \boldsymbol{\phi}_t. \end{aligned} \qquad (3.42) $$
Denoting by **t** the unit slip vector, we have

$$ \phi_t = \mu \phi_n \mathbf{t}, \quad \frac{\phi_t - \gamma_t \mathbf{w}_t^\alpha}{\|\phi_t - \gamma_t \mathbf{w}_t^\alpha\|} = \mathbf{t}, \quad \delta\phi_t - \gamma_t \mathbf{w}_t^\alpha = \eta \mathbf{t}. \qquad (3.43) $$

Combining (3.42)–(3.43), we get

$$ \mathcal{D}_{C_t^\phi}^{\text{slip}}(\mathbf{w}_t^\alpha, \phi_t)(\delta\mathbf{w}_t^\alpha, \delta\phi_t) = \mu\phi_n\eta(\mathbf{t}\mathbf{t}^T - \mathbf{I}_2)t + \|\phi_t - \gamma_t\mathbf{w}_t^\alpha\|\delta\phi_t. \qquad (3.44) $$

Since $\mathbf{t}\mathbf{t}^T + \mathbf{n}\mathbf{n}^T = \mathbf{I}_2$ in the 2D case, we have $(\mathbf{t}\mathbf{t}^T - \mathbf{I}_2)\mathbf{t} = \mathbf{n}\mathbf{n}^T\mathbf{t} = \mathbf{0}$. Using (3.26), we obtain

$$ \|\boldsymbol{\phi}_t^k - \gamma_t \mathbf{w}_t^k\| (\boldsymbol{\phi}_t^{k+1} - \boldsymbol{\phi}_t^k) = - \|\boldsymbol{\phi}_t^k - \gamma_t \mathbf{w}_t^k\| \boldsymbol{\phi}_t^k + \mu \phi_n^k (\boldsymbol{\phi}_t^k - \gamma_t \mathbf{w}_t^k). \qquad (3.45) $$

Therefore, (3.45) becomes

$$ \boldsymbol{\phi}_t^{k+1} = \mu \phi_n^k \frac{\boldsymbol{\phi}_t^k - \gamma_t \mathbf{w}_t^k}{\|\boldsymbol{\phi}_t^k - \gamma_t \mathbf{w}_t^k\|} = \mu \phi_n^k \mathbf{t}. \qquad (3.46) $$

Finally, the slip 2D conditions of the semi-smooth Newton formalism are as follows:

$$ \mathbf{w}_n^{k+1} = 0, \qquad (3.47) $$

$$ \boldsymbol{\phi}_t^{k+1} = \mu \phi_n^k \frac{\boldsymbol{\phi}_t^k - \gamma_t \mathbf{w}_t^k}{\|\boldsymbol{\phi}_t^k - \gamma_t \mathbf{w}_t^k\|}. \qquad (3.48) $$

# 4 Primal-dual active set methods

This section is devoted to the numerical treatment of the contact conditions using a primal-dual active set method within the framework of dynamic contact problems. After defining the active and inactive subsets of all nodes that are currently in contact, we compute the contact conditions on each subset only in terms of contact reaction or impulses, using the local general equations of motion in the form of (3.3) and (3.4).

## 4.1 "Exact" primal-dual active set method (EPDAS)

Let us denote by $\mathcal{S}$ the set of potential contact and by $\alpha$ a potential contact node belonging to $\mathcal{S}$. The frictional contact conditions (2.15)–(2.17) are realised by applying an active set strategy which derives directly from the computation of the fixed point on the nonlinear complementary functions $C_n^\phi$ and $C_t^\phi$ based on the Newton semi-smooth scheme derived in part 3. The active and inactive sets are defined as follows:

$$ \mathcal{A}_n^{k+1} = \{\alpha \in \mathcal{S} : \phi_n^{\alpha,k} - \gamma_n \mathbf{w}_n^{\alpha,k} \geq 0\}, $$

$$ \mathcal{I}_n^{k+1} = \{\alpha \in \mathcal{S} : \phi_n^{\alpha,k} - \gamma_n \mathbf{w}_n^{\alpha,k} < 0\}, $$

$$ \mathcal{A}_t^{k+1} = \{\alpha \in \mathcal{S} : \|\boldsymbol{\phi}_t^{\alpha,k} - \gamma_t \mathbf{w}_t^{\alpha,k}\| - \mu \phi_n^{\alpha,k} \geq 0\}, $$

$$ \mathcal{I}_t^{k+1} = \{\alpha \in \mathcal{S} : \|\boldsymbol{\phi}_t^{\alpha,k} - \gamma_t \mathbf{w}_t^{\alpha,k}\| - \mu \phi_n^{\alpha,k} \leq 0\}. $$

The status of a given potential $\alpha$ at the nonlinear iteration $j$ depends on the set it belongs to. It can be either in the noncontact, slip or stick status.
Denoting by **t** the unit slip vector, we have

$$ \phi_t = \mu \phi_n \mathbf{t}, \quad \frac{\phi_t - \gamma_t \mathbf{w}_t^\alpha}{\|\phi_t - \gamma_t \mathbf{w}_t^\alpha\|} = \mathbf{t}, \quad \delta \phi_t - \gamma_t \mathbf{w}_t^\alpha = \eta \mathbf{t}. \qquad (3.43) $$

Combining (3.42)–(3.43), we get

$$ \mathcal{D}_{C_t^\phi}^{\text{slip}}(\mathbf{w}_t^\alpha, \phi_t)(\delta \mathbf{w}_t^\alpha, \delta \phi_t) = \mu \phi_n \eta (\mathbf{t}\mathbf{t}^T - \mathbf{I}_2)t + \|\phi_t - \gamma_t \mathbf{w}_t^\alpha\| \delta \phi_t. \qquad (3.44) $$

Since $\mathbf{t}\mathbf{t}^T + \mathbf{n}\mathbf{n}^T = \mathbf{I}_2$ in the 2D case, we have $(\mathbf{t}\mathbf{t}^T - \mathbf{I}_2)\mathbf{t} = \mathbf{n}\mathbf{n}^T\mathbf{t} = \mathbf{0}$. Using (3.26), we obtain

$$ \|\boldsymbol{\phi}_t^k - \gamma_t \mathbf{w}_t^k\| (\boldsymbol{\phi}_t^{k+1} - \boldsymbol{\phi}_t^k) = - \|\boldsymbol{\phi}_t^k - \gamma_t \mathbf{w}_t^k\| \boldsymbol{\phi}_t^k + \mu \phi_n^k (\boldsymbol{\phi}_t^k - \gamma_t \mathbf{w}_t^k). \qquad (3.45) $$

Therefore, (3.45) becomes

$$ \boldsymbol{\phi}_t^{k+1} = \mu \phi_n^k \frac{\boldsymbol{\phi}_t^k - \gamma_t \mathbf{w}_t^k}{\|\boldsymbol{\phi}_t^k - \gamma_t \mathbf{w}_t^k\|} = \mu \phi_n^k \mathbf{t}. \qquad (3.46) $$

Finally, the slip 2D conditions of the semi-smooth Newton formalism are as follows:

$$ \mathbf{w}_n^{k+1} = 0, \qquad (3.47) $$

$$ \boldsymbol{\phi}_t^{k+1} = \mu \phi_n^k \frac{\boldsymbol{\phi}_t^k - \gamma_t \mathbf{w}_t^k}{\|\boldsymbol{\phi}_t^k - \gamma_t \mathbf{w}_t^k\|}. \qquad (3.48) $$

# 4 Primal-dual active set methods

This section is devoted to the numerical treatment of the contact conditions using a primal-dual active set method within the framework of dynamic contact problems. After defining the active and inactive subsets of all nodes that are currently in contact, we compute the contact conditions on each subset only in terms of contact reaction or impulses, using the local general equations of motion in the form of (3.3) and (3.4).

## 4.1 "Exact" primal-dual active set method (EPDAS)

Let us denote by $\mathcal{S}$ the set of potential contact and by $\alpha$ a potential contact node belonging to $\mathcal{S}$. The frictional contact conditions (2.15)–(2.17) are realised by applying an active set strategy which derives directly from the computation of the fixed point on the nonlinear complementary functions $C_n^\phi$ and $C_t^\phi$ based on the Newton semi-smooth scheme derived in part 3. The active and inactive sets are defined as follows:

$$ \mathcal{A}_n^{k+1} = \{\alpha \in \mathcal{S} : \phi_n^{\alpha,k} - \gamma_n \mathbf{w}_n^{\alpha,k} \geq 0\}, $$

$$ \mathcal{I}_n^{k+1} = \{\alpha \in \mathcal{S} : \phi_n^{\alpha,k} - \gamma_n \mathbf{w}_n^{\alpha,k} < 0\}, $$

$$ \mathcal{A}_t^{k+1} = \{\alpha \in \mathcal{S} : \|\boldsymbol{\phi}_t^{\alpha,k} - \gamma_t \mathbf{w}_t^{\alpha,k}\| - \mu \phi_n^{\alpha,k} \geq 0\}, $$

$$ \mathcal{I}_t^{k+1} = \{\alpha \in \mathcal{S} : \|\boldsymbol{\phi}_t^{\alpha,k} - \gamma_t \mathbf{w}_t^{\alpha,k}\| - \mu \phi_n^{\alpha,k} \leq 0\}. $$

The status of a given potential $\alpha$ at the nonlinear iteration $j$ depends on the set it belongs to. It can be either in the noncontact, slip or stick status.
Remark For rigid problems, according to Sect. 2.2.3, the numerical computation of the local contact step inside the NLGS iteration loop of index j leads to the following primal-dual active set algorithm indexed by k. That is to say, we make only one iteration of the PDAS algorithm at each NLGS iteration, because the convergence of the NLGS method is slower than the convergence of the PDAS method. Therefore k = j in the following algorithm.

(i) Choose $(\mathbf{w}^{(0)}, \boldsymbol{\phi}^{(0)}), \gamma_n > 0, \gamma_t > 0$ and set $k = 0$.

(ii) Compute: $\tau_n^{\alpha,k+1} = \phi_n^{\alpha,k} - \gamma_n w_n^{\alpha,k}$ for each $\alpha \in S$.

(iii) Compute: $\tau_t^{\alpha,k+1} = \|\boldsymbol{\phi}_t^{\alpha,k} - \gamma_t \mathbf{w}_t^{\alpha,k}\| - \mu \phi_n^{\alpha,k}$ for each $\alpha \in S$.

(iv) Set the active and inactive sets:

$$A_n^{k+1} = \{\alpha \in S : \tau_n^{\alpha,k+1} \geq 0\},$$

$$\mathcal{I}_n^{k+1} = S \setminus A_n^{k+1},$$

$$A_t^{k+1} = \{\alpha \in S : \tau_t^{\alpha,k+1} \geq 0\},$$

$$\mathcal{I}_t^{k+1} = S \setminus A_t^{k+1}.$$

(v) Find $(\mathbf{w}^{\alpha,k+1}, \boldsymbol{\phi}^{\alpha,k+1})$ such that

$$\phi_n^{\alpha,k+1} = 0, \quad \boldsymbol{\phi}_t^{\alpha,k+1} = \mathbf{0} \quad \text{for all } \alpha \in \mathcal{I}_n^{k+1}, \qquad (4.1)$$

$$w_n^{\alpha,k+1} = 0 \quad \text{for all } \alpha \in A_n^{k+1}, \qquad (4.2)$$

$$\boldsymbol{\phi}_t^{\alpha,k+1} = \mu \phi_n^{\alpha,k} \frac{(\boldsymbol{\phi}_t^{\alpha,k} - \gamma_t \mathbf{w}_t^{\alpha,k})}{\|\boldsymbol{\phi}_t^{\alpha,k} - \gamma_t \mathbf{w}_t^{\alpha,k}\|} \quad \text{for all } \alpha \in A_t^{k+1} \cap A_n^{k+1}, \qquad (4.3)$$

$$\mathbf{w}_t^{\alpha,k+1} - \mathbf{w}_t^{\alpha,k} = -\mathbf{w}_t^{\alpha,k} \frac{\boldsymbol{\phi}_n^{\alpha,k+1}}{\boldsymbol{\phi}_n^{\alpha,k}}. \qquad (4.4)$$

(vi) If $\|(\mathbf{w}^{\alpha,k+1}, \boldsymbol{\phi}^{\alpha,k+1}) - (\mathbf{w}^{\alpha,j}, \boldsymbol{\phi}^{\alpha,k})\| \le \epsilon$, $A_n^{k+1} = A_n^k$ and $A_t^{k+1} = A_t^k$ stop, else go to (ii).

## 4.2 "Iterative" primal-dual active set method (IPDAS)

Since we have noticed that (4.4) can lead to numerical instabilities when $\phi_n^{\alpha,k}$ is small, we introduce a variant of the previous method. For that purpose, we have to introduce the couple $(\mathbf{w}^{i,j+1}, \boldsymbol{\phi}^{i,k+1})$, which represents the value of $(\mathbf{w}, \boldsymbol{\phi})$ at the $i$th fixed point iteration and at the $k$th active set iteration. Moreover, $(\mathbf{w}^{(i-1,j)}, \boldsymbol{\phi}^{(i-1,j)})$ represents the values of $(\mathbf{w}, \boldsymbol{\phi})$ obtained at the convergence of the active set method for the fixed point iteration $i-1$. Therefore, for each active set iteration, we fix $\phi_n^{\alpha,(i,k+1)}$ to the value $\phi_n^{\alpha,(i-1,j)}$ obtained at the previous fixed point iteration. Then, condition (4.4) leads to $\mathbf{w}_t^{\alpha,k} = \mathbf{0}$.

(i) Choose $(\mathbf{w}^{(0)}, \boldsymbol{\phi}^{(0)}), \gamma_n > 0, \gamma_t > 0$ and set $k = 0$.

(ii) Compute: $\tau_n^{\alpha,(i,k+1)} = \phi_n^{\alpha,(i,k)} - \gamma_n w_n^{\alpha,(i,k)}$ for each $\alpha \in S$.

(iii) Compute: $\tau_t^{\alpha,(i,k+1)} = \|\boldsymbol{\phi}_t^{\alpha,(i,k)} - \gamma_t \mathbf{w}_t^{\alpha,(i,k)}\| - \mu \phi_n^{\alpha,(i,k)}$ for each $\alpha \in S$.

(iv) Set the active and inactive sets:

$$A_n^{k+1} = \{\alpha \in S : \tau_n^{\alpha,(i,k+1)} \geq 0\},$$

$$\mathcal{I}_n^{k+1} = S \setminus A_n^{k+1},$$

$$A_t^{k+1} = \{\alpha \in S : \tau_t^{\alpha,(i,k+1)} \geq 0\},$$

$$\mathcal{I}_t^{k+1} = S \setminus A_t^{k+1}.$$
Remark For rigid problems, according to Sect. 2.2.3, the numerical computation of the local contact step inside the NLGS iteration loop of index j leads to the following primal-dual active set algorithm indexed by k. That is to say, we make only one iteration of the PDAS algorithm at each NLGS iteration, because the convergence of the NLGS method is slower than the convergence of the PDAS method. Therefore k = j in the following algorithm.

(i) Choose $(\mathbf{w}^{(0)}, \boldsymbol{\phi}^{(0)}), \gamma_n > 0, \gamma_t > 0$ and set $k = 0$.

(ii) Compute: $\tau_n^{\alpha,k+1} = \phi_n^{\alpha,k} - \gamma_n w_n^{\alpha,k}$ for each $\alpha \in S$.

(iii) Compute: $\tau_t^{\alpha,k+1} = \|\boldsymbol{\phi}_t^{\alpha,k} - \gamma_t \mathbf{w}_t^{\alpha,k}\| - \mu \phi_n^{\alpha,k}$ for each $\alpha \in S$.

(iv) Set the active and inactive sets:

$$A_n^{k+1} = \{\alpha \in S : \tau_n^{\alpha,k+1} \ge 0\},$$

$$\mathcal{I}_n^{k+1} = S \setminus A_n^{k+1},$$

$$A_t^{k+1} = \{\alpha \in S : \tau_t^{\alpha,k+1} \ge 0\},$$

$$\mathcal{I}_t^{k+1} = S \setminus A_t^{k+1}.$$

(v) Find $(\mathbf{w}^{\alpha,k+1}, \boldsymbol{\phi}^{\alpha,k+1})$ such that

$$\phi_n^{\alpha,k+1} = 0, \quad \boldsymbol{\phi}_t^{\alpha,k+1} = \mathbf{0} \quad \text{for all } \alpha \in \mathcal{I}_n^{k+1}, \qquad (4.1)$$

$$w_n^{\alpha,k+1} = 0 \quad \text{for all } \alpha \in A_n^{k+1}, \qquad (4.2)$$

$$\boldsymbol{\phi}_t^{\alpha,k+1} = \mu \phi_n^{\alpha,k} \frac{(\boldsymbol{\phi}_t^{\alpha,k} - \gamma_t \mathbf{w}_t^{\alpha,k})}{\|\boldsymbol{\phi}_t^{\alpha,k} - \gamma_t \mathbf{w}_t^{\alpha,k}\|} \quad \text{for all } \alpha \in A_t^{k+1} \cap A_n^{k+1}, \qquad (4.3)$$

$$\mathbf{w}_t^{\alpha,k+1} - \mathbf{w}_t^{\alpha,k} = -\mathbf{w}_t^{\alpha,k} \frac{\boldsymbol{\phi}_n^{\alpha,k+1}}{\boldsymbol{\phi}_n^{\alpha,k}}. \qquad (4.4)$$

(vi) If $\|(\mathbf{w}^{\alpha,k+1}, \boldsymbol{\phi}^{\alpha,k+1}) - (\mathbf{w}^{\alpha,j}, \boldsymbol{\phi}^{\alpha,k})\| \le \epsilon$, $A_n^{k+1} = A_n^k$ and $A_t^{k+1} = A_t^k$ stop, else go to (ii).

## 4.2 "Iterative" primal-dual active set method (IPDAS)

Since we have noticed that (4.4) can lead to numerical instabilities when $\phi_n^{\alpha,k}$ is small, we introduce a variant of the previous method. For that purpose, we have to introduce the couple $(\mathbf{w}^{i,j+1}, \boldsymbol{\phi}^{i,k+1})$, which represents the value of $(\mathbf{w}, \boldsymbol{\phi})$ at the $i$th fixed point iteration and at the $k$th active set iteration. Moreover, $(\mathbf{w}^{(i-1,j)}, \boldsymbol{\phi}^{(i-1,j)})$ represents the values of $(\mathbf{w}, \boldsymbol{\phi})$ obtained at the convergence of the active set method for the fixed point iteration $i-1$. Therefore, for each active set iteration, we fix $\phi_n^{\alpha,(i,k+1)}$ to the value $\phi_n^{\alpha,(i-1,j)}$ obtained at the previous fixed point iteration. Then, condition (4.4) leads to $\mathbf{w}_t^{\alpha,k} = \mathbf{0}$.

(i) Choose $(\mathbf{w}^{(0)}, \boldsymbol{\phi}^{(0)}), \gamma_n > 0, \gamma_t > 0$ and set $k = 0$.

(ii) Compute: $\tau_n^{\alpha,(i,k+1)} = \phi_n^{\alpha,(i,k)} - \gamma_n w_n^{\alpha,(i,k)}$ for each $\alpha \in S$.

(iii) Compute: $\tau_t^{\alpha,(i,k+1)} = \|\boldsymbol{\phi}_t^{\alpha,(i,k)} - \gamma_t \mathbf{w}_t^{\alpha,(i,k)}\| - \mu \phi_n^{\alpha,(i,k)}$ for each $\alpha \in S$.

(iv) Set the active and inactive sets:

$$A_n^{k+1} = \{\alpha \in S : \tau_n^{\alpha,(i,k+1)} \ge 0\},$$

$$\mathcal{I}_n^{k+1} = S \setminus A_n^{k+1},$$

$$A_t^{k+1} = \{\alpha \in S : \tau_t^{\alpha,(i,k+1)} \ge 0\},$$

$$\mathcal{I}_t^{k+1} = S \setminus A_t^{k+1}.$$
(v) Find $(\mathbf{w}^{\alpha(i,k+1)}, \boldsymbol{\phi}^{\alpha(i,k+1)})$ such that

$$p_n^{\alpha(i,k+1)} = 0, \quad \boldsymbol{\phi}_t^{\alpha(i,k+1)} = \mathbf{0} \quad \text{for all } \alpha \in \mathcal{I}_n^{k+1}, \qquad (4.5)$$

$$\mathbf{w}_n^{\alpha(i,k+1)} = 0 \quad \text{for all } \alpha \in \mathcal{A}_n^{k+1}, \qquad (4.6)$$

$$\boldsymbol{\phi}_t^{\alpha(i,k+1)} = \mu \frac{\boldsymbol{\phi}_t^{\alpha(i,k)} - \gamma_t \mathbf{w}_t^{\alpha(i,k)}}{\|\boldsymbol{\phi}_t^{\alpha(i,k)} - \gamma_t \mathbf{w}_t^{\alpha(i,k)}\|} \quad \text{for all } \alpha \in \mathcal{A}_t^{k+1} \cap \mathcal{A}_n^{k+1}, \qquad (4.7)$$

$$\mathbf{w}_t^{\alpha(i,k)} = \mathbf{0} \quad \text{for all } \alpha \in \mathcal{I}_t^{k+1} \cap \mathcal{A}_n^{k+1}. \qquad (4.8)$$

(vi) If $\| (\mathbf{w}^{\alpha(i,k+1)}, \boldsymbol{\phi}^{\alpha(i,k+1)}) - (\mathbf{w}^{\alpha(i,k)}, \boldsymbol{\phi}^{\alpha(i,k)}) \| \le \epsilon$, $\mathcal{A}_n^{k+1} = \mathcal{A}_n^k$ and $\mathcal{A}_t^{k+1} = \mathcal{A}_t^k$ stop, else go to (ii).

*Remark* For rigid problems, we do not use this fixed point iteration *i* because the number of iterations of NLGS is very high, it suffices to replace relation (4.4) by

$$\boldsymbol{\phi}_t^{\alpha(i,k+1)} + \frac{\mathbf{w}_t^{\alpha(i,k)}}{W_{tt}} = \boldsymbol{\phi}_t^{\alpha(i,k)} \quad \text{for all } \alpha \in \mathcal{I}_t^{k+1} \cap \mathcal{A}_n^{k+1}. \qquad (4.9)$$

*Remark* For deformable problems, at the beginning to step (v), we must solve the nonlinear systems (2.8).

*Remark* In the case of rigid problem, the stopping criterion (vi) can be enriched with physical considerations, such as the properties of the networks of forces in the granular medium.

# 5 Numerical experiments

A sweeping process based formulation for the frictional contact problem has been drawn in the previous section. Combining the semi-smooth Newton approach with the primal-dual active set method leads to a specific numerical method to tackle contact mechanics problems. For illustration purpose, some simulations are presented hereafter. This intends to outline the ability of the present numerical scheme to undergo frictional contact problems.

## 5.1 Bounce of a hyper-elastic ring against a rigid foundation

We first consider a representative application, a hyper-elastic ring bouncing on a perfectly rigid foundation, to assess the performances of the active set type in a large deformation framework. The ring is thrown toward a rigid foundation with an initial velocity at an angle of 45°, as depicted in Fig. 2. We provide below a description of the numerical settings:

$$\rho = 1000 \text{ kg/m}^3, \quad T = 5 \text{ s}, \quad k = \frac{1}{500},$$

$$\mathbf{u}_0 = (0, 0) \text{ m}, \quad \mathbf{u}_1 = (10, -10) \text{ m/s},$$

$$c_1 = 0.5 \text{ MPa}, \quad c_2 = 0.05 \text{ MPa}, \quad a = 0.5 \times 10^{-4} \text{ MPa},$$

$$c_v = 10, \quad c_\tau = 10, \quad r_{\text{lagrangian}} = 0.1,$$

$$g = 50 \text{ m}, \quad \mu = 0.2,$$
(v) Find $(\mathbf{w}^{\alpha(i,k+1)}, \boldsymbol{\phi}^{\alpha(i,k+1)})$ such that

$$p_n^{\alpha(i,k+1)} = 0, \quad \boldsymbol{\phi}_t^{\alpha(i,k+1)} = \mathbf{0} \quad \text{for all } \alpha \in \mathcal{I}_n^{k+1}, \qquad (4.5)$$

$$\mathbf{w}_n^{\alpha(i,k+1)} = 0 \quad \text{for all } \alpha \in \mathcal{A}_n^{k+1}, \qquad (4.6)$$

$$\boldsymbol{\phi}_t^{\alpha(i,k+1)} = \mu \frac{\boldsymbol{\phi}_t^{\alpha(i,k)} - \gamma_t \mathbf{w}_t^{\alpha(i,k)}}{\|\boldsymbol{\phi}_t^{\alpha(i,k)} - \gamma_t \mathbf{w}_t^{\alpha(i,k)}\|} \quad \text{for all } \alpha \in \mathcal{A}_t^{k+1} \cap \mathcal{A}_n^{k+1}, \qquad (4.7)$$

$$\mathbf{w}_t^{\alpha(i,k)} = \mathbf{0} \quad \text{for all } \alpha \in \mathcal{I}_t^{k+1} \cap \mathcal{A}_n^{k+1}. \qquad (4.8)$$

(vi) If $\|(\mathbf{w}^{\alpha(i,k+1)}, \boldsymbol{\phi}^{\alpha(i,k+1)}) - (\mathbf{w}^{\alpha(i,k)}, \boldsymbol{\phi}^{\alpha(i,k)})\| \le \epsilon$, $\mathcal{A}_n^{k+1} = \mathcal{A}_n^k$ and $\mathcal{A}_t^{k+1} = \mathcal{A}_t^k$ stop, else go to (ii).

*Remark* For rigid problems, we do not use this fixed point iteration *i* because the number of iterations of NLGS is very high, it suffices to replace relation (4.4) by

$$\boldsymbol{\phi}_t^{\alpha(i,k+1)} + \frac{\mathbf{w}_t^{\alpha(i,k)}}{W_{tt}} = \boldsymbol{\phi}_t^{\alpha(i,k)} \quad \text{for all } \alpha \in \mathcal{I}_t^{k+1} \cap \mathcal{A}_n^{k+1}. \qquad (4.9)$$

*Remark* For deformable problems, at the beginning to step (v), we must solve the nonlinear systems (2.8).

*Remark* In the case of rigid problem, the stopping criterion (vi) can be enriched with physical considerations, such as the properties of the networks of forces in the granular medium.

# 5 Numerical experiments

A sweeping process based formulation for the frictional contact problem has been drawn in the previous section. Combining the semi-smooth Newton approach with the primal-dual active set method leads to a specific numerical method to tackle contact mechanics problems. For illustration purpose, some simulations are presented hereafter. This intends to outline the ability of the present numerical scheme to undergo frictional contact problems.

## 5.1 Bounce of a hyper-elastic ring against a rigid foundation

We first consider a representative application, a hyper-elastic ring bouncing on a perfectly rigid foundation, to assess the performances of the active set type in a large deformation framework. The ring is thrown toward a rigid foundation with an initial velocity at an angle of 45°, as depicted in Fig. 2. We provide below a description of the numerical settings:

$$\rho = 1000 \text{ kg/m}^3, \quad T = 5 \text{ s}, \quad k = \frac{1}{500},$$

$$\mathbf{u}_0 = (0, 0) \text{ m}, \quad \mathbf{u}_1 = (10, -10) \text{ m/s},$$

$$c_1 = 0.5 \text{ MPa}, \quad c_2 = 0.05 \text{ MPa}, \quad a = 0.5 \times 10^{-4} \text{ MPa},$$

$$c_v = 10, \quad c_\tau = 10, \quad r_{\text{lagrangian}} = 0.1,$$

$$g = 50 \text{ m}, \quad \mu = 0.2,$$
<figure><img src="image_3.png" /><figcaption>Figure 2 Sequence of the deformed hyper-elastic ring before, during and after impact</figcaption></figure>

stopping criterion: $\epsilon = 10^{-6}$.

Concerning the discretization size, it is described by a parameter nbc, which denotes the number of contact nodes on $\Gamma$. The number of nodes on the outer contour of the ring is thus taken in order to obtain sufficient regular elements. The boundary $\Gamma_3$ is divided into 128 equal parts, 1664 hyper-elastic finite elements were used for a total number of degrees of freedom equal to 2048.

The compressible material response illustrated here is described by a variant of the Ogden constitutive law (cf. [40]) defined by this energy density function

$$W(\mathbf{C}) = c_1(I_1 - 3) + c_2(I_2 - 3) + a(I_3 - 1) - (c_1 + 2c_2 + a) \ln I_3$$

with $I_1, I_2$ and $I_3$, the three invariants of the tensor $\mathbf{C}$, with $\mathbf{C} = \mathbf{F}^T\mathbf{F}$.

Details on the problem's physical settings are given as follows:

$$\Omega = \{(x_1, x_2) \in \mathbb{R}^2 : 81 < (x_1 - 100)^2 + (x_2 - 100)^2 < 100\},$$

$$\Gamma_1 = \emptyset, \quad \Gamma_2 = \emptyset,$$

$$\Gamma_3 = \{(x_1, x_2) \in \mathbb{R}^2 : (x_1 - 100)^2 + (x_2 - 100)^2 = 100\}.$$

$\Omega$ represents the section of a three-dimensional deformable body under the plane constraint assumption. The foundation is given by $\{(x_1, x_2) \in \mathbb{R}^2 : x_2 \le -50\}$. For the discretization, we use 1664 elastic nodes and 128 contact nodes.

As previously announced, we will compare the results obtained by the active set method and the augmented quasi-Lagrangian method with regard to the convergence and the performances of the methods (see Tables 1 and 2).

The results obtained above show that the performances of the active set method are slightly better than those of the augmented Lagrangian one. According to Fig. 3, the higher number of nodes in contact we have, the more the gap in CPU time between the methods widens, which confirms the robustness of the active set type methods on complex and relevant configurations. Moreover, even if the convergence of the active set method is ensured only for small values of the numerical parameters $c_v$ and $c_\tau$ (see [26, 27]), we do
<figure><img src="image_4.png" /><figcaption>Figure 2 Sequence of the deformed hyper-elastic ring before, during and after impact</figcaption></figure>

stopping criterion: $\epsilon = 10^{-6}$.

Concerning the discretization size, it is described by a parameter nbc, which denotes the number of contact nodes on $\Gamma$. The number of nodes on the outer contour of the ring is thus taken in order to obtain sufficient regular elements. The boundary $\Gamma_3$ is divided into 128 equal parts, 1664 hyper-elastic finite elements were used for a total number of degrees of freedom equal to 2048.

The compressible material response illustrated here is described by a variant of the Ogden constitutive law (cf. [40]) defined by this energy density function

$$W(\mathbf{C}) = c_1(I_1 - 3) + c_2(I_2 - 3) + a(I_3 - 1) - (c_1 + 2c_2 + a) \ln I_3$$

with $I_1, I_2$ and $I_3$, the three invariants of the tensor $\mathbf{C}$, with $\mathbf{C} = \mathbf{F}^T\mathbf{F}$.

Details on the problem's physical settings are given as follows:

$$\Omega = \{(x_1, x_2) \in \mathbb{R}^2 : 81 < (x_1 - 100)^2 + (x_2 - 100)^2 < 100\},$$

$$\Gamma_1 = \emptyset, \quad \Gamma_2 = \emptyset,$$

$$\Gamma_3 = \{(x_1, x_2) \in \mathbb{R}^2 : (x_1 - 100)^2 + (x_2 - 100)^2 = 100\}.$$

$\Omega$ represents the section of a three-dimensional deformable body under the plane constraint assumption. The foundation is given by $\{(x_1, x_2) \in \mathbb{R}^2 : x_2 \le -50\}$. For the discretization, we use 1664 elastic nodes and 128 contact nodes.

As previously announced, we will compare the results obtained by the active set method and the augmented quasi-Lagrangian method with regard to the convergence and the performances of the methods (see Tables 1 and 2).

The results obtained above show that the performances of the active set method are slightly better than those of the augmented Lagrangian one. According to Fig. 3, the higher number of nodes in contact we have, the more the gap in CPU time between the methods widens, which confirms the robustness of the active set type methods on complex and relevant configurations. Moreover, even if the convergence of the active set method is ensured only for small values of the numerical parameters $c_v$ and $c_\tau$ (see [26, 27]), we do
<table><caption>**Table 1** Results of the Active set method for the unilateral contact and Coulomb friction laws in comparison with the number of degrees of freedom (dof), the number of contact nodes (nbc), the total numbers of Newton iterations (Ntit), the average fixed point iterations (afpit) and the total CPU time (CPU) in seconds</caption><tr><td>nbc</td><td>32</td><td>64</td><td>128</td><td>256</td><td>512</td></tr><tr><td>dof</td><td>192</td><td>384</td><td>1792</td><td>4608</td><td>15,360</td></tr><tr><td>Ntit</td><td>1876</td><td>1904</td><td>1944</td><td>1961</td><td>2022</td></tr><tr><td>afpit</td><td>4</td><td>4</td><td>4</td><td>4</td><td>4</td></tr><tr><td>CPU</td><td>6.32</td><td>12.39</td><td>99.22</td><td>304.62</td><td>1841.23</td></tr></table>

<table><caption>**Table 2** Results of the Augmented Lagrangian method for the unilateral contact and Coulomb friction laws in comparison with the number of degrees of freedom (dof), the number of contact nodes (nbc), the total numbers of Newton iterations (Ntit) and the total CPU time (CPU) in seconds</caption><tr><td>nbc</td><td>32</td><td>64</td><td>128</td><td>256</td><td>512</td></tr><tr><td>dof</td><td>256</td><td>512</td><td>2048</td><td>5120</td><td>16,384</td></tr><tr><td>Ntit</td><td>1462</td><td>1568</td><td>1705</td><td>1782</td><td>1866</td></tr><tr><td>CPU</td><td>8.59</td><td>18.03</td><td>139.69</td><td>428.80</td><td>2693.86</td></tr></table>

<figure><img src="image_4.png" /><figcaption>**Figure 3** Active set and Augmented Lagrangian CPU time with respect to nbc</figcaption></figure>

not observe any influence on the results, both in terms of convergence and the quality of the solution.

The robustness of PDAS methods lies in the fact that these methods do not use the Lagrange multipliers, unlike the classical augmented Lagrangian method. Avoiding the Lagrange multipliers is a relevant and strong argument in favour of the active set type methods.

## 5.2 Granular materials in a 2D rotating drum

The second numerical example concerns the rotating drum in which we evaluate the behaviour of a 2D granular flow. The choice of this test case lies in the fact that it is representative of multi-rigid body contact problems due to the large number of simultaneous contacts to deal with and the large amount of CPU time needed to compute the solution. Besides, this test case has been studied in several works using both experiments and simulations, motivating the present choice. Specifically, a collection of 400 particles is consid-
<table><caption>**Table 1** Results of the Active set method for the unilateral contact and Coulomb friction laws in comparison with the number of degrees of freedom (dof), the number of contact nodes (nbc), the total numbers of Newton iterations (Ntit), the average fixed point iterations (afpit) and the total CPU time (CPU) in seconds</caption><tr><td>nbc</td><td>32</td><td>64</td><td>128</td><td>256</td><td>512</td></tr><tr><td>dof</td><td>192</td><td>384</td><td>1792</td><td>4608</td><td>15,360</td></tr><tr><td>Ntit</td><td>1876</td><td>1904</td><td>1944</td><td>1961</td><td>2022</td></tr><tr><td>afpit</td><td>4</td><td>4</td><td>4</td><td>4</td><td>4</td></tr><tr><td>CPU</td><td>6.32</td><td>12.39</td><td>99.22</td><td>304.62</td><td>1841.23</td></tr></table>

<table><caption>**Table 2** Results of the Augmented Lagrangian method for the unilateral contact and Coulomb friction laws in comparison with the number of degrees of freedom (dof), the number of contact nodes (nbc), the total numbers of Newton iterations (Ntit) and the total CPU time (CPU) in seconds</caption><tr><td>nbc</td><td>32</td><td>64</td><td>128</td><td>256</td><td>512</td></tr><tr><td>dof</td><td>256</td><td>512</td><td>2048</td><td>5120</td><td>16,384</td></tr><tr><td>Ntit</td><td>1462</td><td>1568</td><td>1705</td><td>1782</td><td>1866</td></tr><tr><td>CPU</td><td>8.59</td><td>18.03</td><td>139.69</td><td>428.80</td><td>2693.86</td></tr></table>

<figure><img src="image_5.png" /><figcaption>**Figure 3** Active set and Augmented Lagrangian CPU time with respect to nbc</figcaption></figure>

not observe any influence on the results, both in terms of convergence and the quality of the solution.

The robustness of PDAS methods lies in the fact that these methods do not use the Lagrange multipliers, unlike the classical augmented Lagrangian method. Avoiding the Lagrange multipliers is a relevant and strong argument in favour of the active set type methods.

## 5.2 Granular materials in a 2D rotating drum

The second numerical example concerns the rotating drum in which we evaluate the behaviour of a 2D granular flow. The choice of this test case lies in the fact that it is representative of multi-rigid body contact problems due to the large number of simultaneous contacts to deal with and the large amount of CPU time needed to compute the solution. Besides, this test case has been studied in several works using both experiments and simulations, motivating the present choice. Specifically, a collection of 400 particles is consid-
<table><caption>**Table 3** Physical parameters of the flow in a rotating drum</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Particles radius (mm)</td><td>r = [3.5–4]</td></tr><tr><td>Density (kg · m<sup>-3</sup>)</td><td>ρ = 2700</td></tr><tr><td>Drum diameter (mm)</td><td>D<sub>drum</sub> = 450</td></tr><tr><td>rotation rate (rpm)</td><td>Ω<sub>drum</sub> = 6</td></tr><tr><td>Gravitational acceleration (m · s<sup>-2</sup>)</td><td>g = -9.80665</td></tr><tr><td>Mass particle (kg)</td><td>m<sub>p</sub> = ρπr<sub>p</sub><sup>2</sup></td></tr></tbody></table>

<table><caption>**Table 4** Numerical parameters for simulation in a rotating drum for each contact approach</caption><thead><tr><td>Numerical parameter DEM</td><td>Numerical parameter NSCD</td><td>Other</td></tr></thead><tbody><tr><td>e<sub>n</sub> = 0.59</td><td>e<sub>n</sub> = 0.46</td><td>T = 5 s</td></tr><tr><td>k<sub>n</sub> = 10<sup>3</sup> N · m<sup>-1</sup></td><td>e<sub>t</sub> = 1</td><td>–</td></tr><tr><td>k<sub>t</sub> = <sup>2</sup>&frasl;<sub>7</sub>k<sub>n</sub> N · m<sup>-1</sup></td><td>γ<sub>n</sub> = 10<sup>-8</sup></td><td>μ = 0.2</td></tr><tr><td>η<sub>ratio</sub> = 0.5</td><td>γ<sub>t</sub> = 10<sup>-8</sup></td><td>–</td></tr><tr><td>–</td><td>ε<sub>NLGS</sub> = 10<sup>-8</sup></td><td>–</td></tr><tr><td>dt ≈ 3.810<sup>-6</sup> s</td><td>dt = 1.10<sup>-5</sup> s</td><td>–</td></tr></tbody></table>

ered inside a 450 mm diameter drum with a rotation speed equal to 6 rpm. The uniform random distribution chosen between 7 and 8 mm avoids the blockages during the drum rotation process and the particles are organised in two layers: one in a sticked contact with the inner wall of the drum and the other as a bed.

To assess the correctness of the simulation carried out using NSCD-PDAS, we compare the numerical results with the ones obtained with the well-known discrete element method (DEM) initiated and developed by Cundall and Strack [41], which consists of an explicit treatment of particle-particle collisions using soft-sphere approach (based on the “Spring–Dashpot” model) as detailed in [42], and opposed to the nonsmooth approach. To ease the comparisons with the DEM-CUNDALL approach and to verify the consistency of the numerical results, the NSCD-PDAS was implemented with the MFiX-EXA framework (see [42] for more details). In this way, one can take advantage from the parallel computing to speed up the computations. Thus, all the simulations were performed using 96 processors. A summary of the physical model and numerical parameters of the bidimensional rotating drum test case is reported in Tables 3 and 4. Let us notice that in this example the active set numerical parameters $c_v$ and $c_t$ can be chosen in a large interval, ranging from $10^{-8}$ to $10^2$, in order to obtain a convergence of the fixed point in the active set method.

We first evaluate the surface flow obtained by DEM-CUNDALL and NSCD-PDAS in the steady state. For this purpose, we first make a simulation over a time integration of 5 s. Therefore, one can observe according to Fig. 4 that the steady surface flow reached after 5 s of simulations is quite similar for both approaches in parallel, with $dt_{NSCD} > dt_{DEM}$.

With regard to the interpenetrations generated during a normal persistent contact between two particles, we measure them for both approaches, as depicted in Fig. 5.

The nonsmooth formulation of the frictional contact conditions in the NSCD framework ensures a non-interpenetrability between the bodies (average penetration of the order of $10^{-6}$ m on some contacts)(see Fig. 5(b)), while the DEM-CUNDALL approach is characterised by a more important interpenetrability (average penetration of the order of $10^{-6}$ m on a good part of the contacts) for a normal spring constant $k_H = 10^3$ Nm$^{-1}$, as
<table><caption>**Table 3** Physical parameters of the flow in a rotating drum</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Particles radius (mm)</td><td>r = [3.5–4]</td></tr><tr><td>Density (kg · m<sup>-3</sup>)</td><td>ρ = 2700</td></tr><tr><td>Drum diameter (mm)</td><td>D<sub>drum</sub> = 450</td></tr><tr><td>rotation rate (rpm)</td><td>Ω<sub>drum</sub> = 6</td></tr><tr><td>Gravitational acceleration (m · s<sup>-2</sup>)</td><td>g = -9.80665</td></tr><tr><td>Mass particle (kg)</td><td>m<sub>p</sub> = ρπr<sub>p</sub><sup>2</sup></td></tr></tbody></table>

<table><caption>**Table 4** Numerical parameters for simulation in a rotating drum for each contact approach</caption><thead><tr><td>Numerical parameter DEM</td><td>Numerical parameter NSCD</td><td>Other</td></tr></thead><tbody><tr><td>e<sub>n</sub> = 0.59</td><td>e<sub>n</sub> = 0.46</td><td>T = 5 s</td></tr><tr><td>k<sub>n</sub> = 10<sup>3</sup> N · m<sup>-1</sup></td><td>e<sub>t</sub> = 1</td><td>–</td></tr><tr><td>k<sub>t</sub> = <sup>2</sup>&frasl;<sub>7</sub>k<sub>n</sub> N · m<sup>-1</sup></td><td>γ<sub>n</sub> = 10<sup>-8</sup></td><td>μ = 0.2</td></tr><tr><td>η<sub>ratio</sub> = 0.5</td><td>γ<sub>t</sub> = 10<sup>-8</sup></td><td>–</td></tr><tr><td>–</td><td>ε<sub>NLGS</sub> = 10<sup>-8</sup></td><td>–</td></tr><tr><td>dt ≈ 3.8.10<sup>-6</sup> s</td><td>dt = 1.10<sup>-5</sup> s</td><td>–</td></tr></tbody></table>

ered inside a 450 mm diameter drum with a rotation speed equal to 6 rpm. The uniform random distribution chosen between 7 and 8 mm avoids the blockages during the drum rotation process and the particles are organised in two layers: one in a sticked contact with the inner wall of the drum and the other as a bed.

To assess the correctness of the simulation carried out using NSCD-PDAS, we compare the numerical results with the ones obtained with the well-known discrete element method (DEM) initiated and developed by Cundall and Strack [41], which consists of an explicit treatment of particle-particle collisions using soft-sphere approach (based on the “Spring–Dashpot” model) as detailed in [42], and opposed to the nonsmooth approach. To ease the comparisons with the DEM-CUNDALL approach and to verify the consistency of the numerical results, the NSCD-PDAS was implemented with the MFiX-EXA framework (see [42] for more details). In this way, one can take advantage from the parallel computing to speed up the computations. Thus, all the simulations were performed using 96 processors. A summary of the physical model and numerical parameters of the bidimensional rotating drum test case is reported in Tables 3 and 4. Let us notice that in this example the active set numerical parameters $c_v$ and $c_τ$ can be chosen in a large interval, ranging from $10^{-8}$ to $10^2$, in order to obtain a convergence of the fixed point in the active set method.

We first evaluate the surface flow obtained by DEM-CUNDALL and NSCD-PDAS in the steady state. For this purpose, we first make a simulation over a time integration of 5 s. Therefore, one can observe according to Fig. 4 that the steady surface flow reached after 5 s of simulations is quite similar for both approaches in parallel, with $dt_{NSCD} > dt_{DEM}$.

With regard to the interpenetrations generated during a normal persistent contact between two particles, we measure them for both approaches, as depicted in Fig. 5.

The nonsmooth formulation of the frictional contact conditions in the NSCD framework ensures a non-interpenetrability between the bodies (average penetration of the order of $10^{-6}$ m on some contacts)(see Fig. 5(b)), while the DEM-CUNDALL approach is characterised by a more important interpenetrability (average penetration of the order of $10^{-6}$ m on a good part of the contacts) for a normal spring constant $k_H = 10^3$ Nm$^{-1}$, as
<figure><img src="image_5.png" /><figcaption>**Figure 4** Perspective views of granular flow profiles in a rotating drum obtained by DEM-CUNDALL [(a) particle diameter, (b) processor rank] and NSCD-PDAS [(c) particle diameter, (d) processor rank] when the steady surface flow is reached</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5** Visualization of interpenetrations measured (PN) during the simulations of granular flows in a rotating drum obtained by (a) DEM-CUNDALL and (b) NSCD-PDAS when the surface flow is reached</figcaption></figure>

shown in Fig. 5(a). This is an expected result since the CUNDALL model allows interpenetratión.

To highlight the conservation property of the NSCD method, we compare the evolution of the kinetic energy over time between the two approaches. Since the DEM-CUNDALL approach is dependant on the normal spring constant $k_n$, a parametric analysis is performed by taking two values ($k_n = 10^3 \text{ N} \cdot \text{m}^{-1}$ and $k_n = 10^5 \text{ N} \cdot \text{m}^{-1}$). Figure 6 shows that the system's kinetic energy is conserved with the NSCD-PDAS approach (red curve), whereas the influence of the spring constant is noticeable. The DEM-CUNDALL approach shows a rise of the kinetic energy with respect to increasing normal spring constant $k_n$. One can expect that for higher values of $k_n$ the kinetic energy level will converge to the NSCD-PDAS, but using infinitely small time steps to ensure the numerical stability.

# 6 Conclusion

A unified modelling and numerical treatment of contact problems for both deformable and rigid bodies has been developed in this work. In this context, the complementary
<figure><img src="image_6.png" /><figcaption>**Figure 4** Perspective views of granular flow profiles in a rotating drum obtained by DEM-CUNDALL [(a) particle diameter, (b) processor rank] and NSCD-PDAS [(c) particle diameter, (d) processor rank] when the steady surface flow is reached</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5** Visualization of interpenetrations measured (PN) during the simulations of granular flows in a rotating drum obtained by (a) DEM-CUNDALL and (b) NSCD-PDAS when the surface flow is reached</figcaption></figure>

shown in Fig. 5(a). This is an expected result since the CUNDALL model allows interpenetratión.

To highlight the conservation property of the NSCD method, we compare the evolution of the kinetic energy over time between the two approaches. Since the DEM-CUNDALL approach is dependant on the normal spring constant $k_n$, a parametric analysis is performed by taking two values ($k_n = 10^3 \text{ N} \cdot \text{m}^{-1}$ and $k_n = 10^5 \text{ N} \cdot \text{m}^{-1}$). Figure 6 shows that the system's kinetic energy is conserved with the NSCD-PDAS approach (red curve), whereas the influence of the spring constant is noticeable. The DEM-CUNDALL approach shows a rise of the kinetic energy with respect to increasing normal spring constant $k_n$. One can expect that for higher values of $k_n$ the kinetic energy level will converge to the NSCD-PDAS, but using infinitely small time steps to ensure the numerical stability.

# 6 Conclusion

A unified modelling and numerical treatment of contact problems for both deformable and rigid bodies has been developed in this work. In this context, the complementary
<figure><img src="image_7.png" /><figcaption>**Figure 6** Evolution of the kinetic energy in DEM-CUNDALL (for different normal stiffness constants *k*<sub>n</sub>) and NSCD-PDAS simulations of granular flows in a rotating drum over time</figcaption></figure>

functions for both contact and friction formulated using a semi-smooth Newton approach lead to a fixed point problem. Therefore, the numerical treatment by the PDAS method leads for these two problems to an algorithm easy to implement. Numerical experiments of both hyper-elastic ring bouncing and granular materials in a rotating drum reported above have been performed and highlight the relevance of the proposed method and its efficiency in terms of CPU time, in comparison with the well-established DEM approach. Thus, the present method provides a unified modelling framework for the contact mechanics of contact discussed above and will be used for further developments.

**Acknowledgements**

Not applicable.

**Funding**

The project has received funding from the European Commission Horizon 2020 Research and Innovation Programme under the Marie Skłodowska-Curie grant agreement No. 823731 - CONMECH. This work was also realized with the support of MESO@LR, a Center of Competence in High-Performance Computing from the Languedoc-Roussillon region.

**Availability of data and materials**

Not applicable.

**Declarations**

**Competing interests**

The authors declare that they have no competing interests.

**Author contribution**

All authors contributed equally in writing this article. All authors read and approved the final manuscript.

**Author details**

¹Laboratoire de Mathématiques et Physique, Université de Perpignan Via Domitia, 52 Avenue Paul Alduy, 66860 Perpignan, France. ²Institut Montpelliérain Alexander Grothendieck, Université de Nîmes, Site des Carmes, Place Gabriel Péri, 30000 Nîmes, France.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
<figure><img src="image_8.png" /><figcaption>**Figure 6** Evolution of the kinetic energy in DEM-CUNDALL (for different normal stiffness constants *k*<sub>n</sub>) and NSCD-PDAS simulations of granular flows in a rotating drum over time</figcaption></figure>

functions for both contact and friction formulated using a semi-smooth Newton approach lead to a fixed point problem. Therefore, the numerical treatment by the PDAS method leads for these two problems to an algorithm easy to implement. Numerical experiments of both hyper-elastic ring bouncing and granular materials in a rotating drum reported above have been performed and highlight the relevance of the proposed method and its efficiency in terms of CPU time, in comparison with the well-established DEM approach. Thus, the present method provides a unified modelling framework for the contact mechanics of contact discussed above and will be used for further developments.

**Acknowledgements**

Not applicable.

**Funding**

The project has received funding from the European Commission Horizon 2020 Research and Innovation Programme under the Marie Skłodowska-Curie grant agreement No. 823731 - CONMECH. This work was also realized with the support of MESO@LR, a Center of Competence in High-Performance Computing from the Languedoc-Roussillon region.

**Availability of data and materials**

Not applicable.

**Declarations**

**Competing interests**

The authors declare that they have no competing interests.

**Author contribution**

All authors contributed equally in writing this article. All authors read and approved the final manuscript.

**Author details**

¹Laboratoire de Mathématiques et Physique, Université de Perpignan Via Domitia, 52 Avenue Paul Alduy, 66860 Perpignan, France. ²Institut Montpelliérain Alexander Grothendieck, Université de Nîmes, Site des Carmes, Place Gabriel Péri, 30000 Nîmes, France.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
References

1. Laursen, T.A.: Computational Contact and Impact Mechanics: Fundamentals of Modeling Interfacial Phenomena in Nonlinear Finite Element Analysis. Springer, Berlin (2013)

2. Wriggers, P., Zavarise, G.: Computational Contact Mechanics. Encyclopedia of Computational Mechanics (2004)

3. Brogliato, B.: Impacts in Mechanical Systems: Analysis and Modelling. Springer, Berlin (2000)

4. Acary, V., Brogliato, B.: Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics. Springer, Berlin (2008)

5. Moreau, J.J.: Application of convex analysis to some problems of dry friction. In: Trends in Applications of Pure Mathematics to Mechanics, pp. 263–280. Pitman, London (1977)

6. Moreau, J.J.: Unilateral contact and dry friction in finite freedom dynamics. In: Nonsmooth Mechanics and Applications, pp. 1–82. Springer, Berlin (1988)

7. Monteiro Marques, M.D.P.: Differential Inclusions in Nonsmooth Mechanical Problem: Shocks and Dry Friction. Springer, Basel (1993)

8. Jean, M.: The non-smooth contact dynamics method. Comput. Methods Appl. Mech. Eng. 177(3–4), 235–257 (1999)

9. Moreau, J.J.: Numerical aspects of the sweeping process. Comput. Methods Appl. Mech. Eng. 177(3–4), 329–349 (1999)

10. Jean, M., Moreau, J.J.: Unilaterality and dry friction in the dynamics of rigid body collections. In: 1st Contact Mechanics International Symposium, pp. 31–48 (1992)

11. Dubois, F., Acary, V., Jean, M.: The contact dynamics method: a nonsmooth story. C. R., Méc. 346(3), 247–262 (2018)

12. Moreau, J.J.: Some numerical methods in multibody dynamics: application to granular materials. Eur. J. Mech. A, Solids 13(4–suppl), 93–114 (1994)

13. de Saxcé, G., Feng, Z.-Q.: New inequality and functional for contact with friction: the implicit standard material approach. J. Struct. Mech. 19(3), 301–325 (1991)

14. Fortin, J., Millet, O., de Saxcé, G.: Numerical simulation of granular materials by an improved discrete element method. Int. J. Numer. Methods Eng. 62(5), 639–663 (2005)

15. Feng, Z.-Q., Joli, P., Cros, J.-M., Magnain, B.: The bi-potential method applied to the modeling of dynamic problems with friction. Comput. Mech. 36(5), 375–383 (2005)

16. Joli, P., Feng, Z.-Q.: Uzawa and Newton algorithms to solve frictional contact problems within the bi-potential framework. Int. J. Numer. Methods Eng. 73(3), 317–330 (2008)

17. Dumont, S.: On enhanced descent algorithms for solving frictional multicontact problems: application to the discrete element method. Int. J. Numer. Methods Eng. 93(11), 1170–1190 (2013)

18. Kikuchi, N., Oden, J.T.: Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods. SIAM Studies in Applied Mathematics, vol. 8. Philadelphia (1988)

19. Oden, J.T., Kim, S.J.: Interior penalty methods for finite element approximations of the Signorini problem in elastostatics. Comput. Math. Appl. 8(1), 35–56 (1986)

20. Alart, P., Curnier, A.: A mixed formulation for frictional contact problems prone to Newton like solution methods. Comput. Methods Appl. Mech. Eng. 92(3), 353–375 (1991)

21. Raous, M., Chabrand, P., Lebon, F.: Numerical methods for solving unilateral contact problem with friction. J. Theor. Appl. Mech. 7, 111–128 (1988)

22. Joli, P., Feng, Z.-Q.: Uzawa and Newton algorithms to solve frictional contact problems within the bi-potential framework. Int. J. Numer. Methods Eng. 73(3), 317–330 (2008)

23. Houly, F.: An adaptation of Nitsche's method to the Tresca friction problem. J. Math. Anal. Appl. 411(1), 329–339 (2014)

24. Houly, F., Hild, P., Renard, Y.: A Nitsche finite element method for dynamic contact: 1. Space semi-discretization and time-marching schemes. ESAIM: Math. Model. Numer. Anal. 49(2), 481–502 (2015)

25. Houly, F., Fabre, M., Hild, P., Mlika, R., Pousin, J., Renard, Y.: An overview of recent results on Nitsche's method for contact problems. In: Geometrically Unfitted Finite Element Methods and Applications, pp. 93–141 (2017)

26. Hintermüller, M., Kovtunenko, V.A., Kunisch, K.: Semismooth Newton methods for a class of unilaterally constrained variational problems. In: Universität Graz/Technische Universität Graz. SFB F003-Optimierung und Kontrolle (2003)

27. Hintermüller, M., Ito, K., Kunisch, K.: The primal-dual active set strategy as a semismooth Newton method. SIAM J. Optim. 13(3), 865–888 (2002)

28. Hintermüller, M., Kovtunenko, V.A., Kunisch, K.: Obstacle problems with cohesion: a hemivariational inequality approach and its efficient numerical solution. SIAM J. Optim. 21(2), 491–516 (2011)

29. Hüeber, S., Stadler, G., Wohlmuth, B.I.: A primal-dual active set algorithm for three-dimensional contact problems with Coulomb friction. SIAM J. Sci. Comput. 30(2), 572–596 (2008)

30. Hüeber, S., Wohlmuth, B.I.: A primal-dual active set strategy for non-linear multibody contact problems. Comput. Methods Appl. Mech. Eng. 194(27–29), 3147–3166 (2005)

31. Koziara, T., Bićanić, N.: Semismooth Newton method for frictional contact between pseudo-rigid bodies. Comput. Methods Appl. Mech. Eng. 197(33–40), 2763–2777 (2008)

32. Sharaf, I.M.: An active set algorithm for a class of linear complementarity problems arising from rigid body dynamics. Pak. J. Stat. Oper. Res. 339–352 (2016)

33. Barboteu, M., Dumont, S.: A primal-dual active set method for solving multi-rigid-body dynamic contact problems. Math. Mech. Solids 23(3), 489–503 (2018)

34. Abide, S., Barboteu, M., Cherkaoui, S., Danan, D., Dumont, S.: Inexact primal-dual active set method for solving elastodynamic frictional contact problems. Comput. Math. Appl. 82, 36–59 (2021)

35. Abide, S., Barboteu, M., Danan, D.: Analysis of two active set type methods to solve unilateral contact problems. Appl. Math. Comput. 284, 286–307 (2016)

36. Abide, S., Barboteu, M., Cherkaoui, S., Dumont, S.: A semi-smooth Newton and primal-dual active set method for non-smooth contact dynamics. Comput. Methods Appl. Mech. Eng. 387, 114153 (2021)

37. Moreau, J.J.: Standard inelastic shocks and the dynamics of unilateral constraints. In: Del Piero, G., Maceri, F. (eds.) Unilateral Problems in Structural Analysis. CISM Course and Lectures, vol. 288. pp. 173–221. Springer, Berlin (1985)

38. Jourdan, F., Alart, P., Jean, M.: A Gauss-Seidel like algorithm to solve frictional contact problems. Comput. Methods Appl. Mech. Eng. 155(1–2), 31–47 (1998)
## References

1. Laursen, T.A.: Computational Contact and Impact Mechanics: Fundamentals of Modeling Interfacial Phenomena in Nonlinear Finite Element Analysis. Springer, Berlin (2013)

2. Wriggers, P., Zavarise, G.: Computational Contact Mechanics. Encyclopedia of Computational Mechanics (2004)

3. Brogliato, B.: Impacts in Mechanical Systems: Analysis and Modelling. Springer, Berlin (2000)

4. Acary, V., Brogliato, B.: Numerical Methods for Nonsmooth Dynamical Systems: Applications in Mechanics and Electronics. Springer, Berlin (2008)

5. Moreau, J.J.: Application of convex analysis to some problems of dry friction. In: Trends in Applications of Pure Mathematics to Mechanics, pp. 263–280. Pitman, London (1977)

6. Moreau, J.J.: Unilateral contact and dry friction in finite freedom dynamics. In: Nonsmooth Mechanics and Applications, pp. 1–82. Springer, Berlin (1988)

7. Monteiro Marques, M.D.P.: Differential Inclusions in Nonsmooth Mechanical Problem: Shocks and Dry Friction. Springer, Basel (1993)

8. Jean, M.: The non-smooth contact dynamics method. Comput. Methods Appl. Mech. Eng. 177(3–4), 235–257 (1999)

9. Moreau, J.J.: Numerical aspects of the sweeping process. Comput. Methods Appl. Mech. Eng. 177(3–4), 329–349 (1999)

10. Jean, M., Moreau, J.J.: Unilaterality and dry friction in the dynamics of rigid body collections. In: 1st Contact Mechanics International Symposium, pp. 31–48 (1992)

11. Dubois, F., Acary, V., Jean, M.: The contact dynamics method: a nonsmooth story. C. R., Méc. 346(3), 247–262 (2018)

12. Moreau, J.J.: Some numerical methods in multibody dynamics: application to granular materials. Eur. J. Mech. A, Solids 13(4–suppl), 93–114 (1994)

13. de Saxcé, G., Feng, Z.-Q.: New inequality and functional for contact with friction: the implicit standard material approach. J. Struct. Mech. 19(3), 301–325 (1991)

14. Fortin, J., Millet, O., de Saxcé, G.: Numerical simulation of granular materials by an improved discrete element method. Int. J. Numer. Methods Eng. 62(5), 639–663 (2005)

15. Feng, Z.-Q., Joli, P., Cros, J.-M., Magnain, B.: The bi-potential method applied to the modeling of dynamic problems with friction. Comput. Mech. 36(5), 375–383 (2005)

16. Joli, P., Feng, Z.-Q.: Uzawa and Newton algorithms to solve frictional contact problems within the bi-potential framework. Int. J. Numer. Methods Eng. 73(3), 317–330 (2008)

17. Dumont, S.: On enhanced descent algorithms for solving frictional multicontact problems: application to the discrete element method. Int. J. Numer. Methods Eng. 93(11), 1170–1190 (2013)

18. Kikuchi, N., Oden, J.T.: Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods. SIAM Studies in Applied Mathematics, vol. 8. Philadelphia (1988)

19. Oden, J.T., Kim, S.J.: Interior penalty methods for finite element approximations of the Signorini problem in elastostatics. Comput. Math. Appl. 8(1), 35–56 (1986)

20. Alart, P., Curnier, A.: A mixed formulation for frictional contact problems prone to Newton like solution methods. Comput. Methods Appl. Mech. Eng. 92(3), 353–375 (1991)

21. Raous, M., Chabrand, P., Lebon, F.: Numerical methods for solving unilateral contact problem with friction. J. Theor. Appl. Mech. 7, 111–128 (1988)

22. Joli, P., Feng, Z.-Q.: Uzawa and Newton algorithms to solve frictional contact problems within the bi-potential framework. Int. J. Numer. Methods Eng. 73(3), 317–330 (2008)

23. Houly, F.: An adaptation of Nitsche's method to the Tresca friction problem. J. Math. Anal. Appl. 411(1), 329–339 (2014)

24. Houly, F., Hild, P., Renard, Y.: A Nitsche finite element method for dynamic contact: 1. Space semi-discretization and time-marching schemes. ESAIM: Math. Model. Numer. Anal. 49(2), 481–502 (2015)

25. Houly, F., Fabre, M., Hild, P., Mlika, R., Pousin, J., Renard, Y.: An overview of recent results on Nitsche's method for contact problems. In: Geometrically Unfitted Finite Element Methods and Applications, pp. 93–141 (2017)

26. Hintermüller, M., Kovtunenko, V.A., Kunisch, K.: Semismooth Newton methods for a class of unilaterally constrained variational problems. In: Universität Graz/Technische Universität Graz. SFB F003-Optimierung und Kontrolle (2003)

27. Hintermüller, M., Ito, K., Kunisch, K.: The primal-dual active set strategy as a semismooth Newton method. SIAM J. Optim. 13(3), 865–888 (2002)

28. Hintermüller, M., Kovtunenko, V.A., Kunisch, K.: Obstacle problems with cohesion: a hemivariational inequality approach and its efficient numerical solution. SIAM J. Optim. 21(2), 491–516 (2011)

29. Hüeber, S., Stadler, G., Wohlmuth, B.I.: A primal-dual active set algorithm for three-dimensional contact problems with Coulomb friction. SIAM J. Sci. Comput. 30(2), 572–596 (2008)

30. Hüeber, S., Wohlmuth, B.I.: A primal-dual active set strategy for non-linear multibody contact problems. Comput. Methods Appl. Mech. Eng. 194(27–29), 3147–3166 (2005)

31. Koziara, T., Bićanić, N.: Semismooth Newton method for frictional contact between pseudo-rigid bodies. Comput. Methods Appl. Mech. Eng. 197(33–40), 2763–2777 (2008)

32. Sharaf, I.M.: An active set algorithm for a class of linear complementarity problems arising from rigid body dynamics. Pak. J. Stat. Oper. Res. 339–352 (2016)

33. Barboteu, M., Dumont, S.: A primal-dual active set method for solving multi-rigid-body dynamic contact problems. Math. Mech. Solids 23(3), 489–503 (2018)

34. Abide, S., Barboteu, M., Cherkaoui, S., Danan, D., Dumont, S.: Inexact primal-dual active set method for solving elastodynamic frictional contact problems. Comput. Math. Appl. 82, 36–59 (2021)

35. Abide, S., Barboteu, M., Danan, D.: Analysis of two active set type methods to solve unilateral contact problems. Appl. Math. Comput. 284, 286–307 (2016)

36. Abide, S., Barboteu, M., Cherkaoui, S., Dumont, S.: A semi-smooth Newton and primal-dual active set method for non-smooth contact dynamics. Comput. Methods Appl. Mech. Eng. 387, 114153 (2021)

37. Moreau, J.J.: Standard inelastic shocks and the dynamics of unilateral constraints. In: Del Piero, G., Maceri, F. (eds.) Unilateral Problems in Structural Analysis. CISM Course and Lectures, vol. 288. pp. 173–221. Springer, Berlin (1985)

38. Jourdan, F., Alart, P., Jean, M.: A Gauss-Seidel like algorithm to solve frictional contact problems. Comput. Methods Appl. Mech. Eng. 155(1–2), 31–47 (1998)
39. Renouf, M., Alart, P.: Conjugate gradient type algorithms for frictional multi-contact problems: applications to granular materials. Comput. Methods Appl. Mech. Eng. 194(18–20), 2019–2041 (2005)

40. Ciarlet, P.G., Geymonat, G.: Sur les lois de comportement en élasticité non-linéaire compressible. C. R. Acad. Sci. 295, 423–426 (1982)

41. Cundall, P.A., Strack, O.D.L.: A discrete numerical model for granular assemblies. Geotechnique 29(1), 47–65 (1979)

42. Garg, R., Galvin, J., Li, T., Pannala, S.: Documentation of open-source MFIX-DEM software for gas-solids flows. From https://mfix.netl.doe.gov/download/mfix/mfix_current_documentation/dem_doc_2012-1.pdf (2012)

Submit your manuscript to a SpringerOpen journal and benefit from:

► Convenient online submission

► Rigorous peer review

► Open access: articles freely available online

► High visibility within the field

► Retaining the copyright to your article

Submit your next manuscript at ► springeropen.com
39. Renouf, M., Alart, P.: Conjugate gradient type algorithms for frictional multi-contact problems: applications to granular materials. Comput. Methods Appl. Mech. Eng. 194(18–20), 2019–2041 (2005)

40. Ciarlet, P.G., Geymonat, G.: Sur les lois de comportement en élasticité non-linéaire compressible. C. R. Acad. Sci. 295, 423–426 (1982)

41. Cundall, P.A., Strack, O.D.L.: A discrete numerical model for granular assemblies. Geotechnique 29(1), 47–65 (1979)

42. Garg, R., Galvin, J., Li, T., Pannala, S.: Documentation of open-source MFIX-DEM software for gas-solids flows. From https://mfix.netl.doe.gov/download/mfix/mfix_current_documentation/dem_doc_2012-1.pdf (2012)

Submit your manuscript to a SpringerOpen journal and benefit from:

► Convenient online submission

► Rigorous peer review

► Open access: articles freely available online

► High visibility within the field

► Retaining the copyright to your article

Submit your next manuscript at ► springeropen.com