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

HAL
open science

# Various numerical methods for solving unilateral contact problems with friction

Patrick Chabrand, Frédéric Dubois, Michel Raous

► To cite this version:

Patrick Chabrand, Frédéric Dubois, Michel Raous. Various numerical methods for solving unilateral contact problems with friction. Mathematical and Computer Modelling, 1998, 28 (4-8), pp.97 - 108. 10.1016/S0895-7177(98)00111-3. hal-01393691

HAL Id: hal-01393691

https://hal.science/hal-01393691

Submitted on 8 Nov 2016

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Various Numerical Methods for Solving Unilateral Contact Problems with Friction

P. CHABRAND, F. DUBOIS AND M. RAOUS

Laboratoire de Mécanique et d'Acoustique-CNRS
31, chemin Joseph Aiguier, 13402 Marseille Cedex 20, France

**Abstract**—A number of numerical methods developed in the group “Mechanics and Modelling of Contact” for solving frictional problems, are presented. These methods are based on different formulations. The presentation is done in the context of finite elastoplastic strains. The efficiency of the methods is compared on two tests and benchmarks related to metal forming processes. © 1998 Elsevier Science Ltd. All rights reserved.

**Keywords**—Contact, Friction, Finite elastoplasticity, Numerical methods.

## 1. INTRODUCTION

Various numerical methods for solving unilateral contact problems with friction have been developed in the group “Mechanic and Modelling of Contact” during the last few years. They are based on different formulations: quasi-variational inequalities, complementarity problems, and mixed formulations. In this paper, four numerical methods together with the appropriate formulations are presented: minimization under constraints (Gauss-Seidel with acceleration procedures), mathematical programming (Lemke method), and Lagrangian methods (augmented Lagrangian, penalization). Difficulties of implementation as well as efficiency are discussed and compared. These methods were first developed for elasticity problems under small deformation hypothesis in the finite element codes Protis and Gyptis (see [1-4]). In this paper, their implementation for finite elastoplastic strain and large displacement is presented. Thus, the efficiency of the methods is tested on the coupling of the kinematic and constitutive nonlinearities. Implementation is done in finite element code Simem3 with Renault [5].

Section 2 is devoted to the formulation of the finite elastoplasticity problem and Section 3 presents the associated numerical treatment. The various formulations and numerical methods for the unilateral contact with friction are given in Section 4. These methods are compared on two examples in Section 5.

## 2. FINITE STRAIN ELASTOPLASTICITY

We characterize the behaviour of the material using large strains and displacement elastoplastic formulation, assuming isotropic hardening [6,7]. The model adopted is based on two basic hypotheses:

* a multiplicative decomposition of the deformation gradient $F$ into its elastic $F^e$ and plastic $F^p$ parts, as proposed by Lee [8]:
$$F = F^e \cdot F^p, \qquad (1)$$

where $F^e$ defines the local, stress-free unloaded configuration;

* Clausius-Duhem inequality:

$$D_0 = \tau : d - \bar{\Psi}_0 \geq 0, \qquad (2)$$

where $D_0$ denotes the dissipation, $\tau$ the Kirchhoff stress tensor, $d$ the rate of deformation, and $\bar{\Psi}_0$ the free energy density.

The behaviour of the material is determined by the free decoupled energy density function

$$\Psi_0(b^e, \bar{\epsilon}^p) = \Psi_0^e(b^e) + \Psi_0^p(\bar{\epsilon}^p), \qquad (3)$$

where its elastic part $\Psi_0^e$ is an isotropic function of the left Cauchy-Green strain tensor $b^e$, while its plastic part $\Psi_0^p$ depends on the equivalent plastic strain $\bar{\epsilon}^p$.

The behaviour between the intermediate and current configurations is characterized by a hyperelastic response. The Kirchhoff strain tensor is then given by the constitutive equation [9]

$$\tau = 2b^e \frac{\partial \Psi_0^e}{\partial b^e}. \qquad (4)$$

The von-Mises yield function is defined as

$$f(\tau, \bar{\epsilon}^p) = \sqrt{\frac{3}{2} \operatorname{dev}(\tau) : \operatorname{dev}(\tau) - A_0(\bar{\epsilon}^p)} \le 0, \quad \text{with } A_0 = \frac{\partial \Psi_0^p}{\partial \bar{\epsilon}^p}, \qquad (5)$$

where $\operatorname{dev}(\tau)$ denotes the deviatoric part of the Kirchhoff stress.

The intermediate configuration is updated by integrating plastic evolutionary laws [10]:

$$\begin{align}
d^p &= \lambda \frac{\partial f}{\partial \tau}, & \lambda \ge 0, & f(\tau, \bar{\epsilon}^p) \le 0, & \lambda f(\tau, \bar{\epsilon}^p) = 0, \nonumber \\
\bar{\epsilon}^p &= -\lambda \frac{\partial f}{\partial A_0}, & & & \nonumber \\
d^p &= F^{e-T} (\dot{F}^p F^{p-1})^s \mathbb{F}^e T = -\frac{1}{2} b^{e-1} \dot{b}^e, \tag{6}
\end{align}$$

where $d^p$ is the plastic strain rate, $\lambda$ is the plastic multiplier, and $\dot{b}^c$ denotes the contravariant derivative of $b$.

Last, we assume that the deviatoric and spherical parts of the density of energy function are uncoupled, having the following form [11]:

$$\Psi_0^e(b^e) = U(J^e) + \bar{\Psi}_0^e(\bar{b}^e), \qquad (7)$$

$$= \frac{K}{2} \left[ \frac{1}{2} (J^{e2} - 1) - \ln J^e \right] + \frac{1}{2} \mu (I_{\bar{b}^e} - 3), \qquad (8)$$

where $J^e = \det(F^e) = \det(F), \bar{b}^e = J^{e-1/3}b^e, I_{\bar{b}^e} = \text{trace}(\bar{b}^e)$. K and $\mu$ are material constants.

From the computational point of view, the above constitutive equations are integrated using an elastic predictor/plastic corrector algorithm [6]. Assume that at time $t_n$, the configuration $\Omega_n$ of the body and state variables $(\tau, C_n^{p-1} = F_n^{p-1}F_n^{p-T}, \bar{\epsilon}_n^p)$ are known. Let $\Delta u_{n+1}$ be the (known) incremental displacement between the configurations $\Omega_n$ and $\Omega_{n+1}$. The total deformation gradient at time $t_{n+1}$ can be computed.

We then make an elastic prediction, assuming that the plastic variables remain frozen ($F_{n+1}^p = F_n^p$, i.e., $C_{n+1}^{p-1} = C_n^{p-1}, \bar{\epsilon}_{n+1}^p = \bar{\epsilon}_n^p$). Using the notion of operator splitting, the trial elastic part of the deformation gradient, as well as the elastic trial Kirchhoff stress can then be computed:

$$\tilde{b}_{n+1}^{e \mathrm{Tr}} = J^{-2/3} \mathbb{F} C_n^{p-1} \mathbb{F}^T, \qquad (9)$$

$$\tau_{n+1}^{\mathrm{Tr}} = JU'(J)\mathbb{1} + \mu \operatorname{dev}(\tilde{b}_{n+1}^{\mathrm{Tr}}). \qquad (10)$$
We check whether the stress lies within the elastic region, i.e.,

$$f(\tau_{n+1}^{\text{Tr}}, \bar{\epsilon}_n^p) \le 0. \quad (11)$$

If not, a plastic correction step is then carried out using a classical radial return algorithm [12]. After computing the plastic multiplier correction ($\Delta\lambda$), the state variables are updated

$$\text{dev}(\tau) = \text{dev}(\tau^{\text{Tr}}) - \frac{2}{3}\mu \text{trace} (\tilde{b}^e) \sqrt{\frac{3}{2}} \Delta\lambda \frac{\text{dev}(\tau^{\text{Tr}})}{\|\text{dev}(\tau^{\text{Tr}})\|}, \quad (12)$$

$$C_{n+1}^{p-1} = J^{2/3} \mathbb{F}^{-1} \left( \frac{\text{dev}(\tau)}{m} u + \frac{1}{3} I_{\tilde{b}^e} \right) \mathbb{F}^{-T}, \quad (13)$$

$$\bar{\epsilon}_{n+1}^p = \bar{\epsilon}_n^p + \Delta\lambda. \quad (14)$$

To improve the treatment of the plastic incompressibility (det $F^p = 1$), a special procedure is used based on a three-field formulation (displacement, pressure, dilatation). It leads to the construction of a class of mixed finite elements. In our numerical applications, we will use a Q4/P0, which is one of these mixed elements (bilinear in displacements and constant with respect to pressure and dilatation).

### 3. FINITE-ELEMENT MODELLING

In the displacement based finite-element method, the discretized form of the equilibrium equations is used to calculate an estimated incremental displacement. A modified Newton-Raphson method is used to deal with the nonlinear equations arising from the constitutive equations and the large deformations. The algorithm is organized as follows: loading is given as a sequence of loading steps. At a given loading step, an iterative process is then performed to solve the set of nonlinear equations associated with successive intermediate configurations. We take $\Omega_n$ as the reference configuration at the loading step $n + 1$. Configuration $\Omega_{n+1}$ is computed iteratively from $\Omega_n$ (which is known from the previous loading step) by iterations $i$ of the Newton-Raphson method with corresponding intermediate configurations $\Omega_{n+1}^i$.

Let $u_n$ be the nodal displacement vector at the end of loading step $n$. We denote $\Delta u_{n+1}^i$, the incremental displacement between configurations $\Omega_n$ and $\Omega_{n+1}^i$, and $du^i$, the displacement between $\Omega_{n+1}^i$ and $\Omega_{n+1}^{i+1}$. Hence,

$$\Delta u_{n+1}^{i+1} = \Delta u_{n+1}^i + du^i, \quad (15)$$

$$u_{n+1}^{i+1} = u_n + \Delta u_{n+1}^{i+1}. \quad (16)$$

Omitting the contact terms, the equilibrium equation of the configuration $\Omega_{n+1}^i$ is

$$\{\text{Res}(u_{n+1}^i)\} = \{F^{\text{int}}(u_{n+1}^i)\} - \{F_{n+1}^{\text{ext}}\} = 0, \quad (17)$$

where $\{\text{Res}(u_{n+1}^i)\}$ denotes the equilibrium residual vector obtained by assembling the internal and external forces. $\{F^{\text{int}}(u_{n+1}^i)\}$ is the discrete load vector corresponding to the internal stresses, $\{F_{n+1}^{\text{ext}}\}$ is the discrete load vector corresponding to the external forces excluding the contact forces. A linearized form of equation (17) can be obtained

$$[K_T]_{n+1}^i \{du^{i+1}\} = - \{F^{\text{int}}(u_{n+1}^i)\} + \{F_{n+1}^{\text{ext}}\}, \quad (18)$$

where $[K_T]_{n+1}^i$ is the tangent consistent stiffness matrix. From now on, we omit the subscript $n + 1$ denoting the dependence on the load increment.
<figure><img src="image_3.png" /><figcaption>Figure 1. Contact definitions.</figcaption></figure>

4. FRICTIONAL CONTACT

We restrict this part of the study to the contact between a deformable body and a rigid obstacle (Figure 1). At each node P on the surface of the deformable body that is close to being in contact, we denote by $P_c$ the projection of P on the rigid obstacle surface [6] along its outward normal $n$, which is unknown a priori. Let $w$ denote the relative position between the points P and $P_c$.

We perform the following decomposition into the normal and tangential components of the displacements ($u$), the relative position ($w$), and the contact stress vector ($r$):

$$u_N = u \cdot n, \quad u_T = u - u_N n, \quad w = w_N n + w_T, \quad r = r_N n + r_T. \quad (19)$$

The unilateral contact conditions can then be written in terms of the relative position ($w$) as follows:

$$w_N \ge 0, \quad r_N \ge 0, \quad w_N r_N = 0. \quad (20)$$

As a contact condition for the tangential direction, we take Coulomb's law of friction, written below in an incremental form in which $\Delta w_T$ denotes the relative tangential displacement

$$\|r_T\| \le \mu r_N \begin{cases} \|r_T\| < \mu r_N \Rightarrow \Delta w_T = 0, \\ \|r_T\| = \mu r_N \Rightarrow \Delta w_T = -\alpha r_T, \quad \alpha \ge 0. \end{cases} \quad (21)$$

After finite element discretization, one can obtain the following mixed discretized system [13]:

$$\begin{align} \{\operatorname{Res}(u^i, R_N^i, R_T^i)\} &= \{F^{\operatorname{int}}(u^i)\} - \{F^{\operatorname{ext}}\} - [H^i]^T \{R^i\} = 0, \nonumber \\ W_N^i T (R_N^* - R_N^i) &\ge 0, \quad \forall R_N^* \in C_N^h, \nonumber \\ \Delta W_T^i T (R_T^* - R_T^i) &\ge 0, \quad \forall R_T^* \in C_T^h (R_N^i), \tag{22} \end{align}$$

with $nc$ the number of contact nodes, $\dim$ the dimension of the problem, $R_N$ the normal reaction vector ($nc$ terms), $R_T$ the tangential reaction vector ($nc \cdot (\dim - 1)$ terms), $\{R\}$ the local contact reaction vector ($nc \cdot \dim$ terms), $W_N$ the normal relative position vector, and $\Delta W_T$ the tangential relative displacement vector. The matrix $H$ can be used to transform the expressions from the global frame to the local one [14]. Furthermore, we have introduced the subsets

$$\begin{align} C_N^h &= \{P_N \in \mathbb{R}^{nc}; P_{N_j} \ge 0; \forall j = 1, \dots, n_c\}, \nonumber \\ C_T^h(R_{N_j}) &= \{P_{T_j} \in \mathbb{R}^{\dim - 1}; \|P_{T_j}\| \le \mu R_{N_j}\}, \nonumber \\ C_T^h(R_n) &= \{P_T \in \mathbb{R}^{(\dim - 1)*nc}; P_{T_j} \in C_T^h(R_{N_j}); \forall j = 1, \dots, nc\}. \tag{23} \end{align}$$

It is important to note that all these subsets depend implicitly on the solution (the contact nodes, local frame, etc., are unknown a priori). In the formulation proposed here, we assume that
the vector of the contact forces {$R^i$} is an unknown, like the displacements. This formulation is well suited for using minimization methods under constraints of quadratic functional (Gauss-Seidel methods with projection), and Lemke's method, which are presented below. On the other hand, if the numerical methods used are the penalty/Lagrangian methods, the regularization will eliminate the unknown vector {$R$}.

Using an extension of Newton's method to these generalized equations, as proposed in [13], leads to a linearized system

$$[\tilde{K}_T]^i \{du^{i+1}\} = - \{F^{\text{int}}(u^i)\} + \{F^{\text{ext}}\} + [H^i]^T \{R^{i+1}\}, \quad (24)$$

with

$$(W_N + H_N(u^i) du^{i+1})^T (R_N^* - R_N^{i+1}) \ge 0, \quad \forall R_N^* \in C_N^h, \\ (\Delta W_T + H_T(u^i) du^{i+1})^T (R_T^* - R_T^{i+1}) \ge 0, \quad \forall R_T^* \in C_T^h (R_N^{i+1}), \quad (25)$$

$$[\tilde{K}_T]^i = [K_T]^i + [R^i \nabla H^i T]. \quad (26)$$

One can see that the contact conditions are expressed in a local frame which is known at the beginning of each equilibrium equation $i + 1$, and that the rigid obstacle surface is locally linearized by its tangent plane (Figure 1).

From the numerical point of view, we have made the following hypotheses.

* The second term in the expression (26) of $[\tilde{K}_T]^i$ is neglected. This term gives the variation of the local frame.

* In order to decouple the nonlinearities due to the volumic and surface behaviour, we add a contact loop to our resolution (Table 1). The basic idea is that we assume the local frame determined at the beginning of a contact step to be “frozen” during the equilibrium loop. After convergence, if necessary, the local frame is updated and a new equilibrium loop is performed. This strategy, along with the prediction of the displacement at the beginning of the equilibrium loop [6], is well suited to our problems. Furthermore, it allows us to express the unilateral contact conditions in terms of incremental values. Denoting $\Delta w_N$, the new expression for the normal relative position, we obtain the following relations:

$$\Delta w_N = \Delta u_{n+1} \cdot n - \Delta g_N \ge 0, \quad r_N \ge 0, \quad \Delta w_N r_N = 0, \quad (27)$$

where $\Delta g_N$ is the gap function computed at each contact step.

<table><caption>Table 1. Solution strategy.</caption><tr><td>Loading loop (n), do while n &le; max_n</td></tr><tr><td>Initialization of contact conditions</td></tr><tr><td>Contact loop (ic), do while ic &le; max_ic</td></tr><tr><td>Equilibrium loop (i), do while i &le; max_i</td></tr><tr><td>Updating &sigma;<sup>i</sup><sub>n</sub>, K<sup>i</sup><sub>n</sub>, F<sup>i</sup><sub>n</sub></td></tr><tr><td>If equilibrium convergence go to (jj)</td></tr><tr><td>Solve K<sup>i</sup><sub>n</sub> du<sup>i+1</sup> = F<sup>i</sup><sub>n</sub> + H<sup>T</sup>R<sup>i+1</sup><sub>n</sub></td></tr><tr><td>&Delta;u<sup>i+1</sup><sub>n</sub> = &Delta;u<sup>i+1</sup><sub>n</sub> + du<sup>i+1</sup>, i = i + 1</td></tr><tr><td>End do while i &le; max_i</td></tr><tr><td>(jj) If contact convergence go to (j)</td></tr><tr><td>Contact updating</td></tr><tr><td>ic = ic + 1, i = 0</td></tr><tr><td>End do while ic &le; max_ic</td></tr><tr><td>(j) If program ends, STOP</td></tr><tr><td>n = n + 1, i = 0, ic = 0</td></tr><tr><td>End do while n &le; max_n</td></tr></table>
The numerical methods available for dealing with the discretized frictional contact problem set in equations (24) and (25) are the penalty/Lagrangian methods, the minimization under constraints of a quadratic growth functional (Gauss-Seidel methods with projection) methods, and Lemke's method, which can be used to solve a linear complementarity problem associated with the contact and friction conditions.

## 4.1. Augmented Lagrangian Formulation

The method used in this study is that described by Simo and Laursen [15]. Formulating Coulomb's friction law by analogy with the theory of plasticity leads to the following definition of the slip surface:

$$f_s = \|R_T\| - \mu R_N = 0. \quad (28)$$

If $f_s < 0$, no slip occurs (stick contact) and when $f_s = 0$, slip is impending. The condition $f_s > 0$ is not admissible. The use of these notations makes it possible to write the Coulomb friction law as the following Kuhn-Tucker conditions:

$$f_s \le 0, \quad \Delta W_T = -\Delta \xi \frac{\partial f_s}{\partial R_T}, \quad \Delta \xi \ge 0, \quad \Delta \xi \cdot f_s = 0. \quad (29)$$

Additively decomposing $R_N$ and $R_T$ into their penalty and Lagrange multiplier parts yields:

$$R_N = (\lambda_N - \epsilon_N W_N), \quad (30)$$

$$\Delta R_T = \Delta \lambda_T - \epsilon_T \left( \Delta w_T + \Delta \xi \frac{\partial f_s}{\partial R_T} \right). \quad (31)$$

Introducing the unknown Lagrange terms into the problem leads to an additive loop (indicated below by superscript $k$) in which the Lagrange multipliers are updated. At each of these iterations, the Newton-Raphson method (iteration indicated below by superscript $i$) is used to determine the equilibrium state. The normal component ($R_N$) of the contact forces is computed using

$$^k R_{N+1}^i = \langle {}^k \lambda_{N+1} - \epsilon_N \langle {}^k \Delta W_{N+1}^i \rangle \rangle. \quad (32)$$

When a backward Euler scheme associated with a radial return algorithm is used to determine the tangential component ($R_T$) of the contact forces, a trial value of ($R_T$) is first calculated assuming a sticking state:

$$^k R_{T+1}^{i(\text{trial})} = R_{T_n} + ^k \Delta \lambda_T - \epsilon_T \langle {}^k \Delta W_{T+1}^i \rangle. \quad (33)$$

The corrected value of ($R_T$) is then evaluated using a radial return algorithm, which is a straightforward extension of those used in plasticity (see [16,17])

$$^k R_{T+1}^i = {}^k R_{T+1}^{i(\text{trial})} - \Delta \xi \frac{^k R_{T+1}^{i(\text{trial})}}{\| {}^k R_{T+1}^{i(\text{trial})} \|}, \quad (34)$$

where $\Delta \xi = 0$, if $f_{s+1}^{\text{trial}} \le 0$ and $\Delta \xi = (f_{s+1}^{\text{trial}})/\epsilon_T$, if $f_{s+1}^{\text{trial}} > 0$.

After convergence of the Newton-Raphson iterations, the updated forms of the Lagrange multipliers are then given by the following expressions:

$${}^{k+1}\lambda_{N+1} = \langle {}^k \lambda_{N+1} - \epsilon_N \langle {}^k \Delta W_{N+1} \rangle \rangle, \quad (35)$$

$${}^{k+1}\Delta\lambda_{T+1} = {}^k\Delta\lambda_{T+1} - \epsilon_T \left( \Delta W_T^k + \Delta\xi \frac{^k R_{T+1}^{(\text{trial})}}{\|^k R_{T+1}^{(\text{trial})}\|} \right). \quad (36)$$

With this method, it is possible to check exactly the unilateral contact and the Coulomb friction law without using a large penalty coefficient. However, a good compromise between the value of the penalty coefficients $\epsilon_N$ and $\epsilon_T$, the value of the Lagrange multipliers $\lambda_N$ and $\lambda_T$ and the number of augmentations has to be found. Precautions have to be taken in choosing these parameters, otherwise it can be a very CPU time consuming method.
## 4.2. Minimization with Projection Method

Here we use a straightforward extension of the method presented by Raous et al. [1] and Lebon-Raous [2] for dealing with unilateral contact and friction in elasticity with infinitesimal deformations. Let us recall first the contactless case ($R^{i+1}$ vanishes in (24)). At each equilibrium iteration (as $K_T$ is a symmetric positive definite matrix), it can easily be proved that the problem can be set as the following minimization one, where for the sake of simplicity, the index $n$ is not shown.

$$ \text{Find the vector } du^{i+1} \text{ in } R^{2\text{no}} \text{ (no being the number of nodes),} \\ \text{such that } J(du^{i+1}) \le J(v), \forall v \in R^{2\text{no}} \text{ with} \\ J(v) = \frac{1}{2}v^t K_T v - v^t F, \text{ with } F = -\{F^{\text{int}}(u_{n+1})^i\} + \{F^{\text{ext}}_{n+1}\}. \qquad (37) $$

As shown in [18], the frictionless contact problem can be written as a constrained minimization problem, in which the constraints are the discretized form of the unilateral conditions and where the unknown vector $du^{i+1}$ is sought in a convex set $K$ defined by

$$ K = \{v \in R^{2\text{no}}, \text{ such that } \Delta W_{N_j}^{i+1} = \Delta W_{N_j}^i + v_{N_j} \ge 0, \forall j \in I\}, $$

where $I$ is the set of contact nodes. This problem is then solved using a Gauss-Seidel with projection method. This technique has been generalized to the case involving friction [1,2]. A fixed-point algorithm on the sliding threshold leads to solving a sequence of constrained minimization problems where the constraints affect only the normal components of the contact nodes, as in frictionless contact.

These problems include an additive nondifferentiable term for the tangential components of the contact nodes corresponding to friction, which is updated at each fixed-point iteration as shown below.

$$ \text{Find a fixed point } b \text{ of the application } s^{i+1} \to \mu R_N (du^{i+1}(s)), \\ \text{with } du^{i+1}, \text{ such that } Q(du^{i+1}) \le Q(v), \forall v \in K, \qquad (38) \\ Q(v) = \frac{1}{2} v^t K_T v - v^t (F + \epsilon_T b^{i+1}), $$

where $w_N$ and $w_T$ are the normal and tangential components of the vector $w$ defined in (20) and (21).

$\epsilon_T$ concerns the friction and is defined as follows:

$$ \begin{align*} \epsilon_T &= 1, & \text{if } \Delta W_T(v) < 0, \\ \epsilon_T &= -1, & \text{if } \Delta W_T(v) > 0, \\ \epsilon_T &= 0, & \text{otherwise.} \end{align*} $$

This is accelerated by using an Aitken procedure [4] which does not need the determination of an optimal parameter as over-relaxation does. This method is very robust and works even on ill-conditioned problems. Nevertheless, in some cases, it would require many iterations.

## 4.3. Linear Complementarity Problem

In dealing with the three-dimensional case, Klarbring and Bjorkman [13] have introduced a piecewise linear friction law, approximating Coulomb's friction law. This discretization procedure then makes it possible to write the friction relations as complementarity conditions and then to set the problem as a linear complementarity one. In the present study, we shall restrict ourselves to a two-dimensional analysis. In this situation, their approach leads to introducing two new variables,
$\lambda$ and $\phi$, which define the boundary of the Coulomb's cone. The Kuhn-Tucker conditions for the frictional Coulomb problem can then be written in the following form:

$$
\begin{align}
R_T \in C(R_N) &= \{P_T, \phi_m(P_T, R_N) \ge 0, m = 1, 2\}, \nonumber \\
\phi_1(P_T, R_N) &= -P_T + \mu R_N, \nonumber \\
\phi_2(P_T, R_N) &= P_T + \mu R_N, \nonumber \\
\Delta W_T &= -\sum_{m=1}^{2} \lambda_m \frac{\partial \phi_m}{\partial R_T}, \tag{39} \\
\lambda_m \ge 0, \quad \phi_m \ge 0, \quad \lambda_m \phi_m = 0, \quad \text{for } m = 1, 2. \nonumber
\end{align}
$$

Using a condensation procedure, two connected systems can be written. The first one deals only with the normal and tangential components of the contact nodes which, upon introducing variables $\lambda$ and $\phi$, are all constrained by complementarity conditions. This linear complementarity problem with a $2*nc$ by $2*nc$ square singular matrix (nc being the number of contact nodes) can then be straightforwardly solved using a pivot algorithm such as Lemke's method. The second system deals with the nodes which are not involved in the contact. This is a nonconstrained problem in which the only unknown vector is the displacement one. Its solution obviously depends on the solution of the previous system, and can be obtained using a more classical algorithm.

## Lemke's Method Principle

To describe Lemke's method, let us consider the following linear complementary problem.

PROBLEM P1. Let $q$ be a given vector $\in \mathbb{R}^n$, find $W$ and $Z \in \mathbb{R}^n$ such that

$$
\begin{gather}
\bar{I}W - \bar{M}Z = q, \tag{40} \\
W_i \ge 0, \quad Z_i \ge 0, \quad W_i \cdot Z_i = 0, \quad i = 1, \dots, n.
\end{gather}
$$

Here $\bar{I}$ is the $n \times n$ identity matrix and $\bar{M}$ a regular $n \times n$ matrix. Lemke's method is based on the following remark: if $q_i \ge 0, \forall i = 1, \dots, n$, the solution of Problem P1 is given by

$$ W_i = q_i, \quad Z_i = 0, \quad i = 1, \dots, n. \tag{41} $$

A pivoting algorithm is used to construct by linear combinations, a sequence of nonnegative vectors $q^k$, which will be defined below. If there is an index $j$ such that $q_j \le 0$, the pivoting algorithm is initialized by introducing supplementary variable $z_0 \in \mathbb{R}$

PROBLEM P2. Let $q \in \mathbb{R}^n$ be given, find $W$ and $Z \in \mathbb{R}^n$ and $z_0 \in \mathbb{R}$ such that

$$
\begin{gather}
\bar{I}W - \bar{M}Z - \bar{I}z_0 = q, \tag{42} \\
W_i \ge 0, \quad Z_i \ge 0, \quad W_i \cdot Z_i = 0, \quad z_0 \ge 0, \quad i = 1, \dots, n.
\end{gather}
$$

In Problem P2, $\bar{I}$ is a $n$-vector of ones. $z_0$ is introduced to obtain a nonnegative right-hand side member and to keep it so by pivoting. Let $\bar{A}^0 = [\bar{I}, \bar{M}, \bar{I}]$, the $n*(2n+1)$ matrix, $X^0 \in \mathbb{R}^{2n+1}$ and $q^0 \in \mathbb{R}^n$ be defined by $X^{0T} = (W_1, \dots, W_n, Z_1, \dots, Z_n, z_0)$ and $q^0 = q$.

The components of vector $X^i$ which are associated with $\bar{I}$ in $\bar{A}^i$ are usually called the basic variables. This means that in the initialization step defined above, the basic variables are $(W_1, \dots, W_n)$. The pairs $(W_i, Z_i)$ are called complementary variables.

Using the notation introduced above, Lemke's algorithm can be written as follows.
## Lemke's Algorithm

### 1/INITIALIZATION. Let s be the index such that

$$ q_s = \min_{i=1,...,n} q_i. \quad (43) $$

Matrix $\bar{A}^1$ and right-hand side member $q^1$ are then defined by

$$ \begin{aligned} A_{i,j}^1 &= A_{i,j}^0 - A_{s,j}^0, & j &= 1, \dots, 2n+1 \text{ and } i \neq s, \\ A_{s,j}^1 &= -A_{s,j}^0, & j &= 1, \dots, 2n+1, \end{aligned} \quad (44) $$

$$ q_i^1 = q_i^0 - q_s^0, \quad \text{for } i \neq s \text{ and } q_s^1 = -q_s^0. \quad (45) $$

The system to be solved is then $\bar{A}^1 \cdot X = q^1$ and its right-hand side member $q^1$ is nonnegative. $z_0$ then becomes the $i^{th}$ basic variable ($A_{i,2n+1}^1 = \delta_{2n+1}^i, i = 1, \dots, n$) instead of $W_s$ ($\exists j$ such that $A_{j,s}^1 \neq \delta_j^i$). One solution of the problem is then $W_1 = q_1^1, W_2 = q_2^1, \dots, W_{s-1} = q_{s-1}^1$, $z_0 = q_s^1, W_{s+1} = q_{s+1}^1, \dots, W_n = q_n^1$ and $Z_1 = 0, Z_2 = 0, \dots, Z_{s-1} = 0, Z_s = W_s = 0, Z_{s+1} = 0, \dots, Z_n = 0$.

### 2/STEP $k+1$. From the initializing step, the pivoting is carried out until $z_0$ leaves the basis. The $n$ basic variables will be the $n$ components of $q^k$ and their complementary variables and $z_0$ will be zero.

Let $t$ be the row of complementary variables which left the basis in the previous step and $r$ the row index such that

$$ \frac{q_r^k}{A_{r,t}^k} = \min_{i=1,...,n} \frac{q_i^k}{A_{i,t}^k}, \quad A_{i,t}^k > 0. \quad (46) $$

In the next step, $y_t$ is incorporated in the basis instead of the $r^{th}$ variable by

$$ \begin{aligned} A_{i,j}^{k+1} &= A_{i,j}^k - \frac{A_{i,t}^k}{A_{r,t}^k} A_{r,j}^k, & j &= 1, \dots, 2n+1 \text{ and } i \neq r, \\ A_{r,j}^{k+1} &= \frac{A_{r,j}^k}{A_{r,t}^k}, & j &= 1, \dots, n, \end{aligned} \quad (47) $$

$$ \begin{aligned} q_i^{k+1} &= q_i^k - \frac{A_{i,t}^k}{A_{r,t}^k} q_r^k, & i &\neq r, \\ q_r^{k+1} &= \frac{q_r^k}{A_{r,t}^k}. \end{aligned} \quad (48) $$

Thus $A_{i,s}^{k+1} = \delta_i^t, i = 1, \dots, n$ and $(\exists j \neq t) A_{j,r}^{k+1} \neq \delta_r^j)$. If $r = s, z_0$ leaves the basis and the solution is reached, if not we restart in 2.

### REMARK. If at step $(k)$ for each $i = 1, \dots, n$, $A_{i,s}^{k+1} \le 0$, Lemke's algorithm has failed. One obtains, an infinite number of solutions characterized by, if $B_i, i = 1, \dots, n$ are the basic variables and $H_i, i = 1, \dots, n+1$ the variables outside the basis

$$ \begin{aligned} B_i + \lambda A_{i,s}^{k+1} &= q_i, & \forall \lambda \ge 0, & i = 1, \dots, n. \quad (49) \\ H_i + \lambda A_{i,s}^{k+1} &= 0, \end{aligned} $$

Convergence of the direct Lemke's method is based on the copositivity of matrix $M$ which, because of the friction, is nonsymmetric. If the friction coefficient is too large, this property is not conserved. This corresponds to the well-known loss of uniqueness which occurs in the frictional problem when the friction coefficient is too large.
# 5. NUMERICAL EXAMPLES

In this section, the algorithms presented above are tested and compared on two simple frictional contact problems. The methods have been implemented in the finite-element code developed by Renault and several University partners in the context of deep drawing [7,19]. The first example is that on the indentation of an elastoplastic block by a rigid cylinder. The second example is that of extrusion of an aluminum cylinder. For each case, it has been verified that the mechanical solutions obtained in using the various numerical methods were effectively the same.

## 5.1. Indentation Test

Here an elastoplastic block is indented by a rigid cylinder. The cylinder has a diameter of 60 mm and the block measures 40 mm by 80 mm. The block is discretized using the four-node elastoplastic element (Q4/P0) described in [6] and has the following material properties: $E = 2 \times 10^5$ MPa, $\nu = 0.3$. The hardening law is given by $\sigma_0 = 593 (0.00384 + \bar{\epsilon}_p)^{0.202}$. The coefficient of friction between the die and the workpiece is $\mu = 0.2$. The mesh and the deformed geometry are shown in Figure 2. The arrows correspond to the "contact" reactions.

<figure><img src="image_4.png" /><figcaption>Figure 2. Undeformed and deformed geometry of indentation.</figcaption></figure>

<table><caption>Table 2. Computational times.</caption><thead><tr><td>Method</td><td>Gauss-Seidel</td><td>Lemke</td><td>Augm. Lag. (5)</td><td>Penalty</td></tr></thead><tbody><tr><td>CPU Times</td><td>1 h 17'</td><td>9'</td><td>28'</td><td>11'</td></tr></tbody></table>

## 5.2. Extrusion of an Aluminum Cylinder

We consider here the example of the frictional extrusion of an aluminum cylinder into a rigid, conical die [15]. The aluminum cylinder has a radius of 5.08 cm and an initial length of 25.4 cm. The billet is pushed 17.78 cm into a conical die, with a wall angle of 5°. In this axisymmetric problem, the discretization of the cylinder is performed using 80 Q4/P0 elements. The material properties are as follows: $E = 6.8955 \times 10^4$ MPa, $\nu = 0.32$ and the hardening law $\sigma_0 = 261.2(0.11868 + \bar{\epsilon}_p)$. The friction coefficient is $\mu = 0.1$.

The mesh and the deformed geometry are shown in Figure 3, where the arrows correspond to the contact forces. The computational times (total CPU time) are given in Table 3.

<table><caption>Table 3. Computational times.</caption><thead><tr><td>Method</td><td>Gauss-Seidel</td><td>Lemke</td><td>Augm. Lag.</td><td>Penalty</td></tr></thead><tbody><tr><td>CPU Times</td><td>8'30"</td><td>2'</td><td>6'30"</td><td>1'45"</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>Figure 3. Extrusion of an aluminum cylinder: initial and deformed meshes.</figcaption></figure>

## 6. CONCLUSION

Our aim is not to draw general conclusions based on the comparison of the four numerical methods under consideration in this paper. We will only make several remarks inspired not only by the examples presented here, but also by our experience on the treatment of contact problems.

The efficiency of each method depends on various characteristics (type of geometries, of loading, total number of degrees of freedom, ratio of the total number of contact nodes, constitutive laws, ...) and it cannot be said that one special method is the best one.

The Gauss-Seidel method is easy to implement. The Aitken acceleration is better than the overrelaxation because there is no convergence parameter to determine. Usually, the computational time is not as large as it would be expected, however, in some cases the convergence can be very slow.

The Lemke method is more complicated to implement, it needs the condensation of the problem. It is always the fastest method, but nonconvergence may occur and it will fail for very large numbers of degrees of freedom (as direct methods usually do).

The Penalty method we use here is a particular case of the Augmented Lagrangian method, it consists of only one augmentation and the choice of a large value for $ε_T$. It is of course less accurate than the Augmented Lagrangian method.

In the case of the Augmented Lagrangian method, a compromise between accuracy and computational time is necessary, because the number of augmentations and the values of parameter $ε_T$ can affect greatly the computational time.

## REFERENCES

1. M. Raous, P. Chabrand and F. Lebon, Numerical methods for frictional contact problems and applications, *App. Th. Mech. Jour., Supp. 1 7*, 111-128 (1988).

2. F. Lebon and M. Raous, Multibody contact problems including friction in structure assembly, *Computers & Structures* **43** (5), 925-934 (1992).

3. M. Raous and S. Barbarin, Conjugate gradient for frictional contact, *Contact Mechanics, Presses Polytechniques et Universitaires Romandes* (1992).

4. M. Raous, H. Barbosa and J.C. Latil, Aitken acceleration of Gauss-Seidel with projection algorithm for unilateral contact problems with friction, *Eng. Comp.* (submitted) (1995).

5. P. Chabrand and F. Dubois, A study on the friction in the deep drawing context, NUMIFORM'95, Balkema, Rotterdam, 395-400 (1995).

6. F. Dubois, Contact, frottement, grandes déformations élastoplastiques. Application à l'emboutissage, Ph.D. Thesis, Université d'Aix-Marseille II. (1994).

7. P. Chabrand, F. Dubois and J.C. Gelin, Modelling drawbeads in sheet metal forming, *Int. J. Mech. Sci.* (to appear).

8. E.H. Lee, Elastic plastic deformation at finite strain, *ASME Trans. J. Appl. Mech.* **36**, 1-6 (1969).

9. F. Sidoroff, Incremental constitutive equation for large strain elasto plasticity, *Int. J. Engng. Sci.* **20** (1), 19-26 (1982).

10. J.C. Simo, Algorithms for static and dynamic multiplicative plasticity that preserve the classical return mapping scheme of the infinitesimal theory, *Comp. Meth. in Appl. Mech. and Eng.* **99**, 61-112 (1992).

11. J.C. Simo and C. Miehe, Associative coupled thermoplasticity at finite strains: Formulation, numerical analysis and implementations, *Comp. Meth. in Appl. Mech. and Eng.* **98**, 41-104 (1992).

12. J.C. Simo and R.L. Tayor, Consistent tangent operators for rate-independent elastoplasticity, *Comp. Meth. in Appl. Mech. and Eng.* **48**, 101-118 (1985).
13. A. Klarbring and G. Bjorkman, Solution of large displacement contact problems with friction using Newton's method for generalized equations, *Int. J. Num. Meth. Eng.* **34**, 249-269 (1992).

14. M. Jean and G. Touzot, Implementation of unilateral contact and dry friction in computer codes dealing with large deformation problems, *App. Th. Mech. Jour., Supp. 1 7*, 145-160 (1988).

15. J.C. Simo and T.A. Laursen, An augmented Lagrangian treatment of contact problems involving friction, *Computers & Structures* **42**, 97-116 (1992).

16. A.E. Giannakopoulos, The return mapping method for the integration of friction constitutive relations, *Computers & Structures* **32**, 157-167 (1989).

17. P. Wriggers, T. Vu Van and E. Stein, Finite element formulation of large deformation impact-contact problems with friction, *Computers & Structures* **37**, 319-331 (1990).

18. R. Glowinski, J.L. Lions and R. Trémolières, *Analyse Numérique des Inéquations Variationnelles*, Dunod, Paris, (1976).

19. B. Tathi, A joint project for the numerical simulation of 3D sheet metal forming processes with quasi-static and dynamic approaches, *NUMIFORM'95*, Balkema, Rotterdam, 779-784 (1995).