# Uzawa and Newton algorithms to solve frictional contact problems within the bi-potential framework

Pierre Joli, Zhi-Qiang Feng

► To cite this version:

Pierre Joli, Zhi-Qiang Feng. Uzawa and Newton algorithms to solve frictional contact problems within the bi-potential framework. International Journal for Numerical Methods in Engineering, 2008, 73 (3), pp.317-330. 10.1002/nme.2073. hal-00342939

HAL Id: hal-00342939

https://hal.science/hal-00342939v1

Submitted on 29 Oct 2023

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Uzawa and Newton algorithms to solve frictional contact problems within the bi-potential framework

P.Joli<fnref n="1" /> and Z.-Q. Feng<fnref n="2" />, *, †

This paper is concerned with the numerical modeling of three-dimensional unilateral contact problems in elastostatics with Coulomb friction laws. We propose a Newton-like algorithm to solve the local contact non-linear equations within the bi-potential framework. The piecewise continuous contact tangent matrices are explicitly derived. A comparative study is made between the Newton algorithm and the previously developed Uzawa algorithm. A test example is included to demonstrate the developed algorithms and to highlight their performance.

KEY WORDS: contact and friction; Uzawa algorithm; Newton algorithm; bi-potential method

## 1. INTRODUCTION

Problems involving contact and friction are among the most difficult ones in mechanics and at the same time are of crucial practical importance in many engineering branches. The main mathematical difficulty lies in the severe contact non-linearities because the natural first-order constitutive laws of contact and friction phenomena are expressed by non-smooth multivalued force-displacement or force-velocity relations. In the last decade, substantial progress has been made in the analysis of contact problems using finite element procedures. A large number of algorithms for the numerical solution of the related finite element equations and inequalities have been presented in the literature. Review papers may be consulted for an extensive list of References [1, 2]. See also the monographs by Kikuchi and Oden [3], Zhong [4], Wriggers [5] and Laursen [6]. The popular penalty approximation and ‘mixed’ or ‘trial-and-error’ methods [7, 8] appear, at

Contract/grant sponsor: Publishing Arts Research Council; contract/grant number: 98-1846389

<footnote n="1">¹Laboratoire IBISC, Université d'Évry, 40 rue du Pelvoux, Évry 91020, France</footnote>

<footnote n="2">²Laboratoire de Mécanique d'Évry, Université d'Évry, 40 rue du Pelvoux, Évry 91020, France</footnote>

<footnote>*Correspondence to: Z.-Q. Feng, Université d'Évry, LMEE, 40 rue du Pelvoux, Évry 91020, France.</footnote>

<footnote>†E-mail: feng@iup.univ-evry.fr</footnote>
first glance, suitable for many applications. But in these kinds of methods, the contact boundary conditions and friction laws are not satisfied accurately and it is difficult for the users to choose appropriate penalty factors. They may fail for stiff problems because of unpleasant numerical oscillations between contact statuses. The augmented Lagrangian method first appeared to deal with constrained minimization problems. Since friction problems are not minimization problems, the formulation needs to be extended. Alart and Curnier [9], Simo and Laursen [10] and De Saxcé and Feng [11] have obtained some extensions in mutually independent works. Heegaard and Curnier [12] presented a generalized Newton method for large slip frictionless contact problems. From a simple spring-wall contact problem, they have also provided some interesting discussion on convergence properties of the generalized Newton method compared to a global Uzawa algorithm. Recently, Dostál et al. have solved three-dimensional frictional contact problems by using a FETI-based domain decomposition method [13]. The contact problem is formulated as a constrained quadratic programming problem and the contact conditions are enforced by penalty and augmented Lagrangian methods. The bi-potential method proposed by De Saxcé and Feng provides a powerful tool to model dissipative constitutive laws such as Coulomb friction laws. The application of the augmented Lagrangian method to the contact laws leads to implicit equations of projection onto the Coulomb friction cone, strictly equivalent to the original contact inequality [14]. An iterative Uzawa algorithm can be used to solve the non-linear implicit equations and this algorithm has been successfully applied by Feng [15] and Feng et al. [16] to simulate large deformation contact problems.

The aim of the present paper is to develop a Newton-like algorithm to solve the local contact non-linear equations within the bi-potential framework. Characteristics of Uzawa and Newton algorithms are discussed. A test numerical example is performed in this study to show the validity of the developed algorithms.

## 2. PROBLEM SETTING

### 2.1. Governing equations

The finite element method is often used in computational mechanics. Without going into details, quasi-static non-linear problems involving contact are governed by the following discretized equation:

$$ \mathbf{F}_{\text{int}} + \mathbf{F}_{\text{ext}} + \mathbf{R} = 0 \qquad (1) $$

where $\mathbf{F}_{\text{int}}$ is the vector of internal forces, $\mathbf{F}_{\text{ext}}$ denotes the vector of external loads and $\mathbf{R}$ the vector of contact reaction forces.

This equation is strongly non-linear with respect to the nodal displacements $\mathbf{U}$, because of finite strains and large displacements of solid. Moreover, the constitutive laws of contact with friction are usually represented by inequalities and the contact potential is even non-differentiable as we will see in Sections 2.3 and 2.4. A typical solution procedure for this type of non-linear analysis is obtained by using the Newton-Raphson iterative procedure [17]:

$$ \begin{aligned} \mathbf{K}_{\text{T}}^{i} \Delta \mathbf{U} &= \mathbf{F}_{\text{int}}^{i} + \mathbf{F}_{\text{ext}} + \mathbf{R} \\ \mathbf{U}^{i+1} &= \mathbf{U}^{i} + \Delta \mathbf{U} \end{aligned} \qquad (2) $$
where $i$ and $i+1$ are the iteration numbers at which the equations are computed. $\mathbf{K}_T$ is the tangent stiffness matrix and $\Delta\mathbf{U}$ the vector of nodal displacements correction. Taking the derivative of $\mathbf{F}_{\text{int}}$ with respect to the nodal displacements $\mathbf{U}$ gives the tangent stiffness matrix as

$$ \mathbf{K}_{\text{T}}^{i} = - \frac{\partial \mathbf{F}_{\text{int}}^{i}}{\partial \mathbf{U}^{i}} \qquad (3) $$

It is noted that Equation (2) cannot be solved directly because $\Delta\mathbf{U}$ and $\mathbf{R}$ are both unknown. The key idea is to determine first the reaction vector $\mathbf{R}$ in a reduced system which only concerns the contact nodes. Then, the displacement increments $\Delta\mathbf{U}$ can be computed in the whole structure, using contact reactions as external loading. In the following, we focus our attention on describing how to determine the contact forces. Let us begin with the general description of contact kinematics.

## 2.2. Contact kinematics

First of all, basic definitions and notations used are described. For the sake of simplicity, we consider two deformable bodies $\Omega^a$ (Figure 1), $a = 1, 2$, coming into contact. Each body is decomposed by finite elements and the nodal positions in the global coordinate frame are represented by the vector $\mathbf{X}_a$. The boundary $\Gamma^a$ of each body is assumed to be sufficiently smooth everywhere such that an outward unit normal vector, denoted by $\mathbf{N}_a$, can be defined at any point $P_a$ on $\Gamma^a$. Moreover, it is possible to determine position vectors of each contact point $P_1$ (resp. $P_2$) from an interpolation matrix $\mathbf{B}_1$ (resp. $\mathbf{B}_2$) as follows:

$$ \mathbf{X}(P_1) = \mathbf{B}_1 \mathbf{X}_1, \quad \mathbf{X}(P_2) = \mathbf{B}_2 \mathbf{X}_2 \qquad (4) $$

We consider only the case with $N_c$ contact nodes $P_1^\alpha$ ($\alpha = 1, N_c$) defined on $\Gamma^1$ and $P_2^\alpha$ are target points defined by the normal projection of $P_1^\alpha$ onto $\Gamma^2$. We can build the relative position between $P_1^\alpha$ and $P_2^\alpha$ by

$$ \mathbf{X}^{\alpha} = \mathbf{X}(P_{2}^{\alpha}) - \mathbf{X}(P_{1}^{\alpha}) \qquad (5) $$

We consider a local orthogonal reference frame by means of three vectors (in algebraic form) $\mathbf{T}_1, \mathbf{T}_2$ and $\mathbf{N}$ which are defined with respect to the global reference frame. We set the following notation:

$$ \mathbf{X}^{\alpha} = x_{t_1}^{\alpha} \mathbf{T}_1 + x_{t_2}^{\alpha} \mathbf{T}_2 + x_{n}^{\alpha} \mathbf{N} \Leftrightarrow \mathbf{x}^{\alpha} = \begin{cases} x_{t_1}^{\alpha} = \mathbf{T}_1^{\text{T}} \mathbf{X}^{\alpha} \\ x_{t_2}^{\alpha} = \mathbf{T}_2^{\text{T}} \mathbf{X}^{\alpha} \\ x_{n}^{\alpha} = \mathbf{N}^{\text{T}} \mathbf{X}^{\alpha} \end{cases} \qquad (6) $$

<figure><img src="image_2.png" /><figcaption>Figure 1. Contact kinematics.</figcaption></figure>
and combining Equations (4)–(6), we have

$$ \mathbf{x}^{\alpha} = \mathbf{H}_{\alpha} \mathbf{X} \qquad (7) $$

with

$$ \mathbf{H}_{\alpha} = [-\mathbf{H}_{\alpha}^{1} | \mathbf{H}_{\alpha}^{2}], \quad \mathbf{H}_{\alpha}^{k} = \begin{bmatrix} \mathbf{T}_{1}^{\mathrm{T}} \mathbf{B}_{k} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{T}_{2}^{\mathrm{T}} \mathbf{B}_{k} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{N}^{\mathrm{T}} \mathbf{B}_{k} \end{bmatrix}, \quad \mathbf{X} = \begin{Bmatrix} \mathbf{X}_{1} \\ \mathbf{X}_{2} \end{Bmatrix} \qquad (8) $$

The incremental form of Equation (7) gives the gap vector between $P_1^\alpha$ and $P_2^\alpha$:

$$ \mathbf{x}_{i+1}^{\alpha} = \mathbf{H}_{\alpha} \Delta \mathbf{X}_{i} + \mathbf{g}^{\alpha} \qquad (9) $$

where $\mathbf{g}^\alpha = (0 \ 0 \ g)$ represents the initial gap vector which is determined by a contact collision detector. We have opted to carry out this operation at the beginning of each load step.

In Equation (9), we neglect the variation of the normal at the contact point during one load step. If this assumption is not satisfied, it is possible to reduce the load step or to perform the collision detection at each iteration. In the case of discontinuous curvature of contact surfaces, special smoothing techniques can be used as in [12, 18].

In the local reference frame the contact force $\mathbf{r}^\alpha$ and the gap vector $\mathbf{x}^\alpha$ can be defined by

$$ \mathbf{r}^{\alpha} = \mathbf{r}_{i}^{\alpha} + r_{n}^{\alpha} \mathbf{n} = r_{i}^{\alpha} \mathbf{t}_{1} + r_{i}^{\alpha} \mathbf{t}_{2} + r_{n}^{\alpha} \mathbf{n} \qquad (10) $$

$$ \mathbf{x}^{\alpha} = \mathbf{x}_{i}^{\alpha} + x_{n}^{\alpha} \mathbf{n} = x_{i}^{\alpha} \mathbf{t}_{1} + x_{i}^{\alpha} \mathbf{t}_{2} + x_{n}^{\alpha} \mathbf{n} \qquad (11) $$

with

$$ \mathbf{t}_{1}^{\mathrm{T}} = (100), \quad \mathbf{t}_{2}^{\mathrm{T}} = (010), \quad \mathbf{n}^{\mathrm{T}} = (001) \qquad (12) $$

Application of the contact virtual work results in

$$ \mathbf{r}^{\alpha \mathrm{T}} \delta \mathbf{x}^{\alpha} = \mathbf{R}^{\alpha \mathrm{T}} \delta \mathbf{X}^{\alpha} \qquad (13) $$

We obtain from Equation (7) then

$$ \mathbf{R}^{\alpha} = \mathbf{H}_{\alpha}^{\mathrm{T}} \mathbf{r}^{\alpha} \qquad (14) $$

From now on, we omit deliberately the underscript $i$ (or $i+1$) because the principle of computation of contact force presented later is the same at each iteration. In order to set problems with multiple contact points, we introduce the following complementary notations:

$$ \mathbf{H} = \begin{bmatrix} \mathbf{H}_1 \\ \vdots \\ \mathbf{H}_\alpha \end{bmatrix}, \quad \mathbf{x} = \begin{Bmatrix} \mathbf{x}^1 \\ \vdots \\ \mathbf{x}^{N_c} \end{Bmatrix}, \quad \mathbf{r} = \begin{Bmatrix} \mathbf{r}^1 \\ \vdots \\ \mathbf{r}^{N_c} \end{Bmatrix}, \quad \mathbf{g} = \begin{Bmatrix} \mathbf{g}^1 \\ \vdots \\ \mathbf{g}^{N_c} \end{Bmatrix} \qquad (15) $$

with these notations in hand, we have the following equations:

$$ \mathbf{R} = \sum_{\alpha=1}^{N_c} \mathbf{R}^{\alpha} = \mathbf{H}^{\mathrm{T}} \mathbf{r}, \quad \mathbf{x} = \mathbf{H} \Delta \mathbf{X} + \mathbf{g} \qquad (16) $$
### 2.3. Signorini conditions and Coulomb friction laws

To simplify the notations, the superscript $\alpha$ is omitted in the description of the contact laws. The unilateral contact law is characterized by a geometric condition of non-penetration, a static condition of no-adhesion and a mechanical complementary condition. These three conditions known as Signorini conditions are expressed, for each contact point, in terms of the signed contact distance $x_n$ and the normal contact force $r_n$ by

$$ \text{Signor}(x_n, r_n) \Leftrightarrow x_n \ge 0, \quad r_n \ge 0 \text{ and } x_n r_n = 0 \quad (17) $$

Classically, a rate-independent dry-friction law is characterized by a kinematic slip rule. In this work, the classic Coulomb friction rule is used and defined by

$$ \text{Coul}(\dot{\mathbf{x}}_t, \mathbf{r}_t) \Leftrightarrow \phi = \|\mathbf{r}_t\| - \mu r_n \le 0, \quad \dot{\mathbf{x}}_t = -\lambda \frac{\mathbf{r}_t}{\|\mathbf{r}_t\|}, \quad \lambda \ge 0 \text{ and } \phi \lambda = 0 \quad (18) $$

or equivalently

$$ \text{Coul}(\dot{\mathbf{x}}_t, \mathbf{r}_t) \Leftrightarrow \text{if } \|\dot{\mathbf{x}}_t\| = 0 \quad \text{then} \quad \|\mathbf{r}_t\| \le \mu r_n \quad \text{else} \quad \mathbf{r}_t = -\mu r_n \frac{\dot{\mathbf{x}}_t}{\|\dot{\mathbf{x}}_t\|} \quad (19) $$

where $\mu$ is the coefficient of friction and the superposed dot denotes time derivative. The set of admissible forces, denoted by $K_\mu$, is defined by

$$ K_\mu = \{ \mathbf{r} \in \mathbb{R}^3 \text{ such that } \|\mathbf{r}_t\| - \mu r_n \le 0 \} \quad (20) $$

$K_\mu$ is the so-called Coulomb cone and is convex.

In this work, we deal with the quasi-static contact problem with friction. For a given time (or loading) history $\theta \in [0, T]$, where $[0, T]$ is a time interval which can be partitioned into $N$ sub-intervals of size $\Delta\theta$, we adopt a backward-Euler time discretization of the time derivative $\dot{\mathbf{x}}_t$ as follows:

$$ \dot{\mathbf{x}}_t \approx \frac{\mathbf{x}_{t_i} - \mathbf{x}_{t_0}}{\Delta\theta} \quad (21) $$

In quasi-static cases, $N$ is the total number of load steps and we can set $\Delta\theta = 1$. As we have discussed above, at each load step, a contact detection is performed. In this way, we have $\mathbf{x}_{t_0} = \mathbf{0}$. Then, Equation (21) reduces to $\dot{\mathbf{x}}_t \approx \mathbf{x}_t$ by omitting the underscript iteration number $i$.

The complete contact law (Signorini conditions + Coulomb friction laws) is thus a complex non-smooth dissipative law including three statuses:

$$ \text{No contact: } x_n > 0 \text{ and } \mathbf{r}=0 $$

$$ \text{Contact with sticking: } \|\mathbf{x}_t\| = 0 \text{ and } \mathbf{r} \in \text{int}(K_\mu) \quad (22) $$

$$ \text{Contact with sliding: } \|\mathbf{x}_t\| \neq 0 \text{ and } \mathbf{r} \in \text{bd}(K_\mu) \text{ with } \mathbf{r}_t = -\mu r_n \frac{\mathbf{x}_t}{\|\mathbf{x}_t\|} $$

where ‘int($K_\mu$)’ and ‘bd($K_\mu$)’ denote the interior and the boundary of $K_\mu$, respectively. The multivalued character of the law lies in the first and the second part of the statement. If $r_n$ is null then $\mathbf{x}$ is arbitrary but its normal component $x_n$ should be positive. In other words, one single element of $\mathbb{R}^3$ ($\mathbf{r}=\mathbf{0}$) is associated with an infinite number of gap vectors $\mathbf{x} \in \mathbb{R}^3$. The same arguments can be developed for the second part of the statement. For this reason, the contact forces cannot be derived from a potential function of gap vectors.
## 2.4. The bi-potential method

De Saxcé and Feng [14] have proposed a contact bi-potential as follows:

$$b_c(-\mathbf{x}, \mathbf{r}) = \bigcup_{\mathbb{R}_-} (-x_n) + \bigcup_{K_\mu} (\mathbf{r}) + \mu r_n \|\mathbf{x}_t\| \quad (23)$$

where $\mathbb{R}_- = ]-\infty, 0]$ is the set of the negative and null real numbers. $\bigcup_{K_\mu} \mathbf{r}$ denotes the so-called indicator function of the closed convex set $K_\mu$:

$$\bigcup_{K_\mu}(\mathbf{r}) = \begin{cases} 0 & \text{if } \mathbf{r} \in K_\mu \\ +\infty & \text{otherwise} \end{cases} \quad (24)$$

Then, the complete contact laws can be written in a compact form of implicit subnormality rules

$$-\mathbf{x} \in \partial_r b_c(-\mathbf{x}, \mathbf{r}), \quad \mathbf{r} \in \partial_{-\mathbf{x}} b_c(-\mathbf{x}, \mathbf{r}) \quad (25)$$

In order to avoid non-differentiable potentials, it is convenient to use the Augmented Lagrangian Method [9–11, 14, 19]. Thus, the above inclusion is equivalent to the following projection operation:

$$\mathbf{r} = \text{Proj}_{K_\mu}(\mathbf{r}^*) \quad (26)$$

where $\mathbf{r}^*$ is the so-called augmented contact forces vector and is given by

$$\mathbf{r}^* = \mathbf{r} - \rho \mathbf{x}^* \quad \text{with} \quad \mathbf{x}^* = \mathbf{x} + \mu \|\mathbf{x}_t\| \mathbf{n} \quad (27)$$

where $\rho$ is an arbitrary positive parameter.

The three possible contact statuses as mentioned in Equation (22) are illustrated in Figure 2. Within the bi-potential framework, these statuses can be stated as: $\mathbf{r}^* \in K_\mu$ (contact with sticking), $\mathbf{r}^* \in K_\mu^*$ (separating) and $\mathbf{r}^* \in \mathbb{R}^3 - (K_\mu \cup K_\mu^*)$ (contact with sliding). $K_\mu^*$ is the polar cone of $K_\mu$.

<figure><img src="image_3.png" /><figcaption>Figure 2. The Coulomb cone and contact projection operators.</figcaption></figure>
Consequently, the projection operation can be explicitly defined by

$$
\begin{align}
\operatorname{Proj}_{\mathbf{K}_{\mu}}(\mathbf{r}^*) &= \mathbf{r}^* && \text{if } \|\mathbf{r}_t^*\| < \mu r_n^* \nonumber \\
\operatorname{Proj}_{\mathbf{K}_{\mu}}(\mathbf{r}^*) &= \mathbf{0} && \text{if } \mu \|\mathbf{r}_t^*\| < -r_n^* \tag{28} \\
\operatorname{Proj}_{\mathbf{K}_{\mu}}(\mathbf{r}^*) &= \mathbf{r}^* && - \left( \frac{\|\mathbf{r}_t^*\| - \mu r_n^*}{1 + \mu^2} \right) \left( \frac{\mathbf{r}_t^*}{\|\mathbf{r}_t^*\|} - \mu \mathbf{n} \right) && \text{otherwise} \nonumber
\end{align}
$$

3. SOLUTION METHODS

3.1. Equilibrium equations of contact points

We have defined above the governing equations for the contact problems. In short, the system of
equations to be solved can be summarized as follows:

$$
\begin{align}
\mathbf{K}_T \Delta \mathbf{U} &= \mathbf{F}_{\text{int}} + \mathbf{F}_{\text{ext}} + \mathbf{H}^T \mathbf{r} \nonumber \\
\mathbf{x} &= \mathbf{H} \Delta \mathbf{U} + \mathbf{g} \tag{29} \\
\mathbf{r}^{\alpha} &= \operatorname{Proj}_{\mathbf{K}_{\mu}} (\mathbf{r}^{*\alpha}) \quad (\alpha = 1, N_c) \nonumber
\end{align}
$$

Eliminating ΔU leads to the following reduced system of equations:

$$
\mathbf{x} = \mathbf{W}\mathbf{r} + \tilde{\mathbf{x}} \tag{30}
$$

with

$$
\mathbf{W} = \mathbf{H} \mathbf{K}_{\mathrm{T}}^{-1} \mathbf{H}^{\mathrm{T}} \quad \text{and} \quad \tilde{\mathbf{x}} = \mathbf{H} \mathbf{K}_{\mathrm{T}}^{-1} (\mathbf{F}_{\mathrm{int}} + \mathbf{F}_{\mathrm{ext}}) + \mathbf{g} \tag{31}
$$

Solving the contact problem leads thus to the following problem:

$$
\text{Find } \chi \text{ such that } \mathbf{f}(\chi) = \mathbf{0} \tag{32}
$$

with

$$
\chi = \begin{cases} \chi_1 \\ \vdots \\ \chi_{N_c} \end{cases} \quad \text{and} \quad \chi_\alpha = \begin{cases} \mathbf{r}^\alpha \\ \mathbf{x}^\alpha \end{cases} \tag{33}
$$

$$
\mathbf{f}(\boldsymbol{\chi}) = \begin{Bmatrix} \mathbf{f}_1(\boldsymbol{\chi}) \\ \vdots \\ \mathbf{f}_{N_c}(\boldsymbol{\chi}) \end{Bmatrix} \quad \text{and} \quad \mathbf{f}_{\alpha}(\boldsymbol{\chi}) = \begin{Bmatrix} \boldsymbol{x}^{\alpha} - \sum_{\beta=1}^{N_c} \mathbf{W}_{\alpha\beta} \boldsymbol{r}^{\beta} - \tilde{\boldsymbol{x}}^{\alpha} \\ \mathbf{Z}_{\alpha} \end{Bmatrix} \quad (34)
$$

$$
\mathbf{Z}_{\alpha} = \mathbf{r}^{\alpha} - \operatorname{Proj}_{\mathbf{K}_{\mu}} (\mathbf{r}^{*\alpha}) \quad (35)
$$

$$
\tilde{\mathbf{x}}^{\alpha} = \mathbf{H}_{\alpha} \mathbf{K}_{\mathrm{T}}^{-1} (\mathbf{F}_{\mathrm{int}} + \mathbf{F}_{\mathrm{ext}}) + \mathbf{g}^{\alpha} \quad (36)
$$

$$
\mathbf{W}_{\beta\alpha\beta} = \mathbf{H}_{\alpha}\mathbf{K}_{\mathrm{T}}^{-1}\mathbf{H}^{\mathrm{T}} \quad (37)
$$
## 3.2. Global solution: non-linear Gauss-Seidel-like algorithm

Jourdan et al. [20] have applied a non-linear Gauss-Seidel-like algorithm to simulate deep drawing problems. Signorini conditions and Coulomb friction laws are derived from two distinct pseudo-potentials. As we will see below, this algorithm can be readily extended to be applied to the case of the bi-potential formulation.

The principle of this algorithm is to decompose the global solution of the $(6 \times N_c)$ equations (32) into $N_c$ successive local solutions of the six following equations:

$$ \mathbf{f}_{\alpha}(\boldsymbol{x}) = \left\{ \begin{array}{l} \mathbf{x}^{\alpha} - \mathbf{W}_{\alpha\alpha}\mathbf{r}^{\alpha} - \mathbf{x}^{\alpha\beta} \\ \mathbf{Z}_{\alpha} \end{array} \right\} = \mathbf{0} \qquad (38) $$

with

$$ \mathbf{x}^{\alpha\beta} = \sum_{\beta=1, \beta \neq \alpha}^{N_c} \mathbf{W}_{\alpha\beta}\mathbf{r}^{\beta} + \tilde{\mathbf{x}}^{\alpha} \qquad (39) $$

where $\mathbf{x}^{\alpha\beta}$ represents the part of the relative position at the contact point $\alpha$ due to the initial gap, the external forces and contact forces of $N_c - 1$ other contact nodes $\beta$. This contribution is 'frozen' during each local solution. One series of $N_c$ local solutions corresponds to one iteration $k$ of the algorithm. The iterative process is successively applied for each contact point ($\alpha = 1, N_c$) until the convergence of solution. The contact convergence criterion is stated as

$$ \frac{\|\mathbf{r}^{(k+1)} - \mathbf{r}^{(k)}\|}{\|\mathbf{r}^{(k+1)}\|} \le \varepsilon_g \qquad (40) $$

where $\mathbf{r} = \{\mathbf{r}^1 \mathbf{r}^2 \cdots \mathbf{r}^{N_c}\}$ is the vector of contact reactions of all contact nodes and $\varepsilon_g$ is a user-defined tolerance. The initial condition is given by $\mathbf{r}^{(0)} = \mathbf{0}$.

In the bi-potential formulation, the usual approach to solve the local implicit equations (38) is to use a predictor/corrector Uzawa algorithm. Many examples have been successfully treated by Feng [15] and Feng et al. [16]. The advantage of this approach is the simplicity of programming and the numerical robustness, but it needs more iterations when compared with the implicit Newton algorithm. This point will be discussed later. In the following sections, we present both the Uzawa algorithm and the Newton algorithm for the solution of implicit equations (38). It is noted that, for the first time, the Newton algorithm is applied in the context of the bi-potential framework.

## 3.3. Local solution: Uzawa algorithm

Numerical solution of implicit equation (38) can be carried out by means of the Uzawa algorithm, which leads thus to an iterative process involving one predictor-corrector step:

$$ \begin{aligned} \text{Predictor } \mathbf{r}^{\alpha*(k+1)} &= \mathbf{r}^{\alpha(k)} - \rho^{(k)}(\mathbf{x}^{\alpha\beta(k)} + \mu\|\mathbf{x}_t^{\alpha\beta(k)}\|\mathbf{n}) \\ \text{Corrector } \mathbf{r}^{\alpha(k+1)} &= \text{Proj}_{\mathbf{K}_{\mu}}(\mathbf{r}^{\alpha*(k+1)}) \end{aligned} \qquad (41) $$

where $k$ and $k + 1$ are the iteration numbers at which the contact reactions are computed. The corrector step is explicitly given by Equation (28). In view of Equation (38), the gap vector is updated by

$$ \mathbf{x}^{\alpha(k+1)} = \mathbf{W}_{\alpha\alpha}\mathbf{r}^{\alpha(k+1)} + \mathbf{x}^{\alpha\beta(k)} \qquad (42) $$
It is noted that the solution is controlled by a global convergence criterion (iteration $k$) as stated in Equation (40).

### 3.4. Local solution: Newton algorithm

The Newton algorithm applied to implicit equations (38) leads to the following iterative numerical scheme:

Step 1: initialize

$$ \boldsymbol{\chi}_{\alpha}^{0} = \left\{ \begin{array}{l} \mathbf{r}^{\alpha(0)} = \mathbf{0} \\ \mathbf{x}^{\alpha(0)} = \mathbf{x}^{\alpha\beta(k)} \end{array} \right\}, \quad i = 0 $$

Step 2: augment contact forces

$$ \mathbf{r}^{*\alpha(i)} = \mathbf{r}^{\alpha(i)} - \rho \mathbf{x}^{*\alpha(i)} \quad \text{with} \quad \mathbf{x}^{*\alpha(i)} = \mathbf{x}^{\alpha(i)} + \mu \|\mathbf{x}_t^{\alpha(i)}\| \mathbf{n} \quad (43) $$

Step 3: solve

$$ \begin{gathered} \left[ \frac{\partial \mathbf{f}_{\alpha}(\boldsymbol{\chi}^i)}{\partial \boldsymbol{\chi}_{\alpha}} \right] \Delta \boldsymbol{\chi}_{\alpha} = -\mathbf{f}_{\alpha}(\boldsymbol{\chi}^i) \quad (44) \\ \boldsymbol{\chi}_{\alpha}^{i+1} = \boldsymbol{\chi}_{\alpha}^i + \Delta \boldsymbol{\chi}_{\alpha} \end{gathered} $$

Step 4: check convergence

$$ \begin{cases} \text{if } \|\Delta\boldsymbol{\chi}_{\alpha}\| / \|\boldsymbol{\chi}_{\alpha}^i\| > \varepsilon_1 \text{ then set } i = i + 1 \text{ and goto Step 2} \\ \text{else } \boldsymbol{\chi}_{\alpha}^{k+1} = \boldsymbol{\chi}_{\alpha}^{i+1} \end{cases} \quad (45) $$

where $\boldsymbol{\chi}_{\alpha}^{k+1}$ are the variables of the contact point $\alpha$ taken at the $(k+1)$th iteration of the global solution in the Gauss-Seidel algorithm. $\varepsilon_1$ is the convergence tolerance in the local solution of the Newton algorithm.

$[\partial \mathbf{f}_{\alpha}(\boldsymbol{\chi}^i) / \partial \boldsymbol{\chi}_{\alpha}]$ represents the $(6 \times 6)$ tangent matrix of the local equations at the contact point $\alpha$ and has the general form

$$ \left[ \frac{\partial \mathbf{f}_{\alpha}(\boldsymbol{\chi}^i)}{\partial \boldsymbol{\chi}_{\alpha}} \right] = \begin{bmatrix} -\mathbf{W}_{\alpha\alpha} & \mathbf{I}_{3 \times 3} \\ \mathbf{A}_{\alpha} & \mathbf{B}_{\alpha} \end{bmatrix} \quad (46) $$

with

$$ \mathbf{A}_{\alpha} = \left[ \frac{\partial \mathbf{Z}_{\alpha}}{\partial r_{f_1}^{\alpha}} \quad \frac{\partial \mathbf{Z}_{\alpha}}{\partial r_{f_2}^{\alpha}} \quad \frac{\partial \mathbf{Z}_{\alpha}}{\partial r_{f_n}^{\alpha}} \right], \quad \mathbf{B}_{\alpha} = \left[ \frac{\partial \mathbf{Z}_{\alpha}}{\partial x_{f_1}^{\alpha}} \quad \frac{\partial \mathbf{Z}_{\alpha}}{\partial x_{f_2}^{\alpha}} \quad \frac{\partial \mathbf{Z}_{\alpha}}{\partial x_{f_n}^{\alpha}} \right] \quad (47) $$

There are two explicit forms of the matrices $\mathbf{A}$ and $\mathbf{B}$ corresponding to sticking and sliding contact statuses defined by Equation (28). Consequently, the components of these two matrices are explicitly given in two cases.

#### Case A: Contact with sliding

This case occurs if and only if

$$ \mu \|\mathbf{r}_t^{*\alpha}\| \ge -r_n^{*\alpha} \quad \text{and} \quad \|\mathbf{r}_t^{*\alpha}\| \ge \mu r_n^{*\alpha} \quad (48) $$
In view of Equations (28) and (35), we have

$$Z_{\alpha} = \rho \mathbf{x}^{*\alpha} + \left( \frac{\|\mathbf{r}_t^{*\alpha}\| - \mu r_n^{*\alpha}}{1 + \mu^2} \right) \left( \frac{\mathbf{r}_t^{*\alpha}}{\|\mathbf{r}_t^{*\alpha}\|} - \mu \mathbf{n} \right) \quad (49)$$

The derivation of $Z_\alpha$ with respect to the gap vector and contact forces gives

$$\frac{\partial Z_{\alpha}}{\partial x_{n}^{\alpha}} = \rho \mathbf{n} + \frac{\rho \mu}{1 + \mu^2} \left( \frac{\mathbf{r}_{t}^{*\alpha}}{\|\mathbf{r}_{t}^{*\alpha}\|} - \mu \mathbf{n} \right) \quad (50)$$

$$\frac{\partial Z_{\alpha}}{\partial x_{t_1}^{\alpha}} = \rho \left( \mathbf{t}_1 + \frac{\mu x_{t_1}^{\alpha}}{\|\mathbf{x}_t^{*\alpha}\|} \mathbf{n} \right) - \frac{\rho}{1 + \mu^2} \left[ \left( \frac{r_{t_1}^{*\alpha}}{\|\mathbf{r}_t^{*\alpha}\|} - \frac{\mu^2 x_{t_1}^{\alpha}}{\|\mathbf{x}_t^{*\alpha}\|} \right) \left( \frac{\mathbf{r}_t^{*\alpha}}{\|\mathbf{r}_t^{*\alpha}\|} - \mu \mathbf{n} \right) + \left( \|\mathbf{r}_t^{*\alpha}\| - \mu r_n^{*\alpha} \right) \left( \frac{1}{\|\mathbf{r}_t^{*\alpha}\|} \mathbf{t}_1 - \frac{r_{t_1}^{*\alpha}}{\|\mathbf{r}_t^{*\alpha}\|^3} \mathbf{r}_t^{*\alpha} \right) \right] \quad (51)$$

$$\frac{\partial Z_{\alpha}}{\partial x_{t_2}^{\alpha}} = \rho \left( \mathbf{t}_2 + \frac{\mu x_{t_2}^{\alpha}}{\|\mathbf{x}_t^{*\alpha}\|} \mathbf{n} \right) - \frac{\rho}{1 + \mu^2} \left[ \left( \frac{r_{t_2}^{*\alpha}}{\|\mathbf{r}_t^{*\alpha}\|} - \frac{\mu^2 x_{t_2}^{\alpha}}{\|\mathbf{x}_t^{*\alpha}\|} \right) \left( \frac{\mathbf{r}_t^{*\alpha}}{\|\mathbf{r}_t^{*\alpha}\|} - \mu \mathbf{n} \right) + \left( \|\mathbf{r}_t^{*\alpha}\| - \mu r_n^{*\alpha} \right) \left( \frac{1}{\|\mathbf{r}_t^{*\alpha}\|} \mathbf{t}_2 - \frac{r_{t_2}^{*\alpha}}{\|\mathbf{r}_t^{*\alpha}\|^3} \mathbf{r}_t^{*\alpha} \right) \right] \quad (52)$$

$$\frac{\partial Z_{\alpha}}{\partial r_{n}^{\alpha}} = \frac{\mu}{1 + \mu^2} \left( \mu \mathbf{n} - \frac{\mathbf{r}_{t}^{*\alpha}}{\|\mathbf{r}_{t}^{*\alpha}\|} \right) \quad (53)$$

$$\frac{\partial Z_{\alpha}}{\partial r_{t_1}^{\alpha}} = \frac{r_{t_1}^{*\alpha}}{(1 + \mu^2) \|\mathbf{r}_t^{*\alpha}\|} \left( \frac{\mathbf{r}_t^{*\alpha}}{\|\mathbf{r}_t^{*\alpha}\|} - \mu \mathbf{n} \right) + \frac{\|\mathbf{r}_t^{*\alpha}\| - \mu r_n^{*\alpha}}{1 + \mu^2} \left( \frac{1}{\|\mathbf{r}_t^{*\alpha}\|} \mathbf{t}_1 - \frac{r_{t_1}^{*\alpha}}{\|\mathbf{r}_t^{*\alpha}\|^3} \mathbf{r}_t^{*\alpha} \right) \quad (54)$$

$$\frac{\partial Z_{\alpha}}{\partial r_{t_2}^{\alpha}} = \frac{r_{t_2}^{*\alpha}}{(1 + \mu^2) \|\mathbf{r}_t^{*\alpha}\|} \left( \frac{\mathbf{r}_t^{*\alpha}}{\|\mathbf{r}_t^{*\alpha}\|} - \mu \mathbf{n} \right) + \frac{\|\mathbf{r}_t^{*\alpha}\| - \mu r_n^{*\alpha}}{1 + \mu^2} \left( \frac{1}{\|\mathbf{r}_t^{*\alpha}\|} \mathbf{t}_2 - \frac{r_{t_2}^{*\alpha}}{\|\mathbf{r}_t^{*\alpha}\|^3} \mathbf{r}_t^{*\alpha} \right) \quad (55)$$

### Case B: Contact with sticking

This case is true if and only if

$$\mu \|\mathbf{r}_t^{*\alpha}\| \geq r_n^{*\alpha} \quad \text{and} \quad \|\mathbf{r}_t^{*\alpha}\| < \mu r_n^{*\alpha} \quad (56)$$

In view of Equations (28) and (35), we have

$$Z_{\alpha} = \rho \mathbf{x}^{*\alpha} \quad (57)$$

and then

$$A_{\alpha} = \mathbf{0}_{3 \times 3} \quad (58)$$

$$\frac{\partial Z_{\alpha}}{\partial x_{n}^{\alpha}} = \rho \mathbf{n} \quad (59)$$

$$\frac{\partial Z_{\alpha}}{\partial x_{t_1}^{\alpha}} = \rho \mathbf{t}_1 + \rho \mu \frac{x_{t_1}^{\alpha}}{\|\mathbf{x}_t^{\alpha}\|} \mathbf{n} \quad (60)$$
$$ \frac{\partial \mathbf{Z}_{\alpha}}{\partial x_{t_2}^{\alpha}} = \rho \mathbf{t}_2 + \rho \mu \frac{x_{t_2}^{\alpha}}{\|\mathbf{x}_{t_1}^{\alpha}\|} \mathbf{n} \qquad (61) $$

In the particular case of no contact, it is not necessary to construct the tangent matrix and the solution is directly given by

$$ \boldsymbol{x}_{\alpha}^{k+1} = \begin{Bmatrix} \mathbf{0} \\ \mathbf{x}^{\alpha\beta(k)} \end{Bmatrix} \qquad (62) $$

### Remark 1

In the Newton iterative process, the solution can be obtained by a classic condensation technique as follows:

$$ \Delta \boldsymbol{x}_{\alpha}^{i+1} = \begin{cases} \Delta \mathbf{r}^{\alpha} = [\mathbf{A}_{\alpha} + \mathbf{B}_{\alpha} \mathbf{W}_{\alpha\alpha}]^{-1} (-\mathbf{Z}_{\alpha}^{(i)} + \mathbf{B}_{\alpha} \mathbf{x}^{\alpha(i)}) \\ \Delta \mathbf{x}^{\alpha} = -\mathbf{x}^{\alpha(i)} + \mathbf{W}_{\alpha\alpha} \Delta \mathbf{r}^{\alpha} \end{cases} \qquad (63) $$

### Remark 2

Compared to the Uzawa algorithm, in the Newton algorithm, the solution is controlled by both local (iteration i) and global convergence criteria (iteration k).

## 4. A TEST EXAMPLE

Many application examples, academic or industrial, have been carried out using the bi-potential method. The constitutive law of deformable bodies can be linear or non-linear with large deformations and large displacements. For the sake of simplicity and clarity and in order to focus our attention on the comparative study between Uzawa and Newton algorithms, only a test example with a linear elastic law and in the presence of large slip contact is considered in this study.

The example concerns the contact between a three-dimensional elastic block ABCDEFGH and a rigid surface Γ (Figure 3). The upper surface ABCD is given a rigid motion described by (a, b, θ) where a and b are, respectively, perpendicular and parallel to the rigid surface. The loading program is designed to apply first a vertical displacement following a and then a horizontal displacement following b. The lower surface EFGH comes into contact with the rigid surface whose normal vector is (0, 0, 1). Each side of the block has a length of 1 mm. The other characteristics of this example are as follows:

* Young's modulus: $E = 210\,000 \text{ N/mm}^2$;

* the Poisson ratio: $\nu = 0.3$;

* coefficient of friction: $\mu = 0.3$;

* global convergence tolerance: $\epsilon_g = 10^{-8}$;

* local convergence tolerance (only used in the Newton algorithm): $\epsilon_l = 10^{-5}$;

* boundary conditions: $\|\mathbf{a}\| = 0.1 \text{ mm}, \|\mathbf{b}\| = 0.4 \text{ mm}, \theta = 60$.

The block is subdivided into eight eight-node brick-like elements as shown in Figure 3. Each element has 27 integration points. Fifty load steps are performed for this problem, so a horizontal or vertical displacement of 0.01 mm is applied to ABCD each step. Figure 4 shows the variation of normal contact forces of points F and H. The horizontal lines correspond to cases where the
<figure><img src="image_4.png" /><figcaption>Figure 3. Contact problem under displacement control.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. Variation of normal contact forces of points F and H.</figcaption></figure>

block slides in a stationary regime. Table I gives detailed numerical results of point F, concerning displacements and contact forces for selected load steps. From these results, we can easily verify the contact conditions such as Coulomb friction laws and Signorini conditions. The analysis was performed by means of both Uzawa and Newton algorithms developed above. The results are almost the same as shown in Figure 4 and Table I. However, the number of iterations in the global contact solution is different, as expected. Figure 5 indicates the evolution of iterations with respect to cumulative iterations in the Newton-Raphson procedure. It shows that, globally, the Uzawa algorithm needs more iterations than the Newton algorithm.
<table><caption>Table I. Contact forces and displacements of point F.</caption><thead><tr><th>Step</th><th>r<sub>t1</sub> (N)</th><th>r<sub>t2</sub> (N)</th><th>r<sub>n</sub> (N)</th><th>U<sub>x</sub> (mm)</th><th>U<sub>y</sub> (mm)</th><th>U<sub>z</sub> (mm)</th></tr></thead><tbody><tr><td>1</td><td>-32.5956</td><td>-32.5956</td><td>153.657</td><td>0.000739176</td><td>0.000739176</td><td>-7.15743e-011</td></tr><tr><td>3</td><td>-94.6401</td><td>-94.6401</td><td>446.138</td><td>0.00221233</td><td>0.00221233</td><td>-2.38517e-010</td></tr><tr><td>5</td><td>-152.573</td><td>-152.573</td><td>719.235</td><td>0.00368016</td><td>0.00368016</td><td>-4.05892e-010</td></tr><tr><td>7</td><td>-206.493</td><td>-206.493</td><td>973.417</td><td>0.00514484</td><td>0.00514484</td><td>-6.35787e-010</td></tr><tr><td>10</td><td>-280.066</td><td>-280.066</td><td>1320.24</td><td>0.00734143</td><td>0.00734143</td><td>-1.11922e-009</td></tr><tr><td>11</td><td>-296.295</td><td>-302.973</td><td>1412.58</td><td>0.00793292</td><td>0.00794625</td><td>3.1167e-010</td></tr><tr><td>12</td><td>-312.422</td><td>-325.207</td><td>1503.21</td><td>0.00847988</td><td>0.0085156</td><td>3.90327e-010</td></tr><tr><td>14</td><td>-344.324</td><td>-367.835</td><td>1679.49</td><td>0.00944491</td><td>0.00954032</td><td>5.39157e-010</td></tr><tr><td>16</td><td>-375.738</td><td>-408.179</td><td>1849.29</td><td>0.0102413</td><td>0.0104021</td><td>6.56127e-010</td></tr><tr><td>18</td><td>-406.839</td><td>-446.38</td><td>2013.21</td><td>0.0108732</td><td>0.011094</td><td>-1.32709e-009</td></tr><tr><td>20</td><td>-438.23</td><td>-483.163</td><td>2174.33</td><td>0.0113587</td><td>0.0116285</td><td>-1.74936e-009</td></tr><tr><td>22</td><td>-447.971</td><td>-524.895</td><td>2300.22</td><td>0.0140293</td><td>0.0147479</td><td>9.97336e-012</td></tr><tr><td>25</td><td>-417.835</td><td>-570.219</td><td>2356.4</td><td>0.0267099</td><td>0.0312728</td><td>1.30087e-010</td></tr><tr><td>30</td><td>-375.928</td><td>-601.664</td><td>2364.84</td><td>0.0524647</td><td>0.0705171</td><td>2.01243e-010</td></tr><tr><td>35</td><td>-362.628</td><td>-609.521</td><td>2364.12</td><td>0.0780964</td><td>0.112958</td><td>2.54078e-010</td></tr><tr><td>40</td><td>-358.261</td><td>-611.94</td><td>2363.66</td><td>0.103402</td><td>0.155961</td><td>2.83243e-010</td></tr><tr><td>50</td><td>-355.656</td><td>-613.308</td><td>2363.23</td><td>0.153658</td><td>0.24235</td><td>3.23905e-010</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Figure 5. Evolution of contact iterations.</figcaption></figure>

## 5. CONCLUSIONS

In this work, three-dimensional contact problems with friction have been theoretically investigated and numerically implemented. Within the bi-potential framework, a new Newton algorithm has been proposed and closed-form tangent matrices have been carefully derived for different contact statuses. A comparative study has been made between the newly proposed Newton algorithm and the previously developed Uzawa algorithm. The characteristics of each algorithm are discussed. The
numerical test indicates that both algorithms give good results. Signorini conditions and Coulomb friction laws are quite well satisfied. For example, the values of penetration of contact points are in the order of $10^{-10}$ as shown in Table I. The numerical test also shows that the Uzawa algorithm needs more iterations than the Newton algorithm. But the Newton algorithm needs an additional local iterative procedure and, at each iteration, a small system of equations should be solved.

The algorithms presented in this paper can be readily extended to dynamic contact problems including more complex frictional models such as orthotropic friction laws with non-associated flow rules [21].

REFERENCES

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

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

3. Kikuchi N, Oden JT. *Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Elements*. SIAM: Philadelphia, 1988.

4. Zhong ZH. *Finite Element Procedures in Contact-impact Problems*. Oxford University Press: Oxford, 1993.

5. Wriggers P. *Computational Contact Mechanics*. Wiley: New York, 2002.

6. Laursen TA. *Computational Contact and Impact Mechanics: Fundamentals of Modeling Interfacial Phenomena in Nonlinear Finite Element Analysis*. Springer: Berlin, 2002.

7. Chaudhary AB, Bathe KJ. A solution method for static and dynamic analysis of three dimensional contact problems with friction. *Computers and Structures* 1986; **24**:855–873.

8. Parisch HA. Consistent tangent stiffness matrix for three-dimensional non-linear contact analysis. *International Journal for Numerical Methods in Engineering* 1989; **28**:1803–1812.

9. Alart P, Curnier A. A mixed formulation for frictional contact problems prone to newton like solution methods. *Computer Methods in Applied Mechanics and Engineering* 1991; **92**:353–375.

10. Simo JC, Laursen TA. An augmented Lagrangian treatment of contact problems involving friction. *Computers and Structures* 1992; **42**:97–116.

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

12. Heegaard JH, Curnier A. An augmented Lagrangian method for discrete large slip contact problems. *International Journal for Numerical Methods in Engineering* 1993; **36**:569–593.

13. Dostál Z, Horák D, Kučera R, Vondrák V, Haslinger J, Dobiáš J, Pták S. FETI based algorithms for contact problems: scalability, large displacements and 3D Coulomb friction. *Computer Methods in Applied Mechanics and Engineering* 2005; **194**:395–409.

14. De Saxcé G, Feng Z-Q. The bi-potential method: a constructive approach to design the complete contact law with friction and improved numerical algorithms. *Mathematical and Computer Modeling* 1998; **28**(4–8):225–245, special issue: Recent Advances in Contact Mechanics.

15. Feng Z-Q. 2D or 3D frictional contact algorithms and applications in a large deformation context. *Communications in Numerical Methods in Engineering* 1995; **11**:409–416.

16. Feng Z-Q, Peyraut F, Labed N. Solution of large deformation contact problems with friction between Blatz–Ko hyperelastic bodies. *International Journal of Engineering Science* 2003; **41**:2213–2225.

17. Simo JC, Hughes TJR. *Computational Inelasticity*. Springer: New York, 1998.

18. Papadopoulos P, Taylor RL. A mixed formulation for the finite element solution of contact problems. *Computer Methods in Applied Mechanics and Engineering* 1992; **94**:373–389.

19. Klarbring A. Mathematical programming and augmented Lagrangian methods for frictional contact problems. In *Contact Mechanics International Symposium*, Curnier A (ed.). Presses Polytechniques et Universtaires Romandes: Lausanne, Switzerland, 1992; 409–422.

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

21. Hjiaj M, Feng Z-Q, De Saxcé G, Mróz Z. Three-dimensional finite element computations for frictional contact problems with non-associated sliding rule. *International Journal for Numerical Methods in Engineering* 2004; **60**:2045–2076.