# Complementarity methods for multibody friction contact problems in finite deformations

Patrick Chabrand, Olivier Chertier, Frédéric Dubois

► To cite this version:

Patrick Chabrand, Olivier Chertier, Frédéric Dubois. Complementarity methods for multibody friction contact problems in finite deformations. International Journal for Numerical Methods in Engineering, 2001, 10.1002/nme.170. [hal-01297284](https://hal.science/hal-01297284)

HAL Id: hal-01297284

https://hal.science/hal-01297284

Submitted on 4 Apr 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_1.png" /></figure>

Distributed under a Creative Commons Attribution - ShareAlike 4.0 International License
# Complementarity methods for multibody friction contact problems in finite deformations

Patrick Chabrand*,†, Olivier Chertier and Frédéric Dubois

Laboratoire de Mécanique et d'Acoustique—C.N.R.S. 31, Chemin Joseph Aiguier, 13402 Marseille Cedex 20, France

## SUMMARY

This paper deals with the frictional contact occurring between deformable elastoplastic bodies subjected to large displacements and finite deformations. Starting from a standard slave/master formulation we have developed a symmetrical formulation with which the unilateral contact conditions and the friction law are satisfied for each body. From the continuum equations, the discretized frictional contact problem is set as a complementarity problem and solved using Lemke's mathematical programming method. The efficiency of the method is illustrated in the case of several examples.

KEY WORDS: multibody frictional contact; finite elastoplasticity; symmetrical slave-master formulation; complementarity problem; Lemke's method

## 1. INTRODUCTION

Methods for solving frictional contact problems are needed in many industrial processes. Friction plays a crucial role, particularly in metal forming processes such as extrusion, rolling and sheet metal forming. Modelling these processes leads to strongly non-linear problems. Non-linearities arise from the finite displacements, the behaviour of the material and the contact conditions. These non-linearities give rise to problems which are extremely difficult to compute and costly to solve. A great deal of effort is still required to develop efficient algorithms for use in this context.

Numerical methods for treating contact problems with finite deformations have been developed by many authors in the case of contact between a rigid obstacle and a deformable body. The main algorithms used so far to solve the contact inequalities have been direct solutions of the inequalities, the penalty method, and the Lagrangian multiplier and augmented Lagrangian methods.

Contact problems can be directly solved using algorithms which are based on a fixed-point method on the sliding limit combined with a minimization problem under constraints,
corresponding to a Tresca problem. To solve the minimization problem under constraints, over-relaxation procedures [1] and conjugate gradient methods [2], both associated with projection techniques, have been developed. The frictional contact problem is regularized before being solved using the conjugate gradient method. These two iterative methods can be used to deal with many friction contact problems, but they can be costly in terms of the computational time required. With the penalty method or the augmented Lagrangian method [3–8] the inequality constraints are only approximately satisfied when there is an insufficiently large number of augmentations. The penalty method has been widely used, but with this method the unilateral conditions and the friction conditions are again only approximately satisfied. The accuracy depends on the penalty coefficient. With Lagrange multiplier methods, the multipliers, which correspond to the physical contact forces, can be exactly computed.

Friction contact problems have also been formulated in the same way as linear complementarity ones [9]. Mathematical programming methods can then be used. The displacements and the contact forces, which are both unknowns, are exactly determined in order to satisfy the constraints. Moreover, with such methods the tangent stiffness matrix remains symmetric when penalty and Lagrange multipliers methods, introduce unsymmetric contributions. We have developed an extension of the Lemke method which gives quite a good description of the contact state. In addition, the stick/slip status is not required to be studied separately, as it is the case when using a regularized approach.

To deal with multibody contact problems, the notion of slave/master bodies was introduced by Hallquist *et al.* [10] and Wriggers and Simo [11]. The use of a continuum mechanics formulation was originally restricted to frictionless contact problem and was recently extended to the frictional problem by Klarbring [12], Laursen and Simo [13] and Wriggers and Miehe [14].

With this approach, the choice as to which body is the master one and which the slave one has to be made a priori. This choice generally depends on the geometry of the antagonist surfaces and on the material characteristics. However, if the solids are not both meshed with a discretization of the same fineness, the master body will be the one meshed with the coarser grid. The contact is then taken to be the contact between the slave nodes and the master surface, and the numerical methods developed constrain the slave nodes so that they do not penetrate into the master body. With this approach, whether or not penetration of the master nodes into the slave body occurs is not checked and the unilateral conditions can be violated. To prevent this from occurring, it is also possible at the same loading increment to reverse the roles played by the bodies: the slave one becomes the master one and *vice versa* [10]. As noted by Laursen [15] and Zavarise and Wriggers [16], with methods of this kind one loses the physical interpretation of the traction forces to be evaluated at each contact node. Pavanachand and Krishnakumar [17] have furthermore observed that there exist some geometric conditions under which some two-pass algorithms fails to converge.

In the present paper, we present a symmetrical treatment for dealing with frictional contact. Both bodies simultaneously play the role of both the slave and the master body. The contact is then taken to be the contact between a slave node on each of the contact surfaces and a master segment on the antagonist surface. With this formulation, no *a priori* choice between the bodies is required. In addition, a coefficient $\theta$ is introduced to weight the role played by each body and to generalize the notion of two-pass master/slave algorithms. Depending on the value of this coefficient, we can choose either the standard formulation, the symmetrical one or a formulation in which one of the bodies can be more *master* than *slave*. This means that in a general presentation of the problem, it is possible to give one of the bodies a
greater influence on the contact depending on the characteristics of the problem in hand. The
advantage of using a coefficient of this kind is that it introduces into the same algorithm
both the standard slave/master approach and a generalized two-pass algorithm. In addition,
introducing this coefficient into the continuum formulation means that each of the bodies
involved in the contact can be simultaneously master and slave, which preserves the physical
meaning of the numerical contact traction forces.

Another aim of this paper is to describe the various steps which are necessary to transform
the initial discrete problem into a suitable one for being solved using Lemke's method. To be
able to use this solver, it is necessary first to condense the problem to the contact variables
involved in the contact alone and then to write the contact with friction as complementarity
conditions. Another important step consists of introducing into the problem to be solved the
distance between a slave node and a master segment and their tangential velocity. Details are
given here about the matrix construction, both in the standard master/slave formulation and in
the symmetrical one. This paper can be summarized as follows. We begin by presenting the
continuum formulation of the frictional contact problem. The impenetrability conditions are
associated with the definition of the gap between a slave particle and the master's surface. We
introduce a convected basis associated with the parametrization of the contact surfaces and we
write Coulomb's friction law with an objective relative velocity. The symmetrical formulation
of the virtual work of contact is then given. The contact and friction constraints, which
are reduced to the two-dimensional case, are both established and the matrix construction is
described in both cases. The last section deals with the examples. First, we describe some very
simple tests used to illustrate the differences between the two formulations. An extrusion test
is then modelled, and lastly, in the context of the study of friction, we analyse the shearing
of two identical deformable microscopic asperities present at the surface of a sheet metal.

## 2. CONTACT KINEMATICS

In this paragraph, we introduce the variables necessary for modelling the kinematics of the
contact between the deformable bodies. Constitutive equations of frictional contact account
for the impenetrability condition at every point on a body into other bodies or itself. These
conditions also have to be consistent with the fact that the relative sliding between bodies
results from a shear stress. There exist many such constitutive equations, which are often
used in the phenomenological approach to contact. The kinematics of the finite deformations
introduce further difficulties. As a matter of fact, the contact conditions are implicit ones
(the contact surfaces are unknown *a priori*) and the friction equations are associated with
a tangential relative velocity between the bodies. This velocity, which has to satisfy the
objectivity principle is based on objective time derivatives such as convective derivatives and
justifies the use of convected bases. Without loss of generality, we restrict our attention to
the case of two bodies. It is easy to extend the solution to cases involving a larger number
of bodies in contact. We parametrize each surface in its reference and current configurations.
Let $\Omega^{(i)}$ denote body $i$ in a reference configuration $C_0$. In the current configuration $C_t$, the
bodies are denoted $\omega^{(i)}$, $i=1,2$. Classically, let $(\varphi^{(i)}, i=1,2)$ be the map of the associated
deformations between $C_0$ and $C_t$ defined by

$$ \varphi^{(i)}: \Omega^{(i)} \times \mathbb{R} \rightarrow \mathbb{R}^n, \quad \Omega^{(i)} \in \mathbb{R}^n, \quad n=2 \text{ or } 3 \quad (1) $$
<figure><img src="image_2.png" /><figcaption>Figure 1. Contact geometry.</figcaption></figure>

and

$$ \omega^{(i)} = \varphi^{(i)}(\Omega^{(i)}) \qquad (2) $$

The potential contact surfaces [12–14] are denoted as $(\Gamma_c^{(i)}, i=1,2)$ in $C_0$. The current position of the surfaces is given by $\gamma_c^{(i)} = \varphi^{(i)}(\Gamma_c^{(i)})$, $i=1,2$. Surface coordinates are introduced, corresponding to $\Gamma_c^{(i)}$. A sufficiently smooth mapping is defined as

$$ \mathbf{\Phi}^{(i)}: A^{(i)} \in \mathbb{R}^{n-1} \longrightarrow \Gamma_c^{(i)}, \quad i=1,2 \qquad (3) $$

here $A^{(i)}$ is an open subset of $\mathbb{R}^{n-1}$, and one component of $A^{(i)}$ is written $\xi = (\xi^1, \dots, \xi^{n-1})$. Let $\mathbf{X}^{(i)} = \mathbf{\Phi}^{(i)}(\xi)$ denote the position of a particle of $\Gamma_c^{(i)}$ and let $\mathbf{x}^{(i)} = \varphi^{(i)}(\xi) = \varphi^{(i)}(\mathbf{\Phi}^{(i)}(\xi))$ denote the position of a particle of $\gamma_c^{(i)}$.

$\mathbf{F}_t^{(i)}$ is the deformation gradient corresponding to $\varphi_t^{(i)}$.

## 2.1. Distance function

In this section, the aim is to introduce the distance between the contact surfaces of the bodies. Each surface will be simultaneously taken to be both a slave and a master surface. At time $t$, at a given point on one surface (taken to be the slave) $\mathbf{x}^{(s)}$ (Figure 1), let us define $\bar{\mathbf{x}}^{(m)}$ and $\bar{\mathbf{x}}^{(m)}$ the position of the closest point on the other surface in the reference and current configurations, respectively, introducing the curvilinear co-ordinate solution of the problem:

$$ \bar{\xi} = \underset{\xi \in A^{(m)}}{\operatorname{argmin}} \frac{1}{2} \|\mathbf{x}^{(s)} - \varphi^{(m)}(\xi)\|^2 \qquad (4) $$

then we have

$$ \begin{aligned} \bar{\mathbf{x}}^{(m)} &= \varphi^{(m)}(\bar{\xi}) \\ \bar{\mathbf{X}}^{(m)} &= \mathbf{\Phi}^{(m)}(\bar{\xi}) \end{aligned} \qquad (5) $$
For a given slave point $\mathbf{x}^{(s)}$, $\bar{\mathbf{x}}^{(m)} = \Phi^{(m)}(\bar{\xi})$ is determined by writing the minimum of the distance function

$$ (\mathbf{x}^{(s)} - \bar{\mathbf{x}}^{(m)}) \cdot \frac{\partial \Phi^{(m)}(\bar{\xi})}{\partial \bar{\xi}^{\beta}} = 0, \quad \beta = 1, \dots, n-1 \qquad (6) $$

## 2.2. Local basis

It is then possible to define a local basis (Figure 1) on $\Gamma_c^{(m)}$ and $\gamma_c^{(m)}$ via partial derivatives with respect to the curvilinear co-ordinates

$$ \mathbf{G}_{\alpha} = \frac{\partial \mathbf{\Phi}^{(m)}(\bar{\xi})}{\partial \bar{\xi}^{\alpha}} \qquad (7) $$

$$ \mathbf{g}_{\alpha} = \frac{\partial \Phi^{(m)}(\bar{\xi})}{\partial \bar{\xi}^{\alpha}} = \mathbf{F}_{t}^{(m)} \mathbf{G}_{\alpha}(\bar{\xi}) \qquad (8) $$

Writing the orthogonality of unit outward normals to the tangential plane to the surface of body (m):

$$ \mathbf{N} = \frac{\mathbf{G}_1 \wedge \mathbf{G}_2}{\|\mathbf{G}_1 \wedge \mathbf{G}_2\|} \qquad (9) $$

$$ \mathbf{n} = \frac{\mathbf{g}_1 \wedge \mathbf{g}_2}{\|\mathbf{g}_1 \wedge \mathbf{g}_2\|} = \frac{(\mathbf{F}_t^{(m)})^{-T} \mathbf{N}}{\|(\mathbf{F}_t^{(m)})^{-T} \mathbf{N}\|} \qquad (10) $$

As proposed by Klarbring [12], it is also possible to treat $\mathbf{x}^{(s)}$ like a translation of $\bar{\mathbf{x}}^{(m)} \in \gamma_c^{(m)}$ in the neighbourhood of $\Phi^{(m)}(\bar{\xi})$ defined by

$$ \mathbf{x}_{\mathrm{tr}}^{(m)}(\bar{\xi}) = \Phi^{(m)}(\bar{\xi}) + d_n \mathbf{n}(\bar{\xi}) \qquad (11) $$

Natural base vectors of the tangent plane in $\mathbf{x}_{\mathrm{tr}}^m(\bar{\xi})$ are

$$ \hat{\mathbf{g}}_{\alpha} = \mathbf{g}_{\alpha} + d_n \frac{\partial \mathbf{n}}{\partial \bar{\xi}_{\alpha}}, \quad \alpha = 1, n-1 \qquad (12) $$

The signed distance between a slave particle occupying position $\mathbf{x}^{(s)}$ and the master surface can be then introduced:

$$ \mathbf{d}(\bar{\xi}) = (\mathbf{x}^{(s)} - \mathbf{x}^{(m)}(\bar{\xi})) = d_n \mathbf{n} \qquad (13) $$

These two tangential vectors are obviously identical when contact is established, i.e. when $d_n = 0$.
## 2.3. Relative velocity

We now turn to defining the relative velocity between $\mathbf{x}^{(s)}$ and $\bar{\mathbf{x}}^{(m)}$. Taking the total time derivative of Equation (13), we obtain with $\dot{\mathbf{n}} = (\partial\mathbf{n}/\partial t + \dot{\xi}^\alpha \partial\mathbf{n}/\partial\xi^\alpha)$:

$$ \mathbf{v}^{(s)} - \bar{\mathbf{v}}^{(m)} - \dot{\xi}^{\alpha} \mathbf{g}_{\alpha} = \dot{\mathbf{d}}_{n} \mathbf{n} + d_{n} \dot{\mathbf{n}} \quad (14) $$

One of the main points now consists of introducing an objective relative velocity between the slave point and the master surface either in the master local reference co-ordinate system or in the master-translated local reference co-ordinate system. Due to the term $\dot{\mathbf{n}}$, the velocities defined above are not objective ones. Using the local basis introduced above, the following relative Eulerian velocities $\mathbf{v}_r$ or $\hat{\mathbf{v}}_r$ can be written as

$$ \mathbf{v}_r = \mathbf{v}^{(s)} - \bar{\mathbf{v}}^{(m)} - d_n \dot{\mathbf{n}} = \dot{\xi}^\alpha \mathbf{g}_\alpha + \dot{d}_n \mathbf{n} \quad (15) $$

$$ \hat{\mathbf{v}}_r = \mathbf{v}^{(s)} - \bar{\mathbf{v}}^{(m)} - d_n \frac{\partial \mathbf{n}}{\partial t} = \dot{\xi}^\alpha \hat{\mathbf{g}}_\alpha + \dot{d}_n \mathbf{n} \quad (16) $$

Equations (15) and (16) show that in both local bases, $\dot{d}_n$ is the normal component of the relative velocities. $\mathbf{v}_g = \dot{\xi}^\alpha \mathbf{g}_\alpha$ and $\hat{\mathbf{v}}_g = \dot{\xi}^\alpha \hat{\mathbf{g}}_\alpha$ are the relative tangential velocities, which are identical when $d_n = 0$.

# 3. UNILATERAL CONTACT AND FRICTION LAW

First, the unilateral contact conditions require the Cauchy contact traction vector onto the slave surface to be decomposed into

$$ \mathbf{t}^{(s)} = t_n \mathbf{n} + \mathbf{t}_g = t_n \mathbf{n} + \mathbf{t}_g \quad (17) $$

with, depending on the contravariant basis used:

$$ \mathbf{t}_n = \mathbf{t}^{(s)} \cdot \mathbf{n} \quad (18) $$

$$ \mathbf{t}_g = \mathbf{t}^{(s)} - t_n \mathbf{n} = \mathbf{t}_g^\alpha \mathbf{g}^\alpha = \mathbf{t}_g^\alpha \hat{\mathbf{g}}^\alpha \quad (19) $$

The unilateral complementarity conditions for normal contact are specified as the following complementarity conditions:

$$ d_n \ge 0, \quad t_n \ge 0, \quad d_n t_n = 0 \quad (20) $$

An alternative equivalent form to Equation (20) is

$$ t_n \in C_n, \quad \forall s_n \in C_n, \quad d_n(s_n - t_n) \ge 0 \quad (21) $$

with

$$ C_n = \{ p_n : \gamma_c^{(s)} \to \mathbb{R}; p_n \ge 0 \} \quad (22) $$

The friction law used here is Coulomb's friction law. Let $C(t_n, \varphi)$ be the Coulomb's cone:

$$ C(t_n, \varphi) = \{ \mathbf{t}_g \in \mathbb{R}^n, \mathbf{t}_g \cdot \mathbf{n} = 0, \|\mathbf{t}_g\| \le \mu t_n \} \quad (23) $$
where $||\mathbf{a}|| = [\mathbf{a}^\alpha m_{\alpha\beta} \mathbf{a}^\beta]^{1/2}$ with $m_{\alpha\beta} = \mathbf{g}_\alpha \cdot \mathbf{g}_\beta$. An equivalent way of writing this can be obtained if $\mathbf{g}$ is replaced by $\hat{\mathbf{g}}$, then $||\mathbf{a}|| = [\hat{\mathbf{a}}^\alpha \hat{m}_{\alpha\beta} \hat{\mathbf{a}}^\beta]^{1/2}$ with $\hat{m}_{\alpha\beta} = \hat{\mathbf{g}}_\alpha \cdot \hat{\mathbf{g}}_\beta$.

Coulomb's friction law can be expressed in a variational form:

$$ \mathbf{t}_g \in C(t_n, \varphi), \quad \forall \mathbf{s}_g \in C(t_n, \varphi), \quad \mathbf{v}_g \cdot (\mathbf{s}_g - \mathbf{t}_g) = \frac{\dot{\xi}^\alpha}{\xi^\alpha} (s_{g\alpha} - t_{g\alpha}) \ge 0 \quad (24) $$

The unilateral conditions (21) and the Coulomb's friction law (24) both can be written in the following integral forms:

$$ \forall s_n \in C_n \quad \int_{\gamma_c^{(s)}} d_n(s_n - t_n) da_c \ge 0 \quad (25) $$

$$ \forall \mathbf{s}_g \in C(t_n, \varphi) \quad \int_{\gamma_c^{(s)}} \mathbf{v}_g \cdot (\mathbf{s}_g - \mathbf{t}_g) da_c \ge 0 \quad (26) $$

### 3.1. Virtual work of contact

We consider test function $\tilde{\varphi}^{(i)}$ for each body (i) in spaces $T_{\varphi}^{(i)}$ of kinematically admissible variations (virtual displacements). The virtual contact work is then given by

$$ (G_c) = - \int_{\gamma_c^{(s)}} \mathbf{t}^{(s)} \cdot \tilde{\varphi}^{(s)} da_c^{(s)} - \int_{\gamma_c^{(m)}} \mathbf{t}^{(m)} \cdot \tilde{\varphi}^{(m)} da_c^{(m)} \quad (27) $$

The law of action and reaction can be expressed as follows:

$$ \mathbf{t}^{(s)} da_c^{(s)} = -\mathbf{t}^{(m)} (\bar{\mathbf{x}}^{(m)}) da_c^{(m)} \quad (28) $$

### 3.2. Standard slave/master formulation

Using Equation (28), we express the virtual contact work with a single integral:

$$ G_c = - \int_{\gamma_c^{(s)}} \mathbf{t}^{(s)} \cdot (\tilde{\varphi}^{(s)}(\mathbf{X}^{(s)}) - \tilde{\varphi}^{(m)}(\bar{\xi}(\mathbf{X}^{(s)}))) da_c^{(s)} \quad (29) $$

The contribution of the part of the master potential surface $\gamma_c^{(m)}$ not involved in the contact ($\mathbf{t}^{(m)} = 0$) has been simply omitted.

In order to write the virtual contact work in terms of normal and tangential contributions, some preliminary results are needed. Following Laursen and Simo [13] for a quantity $\mathbf{F}(\mathbf{X}^{(s)}, \varphi^{(s)}, \varphi^{(m)})$ we introduce its Gateau derivative $\delta\mathbf{F}(\mathbf{X}^{(s)}, \varphi^{(s)}, \varphi^{(m)})$ by

$$ \delta\mathbf{F}(\mathbf{X}^{(s)}, \varphi^{(s)}, \varphi^{(m)}) = \frac{d}{d\epsilon}\mathbf{F}(\mathbf{X}^{(s)}, \varphi_{\epsilon}^{(s)}, \varphi_{\epsilon}^{(m)})|_{\epsilon=0} \quad (30) $$

where the perturbed fields $\varphi_{\epsilon}^{(s)}$ and $\varphi_{\epsilon}^{(m)}$ are defined by

$$ \varphi_{\epsilon}^{(s)} = \varphi^{(s)} + \epsilon \tilde{\varphi}^{(s)}, \quad \varphi_{\epsilon}^{(m)} = \varphi^{(m)} + \epsilon \tilde{\varphi}^{(m)} \quad (31) $$

Using Equation (8), the variation of $\varphi^{(m)}(\bar{\xi})$ is given by

$$ \delta(\varphi^{(m)}(\bar{\xi})) = \delta\bar{\xi}^{\alpha} \mathbf{g}_{\alpha} + \tilde{\varphi}^{(m)}(\bar{\xi}) \quad (32) $$
Furthermore, noting that $t_{g_x}d_n=0$, $\alpha=1, n-1$ and $t_n d_n=0$, which along with (28) constitute the strong law of action and reaction, we write the following variation:

$$ \mathbf{t}^{(s)} \cdot \delta \left( \mathbf{\phi}^{(s)}(\mathbf{X}^{(s)}) - \mathbf{\phi}^{(m)}(\bar{\mathbf{\xi}}) \right) = \mathbf{t}(s) \cdot \delta(d_n \mathbf{n}) = (t_n \mathbf{n} + t g_x \mathbf{g}^\alpha)(\delta d_n \mathbf{n} + d_n \delta \mathbf{n}) = t_n \delta d_n \quad (33) $$

This leads to the following expression:

$$ \mathbf{t}^{(s)} \cdot (\tilde{\mathbf{\phi}}^{(s)} - \tilde{\mathbf{\phi}}^{(m)}(\bar{\mathbf{\xi}})) = t_n \delta d_n + \delta \bar{\mathbf{\xi}}^\alpha t_{g_x} \quad (34) $$

Using these relations one can write:

$$ \int_{\gamma_{\varepsilon}^{(s)}} \mathbf{t}^{(s)} \cdot (\tilde{\mathbf{\phi}}^{(s)}(\mathbf{X}^{(s)}) - \tilde{\mathbf{\phi}}^{(m)}(\bar{\mathbf{\xi}}(\mathbf{X}^{(s)}))) \, da_{\varepsilon}^{(s)} = \int_{\gamma_{\varepsilon}^{(s)}} t_n \delta d_n + t_{g_x} \delta \bar{\mathbf{\xi}}^{\alpha} \, da_{\varepsilon}^{(s)} \quad (35) $$

### 3.3. Symmetrical master/slave contact formulation

In order to introduce our symmetric master-slave formulation, we split the virtual contact work into two parts, which are weighted by a coefficient $(\theta)$

$$ \int_{\gamma_{\varepsilon}^{(1)}} \mathbf{t}^{(1)} \cdot \tilde{\mathbf{\phi}}^{(1)} \, da_{\varepsilon}^{(1)} + \int_{\gamma_{\varepsilon}^{(2)}} \mathbf{t}^{(2)} \cdot \tilde{\mathbf{\phi}}^{(2)} \, da_{\varepsilon}^{(2)} = \theta \int_{\gamma_{\varepsilon}^{(1)}} (t_n^{(1)} \delta d_n^{(1)} + t_{g_x}^{(1)} \delta \bar{\mathbf{\xi}}^{\alpha}) \, da_{\varepsilon}^{(1)} \\ + (1-\theta) \int_{\gamma_{\varepsilon}^{(2)}} (t_n^{(2)} \delta d_n^{(2)} + t_{g_x}^{(2)} \delta \bar{\mathbf{\xi}}^{\alpha}) \, da_{\varepsilon}^{(2)} \quad (36) $$

Taking $\theta=1$ (resp. $\theta=0$) leads to surface 1 being the slave (resp. master) and surface 2 being the master (resp. slave). The choice of parameter $\theta$ will not be fully analysed with the examples presented in this paper. But one can imagine different strategies, such as taking $\theta$ to depend on the material properties or the mesh size [18].

Using the previous results, we will now obtain the discrete system to be solved. Finally, we have to solve

Find $\mathbf{\phi}^{(i)}: \Omega^i \times \mathbb{R} \to \mathbb{R}^n, t_n^{(i)} \in C_n^{(i)}$ and $\mathbf{t}_g^{(i)} \in C^{(i)}(t_n, \mathbf{\phi})$ for $i=1,2$ such that
$$ \forall \tilde{\mathbf{\phi}}^{(i)} \in T_{\phi}^{(i)}: $$

$$ \sum_{i=1,2} \left( \int_{\omega^{(i)}} \mathbf{\sigma}^{(i)} : \nabla^s \tilde{\mathbf{\phi}}^{(i)} \, dv^{(i)} - \int_{\omega^{(i)}} \mathbf{f}_v^{(i)} \cdot \tilde{\mathbf{\phi}}^{(i)} \, dv^{(i)} - \int_{\gamma_f^{(i)}} \mathbf{f}_s^{(i)} \cdot \tilde{\mathbf{\phi}}^{(i)} \, da_f^{(i)} \right) \\ - \theta \int_{\gamma_{\varepsilon}^{(1)}} t_n^{(1)} \delta d_n^{(1)} + t_{g_x}^{(1)} \delta \bar{\mathbf{\xi}}^{\alpha(1)} \, da_{\varepsilon}^{(1)} - (1-\theta) \int_{\gamma_{\varepsilon}^{(2)}} t_n^{(2)} \delta d_n^{(2)} + t_{g_x}^{(2)} \delta \bar{\mathbf{\xi}}^{\alpha(2)} \, da_{\varepsilon}^{(2)} = 0 \quad (37) $$

$$ \forall s_n^{(i)} \in C_n^{(i)} \quad \int_{\gamma_{\varepsilon}^{(i)}} d_n^{(i)} (s_n^{(i)} - t_n^{(i)}) \, da_{\varepsilon}^{(i)} \ge 0 $$

$$ \forall \mathbf{s}_g^{(i)} \in C^{(i)}(t_n^{(i)}, \mathbf{\phi}^{(1)}, \mathbf{\phi}^{(2)}) \quad \int_{\gamma_{\varepsilon}^{(i)}} \mathbf{v}_r^{(i)} \cdot (\mathbf{s}_g^{(i)} - \mathbf{t}_g^{(i)}) \, da_{\varepsilon}^{(i)} \ge 0 $$
<figure><img src="image_3.png" /><figcaption>Figure 2. Contact area.</figcaption></figure>

where $\mathbf{f}_v^{(i)}$ are body forces per unit volume, and $\mathbf{f}_s^{(i)}$ are surface forces per unit area, both of which are applied to solid $\Omega^{(i)}$.

## 4. FINITE ELEMENT DISCRETIZATION

As the discrete complementarity contact with friction problem will be introduced below in the two-dimensional (2D) case, we restrict our presentation of the finite element discretization procedure in this section to the 2D case. In addition, linear triangular or quadrilateral elements will be used. It is now easy to define the variables of the problem in this restricted framework (Figure 2).

For a linear element we obtain

$$ \bar{\mathbf{x}}^{(m)} = (1 - \bar{\xi}^1)\mathbf{x}^m + \bar{\xi}^1 \mathbf{x}^{m-1} \quad (38) $$

To write the discrete problem, it is convenient to introduce here contact elements in the form of the slave node and master segment couple. It is clear that these elements are liable to evolve during the computation. Let the length $l$ of the master segment be defined by

$$ l = \|\mathbf{x}^{m-1} - \mathbf{x}^m\| \quad (39) $$

For a given slave node $\mathbf{x}^{(s)}$, let $\mathbf{x}^m$ and $\mathbf{x}^{m-1}$ be the position of the nodes of the master segment containing the projection of the slave node. Its curvilinear co-ordinate $\bar{\xi}^1$ is given by

$$ \bar{\xi}^1 = \frac{\mathbf{x}^{(s)} - \mathbf{x}^m}{l} \cdot \mathbf{g} \quad (40) $$

The local unit basis vector $\mathbf{g}$ is given by

$$ \mathbf{g} = \frac{\mathbf{x}^{m-1} - \mathbf{x}^m}{l} \quad \text{and} \quad \mathbf{n} = \mathbf{g} \wedge \mathbf{e}_z \quad (41) $$

The choice of a unit basis vector $\mathbf{g}$ is made in the 2D case, to avoid having to deal with the metric in the Coulomb's friction law. If the covariant or contravariant vectors introduced above are used, the metric appears in the expression of the Coulomb's cone and the change
of variable used to write the whole problem as a complementarity one has to include the metric, which is unknown. The unit vector **g** is related to the covariant or contravariant vector by

$$ \mathbf{g} = \frac{\hat{\mathbf{g}}_l}{l} = \hat{\mathbf{g}}^l l \qquad (42) $$

Hence, the tangential contact vector onto the slave surface can be written as

$$ \mathbf{t}_g = \hat{t}_1 \hat{\mathbf{g}}^1 = t_t \mathbf{g} \quad \text{with} \quad t_t = \frac{\hat{t}_1}{l} \qquad (43) $$

The tangential relative velocity is given by

$$ \mathbf{v}_g = \dot{\xi}^1 \hat{\mathbf{g}}_1 = v_t \mathbf{g} \quad \text{with} \quad v_t = l \dot{\xi}^1 \qquad (44) $$

The relative distance vector **d** is expressed in the local basis:

$$ \mathbf{d} = d_n \mathbf{n} = \mathbf{x}^{(s)} - (1 - \bar{\xi}^1) \mathbf{x}^m - \bar{\xi}^1 \mathbf{x}^{m-1} \qquad (45) $$

Let the variations of position in the slave and master configuration be defined by: $\delta\mathbf{x}^{(s)} = \delta\mathbf{u}^{(s)}$ and $\delta\mathbf{x}^m = \delta\mathbf{u}^m$, respectively, where $\delta\mathbf{u}^{(s)}$ and $\delta\mathbf{u}^m$ are displacement variations defined by

$$ \begin{aligned} (\delta\mathbf{u}^{(s)})^T &= (\delta u_x^{(s)}, \delta u_y^{(s)}) \\ (\delta\mathbf{u}^m)^T &= (\delta u_x^m, \delta u_y^m) \end{aligned} \qquad (46) $$

In order to discretize the integrals including the contact related terms in (37), it is necessary to determine the variations of $d_n$ and $\bar{\xi}$:

$$ \delta d_n = (\delta \mathbf{u}^{(s)} - (1 - \bar{\xi}^1) \delta \mathbf{u}^m - \bar{\xi}^1 \delta \mathbf{u}^{m-1}) \cdot \mathbf{n} \qquad (47) $$

$$ \delta \bar{\xi}^1 = \frac{(\delta \mathbf{u}^{(s)} - (1 - \bar{\xi}^1) \delta \mathbf{u}^m - \bar{\xi}^1 \delta \mathbf{u}^{m-1})}{l} \cdot \mathbf{g} + \frac{(\delta \mathbf{u}^{m-1} - \delta \mathbf{u}^m)}{l^2} \cdot d_n \mathbf{n} \qquad (48) $$

These equations can be written in matrix form. Let us first introduce the following vectors for each contact element $e$:

$$ \begin{aligned} \delta \mathbf{u}^e &= \begin{pmatrix} \delta \mathbf{u}^{(s)} \\ \delta \mathbf{u}^m \\ \delta \mathbf{u}^{m-1} \end{pmatrix}, \quad \mathbf{N}_0^e = \begin{pmatrix} \mathbf{0} \\ -\mathbf{n} \\ \mathbf{n} \end{pmatrix}, \quad \mathbf{T}_0^e = \begin{pmatrix} \mathbf{0} \\ -\hat{\mathbf{g}}_l \\ \hat{\mathbf{g}}_l \end{pmatrix} \\ \mathbf{N}_1^e &= \begin{pmatrix} \mathbf{n} \\ -(1 - \bar{\xi}^1)\mathbf{n} \\ -\bar{\xi}^1 \mathbf{n} \end{pmatrix}, \quad \mathbf{T}_1^e = \begin{pmatrix} \hat{\mathbf{g}}_l \\ -(1 - \bar{\xi}^1)\hat{\mathbf{g}}_l \\ -\bar{\xi}^1 \hat{\mathbf{g}}_l \end{pmatrix} \end{aligned} \qquad (49) $$
Using the variations established in Equations (47), (48) and the matrices given in Equation (49), the virtual work of contact $W_c$ can be written using contact elements:

$$ W_c = \theta \int_{\gamma_c^{(1)}} \delta \mathbf{u}^{eT} \left[ t_n \mathbf{N}_1^e + t_t l \left( \frac{\mathbf{T}_1^e}{l} + \frac{d_n \mathbf{N}_0^e}{l^2} \right) \right]^{(1)} \mathrm{d}a_c \\ + (1-\theta) \int_{\gamma_c^{(2)}} \delta \mathbf{u}^{eT} \left[ t_n \mathbf{N}_1^e + t_t l \left( \frac{\mathbf{T}_1^e}{l} + \frac{d_n \mathbf{N}_0^e}{l^2} \right) \right]^{(2)} \mathrm{d}a_c \quad (50) $$

Introducing a standard finite element discretization one can express the previous continuous problem in term of discrete unknowns:

Find $\mathbf{u}^h \in T_{\phi}^h \Psi, \mathbf{R}_n \in C_n^h$ and $\mathbf{R}_g \in C^h(R_n, \mathbf{u}^h)$ such that

$\forall \delta \mathbf{u}^h \in T_{\phi}^h \Psi:$

$$ \delta \mathbf{u}^h \cdot [\mathbf{F}_{\text{int}}(\mathbf{u}^h) - \mathbf{F}_{\text{ext}}(\mathbf{u}^h) - (\mathbf{H}_n(\mathbf{u}^h)\mathbf{R}_n + \mathbf{H}_g(\mathbf{u}^h)\mathbf{R}_g)] = 0 \quad (51) $$

$$ \forall \mathbf{S}_n \in C_n^h \quad \mathbf{d}_n \cdot (\mathbf{S}_n - \mathbf{R}_n) \ge 0 $$

$$ \forall \mathbf{S}_g \in C_g^h(R_n) \quad \mathbf{v}_t \cdot (\mathbf{S}_g - \mathbf{R}_g) \ge 0 $$

For the sake of simplification, we will omit superscript $h$ below. In (51) $\mathbf{R}_n$ and $\mathbf{R}_g$ are the nodal normal and tangential contact force vector, respectively. The dimension of each of these vectors is nc. nc is the number of contact nodes and ndf will be taken to indicate the total number of degrees of freedom. $\mathbf{d}_n$ and $\mathbf{v}_t$ are vectors of $\mathbb{R}^{nc}$ and denote the nodal relative distance and the nodal relative velocity. The nc × ndf matrices $\mathbf{H}_n$ and $\mathbf{H}_g$ arise from the numerical integration of the discrete problem associated with Equation (51).

The following discrete subspaces have been introduced:

$$ T_{\phi}^{h} \Psi := \{ \mathbf{u} \in \mathbb{R}^{\text{ndf}}; u_i|_{\Gamma_u} = 0; i = 1, \dots, \text{ndf} \} \quad (52) $$

$$ C_n^h := \{ \mathbf{R}_n \in \mathbb{R}^{nc} \text{ such that } R_{n_i} \ge 0 \ \forall i = 1, \dots, \text{nc} \} \quad (53) $$

$$ C_g^h(R_n) := \{ \mathbf{R}_g \in \mathbb{R}^{nc} \text{ such that } |R_{g_i}| \le \mu R_{n_i} \ \forall i = 1, \dots, \text{nc} \} \quad (54) $$

## 5. NEWTON'S METHOD FOR GENERALIZED EQUATIONS

As proposed by Klarbring and Bjorkman [19], we write the discretized problem as generalized equations in the sense of Robinson [20], and we use the extension of Newton's method to generalized equations [21]. For this purpose, the problem has to be written on the following form:

$$ \mathbf{0} \in \mathbf{f}(\mathbf{x}) + \mathcal{N}_c(\mathbf{x}) \quad (55) $$
where $\mathbf{f}(\mathbf{x})$ is a non-linear operator and $\mathcal{N}_C(\mathbf{x})$ is the normal cone for a non-empty closed convex set $C$ in $\mathbb{R}^n$:

$$ \mathcal{N}_C(\mathbf{x}) = \begin{cases} \{\mathbf{y} \in \mathbb{R}^n / \mathbf{y}^\mathrm{T}(\mathbf{x}^\star - \mathbf{x}) \le 0 \ \forall \mathbf{x}^\star \in C\} & \text{if } \mathbf{x} \in C \\ \emptyset & \text{if } \mathbf{x} \notin C \end{cases} \quad (56) $$

When $C$ is the whole $\mathbb{R}^n$, we have $\mathcal{N}_{\mathbb{R}^n}(\mathbf{x}) = \{\mathbf{0}\}$ and $\mathbf{f}(\mathbf{x}) = 0$.

Furthermore, Equation (55) is equivalent to the variational inequality:

$$ \mathbf{x} \in C, \quad \mathbf{f}(\mathbf{x})^{\mathrm{T}}(\mathbf{x}^{\star} - \mathbf{x}) \ge 0 \quad \forall \mathbf{x}^{\star} \in C \quad (57) $$

In the framework of an incremental decomposition of the loading, at each increment, at an iteration $(j+1)$ of Newton's method, the problem to be solved is:

$$ \text{Find } \mathbf{x}^{j+1} \text{ such that} \\ 0 \in \mathbf{f}(\mathbf{x}^j) + \nabla_x \mathbf{f}(\mathbf{x}^j)(\mathbf{x}^{j+1} - \mathbf{x}^j) + \mathcal{N}_C(\mathbf{x}^{j+1}) \quad (58) $$

The frictional contact problem previously set can then be written in the identical form to problem (55):

Find $\mathbf{u} \in T_{\varphi} \Psi, \mathbf{R}_n \in C_n, \mathbf{R}_g \in C_g(\mathbf{R}_n)$ such that

$$ \mathbf{0} \in \begin{bmatrix} \mathcal{L}(\mathbf{u}, \mathbf{R}_n, \mathbf{R}_g) \\ \mathbf{d}_n \\ \mathbf{v}_t \end{bmatrix} + \mathcal{N}_C \begin{bmatrix} \mathbf{u} \\ \mathbf{R}_n \\ \mathbf{R}_g \end{bmatrix} \quad (59) $$

with

$$ \mathcal{L}(\mathbf{u}, \mathbf{R}_n, \mathbf{R}_g) = (\mathbf{F}_{\mathrm{int}}(\mathbf{u}) - \mathbf{F}_{\mathrm{ext}}(\mathbf{u}) - (\mathbf{H}_n(\mathbf{u})\mathbf{R}_n + \mathbf{H}_g(\mathbf{u})\mathbf{R}_g)) $$

In Equation (59), $C = \mathbb{R}^{\mathrm{ndf}} \times C_n \times C_g(\mathbf{R}_n)$. At a given loading increment, at an iteration $(j+1)$ of the Newton-Raphson method, the problem to be solved is:

Let $(\mathbf{u}^0, \mathbf{R}_n^0, \mathbf{R}_g^0)$ be given

Find $\mathbf{u}^{j+1} \in T_{\varphi} \Psi, \mathbf{R}_n^{j+1} \in C_n$ and $\mathbf{R}_g^{j+1} \in C_g(\mathbf{R}_n^{j+1})$ such that

$$ 0 \in \begin{bmatrix} \mathcal{L}^j \\ \mathbf{d}_n^j \\ \mathbf{v}_t^j \end{bmatrix} + \begin{bmatrix} \nabla_{u, R_n, R_g} \mathcal{L}^j \\ \nabla_{u, R_n, R_g} \mathbf{d}_n^j \\ \nabla_{u, R_n, R_g} \mathbf{v}_t^j \end{bmatrix}, \quad \begin{bmatrix} \mathbf{d} \mathbf{u}^{j+1} \\ \mathbf{R}_n^{j+1} - \mathbf{R}_n^j \\ \mathbf{R}_g^{j+1} - \mathbf{R}_g^j \end{bmatrix} + \mathcal{N}_C \begin{bmatrix} \mathbf{u}^{j+1} \\ \mathbf{R}_n^{j+1} \\ \mathbf{R}_g^{j+1} \end{bmatrix} \quad (60) $$

with

$$ \begin{align} \nabla_{u, R_n, R_g} \mathcal{L}^j &= [\nabla_u \mathcal{L}^j - \mathbf{H}_n(\mathbf{u}^j) - \mathbf{H}_g(\mathbf{u}^j)] \nonumber \\ \nabla_{u, R_n, R_g} \mathbf{d}_n^j &= [\nabla_u \mathbf{d}_n^j \quad 0 \quad 0] = [\Delta \mathbf{H}_n(\mathbf{u}^j) \quad 0 \quad 0] \nonumber \\ \nabla_{u, R_n, R_g} \mathbf{v}_t^j &= [\nabla_u \mathbf{v}_t^j \quad 0 \quad 0] = [\Delta \mathbf{H}_g(\mathbf{u}^j) \quad 0 \quad 0] \tag{61} \\ \nabla_u \mathcal{L}^j &= \mathbf{K}(\mathbf{u}^j) - \nabla_u \mathbf{H}_n^j \mathbf{R}_n^j - \nabla_u \mathbf{H}_g^j \mathbf{R}_g^j = \hat{\mathbf{K}} \nonumber \\ \mathcal{L}^j &= \mathcal{L}(\mathbf{u}^j, \mathbf{R}_n^j, \mathbf{R}_g^j) \nonumber \\ \mathbf{d} \mathbf{u}^{j+1} &= \mathbf{u}^{j+1} - \mathbf{u}^j \nonumber \end{align} $$
If the previous system is split into three contributions, a more usual presentation is recovered:

* Equilibrium equation

$$ \hat{\mathbf{K}} \, du^{j+1} - \mathbf{H}_n(\mathbf{u}^j) \mathbf{R}_n^{j+1} - \mathbf{H}_g(\mathbf{u}^j) \mathbf{R}_g^{j+1} = - \mathbf{F}_{\text{int}}(\mathbf{u}^j) + \mathbf{F}_{\text{ext}}(\mathbf{u}^j) \quad (62) $$

* Contact variational inequality

$$ \begin{gathered} \mathbf{R}_n^{j+1} \in C_n \quad \forall \mathbf{S}_n \in C_n \\ [\mathbf{d}_n^j + \Delta \mathbf{H}_n^T(\mathbf{u}^j) \, du^{j+1}]_t [\mathbf{S}_n - \mathbf{R}_n^{j+1}]_t \ge 0 \quad \forall t = 1, \dots, nc \end{gathered} \quad (63) $$

* Friction quasi-variational inequality

$$ \begin{gathered} \mathbf{R}_g^{j+1} \in C_g(\mathbf{R}_n^{j+1}) \quad \forall \mathbf{S}_g \in C_g(\mathbf{R}_n^{j+1}) \\ [\mathbf{v}_t^j + \Delta \mathbf{H}_g^T(\mathbf{u}^j) \, du^{j+1}]_t (\mathbf{S}_g - \mathbf{R}_g^{j+1})_t \ge 0 \quad \forall t = 1, \dots, nc \end{gathered} \quad (64) $$

Comments

* The contact condition can be equally well written in the form of complementarity relations:

$$ (\mathbf{R}_n^{j+1})_t \ge 0 $$

$$ [\mathbf{d}_n^j + \Delta \mathbf{H}_n^T(\mathbf{u}^j) \, du^{j+1}]_t \ge 0 $$

and

$$ (\mathbf{R}_n^{j+1})_t \cdot [\mathbf{d}_n^j + \Delta \mathbf{H}_n^T(\mathbf{u}^j) \, du^{j+1}]_t = 0 \quad \forall t = 1, \dots, nc \quad (65) $$

* The friction conditions can also be written as a threshold law:

$$ \| (\mathbf{R}_g^{j+1})_t \| \le \mu (\mathbf{R}_n^{j+1})_t \quad \forall t = 1, \dots, nc $$

with

$$ \| (\mathbf{R}_g^{j+1})_t \| < \mu (\mathbf{R}_n^{j+1})_t \Rightarrow [\mathbf{v}_t^j + \Delta \mathbf{H}_g^T(\mathbf{u}^j) \, du^{j+1}]_t = 0 \quad (66) $$

$$ \| (\mathbf{R}_g^{j+1})_t \| = \mu (\mathbf{R}_n^{j+1})_t \Rightarrow [\mathbf{v}_t^j + \Delta \mathbf{H}_g^T(\mathbf{u}^j) \, du^{j+1}]_t = -\lambda (\mathbf{R}_g^{j+1})_t $$

## 6. MATHEMATICAL PROGRAMMING SOLVER

The frictional contact problem defined above can be solved using a mathematical programming method, such as the Lemke's one. Methods of this kind are appropriate for solving complementarity problems written in the following form:

Find $\mathbf{U} \in \mathbb{R}^p, \mathbf{R} \in \mathbb{R}^p$ such that:

$$ \mathbf{R} - \mathbf{KU} = \mathbf{F} \quad (67) $$

$$ U_j \ge 0, R_j \ge 0, \quad U_j R_j = 0 \quad \forall j = 1, \dots, p $$

where $\mathbf{F}$ is a known vector of $\mathbb{R}^p$ and $\mathbf{K}$ a $\mathbb{R}^p \times \mathbb{R}^p$ matrix.
It can be observed in problem (67) that all the components of the unknowns are constrained by complementarity conditions. In order to use this kind of method to solve the frictional contact problem posed above, several stages are required. First, the constraints have to be set as complementarity ones. This means that the friction conditions have to be rewritten, the contact kinematic variables have to be introduced into the system to be solved and the contact forces have to be described in terms of the local contact referentials. Lastly, a condensation procedure has to be used to reduce the problem to the sole variables involved in the contact. Details of these steps are given in the following sections.

## 6.1. Friction complementarity relations

In dealing with the 3D case, Klarbring and Bjorkman [19] have introduced a piecewise linear friction law approximating Coulomb's friction law. This discretization procedure makes it possible to write the friction relations as complementarity conditions and to set the problem as a linear complementarity one. In the present study, we restrict ourselves to a 2D analysis and their approach leads to introducing new couples of variables, ($\lambda$ and $\phi$), which define the boundary of the Coulomb's cone. The Kuhn Tucker conditions for the frictional contact problem leads to the following expression for a contact node:

$$ 
\begin{aligned}
\phi_1(R_n, R_g) &= -R_g + \mu R_n \\
\phi_2(R_n, R_g) &= +R_g + \mu R_n \\
v_t &= - \left( \lambda_1 \frac{\partial \phi_1}{\partial R_g} + \lambda_2 \frac{\partial \phi_2}{\partial R_g} \right) \\
\lambda_i \ge 0, \quad \phi_i \ge 0, \quad \lambda_i \phi_i = 0, \quad i=1,2
\end{aligned} 
\qquad (68) $$

Thus for each contact node, the unknowns are $d_n$, $\lambda_1$, $\lambda_2$, and their complementarity variable, $R_n$, $\phi_1$ and $\phi_2$. The use of the condensation procedure discussed below makes it possible, upon introducing variables $\lambda$ and $\phi$, to write a system in which all the components are constrained by complementarity conditions. This linear complementarity problem with a 3nc × 3nc square singular matrix (nc being the total number of nodes involved in the contact) can then be straightforwardly solved using a pivoting algorithm such Lemke's method [22].

## 6.2. Matrix construction

If the whole system is translated into the global referential co-ordinate system, first the components of the contact forces have to be accounted for in the local referential co-ordinate system (n,g). In the second step, the kinematic contact variables, namely the contact distance between a slave node and its master segment $d_n$ and the relative tangential velocity, have to be introduced for the slave nodes or in the symmetrical formulation for the slave/master nodes. The unknown variables of the master nodes above are the local components of their incremental displacements. The problem is then written in a similar form to Equation (67) and it can be straightforwardly solved using Lemke's method. These procedures will be described below in greater details in the case of standard and symmetrical slave/master formulation. As a matter of fact, there are large differences between the two cases, which justifies the two presentations.
6.2.1. Standard slave/master formulation. Let us define a contact element by the association between a contact node (denoted in what follows by the superscript (s)) and the nodes which characterize its contacting segment (denoted by the superscript (m) below). Then for each contact element it is possible to write the following linear system schematically as follows:

$$ \begin{bmatrix} \mathbf{K}^{ss} & \mathbf{K}^{sm} & \mathbf{K}^{l1} \\ \mathbf{K}^{ms} & \mathbf{K}^{mm} & \mathbf{K}^{l2} \end{bmatrix}^j \begin{bmatrix} \mathbf{du}^{(s)j+1} \\ \mathbf{du}^{(m)j+1} \\ \mathbf{du}^{(l)j+1} \end{bmatrix} = \begin{bmatrix} \mathbf{F}^{(s)} \\ \mathbf{F}^{(m)} \end{bmatrix} + \begin{bmatrix} \mathbf{Q} & \mathbf{0} \\ \mathbf{H} & \mathbf{I} \end{bmatrix}^j \begin{bmatrix} \mathbf{R}^{(s)j+1} \\ \mathbf{R}^{(m)j+1} \end{bmatrix} \quad (69) $$

where

$$ \begin{align} \left[ \mathbf{du}^{(s)j+1} \right]^T &= \left[ \mathbf{du}_x^{(s)j+1} \mathbf{du}_y^{(s)j+1} \right] \nonumber \\ \left[ \mathbf{du}^{(m)j+1} \right]^T &= \left[ \mathbf{du}_x^{(m)j+1} \mathbf{du}_y^{(m)j+1} \mathbf{du}_x^{(m-1)j+1} \mathbf{du}_y^{(m-1)j+1} \right] \tag{70} \\ \left[ \mathbf{R}^{(s)j+1} \right]^T &= \left[ \mathbf{R}_n^{(s)j+1} \mathbf{R}_g^{(s)j+1} \right] \nonumber \\ \left[ \mathbf{R}^{(m)j+1} \right]^T &= \left[ \mathbf{R}_n^{(m)j+1} \mathbf{R}_g^{(m)j+1} \mathbf{R}_n^{(m-1)j+1} \mathbf{R}_g^{(m-1)j+1} \right] \nonumber \end{align} $$

Superscript *l* indicates nodes which do not belong to the contact element under consideration (i.e. *l* = ndf - 6) and in Equation (69), du<sup>(l)j+1</sup> is the unknown displacements vector of these nodes. It is worth noting that in Equation (70), the slave contact forces and their contribution to the master nodes have been written in terms of the local referential. In Equation (69), matrix **Q** and **H** have the following expressions:

$$ \mathbf{Q} = \begin{bmatrix} n_x & g_x \\ n_y & g_y \end{bmatrix}, \quad \mathbf{H} = \begin{bmatrix} -(1 - \bar{\xi})n_x & -(1 - \bar{\xi})g_x - n_x \frac{d_n}{l} \\ -(1 - \bar{\xi})n_y & -(1 - \bar{\xi})g_y - n_y \frac{d_n}{l} \\ -\bar{\xi}n_x & -\bar{\xi}g_x + n_x \frac{d_n}{l} \\ -\bar{\xi}n_y & -\bar{\xi}g_y + n_y \frac{d_n}{l} \end{bmatrix} \quad (71) $$

These transformation matrices are constructed for each contact elements, and their analytical inverse is given by

$$ \begin{bmatrix} \mathbf{Q} & \mathbf{0} \\ \mathbf{H} & \mathbf{I} \end{bmatrix}^{-1} = \begin{bmatrix} \mathbf{Q}^T & \mathbf{0} \\ -\mathbf{H}\mathbf{Q}^T & \mathbf{I} \end{bmatrix} \quad (72) $$
Substituting from Equation (72), system (69) leads to the following one for each contact element:

$$ \begin{bmatrix} \mathbf{Q}^{\mathrm{T}}\mathbf{K}^{\mathrm{s}s} & \mathbf{Q}^{\mathrm{T}}\mathbf{K}^{\mathrm{s}m} & \mathbf{Q}^{\mathrm{T}}\mathbf{K}^{\mathrm{l}1} \\ -\mathbf{H}\mathbf{Q}^{\mathrm{T}}\mathbf{K}^{\mathrm{s}s} + \mathbf{K}^{\mathrm{m}s} & -\mathbf{H}\mathbf{Q}^{\mathrm{T}}\mathbf{K}^{\mathrm{s}m} + \mathbf{K}^{\mathrm{m}m} & -\mathbf{H}\mathbf{Q}^{\mathrm{T}}\mathbf{K}^{\mathrm{l}2} + \mathbf{K}^{\mathrm{m}m} \end{bmatrix} \begin{bmatrix} \mathbf{du}^{(s)}{}^{j+1} \\ \mathbf{du}^{(m)}{}^{j+1} \end{bmatrix} = \begin{bmatrix} \mathbf{Q}^{\mathrm{T}}\mathbf{F}^{(s)} \\ -\mathbf{H}\mathbf{Q}^{\mathrm{T}}\mathbf{F}^{(s)} + \mathbf{F}^{(m)} \end{bmatrix}^j + \begin{bmatrix} \mathbf{R}^{(s)}{}^{j+1} \\ \mathbf{R}^{(m)}{}^{j+1} \end{bmatrix} \quad (73) $$

In Equation (73), the term $-\mathbf{H}\mathbf{Q}^{\mathrm{T}}\mathbf{F}^{(s)}$ corresponds to the effects of the contact on the master nodes. As a matter of fact, in considering the contact system outlined above, the discretization of the virtual contact work leads to additional external forces being exerted on the master nodes due to the action/reaction principle.

The contact kinematic unknowns with which the unilateral constraints and the Coulomb's friction law are written are the normal distance $d_n$ and the tangential velocity $v_t$. Additional work therefore has to be carried out in order to introduce these unknowns into the system to be solved.

The linearization of $d_n$ used in Equation (63) can also be written as follows:

$$ d_n^{j+1} = d_n^j + (\mathbf{N}_1^j)^\mathrm{T} \, \mathbf{du}^{j+1} \quad (74) $$

The linearized expression at time $t_k$ of the tangential relative velocity $v_t$ is chosen so that it is consistent with the time discretization algorithm used here:

$$ v_t^{j+1} = v_t^j + (\nabla_u^j v_t)^\mathrm{T} \, \mathbf{du}^{j+1} \quad (75) $$

with

$$ v_t^j = l^j \dot{\bar{\xi}}^j \quad (76) $$

where superscript $j$ refers to Newton's iteration. For the sake of simplification, index $k$, which refers to the time, has been omitted and will not be used. $\nabla_u^j v_t$ corresponds to the variation of the tangential velocity at time $t_k$. Using the forward Euler method to integrate $\dot{\bar{\xi}}^j$:

$$ \begin{gathered} \Delta t (\nabla_u^j v_t)^\mathrm{T} \, \mathbf{du}^{j+1} = \Delta \bar{\xi}^j (\nabla_u^j l)^\mathrm{T} \, \mathbf{du}^{j+1} + l^j (\nabla_u^j \bar{\xi})^\mathrm{T} \, \mathbf{du}^{j+1} \quad (77) \\ \text{with } \Delta \bar{\xi}^j = \bar{\xi}^j - \bar{\xi}_{k-1} \quad \text{and} \quad \Delta t = t_k - t_{k-1} \end{gathered} $$

Substituting Equation (77) into Equation (75) leads to

$$ v_t^{j+1} = v_t^j + \left[ \Delta \bar{\xi}^j \mathbf{T}_0^{e'j} + \mathbf{T}_1^{e'j} + \frac{d_n^j \mathbf{N}_0^{e'j}}{lj} \right]^\mathrm{T} \frac{\mathbf{du}^{j+1}}{\Delta t} \quad (78) $$
Thus, for one contact element, we obtain the following system in which all the variables are expressed with components in the local referential:

$$ \begin{bmatrix} d_n^{j+1} \\ \Delta t \cdot v_t^{j+1} \end{bmatrix} = \begin{bmatrix} d_n^j \\ \Delta t \cdot v_t^j \end{bmatrix} + \begin{bmatrix} du_n^{(s)j+1} \\ du_g^{(s)j+1} \end{bmatrix} \qquad (79) $$

The transfer of the displacement vector from the global coordinate system (x, y) to the local one (n, g) is carried out by

$$ \begin{bmatrix} du_n^{(s)j+1} \\ du_g^{(s)j+1} \end{bmatrix} = [\mathbf{Q}^T \ \Delta \mathbf{H}^T]^j \begin{bmatrix} du^{(s)j+1} \\ du^{(m)j+1} \end{bmatrix}_{(x,y)} \qquad (80) $$

where matrix $\Delta \mathbf{H}$ is defined by

$$ \Delta \mathbf{H} = \begin{bmatrix} -(1 - \bar{\xi})n_x & -(1 - \bar{\xi})g_x - n_x \frac{d_n}{l} - \Delta \bar{\xi} g_x \\ -(1 - \bar{\xi})n_y & -(1 - \bar{\xi})g_y - n_y \frac{d_n}{l} - \Delta \bar{\xi} g_y \\ \bar{\xi} n_x & -\bar{\xi} g_x - n_x \frac{d_n}{l} - \Delta \bar{\xi} g_x \\ \bar{\xi} n_y & -\bar{\xi} g_y - n_y \frac{d_n}{l} - \Delta \bar{\xi} g_y \end{bmatrix} \qquad (81) $$

In the standard slave/master approach, these matrices are constructed and reversed analytically for each contact element, so that it is not necessary to construct global matrices and to invert them using numerical methods.

### 6.2.2. Symmetrical slave/master formulation.
In this case, each node on the contact surfaces can simultaneously play the role of a slave and a master node, except for the nodes on the contact surface which have no opposite master segment. The latter nodes will therefore be treated only as master ones. We will use subscript (sm) to indicate the nodes which can be both slave and master nodes and subscript (mm) for nodes which play only the role of master nodes. The total force applied to a (sm) node is obtained by adding the contact force: this node is a slave node submitted to the forces resulting from the action of the other body, which is taken to be the master one. The contact force to which a (mm) node is subjected can obviously only be due to the action of the other body, which is in contact with the segment of the master node.

The local contact element matrices can then be constructed in a similar way to what was done in the previous standard slave/master case. However, due to the coupled action of (sm) nodes, it is not possible here to determine the local inverse of these matrices. These matrices are therefore assembled so as to give two global systems. The one relates to the change of contact variables, and other relates to the contact forces. These two systems can be schematically represented in the following forms for both the kinematic contact unknowns and the contact forces:
* Kinematic contact unknowns: The system used to express the components of the displacement vectors of the contact nodes in the global coordinate systems as a function of the components in the local coordinate systems cannot be written for a single contact element as was done in (80). These relations are related to each other by the following 2nc × 2nc system:

$$ \begin{bmatrix} du_{(n,g)}^{1(sm)} \\ du_{(n,g)}^{1(mm)} \\ du_{(n,g)}^{2(sm)} \\ du_{(n,g)}^{2(mm)} \end{bmatrix} = \begin{bmatrix} Q_{sm}^{1\mathrm{T}} & 0 & \Delta \tilde{H}_{sm}^{2\mathrm{T}} & \Delta \tilde{H}_{mm}^{2\mathrm{T}} \\ 0 & Q_{mm}^{1\mathrm{T}} & 0 & 0 \\ \Delta \tilde{H}_{sm}^{1\mathrm{T}} & \Delta \tilde{H}_{mm}^{1\mathrm{T}} & Q_{sm}^{2\mathrm{T}} & 0 \\ 0 & 0 & 0 & Q_{sm}^{2\mathrm{T}} \end{bmatrix} \begin{bmatrix} du_{(x,y)}^{1(sm)} \\ du_{(x,y)}^{1(mm)} \\ du_{(x,y)}^{2(sm)} \\ du_{(x,y)}^{2(mm)} \end{bmatrix} \quad (82) $$

Once the system (82) has been numerically inverted, the kinematic contact variables can be straightforwardly introduced for the (mm) nodes in the global system, using the analog of relation (79).

* Reactions:

$$ \begin{bmatrix} R_{(x,y)}^{1(sm)} \\ R_{(x,y)}^{1(mm)} \\ R_{(x,y)}^{2(sm)} \\ R_{(x,y)}^{2(mm)} \end{bmatrix} = \begin{bmatrix} \theta Q_{sm}^{1} & 0 & (1-\theta)H_{sm}^{2} & 0 \\ 0 & \theta Q_{mm}^{1} & (1-\theta)H_{sm}^{2} & 0 \\ \theta H_{sm}^{1} & 0 & (1-\theta)Q_{sm}^{2} & 0 \\ \theta H_{sm}^{1} & 0 & 0 & (1-\theta)Q_{mm}^{2} \end{bmatrix} \begin{bmatrix} R_{(n,g)}^{1(sm)} \\ R_{(n,g)}^{1(mm)} \\ R_{(n,g)}^{2(sm)} \\ R_{(n,g)}^{2(mm)} \end{bmatrix} \quad (83) $$

The constraints $R_{ng}^{1(mm)} = R_{ng}^{2(mm)} = 0$ are necessary to be able to write a whole complementarity problem. To deal with these constraints, a new change of variable will be made. A sufficiently large positive number will be added to the kinematic contact unknowns which will ensure the strict positivity of the kinematic variables and consequently, the above constraints will be satisfied, due to the complementarity condition.

## 6.3. Condensation procedure

At each equilibrium equation, the system to be solved can be reformulated in terms of contact components (c) and free components (f); one obtains

$$ \begin{bmatrix} \hat{\mathbf{K}}_{cc} & \hat{\mathbf{K}}_{cf} \\ \hat{\mathbf{K}}_{fc} & \hat{\mathbf{K}}_{ff} \end{bmatrix} \begin{bmatrix} \mathrm{d}\mathbf{u}_c \\ \mathrm{d}\mathbf{u}_f \end{bmatrix} = \begin{bmatrix} \mathbf{F}_c \\ \mathbf{F}_f \end{bmatrix} + \begin{bmatrix} \mathbf{R}_c \\ \mathbf{0} \end{bmatrix} \quad (84) $$

With $\mathbf{F} = \mathbf{F}_{\text{ext-cont}}(\mathbf{u}^j) - \mathbf{F}_{\text{int}}(\mathbf{u}^j)$ and $\mathbf{R}_c = [\mathbf{H}_n(\mathbf{u}^j)\mathbf{R}_n^{j+1} - \mathbf{H}_g(\mathbf{u}^j)\mathbf{R}_g^{j+1}]^T$. The condensation procedure leads to:

$$ \begin{gathered} \hat{\mathbf{K}}^{\star} \mathrm{d}\mathbf{u}_c = \mathbf{F}^{\star} + \mathbf{R}_c \quad (85) \\ \mathrm{d}\mathbf{u}_f = \hat{\mathbf{K}}_{ff}^{-1} \mathbf{F}_f - \hat{\mathbf{K}}_{ff}^{-1} \hat{\mathbf{K}}_{fc} \mathrm{d}\mathbf{u}_c \end{gathered} $$
<figure><img src="image_4.png" /><figcaption>Figure 3. Symmetrical test—geometry.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. Standard approach.</figcaption></figure>

with

$$ \hat{\mathbf{K}}^* = \hat{\mathbf{K}}_{cc} - \hat{\mathbf{K}}_{cf} \hat{\mathbf{K}}_{ff}^{-1} \hat{\mathbf{K}}_{fc} $$

$$ \mathbf{F}^* = \mathbf{F}_c - \hat{\mathbf{K}}_{cf} \mathbf{K}_{ff}^{-1} \mathbf{F}_f $$

Two connected systems have been described. The first one deals only with the contact nodes. The second system deals with the nodes which are not involved in the contact. This is a non-constrained problem in which the only unknown vector is the displacement one. Its solution, which obviously depends on the solution of the previous system, can be obtained using a more standard algorithm.

# 7. NUMERICAL APPLICATIONS

## 7.1. Application tests

### 7.1.1. Symmetrical test.
The standard master/slave approach can lead to large errors in the contact kinematics. A very simple example is presented below to illustrate this. In this example, the geometry of which is given in Figure 3, both bodies have the same material properties, are discretized with an identical mesh and are subjected to symmetrical loads. In the standard approach, only one of these bodies, the upper body here, is taken to be the master one, and the lower one is the slave body.

Figure 4 shows that the contact constraints are satisfied: none of the slave nodes penetrate into the master body. However, as penetration of the master body into the slave one is not checked here, the upper left-hand side node of the master surface is located inside the slave body. The results shown in Figure 5 were obtained with the symmetrical formulation. The contact was checked simultaneously and there was obviously no penetration of one body into the other one.
<figure><img src="image_6.png" /><figcaption>Figure 5. Symmetrical approach.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6. Constant pressure patch test.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Mesh A</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Mesh B</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Mesh C</figcaption></figure>

<caption>Figure 7. Finite element meshes for patch test.</caption>

### 7.1.2. Patch test.
Here we consider the constant pressure patch test for contact problems presented in Reference [23]. A rectangular elastic punch is placed in frictionless contact with an elastic body by means of a uniformly distributed constant pressure (Figure 6). Since the two bodies have the same material properties, the correct solution corresponds to a constant contact pressure on the contact area. The meshes A, B, and C given in Figure 7 presented in Reference [23] were tested. With mesh A, constant contact pressure was found with the single- and two-pass algorithm, whereas for meshes B and C the standard single-pass master/slave algorithm failed. Our $\theta$-pass master/slave algorithm ($\theta \neq 0$) makes it possible to determine in each case a constant pressure in the contact area. In Figure 8, the contact forces obtained using mesh C are presented.

### 7.1.3. Large sliding test.
With the symmetrical formulation, it is possible to improve the treatment of large sliding problems. From the numerical point of view, to deal with face-to-face nodes in the symmetrical formulation, we introduced node neighbourhoods which are sufficiently small domains around the nodes. If the neighbourhood of a node contains a node belonging to the antagonist surface, then one of these nodes will be taken to be the slave one, and the other one will be the master one. Figure 9 shows the ability of the algorithms to deal with cases of this kind, which arise particularly when large sliding movements occur. In this example, the upper body slides onto the lower one, and obviously during the sliding, face-to-face nodes can occur.
<figure><img src="image_11.png" /><figcaption>Figure 8. Patch test: mesh c, contact forces.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 9. Large sliding test.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 10. Geometry of the extrusion test.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 11. Extrusion results—frictionless and friction cases.</figcaption></figure>

## 7.2. Extrusion

Here we consider the extrusion of an aluminium cylinder into an aluminium die. This example, the geometry of which is given in Figure 10, was first mentioned by Laursen [15], who compares the results with and without friction. In this axisymmetric problem, the material properties of both bodies were taken to be $E = 68,956$ MPa, $v = 0.32$. Zero tangential displacements are prescribed on the right-hand side of the piece, as shown in Figure 10. The isotropic hardening is given by a Swift law: $\sigma_0(\bar{e}_p) = 261.2(0.11868 + (\bar{e}_p))$. Both bodies were discretized with Q4/P0 elements. Prescribed displacements were applied at the top of the cylinder. The cylinder was forced into the hole until its top had covered a distance equal to its original length. The final deformed configurations have been plotted in both cases in Figure 11. In the friction case, the cylinder obviously does not penetrate as far into the block as in the frictionless case. In the frictionless case, the reduction in the length of the piece is
<figure><img src="image_15.png" /><figcaption>Figure 12. Asperity shearing—geometry.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 13. Asperity shearing—contact forces.</figcaption></figure>

small. The arrows in Figure 11 indicate the contact nodal forces exerted at a contact node: their orientation is given by that of the contact force, and their length is proportional to the intensity of the contact force. Especially, in the presence of large sliding movements, they depend on the local characteristics of the contact system, namely on the difference between the length of the antagonist segments, the number of nodes and the role played by each of the nodes (master, slave or both) involved. In this work, as contact conditions were written as nodal reactions, contact forces on each contact nodes were plotted in Figure 11 rather than a distribution of these forces on each contact segment. The irregular distribution of these nodal contact forces which can be observed in Figure 11 is caused by the local conditions (contact and geometry). One can see in Figure 11, the effects of the friction on the orientation and the intensity of the contact forces. From the computational point of view, convergence is reached faster in the frictionless case than in the friction one. In the latter case, due to the friction, the plastic strain and the stresses are greater and convergence requires a smaller step size. Our results are very similar to those presented by Laursen. This example illustrates the ability of these methods to accurately describe a severe friction contact problem. As a matter of fact with friction, the sliding of nodes at the corner of the block is more delicate to treat and it can lead to numerical problems, depending on the method used.

### 7.3. Asperity shear off

In this last example, we consider two asperities with identical geometries face to face moving in opposite directions. This part deals with the shear effects located in the contact zone. The contact is assumed to be frictionless. This is a preliminary model for characterizing the friction as a plastic shear process between asperities. The geometry of the test is given in Figure 12. Two different materials are involved. The one was an aluminium with the following material properties: $E = 70,000$ MPa, Poisson's coefficient $\nu = 0.35$. Its hardening law is given by: $\sigma_0(\bar{\epsilon}^p) = 540(0.008057 + \bar{\epsilon}^p)^{0.28}$. The second material is a steel with the following characteristics: $E = 205,900$ MPa, $\nu = 0.3$, and its hardening law is given by:
<figure><img src="image_17.png" /><figcaption>Figure 14. Asperity shearing—friction coefficient.</figcaption></figure>

$\sigma_0(\bar{\epsilon}^p) = 619(0.0078 + \bar{\epsilon}^p)^{0.235}$. The two bodies are meshed with 586 Q4/P0 elements. The height of the asperities is 4.75 µm and their characteristic angle is $\theta = 13.36^\circ$. The meshing between the two asperities is 2.325 µm. The lower aluminium asperity remains fixed. Zero tangential displacements are prescribed on sides AB and CD, and zero vertical displacements are prescribed on BC. The upper steel asperity is moved with prescribed tangential displacements from the left to the right hand side.

Figure 13 shows the evolution of the contact forces with the motion of the asperity and Figure 14 shows the changes in the ratio $F_N/F_T$, where $F_N$ stands for the global normal contact force to which the lower asperity is subjected and $F_T$ stands for the tangential force. Assuming the problem to be a frictionless one, this ratio characterizes an apparent friction coefficient. As can be observed in Figure 14, at the start of the motion, $F_N/F_T$ is equal to the tangent to the characteristic angle of the asperities. As soon as plastic deformations occur, this coefficient decreases, reaching its minimum value when the asperity has been completely sheared off. During the process, a plastic area appears in the contact zones and increases untill all the asperities have become plastic. Lastly, Figure 15 shows various steps in the shearing process. Each arrow indicates a contact point and its length is proportional to the intensity of the contact force. Computations were carried out until the lower asperity was completely sheared. This illustrates the ability of the algorithm to deal with large changes in the contact state from a state with high contact forces to one with very low normal forces.

## 8. CONCLUSION

A symmetrical slave/master formulation for friction contact was developed, and the discrete problem arising from this formulation was set as a complementarity one. This made it possible to use Lemke's method giving the values of the two contact unknowns (distance and tangential velocity) without having to use a penalty coefficient, as required in the standard penalty and augmented Lagrangian methods. Examples are given to illustrate the efficiency of the method. In other studies, we used a formulation of this kind to study the lubricated contact between
<figure><img src="image_18.png" /><figcaption>Figure 15. Asperity shearing—deformed mesh.</figcaption></figure>

a sheet-metal piece and tools. The main advantage of Lemke's method associated with the symmetrical formulation in comparison with more standard methods is that it gives access to the contact state without any approximation, which is of great interest in the case of severe contact situations such as those involved in microscopic surface interactions. As Lemke's method is a pivoting method, it may be time-consuming when there are large numbers of contact nodes. In the examples investigated in our studies, this can occur particularly during the first increments. An adaptive strategy switching between Lemke's method and classical regularized methods depending on the number of contact nodes is studied. In sheet-metal-forming process problems, such as the analysis of a blankholder fitted with a drawbead, Lemke's method gives a good description of the contact state, particularly in the small radius parts of the tools [24]. Various ways of improving the time required by the method are studied, such as the use of predictive techniques taking the previous contact state into account. The $\theta$-pass approach introduced in this paper means that the standard slave/master approach can be associated with a symmetrical one in the same algorithm. Advantages of this formulation are to avoid specific treatments for each case, to avoid the *a priori* choice of the master body and the slave one and, if $\theta \neq 0$ and 1, to allow the determination of the nodal contact forces on each body without any additional computational costs. Algorithm was checked on the patch test and for the large sliding test we verified that results do not depend on the value of $\theta$. In solving other numerical examples in the symmetrical case, when $\theta = \frac{1}{2}$, the validity of the results was checked. The local contact conditions were always satisfied for each contact surface and in posttreatment steps, we checked that Coulomb's sliding conditions were satisfied on every contact node. For any values of $\theta \in \{0, \frac{1}{2}, 1\}$ no significant differences in computational times were observed. Until now, none of the tests carried out led to significant lower computational times with $\theta \notin \{0, \frac{1}{2}, 1\}$ as could be expected when, for example, there are strong differences
between meshes or strong differences between the mechanical properties of each bodies. The subject of this work was to develop the symmetrical complementarity approach presented in the present paper and to show how a such approach could be efficient. Different treatments of certain technical points in the software, when using other values of $\theta$, could perhaps lead to significant differences in terms of computational times.

REFERENCES

1. Raous M, Chabrand P, Lebon F. Numerical methods for frictional contact problems and applications. *Journal of Theoretical and Applied Mechanics* 1988; **7**(1):111-128.

2. Raous M, Barbarin S. Conjugate gradient for frictional contact. In *Contact Mechanics*, Curnier A (ed.). Presses polytechniques et Universitaires Romanes, 1992; 423-432.

3. Bathe KJ, Chaudhary A. A solution method for planar and axisymmetric contact problem. *International Journal for Numerical Methods in Engineering* 1985; **21**:65-88.

4. Wriggers P, Simo JC, Taylor RL. Penalty and augmented Lagrangian formulations for contact problems. In *Proceedings of NUMETA Conference*, Mlddleton J, Pande GN (eds). Balkema: Rotterdam, 1985.

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

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

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

8. Zhong ZH. *Finite Element Procedures for Contact-Impact Problems*. Oxford University Press: Oxford, 1993.

9. Klarbring A, Bjorkman G. A mathematical programming approach to contact problem with friction and varying contact surface. *Computers and Structures* 1988; **30**:1185-1198.

10. Hallquist JO, Goudreau GL, Benson DJ. Sliding interfaces with contact-impact in large-scale lagrangian computations. *Computer Methods in Applied Mechanics and Engineering* 1985; **51**:107-137.

11. Wriggers P, Simo JC. A note on tangent stiffness for fully nonlinear contact problems. *Communications in Applied Numerical Methods* 1985; **1**:199-203.

12. Klarbring A. Large displacement frictional contact: a continuum framework for finite element discretization. *European Journal of Mechanics, A/Solids* 1995; **14**(2):237-253.

13. Laursen TA, Simo JC. A continuum-based finite element formulation for the implicit solution of multibody large deformation frictional contact problems. *International Journal for Numerical Methods in Engineering* 1993; **36**: 3451-3485.

14. Wriggers P, Miehe C. Contact constraints within coupled thermomechanical analysis—a finite element model. *Computer Methods in Applied Mechanics and Engineering* 1994; **113**:301-319.

15. Laursen T. Formulation and treatment of frictional contact problems using finite elements. Ph.D. Thesis, Stanford University, 1992.

16. Zavarise G, Wriggers P. A segment-to-segment contact strategy. *Mathematical and Computer Modelling* 1998; **28**:497-515.

17. Pavanachand Ch, Krishnakumar R. A new one-pass approach for large deformation multibody frictional contact analysis. *Communications in Numerical Methods in Engineering* 1996; **12**:569-579.

18. Cescotto S. Contact quasi-static. In *Modelisation mécanique et numérique du contact et du frottement*. IPSI, 1996.

19. Klarbring A, Bjorkman G. Solution of large displacement contact problems with friction using Newton's method for generalized equations. *International Journal for Numerical Methods in Engineering* 1992; **34**:97-116.

20. Robinson SM. Generalized equations. In *Mathematical Programming: the State of the Art*, Bachem et al. (eds). Springer: Berlin, 1983; 346-367.

21. Josephy NH. Newton's method for generalized equations. And the pies energy model. Ph.D. Thesis, University of Wisconsin, Mathematics Research Center, 1979.

22. Chabrand P, Dubois F, Raous M. Various numerical methods for solving unilateral contact problems with friction. *Mathematical and Computer Modelling* 1997; **28**(4-8):97-108.

23. Taylor RL, Papadopoulos P. On a patch test for contact problems in two dimensions. In *Nonlinear Computational Mechanics: State of the Art*, Wriggers P, Wagner W (eds). Springer: Berlin, 1991.

24. Chabrand P, Dubois F, Gelin JC. Modelling drawbeads in sheet metal forming. *International Journal of Mechanical Sciences* 1996; **38**(1):59-77.

25. Chaudhary A, Bathe KJ. A solution method for static and dynamic analysis of three-dimensional contact problems with friction. *Computers and Structures* 1986; **24**:855-873.
26. Chabrand P, Chertier O. Prise en compte du frottement sur les petits rayons: comparaison numérique expérimental. loi de frottement à coefficient variables. *Rapport de contrat* H5.12.51, CNRS/Renault, 1994.

27. Chabrand P. Modélisation d'effets locaux de surface en emboutissage. In *Tribologie et Ingénierie des surfaces*, Sirpe (ed.). 1996.

28. Chertier O. Contact et frottement en grandes déformations élastoplastiques entre solides déformables. Ph.D. Thesis, Université de la Méditerranée, 1997.

29. He Q-C, Curnier A. Anisotropic dry friction between two orthotropic surfaces undergoing large displacements. *European Journal of Mechanic, A/Solids* 1993; **12**(5):631–666.

30. Dubois F. Contact, frottement, grandes déformations élastoplastiques. Application à l'emboutissage. Ph.D. Thesis, Université de la Méditerranée, 1994.

31. Emmens WC, Montford G. The influence of process conditions and surface characteristics on friction at low pressure. In *Proceedings of the 3rd ICTP*, 1990.

32. Ike H, Makinouchi A. Effect of lateral tension and compression on plane strain flattening asperities lying over a plastically deformable bulk. *Wear* 1990; **140**:17–38.

33. Jeuzette JP, Sonzogni V. A projected conjugate gradient method for structural stability analysis with linear constraints. *Computers and Structures* 1989; **33**(1).

34. Martinet F. Formulation Lagrangienne/Eulérienne pour un problème couplé de fluide structure en emboutissage. Ph.D. Thesis, Université de la Méditerranée, 1998.

35. May HO. The conjugate gradient method for unilateral problems. *Computers and Structures* 1988; **29**(4).