# A weighted residual relationship for the contact problem with Coulomb friction

Anh Le Van, Huynh Tan Tai Nguyen

► To cite this version:

Anh Le Van, Huynh Tan Tai Nguyen. A weighted residual relationship for the contact problem with Coulomb friction. Computers & Structures, 2009, 87 (23-24), pp.1580-1601. 10.1016/j.compstruc.2009.08.013. hal-01004969

HAL Id: hal-01004969

https://hal.science/hal-01004969

Submitted on 5 Oct 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 4.0 International License
# A weighted residual relationship for the contact problem with Coulomb friction

A. Le van *, T.T.H. Nguyen

GeM (Laboratory of Civil and Mechanical Engineering), Faculty of Science, University of Nantes, 2, rue de la Houssinière, BP 92208, 44322 Nantes Cedex 3, France

In this work, a weighted residual relationship is proposed as an extension of the standard virtual work principle to deal with the large deformation contact problem with Coulomb friction. This weak form is a mixed relationship involving the displacements and the multipliers defined on the reference contact surface of the contactor and is shown to be equivalent to the strong form of the initial/boundary value contact problem. The discretization in space by means of the finite element method is carried out on the mixed relationship in a simple way in order to obtain the semi-discrete equation system. The contact tangent stiffness is derived and numerical examples are presented to assess the efficiency of the formulation.

## 1. Introduction

Problems where contact takes place between two or more bodies are widespread and of practical significance in mechanics. Their solution is complex because of the nonlinear and non-smooth laws governing the interface response between the contacting bodies. In the context of large deformations and large slips, the problem becomes even more nonlinear and robust formulations must be designed with a view to efficiently solve it.

Over more than three decades, a great number of formulations and solution methods have been proposed in the literature in order to deal with different kinds of contact problems encountered, from frictionless contact in linear elastostatics up to frictional contact in large deformation with complex material constitutive and interfacial contact laws, either in statics or dynamics. The interested reader is referred to the book of Kikuchi and Oden [1] and the more recent books of Wriggers [2] and Laursen [3] for a review of the numerous references related to the existing methods for the treatment of the contact constraints, in particular the classical Lagrange multiplier, the penalty and the augmented Lagrangian methods. As the formulation proposed in the present work is inspired from methods of the augmented Lagrangian type, we shall limit ourselves to a brief review of these methods in the literature.

The augmented Lagrangian method is often chosen to cope with the contact inequality constraints inasmuch as it combines the regularizing effect of the penalty method and the exact satisfaction of the constraints ensured by the Lagrange multiplier method, without having the ill-conditioning problem inherent to the penalty method. It was applied to different contact problems in many vari-

ants. Kinematically linear problems with frictional contact were solved using the augmented Lagrangian approach by Alart and Curnier [4], Zavarise et al. [5]. Large deformation contact without friction was also successfully addressed by the same method in the early work of Glowinski and Le Tallec [6] considering a plane rigid obstacle, in Heegard and Curnier's paper [7] considering contact between two elastic bodies, and later in Zavarise and Wriggers [8]. The solution of large deformation frictional contact problems between deformable bodies was accomplished by Simo and Laursen [9,10] including the algorithmic symmetrization of the contact tangent stiffness matrix, Wriggers [11], Pietrzak and Curnier [12,13] and more recently by Feng et al. [14] for Blatz-Ko hyperelastic bodies, and Mijar and Arora [15,16]. All previous works considered the Coulomb friction model or an extension of it to incorporate micro-mechanical effects for instance.

In most of the formulations, the tangential contact stresses were computed by locally integrating the contact laws and the contact stiffness matrix is consistently derived from these stresses, following the spirit of the return mapping scheme in elastoplasticity as described by Giannakopoulos in [17]. In some formulations, e.g., [5], the tangential contact law adopted was not of the rate type and the above-mentioned local integration was not necessary. Curnier and co-workers [4,7,13] chose another approach by directly incorporating the contact laws into their augmented Lagrangians, which enabled them to obtain the correct tangential contact stresses without resorting to the local integration.

As regards the global solution of the discrete equation system, Uzawa's method was usually chosen because of its simplicity: the multipliers are updated in an outer augmentation loop and the displacements computed in a standard way inside an inner loop while keeping the multipliers fixed. Otherwise, the discrete equation system can be solved simultaneously for the displacements
and the multipliers, as was done in [13,8], for instance. One can also consult several different global algorithms for solving augmented Lagrangian problems in [6, Chapters 3 and 7].

## 1.1. Motivation

The purpose of this paper is to propose a simple formulation for the contact problem with Coulomb friction which combines the benefits of the aforementioned formulations of the augmented Lagrangian type. The sought formulation is based on a new weak form – similar to the virtual work principle – which should fulfill the following conditions:

(i) The weak form should incorporate the contact laws. More precisely, although it is an equality, it must imply the inequality constraints specific to contact laws. The consequence of this is that the contact laws are automatically satisfied regardless of the constitutive laws of the contacting bodies and no local integration is needed at the contact surface level.

(ii) The weak form may look like an augmented Lagrangian principle in that it involves both the displacements and the multipliers, but, more importantly, it should be stated in the form of a weighted residual relationship, with arbitrary (regular) virtual displacements and multipliers (an overview of the weighted residual methods can be found in [18]). As such, the weak form can be discretized by means of the finite element method in a straightforward way.

## 1.2. Outline of the paper

The paper is organized as follows. The strong form of the initial/boundary value contact problem is summarized in Section 2. In Section 3, a weighted residual relationship is proposed as an extension of the standard virtual work principle with the aim of encompassing the large deformation contact response with Coulomb friction. This mixed weak form which involves both the displacements and the multipliers defined on the reference contact surface of the contactor will be shown to be equivalent to the strong form of the contact problem. Section 4 deals with the semi-discretization in space of the weak form by means of the finite element method, which results in a coupled nonlinear semi-discrete equation system in two unknowns: the nodal displacements and the nodal multipliers. In Section 5, the contact tangent stiffness matrix which is necessary to the Newton-Raphson iterative solution is derived by the exact linearization of the weak form. Different numerical examples assessing the efficiency of the proposed formulation are presented in Section 6.

# 2. Strong form of the contact problem

Consider two bodies whose positions at any time are represented by some regions in a three-dimensional Euclidean affine space. The position at any time of any particle belonging to the bodies is represented by a point in this space. As usually done, a fixed origin being chosen, any point as well as the associated position vector will be identified to a three-component vector in $\mathbb{R}^3$.

We shall describe the mechanical problem using a Lagrangian description. The reference configuration of the two bodies are represented by the regions $\Omega_o^{(1)}$ and $\Omega_o^{(2)}$ in $\mathbb{R}^3$. Let us assume that in the time interval $[0,T]$ of interest, the two bodies may come into contact with each other through the motions denoted $\phi^{(1)}$ and $\phi^{(2)}$:

$$ \forall i \in \{1,2\}, \quad \forall t \in [0,T], $$

$$ \phi^{(i)}(., t) : \quad \Omega_o^{(i)} \ni \mathbf{X}^{(i)} \mapsto \mathbf{X}^{(i)} = \phi^{(i)}(\mathbf{X}^{(i)}, t), \qquad (1) $$

where $\mathbf{X}^{(i)}$ and $\mathbf{x}^{(i)}$ are the reference and current positions of any particle of body $i$, respectively.

The displacement field in body $i$ is defined as $\mathbf{U}^{(i)} = \mathbf{x}^{(i)} - \mathbf{X}^{(i)}$. The strong form of the problem is expressed by the local equations summarized below, which must be satisfied at any time $t \in [0, T]$ and at any point $\mathbf{X}^{(i)} \in \Omega_o^{(i)}, i \in \{1, 2\}$.

The linear momentum balance equation can be written for each body $i$ via:

$$ \mathbf{d}v\Pi^{(i)} + \rho_o^{(i)}\mathbf{f}^{(i)} = \rho_o^{(i)}\ddot{\mathbf{U}}^{(i)}, \qquad (2) $$

where $\Pi^{(i)}$ is the first Piola-Kirchhoff stress tensor, $\rho_o^{(i)}$ is the reference density of body $i$ and $\mathbf{f}^{(i)}$ the prescribed body force per unit mass.

The constitutive relationships relating stress to strain are left unspecified until the numerical examples are considered. This means that the weighted residual relationship proposed hereafter holds for all constitutive materials, as is the case with the standard virtual work principle.

The reference boundary $S_o^{(i)}$ of body $i$ is partitioned into three non-overlapping and pairwise disjoint parts denoted $S_{oU}^{(i)}, S_{oT}^{(i)}$ and $S_{oc}^{(i)}$, where $S_{oU}^{(i)}$ and $S_{oT}^{(i)}$ are the parts where the displacement and the tractions are prescribed, respectively, and $S_{oc}^{(i)}$ is the part where contact potentially takes place between the two bodies during some portion of the time interval $[0,T]$. The boundary conditions on $S_{oU}^{(i)}$ and $S_{oT}^{(i)}$ on each body $i$ can be summarized as

$$ \mathbf{U}^{(i)} = \bar{\mathbf{U}}^{(i)} \text{ in } S_{oU}^{(i)}, \quad \Pi^{(i)}\mathbf{N}^{(i)} = \bar{\mathbf{T}}^{(i)} \text{ in } S_{oT}^{(i)}, \qquad (3) $$

where $\bar{\mathbf{U}}^{(i)}, \bar{\mathbf{T}}^{(i)}$ are the prescribed displacement and nominal traction, respectively, $\mathbf{N}^{(i)}$ the unit outward normal at point $\mathbf{X}^{(i)}$ to surface $S_o^{(i)}$. The boundary conditions on the contact surface will be presented in detail in Section 2.3. The initial conditions for body $i$ are

$$ \forall \mathbf{X}^{(i)} \in \Omega_o^{(i)}, \quad \mathbf{U}^{(i)}(\mathbf{X}^{(i)}, 0) = \mathbf{0}, \quad \dot{\mathbf{U}}^{(i)}(\mathbf{X}^{(i)}, 0) = \bar{\mathbf{V}}_o^{(i)}(\mathbf{X}^{(i)}), \qquad (4) $$

where $\bar{\mathbf{V}}_o^{(i)}$ is the prescribed initial velocity field for body $i$.

## 2.1. Contact kinematics

Let us now focus on the equations governing the response on the contact interface, which are the essential ingredients in the large deformation contact problem. Let the spatial counterparts of $\Omega_o^{(i)}, S_o^{(i)}, S_{oU}^{(i)}, S_{oT}^{(i)}$ and $S_{oc}^{(i)}$ be denoted by $\Omega^{(i)}, S^{(i)}, S_U^{(i)}, S_T^{(i)}$ and $S_c^{(i)}$, respectively, e.g., $S_c^{(i)} = \phi^{(i)}(S_{oc}^{(i)}, t)$. In order to describe the relative motion of the two bodies, we arbitrarily choose one body (the target), say body 2, as the reference and to evaluate the proximity of the other (the contactor), thus body 1, with respect to it.

It is assumed that the reference surface $S_{oc}^{(i)}, i = 1 \text{ or } 2$, is parametrized by the mapping $\theta^{(i)}$ and the current surface $S_c^{(i)}$ at time $t$ by the mapping $\psi^{(i)}(., t) = \phi^{(i)}(., t) \circ \theta^{(i)}$:

$$ \begin{align} \text{A region of } \mathbb{R}^2 &\rightarrow S_{oc}^{(i)} \subset \mathbb{R}^3 \rightarrow S_c^{(i)} \subset \mathbb{R}^3, \nonumber \\ \xi^{(i)} &= (\xi^{(i)1}, \xi^{(i)2}) \quad \mapsto \quad \mathbf{X}^{(i)} = \theta^{(i)}(\xi^{(i)}) \tag{5} \\ &\quad \mapsto \quad \mathbf{x}^{(i)} = \phi^{(i)}(\mathbf{X}^{(i)}, t) = \psi^{(i)}(\xi^{(i)}, t). \nonumber \end{align} $$

Consider a particular point $\mathbf{x} \in S_c^{(1)}$ and define the opposite point $\mathbf{y} \in S_c^{(2)}$ as the closest point to $\mathbf{x}$ in the Euclidean sense (Fig. 1):

$$ \mathbf{y} = \arg \min_{\mathbf{x}^{(2)} \in S_c^{(2)}} \|\mathbf{x} - \mathbf{x}^{(2)}\|. \qquad (6) $$

In writing (6), we assume that point $\mathbf{y}$ is unique without discussing the case of non-uniqueness, which is not of major consequence in numerical computations. The interested readers may consult comprehensive comments on this issue in the literature, e.g., [7,2, p.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Contact kinematics.</figcaption></figure>

48], [3, p. 115], and references quoted therein. If the target surface $S_c^{(2)}$ is smooth at $\mathbf{y}$, then $\mathbf{y}$ is the projection of $\mathbf{x}$ on $S_c^{(2)}$. As may be seen, the notation $\mathbf{x}^{(i)}$, $i \in \{1, 2\}$, in (5) designates any point on $S_c^{(i)}$, whereas the notation $\mathbf{x}$ designates a particular point of interest on $S_c^{(1)}$ and $\mathbf{y}$ the opposite point of $\mathbf{x}$ on $S_c^{(2)}$ defined by (6).

The point $\mathbf{x} \in S_c^{(1)}$ being associated to (at least) one point $\mathbf{y} \in S_c^{(2)}$, the proximity $g$ is defined by

$$g = -\mathbf{v}(\mathbf{x} - \mathbf{y}), \qquad (7)$$

where $\mathbf{v}$ is the normal at point $\mathbf{y}$. The sign convention chosen in (7) implies that the proximity $g$ is positive when there is interpenetration and negative when there is a gap between the two bodies. Thus, symbol $g$ is the opposite of what is referred to as the gap in the literature. It can be checked that $g$ is an objective quantity.

Let us now convert all the spatial quantities defined above into material ones defined on the reference configurations, by means of the motions $\phi^{(1)}$ and $\phi^{(2)}$ in (1). Considering the point $\mathbf{X} \in S_{oc}^{(1)}$ related to the point $\mathbf{x}$ considered by $\mathbf{x} = \phi^{(1)}(\mathbf{X}, t)$, then the point $\mathbf{Y} \in S_{oc}^{(2)}$ related to point $\mathbf{y}$ by $\mathbf{y} = \phi^{(2)}(\mathbf{Y}, t)$ satisfies:

$$\mathbf{Y} = \arg \min_{\mathbf{x}^{(2)} \in S_{oc}^{(2)}} \left\| \phi^{(1)}(\mathbf{X}, t) - \phi^{(2)}(\mathbf{X}^{(2)}, t) \right\|. \qquad (8)$$

A point $\mathbf{Y}$ thus being associated to point $\mathbf{X}$, the proximity $g$ defined in (7) can be rewritten as

$$g = -\mathbf{v}(\phi^{(1)}(\mathbf{X}, t) - \phi^{(2)}(\mathbf{Y}, t)). \qquad (9)$$

Eventually, the inverse image of $\mathbf{X}$ (resp. $\mathbf{Y}$) under the parametrization $\theta^{(1)}$ (resp. $\theta^{(2)}$) will be denoted $\xi$ (resp. $\eta$). By making $\mathbf{x}^{(2)} = \psi^{(2)}(\xi^{(2)}, t)$ in (6), one finds that $\eta$ is characterized by

$$\eta = \arg \min_{\xi^{(2)}} \| \mathbf{x} - \psi^{(2)}(\xi^{(2)}, t) \|. \qquad (10)$$

For notational conveniences, the dependencies on $\phi^{(1)}$, $\phi^{(2)}$ and $\psi^{(2)}$ will be omitted in the material kinematics variables and only the dependency on $(\mathbf{X}, t)$ will be shown. Thus, for instance, one writes $\eta = \eta(\mathbf{X}, t)$ and $g = g(\mathbf{X}, t)$.

Assuming that all the mappings introduced are smooth enough, let us define the natural basis $(\mathbf{a}_1, \mathbf{a}_2)$ at point $\mathbf{y} \in S_c^{(2)}$ as

$$\mathbf{a}_1 = \frac{\partial \psi^{(2)}}{\partial \xi^{(2)1}}(\eta(\mathbf{X}, t), t), \quad \mathbf{a}_2 = \frac{\partial \psi^{(2)}}{\partial \xi^{(2)2}}(\eta(\mathbf{X}, t), t). \qquad (11)$$

The surface $S_c^{(2)}$ is parametrized conveniently so that the normal $\mathbf{n} = \mathbf{a}_1 \times \mathbf{a}_2 / \|\mathbf{a}_1 \times \mathbf{a}_2\|$ at $\mathbf{y}$ to $S_c^{(2)}$ is directed toward the contactor. Use will also be made of the dual basis $(\mathbf{a}^1, \mathbf{a}^2)$ lying in the tangent plane of the current target surface as follows:

$$\forall \alpha \in \{1, 2\}, \quad \mathbf{a}^{\alpha} = \alpha^{\alpha\beta} \mathbf{a}_{\beta}, \qquad (12)$$

where implicit summation from 1 to 2 is assumed over repeated Greek indices and the $2 \times 2$ matrix $[\alpha^{\alpha\beta}]$ is the inverse of matrix $[\alpha_{\alpha\beta}] = [\mathbf{a}_{\alpha} \cdot \mathbf{a}_{\beta}]$.

The contact kinematics with friction is described here by using the following objective relative velocity at point $\mathbf{X}$ and time $t$ [19,20]:

$$\mathbf{V} = \mathbf{V}^{(1)}(\mathbf{X}, t) - \mathbf{V}^{(2)}(\mathbf{Y}(\mathbf{X}, t), t) + g\dot{\mathbf{v}} = -\dot{g}\mathbf{v} + \dot{\eta}^{\alpha}\mathbf{a}_{\alpha}, \qquad (13)$$

where $\mathbf{V}^{(1)}(\mathbf{X}, t)$ and $\mathbf{V}^{(2)}(\mathbf{Y}, t)$ are the velocities of the particles $\mathbf{X}$ and $\mathbf{Y}(\mathbf{X}, t)$, respectively. The rates $\dot{g}$, $\dot{\eta}^{\alpha}$ and $\dot{v}$ are given by [2,3]:

$$\dot{g} = -[\mathbf{V}^{(1)}(\mathbf{X}, t) - \mathbf{V}^{(2)}(\mathbf{Y}(\mathbf{X}, t), t)]\mathbf{v}, \\ \dot{\eta}^{\alpha} = \bar{a}^{\alpha\beta}\mathbf{a}_{\beta} [\mathbf{V}^{(1)}(\mathbf{X}, t) - \mathbf{V}^{(2)}(\mathbf{Y}(\mathbf{X}, t), t)] - g\bar{a}^{\alpha\beta}\mathbf{v}\mathbf{V}_{\beta}^{(2)}(\eta, t), \qquad (14) \\ \dot{v} = -(\kappa_{\alpha\gamma}\bar{a}^{\gamma\beta}\mathbf{a}^{\alpha} \otimes \mathbf{a}_{\beta})(\mathbf{V}^{(1)} - \mathbf{V}^{(2)}) - (\bar{a}^{\alpha\beta}\mathbf{a}_{\beta} \otimes \mathbf{v})\mathbf{V}_{\alpha}^{(2)},$$

where matrix $[\bar{a}^{\alpha\beta}]$ is the inverse of matrix $[\bar{a}_{\alpha\beta}]$, with $\bar{a}_{\alpha\beta} = \bar{a}_{\beta\alpha} + g\kappa_{\alpha\beta}$, $\kappa_{\alpha\beta} = \mathbf{v}\psi_{\alpha\beta}^{(2)}$ is the curvature of the surface $\psi^{(2)}$ at $\eta$, the symbol $\otimes$ designates the tensor product defined as $(\mathbf{a} \otimes \mathbf{b}) \cdot \mathbf{c} = \mathbf{a}(\mathbf{b} \cdot \mathbf{c})$ for all vectors $\mathbf{a}, \mathbf{b}$ and $\mathbf{c}$. Note that tensor $\kappa_{\alpha\gamma}\bar{a}^{\gamma\beta}\mathbf{a}^{\alpha} = \bar{a}^{\alpha\gamma}\kappa_{\gamma\beta}\mathbf{a}_{\alpha} \otimes \mathbf{a}^{\beta}$ in (14c) is symmetric.

On the basis of (13) one can define (i) the relative normal velocity $\mathbf{V}_N$ as the projection of $\mathbf{V}$ on the normal vector $\mathbf{v}$, and (ii) the slip velocity or the relative tangential velocity $\mathbf{V}_T$ as the projection of $\mathbf{V}$ into the tangent plane at $\mathbf{y}$ to $S_{oc}^{(2)}$:

$$\mathbf{V}_N = \mathbf{V} \cdot \mathbf{v} = -\dot{g}\mathbf{v}, \quad \mathbf{V}_T = (\mathbf{I} - \mathbf{v} \otimes \mathbf{v}) \cdot \mathbf{V} = \dot{\eta}^{\alpha}\mathbf{a}_{\alpha}. \qquad (15)$$

It should be noted that the slip rate $\mathbf{V}_T$ associated to the material point $\mathbf{X} \in S_{oc}^{(1)}$ is resolved in terms of the spatial basis $(\mathbf{a}_1, \mathbf{a}_2)$ at point $\mathbf{y} \in S_c^{(2)}$. As done with other kinematical quantities, use will be made of the abbreviated notations $\mathbf{V} = \mathbf{V}(\mathbf{X}, t)$ and $\mathbf{V}_T = \mathbf{V}_T(\mathbf{X}, t)$, omitting the dependencies on motions $\phi^{(1)}, \phi^{(2)}$.

## 2.2. Contact tractions

Let the contact Cauchy traction vector $\mathbf{t}(\mathbf{x}, t)$ at a typical point $\mathbf{x} \in S_c^{(1)}$ be resolved in terms of the local basis $(\mathbf{a}^1, \mathbf{a}^2, \mathbf{v})$ at the projection point $\mathbf{y} \in S_c^{(2)}$ as follows:

$$\mathbf{t} = t_N \mathbf{v} - t_T \mathbf{v} = t_N \mathbf{v} - t_{T\alpha} \mathbf{a}^{\alpha}. \qquad (16)$$

The normal component $t_N$ of traction vector $\mathbf{t}$ in the direction of normal $\mathbf{v}$ is the contact pressure, positive in compression. The tangential component $t_T = t_{T\alpha} \mathbf{a}^{\alpha}$ is the opposite of the frictional traction. It can be verified that the pressure $t_N$ and the frictional traction $t_T$ are objective variables.

Let $\mathbf{T}(\mathbf{X}, t) = \mathbf{\Pi}^{(1)} \cdot \mathbf{N}$ be the nominal Piola-Kirchhoff traction vector at point $\mathbf{X} \in S_{oc}^{(1)}$ (N the normal vector at $\mathbf{X}$ to $S_{oc}^{(1)}$). Like the contact Cauchy traction, the nominal traction $\mathbf{T}$ is resolved in terms of the spatial basis at $\mathbf{y} \in S_c^{(2)}$ as follows:

$$\mathbf{T}(\mathbf{X}, t) = T_N(\mathbf{X}, t)\mathbf{v} - \mathbf{T}_T(\mathbf{X}, t) = T_N(\mathbf{X}, t)\mathbf{v} - T_{T\alpha}\mathbf{a}^{\alpha}. \qquad (17)$$

From relation $t \mathbf{d} S^{(1)} = \mathbf{T} \mathbf{d} S_o^{(1)} (\mathbf{d} S_o^{(1)})$ is a differential reference area, $\mathbf{d} S^{(1)}$ its spatial counterpart), one derives the following relationships between the spatial and material traction components:

$$t_N \mathbf{d} S^{(1)} = T_N \mathbf{d} S_o^{(1)}, \quad t_T \mathbf{d} S^{(1)} = \mathbf{T}_T \mathbf{d} S_o^{(1)}. \qquad (18)$$

## 2.3. Contact laws

As a means of describing the interfacial response at current time $t$, the contact laws are naturally expressed in the spatial form. However, one may use (18) in order to recast them into the mixed form involving the relative material velocity and the nominal contact tractions, which are more suitable for the numerical
implementation within the Lagrangian description framework. The normal contact law then reads:

$$ \forall t \in [0, T], \quad \forall \mathbf{X} \in S_{oc}^{(1)}, \\ g(\mathbf{X}, t) \le 0 \quad \begin{cases} \cdot & \text{if } g(\mathbf{X}, t) < 0, \text{ then } T_N(\mathbf{X}, t) = 0 \\ \cdot & \text{if } g(\mathbf{X}, t) = 0, \text{ then } T_N(\mathbf{X}, t) \ge 0 \end{cases} \tag{19} $$

In this work, the tangential contact response is modeled by the Coulomb friction law:

$$ \forall t \in [0, T], \quad \forall \mathbf{X} \in S_{oc}^{(1)}, \\ (a) g(\mathbf{X}, t) < 0 \Rightarrow \mathbf{T}_T(\mathbf{X}, t) = \mathbf{0} \text{ (no relationship for } \mathbf{V}_T(\mathbf{X}, t)), \\ (b) g(\mathbf{X}, t) = 0 \Rightarrow \|\mathbf{T}_T(\mathbf{X}, t)\| \le \mu T_N(\mathbf{X}, t) \tag{20} $$

$$ \text{where} \quad \begin{cases} \cdot & \text{if } \|\mathbf{T}_T\| < \mu T_N, \text{ then } \mathbf{V}_T = \mathbf{0} \text{ (stick)} \\ \cdot & \text{if } \|\mathbf{T}_T\| = \mu T_N, \text{ then } \mathbf{V}_T \times \mathbf{T}_T = \mathbf{0}, \\ \mathbf{V}_T \cdot \mathbf{T}_T & \ge 0 \text{ (slip)} \end{cases} $$

where $\mu > 0$ is the coefficient of friction.

It can be verified that the contact laws (19) and (20) are objective.

### 3. Weighted residual relationship for the contact problem

The strong form of the contact problem is described by Eqs. (1) and (2), the constitutive law, the boundary conditions (3), the contact laws (19) and (20) and the initial conditions (4). In this section, we convert the strong form into an equivalent weak one which is more convenient for the finite element implementation.

The approach chosen is inspired by the augmented Lagrangian formulation of Pietrzak and Curnier [13] and amounts to rewrite the augmented Lagrangian formulation of Alart and Curnier [4] in a generic weak form. More specifically, a weighted residual relationship is proposed which involves the displacement fields $\mathbf{U}^{(i)}$, $i \in \{1, 2\}$, defined in $Q_o^{(i)}$, and the multiplier field $\lambda_N(\mathbf{X})$ and $\lambda_T(\mathbf{X}) = \lambda_{Tz}a^\alpha$ defined on $S_{oc}^{(1)}$. Accordingly, the weighting functions are the virtual displacements $\mathbf{U}^{(1)}$, $\mathbf{U}^{(2)}$, and the virtual multipliers $\lambda_N^*, \lambda_T^* = \lambda_{Tz}a^\alpha$. Two positive constants $\epsilon_N, \epsilon_T$ being chosen, the weighted residual relationship is stated in the following proposition [21].

#### Proposition 1

* The solution fields of the strong problem – namely $\mathbf{U}^{(1)}, \mathbf{U}^{(2)}$ in $Q_o^{(1)}, Q_o^{(2)}$ and $T_N, \mathbf{T}_T, \mathbf{V}_T$ on $S_{oc}^{(1)}$ – satisfy the following relationship, provided one makes $\lambda_N = T_N$ and $\lambda_T = \mathbf{T}_T$:

$$ \forall t \in [0, T], \forall \mathbf{U}^{(1)*}, \forall \mathbf{U}^{(2)*}, \forall \lambda_N^*, \forall \lambda_{T1}^*, \forall \lambda_{T2}^*, $$

$$ \begin{aligned} & \sum_{i=1}^{2} \left\{ -\int_{\Omega_0^i} \mathbf{\Pi}^{(i)T} : \nabla \mathbf{X}^i \mathbf{U}^{(i)*} d\Omega_0 + \int_{\Omega_0^i} \rho_0^{(i)} \mathbf{f}^{(i)} \mathbf{U}^{(i)*} d\Omega_0 + \int_{S_{oc}^{(i)}} \mathbf{s}_{0i}^{(i)} \mathbf{T}_T^{(i)} \mathbf{U}^{(i)*} dS_0 \right\} \\ & + \int_{S_{oc}^{(i)}} \left[ \langle \lambda_N + \epsilon_N g \rangle \mathbf{v} - \left(1 - \left\langle 1 - \frac{\mu(\lambda_N + \epsilon_N g)}{\|\lambda_T + \epsilon_T \mathbf{v}_T\|} \right\rangle \right) (\lambda_T + \epsilon_T \mathbf{v}_T) \right] (\mathbf{U}^{(1)*}(\mathbf{X}) - \mathbf{U}^{(2)*}(\mathbf{Y}(\mathbf{X}))) dS_0 \\ & + \int_{S_{oc}^{(i)}} \left\{ \langle \lambda_N - (\lambda_N + \epsilon_N g) \rangle \frac{\dot{\lambda}_N}{\epsilon_N} + \left[ \lambda_T - \left(1 - \left\langle 1 - \frac{\mu(\lambda_N + \epsilon_N g)}{\|\lambda_T + \epsilon_T \mathbf{v}_T\|} \right\rangle \right) (\lambda_T + \epsilon_T \mathbf{v}_T) \right] \frac{\dot{\lambda}_T}{\epsilon_T} \right\} dS_0 \\ & = \sum_{i=1}^{2} \int_{\Omega_0^i} \rho_0^{(i)} \mathbf{U}^{(i)*} \mathbf{U}^{(i)} d\Omega_0 \end{aligned} \tag{21} $$

where $\nabla_{\mathbf{X}^i} \mathbf{U}^{(i)*}$ is the gradient tensor of $\mathbf{U}^{(i)*}$ with respect to variables $\mathbf{X}^{(i)}, \langle \cdot \rangle$ is the Macauley bracket: $\langle a \rangle = a$ if $a \ge 0, = 0$ if $a < 0$, and $\langle 1 - \frac{\mu(\lambda_N + \epsilon_N g)}{\|\lambda_T + \epsilon_T \mathbf{v}_T\|} \rangle$ must be replaced by 0 at any point on $S_{oc}^{(1)}$ where $\lambda_T + \epsilon_T \mathbf{v}_T = \mathbf{0}$.

* Conversely, the contact surfaces $S_c^{(1)}$ and $S_c^{(2)}$ being parametrized by $\xi \mapsto \mathbf{x} \in S_c^{(1)}$ and $\eta \mapsto \mathbf{y} \in S_c^{(2)}$, it is assumed that there is a bijective (or piecewise bijective) mapping $S_c^{(1)} \ni \mathbf{x} \mapsto \mathbf{y} \in S_c^{(2)}$ between them. Then, relation (21) implies at any time $t \in [0, T]$ the following local equations:

(1) The momentum balance Eq. (2) for the two bodies 1 and 2.

(2) Relation (3b) on the boundary portion $S_o^{(i)} \setminus S_{oc}^{(i)} = S_{oT}^{(i)} \cup S_{oU}^{(i)}: \mathbf{\Pi}^{(i)} \mathbf{N}^{(i)} = \mathbf{T}^{(i)}$.

(3) The following equalities between the components of the nominal traction vectors and the multipliers on the contactor surface $S_{oc}^{(1)}$.

$$ \begin{array}{l} \forall \mathbf{X} \in S_{oc}^{(1)}, \quad \mathbf{T}^{(1)} = \mathbf{\Pi}^{(1)} \mathbf{N}^{(1)} = \lambda_N \mathbf{v} - \lambda_T \\ \Leftrightarrow \quad T_N = \lambda_N \text{ and } \mathbf{T}_T = \lambda_T. \end{array} \tag{22} $$

(4) The normal and tangential contact laws (19) and (20).

(5) The following relationship which expresses the equilibrium of the traction vectors at the contact interface.

$$ \begin{aligned} & \forall \mathbf{X} \in S_{oc}^{(1)}, \\ & \mathbf{T}^{(2)}(\mathbf{Y}) dS_{oc}^{(2)} = -(\lambda_N \mathbf{v} - \lambda_T) \frac{\|\mathbf{X}_{1i} \times \mathbf{X}_{2i}\|}{\|\mathbf{Y}_{1i} \times \mathbf{Y}_{2i}\|} \frac{D_\xi}{D_\eta} dS_{oc}^{(2)}, \\ & = -\mathbf{T}^{(1)}(\mathbf{X}) dS_{oc}^{(1)} \end{aligned} \tag{23} $$

where point $\mathbf{Y} = \mathbf{Y}(\mathbf{X}, t) \in S_{oc}^{(2)}$ corresponds to point $\mathbf{X}$ through the above-mentioned bijection, $dS_{oc}^{(1)}$ is a differential reference area in $S_{oc}^{(1)}$ and $dS_{oc}^{(2)}$ its counterpart in $S_{oc}^{(2)}$ related to $dS_{oc}^{(1)}$.

#### Remarks

(i) In general, $\lambda_T + \epsilon_T \mathbf{v}_T$ differs from zero since $\lambda_T$ and $\mathbf{v}_T$ do not vanish simultaneously: if the slip rate vanishes, then the tangential traction does not, and vice versa. However, in some problems with a specific symmetry, there may be some points where $\lambda_T + \epsilon_T \mathbf{v}_T$ is zero. In practice, the set of these points is either isolated points or lines on the contact surface $S_{oc}^{(1)}$.

(ii) In dynamics, when shocks and impacts are expected, the acceleration $\dot{\mathbf{U}}^{(i)}$ should be understood as the derivative of the velocity in the sense of the distributions [22].

**Proof.** The first part of the above proposition, which states that the strong form implies the weak one, can be readily checked. Let us prove the reciprocal statement for an arbitrary fixed time $t \in [0, T]$.

* Since (21) holds for all virtual fields $\mathbf{U}^{(1)*}, \mathbf{U}^{(2)*}, \lambda_N^*, \lambda_T^*$ and $\lambda_T^*$, let us pick $\mathbf{U}^{(1)*} = \mathbf{U}^{(2)*} = \mathbf{0}$ and $\lambda_T^* = \mathbf{0}$, while keeping $\lambda_N^*$ arbitrary. Relationship (21) then leads to

$$ \lambda_N = \langle \lambda_N + \epsilon_N g \rangle \quad \text{on } S_{oc}^{(1)}. \tag{24} $$

The first consequence of (24) is $\lambda_N \ge 0$. Next, two mutually exclusive cases arise:

- either $\lambda_N + \epsilon_N g < 0$, then (24) entails:

$$ \lambda_N = 0 \Rightarrow g < 0 \quad (\text{since } \lambda_N + \epsilon_N g < 0) \tag{25} $$

- or $\lambda_N + \epsilon_N g \ge 0$, then (24) entails:

$$ \begin{array}{l} \lambda_N = \lambda_N + \epsilon_N g \Leftrightarrow g = 0, \text{ which gives } \lambda_N \ge 0 \\ \text{since } \lambda_N + \epsilon_N g \ge 0. \end{array} \tag{26} $$

The results for the two cases above can be gathered together as follows:

$$ g \le 0 \text{ where } \begin{cases} \cdot & \text{if } g < 0, \text{ then } \lambda_N = 0 \\ \cdot & \text{if } g = 0, \text{ then } \lambda_N \ge 0 \end{cases} \tag{27} $$

which is the normal contact law (19) provided $\lambda_N = T_N$ which will be shown later.
• Consider (21) again and choose now $\mathbf{U}^{(1)*} = \mathbf{U}^{(2)*} = \mathbf{0}$ and $\hat{\lambda}_N^* = 0$, whereas $\hat{\lambda}_T^*$ is left arbitrary. The localization theorem can then be applied assuming that the coefficient of $\hat{\lambda}_T^*$ is a continuous function in space. At any point on $S_{oc}^{(1)}$ where $\hat{\lambda}_T + \epsilon_T \mathbf{V}_T \neq \mathbf{0}$, one has

$$ \lambda_T = \left( 1 - \left\langle 1 - \frac{\mu(\hat{\lambda}_N + \epsilon_N g)}{\|\hat{\lambda}_T + \epsilon_T \mathbf{V}_T\|} \right\rangle \right) (\hat{\lambda}_T + \epsilon_T \mathbf{V}_T). \quad (28) $$

The first consequence of (28) is that the slip velocity $\mathbf{V}_T$ is parallel to $\hat{\lambda}_T$. Let us show a more precise result, specifically $\mathbf{V}_T$ is oriented in the same direction as $\hat{\lambda}_T$. To do this, let us temporarily denote $\alpha = \frac{\mu(\hat{\lambda}_N + \epsilon_N g)}{\|\hat{\lambda}_T + \epsilon_T \mathbf{V}_T\|} \ge 0$ for brevity. Relation (28) then reads:

$$ \langle 1 - \alpha \rangle \lambda_T = (1 - \langle 1 - \alpha \rangle) \epsilon_T \mathbf{V}_T. \quad (29) $$

As $\alpha \ge 0$ and $1 - \langle 1 - \alpha \rangle \ge 0$, $\forall \alpha$, one has the following mutually exclusive cases depending on $\alpha$'s value:

$$ \left\{ \begin{array}{ll} \cdot \text{ if } \alpha = 0, \text{ then (29) gives } \lambda_T = \mathbf{0}. & \\ \cdot \text{ if } 0 < \alpha < 1, \text{ then } 1 - \alpha > 0, 1 - \langle 1 - \alpha \rangle > 0 & (30) \\ \text{and (29) entails } (1 - \alpha)\lambda_T = \alpha\epsilon_T\mathbf{V}_T. & \\ \cdot \text{ if } \alpha \ge 1, \text{ then (29) gives } \mathbf{V}_T = \mathbf{0}. & \end{array} \right. $$

Thus, in any case $\mathbf{V}_T$ is oriented in the direction of $\hat{\lambda}_T$. Now, let us distinguish two cases:

- If $\hat{\lambda}_N + \epsilon_N g < 0$, then one gets from (25) and (28) $g < 0$ and $\hat{\lambda}_T = \mathbf{0}$.

- If $\hat{\lambda}_N + \epsilon_N g \ge 0$, then one gets from (26) $g = 0$. Two mutually exclusive sub-cases arise here:

- if $\|\hat{\lambda}_T + \epsilon_T \mathbf{V}_T\| < \mu(\hat{\lambda}_N + \epsilon_N g)$, then (28) entails $\mathbf{V}_T = \mathbf{0}$. From $g = 0$ just proven, one finds $\|\hat{\lambda}_T\| < \mu\hat{\lambda}_N$.

- if $\|\hat{\lambda}_T + \epsilon_T \mathbf{V}_T\| \ge \mu(\hat{\lambda}_N + \epsilon_N g)$, then (28) entails $\hat{\lambda}_T = \mu(\hat{\lambda}_N + \epsilon_N g) \frac{\hat{\lambda}_T + \epsilon_T \mathbf{V}_T}{\|\hat{\lambda}_T + \epsilon_T \mathbf{V}_T\|}$.

By taking the norm of both sides of the last equality, recalling here $\hat{\lambda}_N + \epsilon_N g \ge 0$ and $g = 0$, one finds $\|\hat{\lambda}_T\| = \mu\hat{\lambda}_N$. The last results can be summarized as follows:

$$ \left\{ \begin{array}{ll} \cdot \text{ if } g < 0, \text{ then } \lambda_T = \mathbf{0} & \\ \cdot \text{ if } g = 0, \text{ then } \|\hat{\lambda}_T\| \le \mu\hat{\lambda}_N, & (31) \\ \text{and} & \begin{array}{l} \cdot \text{ if } \|\hat{\lambda}_T\| < \mu\hat{\lambda}_N, \text{ then } \mathbf{V}_T = \mathbf{0} \\ \cdot \text{ if } \|\hat{\lambda}_T\| = \mu\hat{\lambda}_N, \text{ then } \mathbf{V}_T \\ \text{is oriented in the direction of } \hat{\lambda}_T \end{array} \end{array} \right. $$

which is merely the tangential contact law (20) provided $\hat{\lambda}_N = T_N$, $\hat{\lambda}_T = T_T$, which will be shown soon.

• Eventually, let us make in (21) $\hat{\lambda}_N^* = 0$, $\hat{\lambda}_T^* = \mathbf{0}$ and keep $\mathbf{U}^{(1)*}$ and $\mathbf{U}^{(2)*}$ arbitrary. According to the divergence theorem, relationship (21) leads to the following equality:

$$ \begin{aligned} & \sum_{i=1}^{2} \left\{ \begin{aligned} & \int_{\Omega_0^{(i)}} \mathbf{U}^{(i)*} (\mathrm{div} \mathbf{\Pi}^{(i)} + \rho_0^{(i)} \mathbf{f}^{(i)}) d\Omega_0 - \int_{S_{oc}^{(i)} \cup S_{0T}^{(i)}} \mathbf{U}^{(i)*} \cdot (\mathbf{\Pi}^{(i)} \cdot \mathbf{N}^{(i)} - \mathbf{T}^{(i)}) dS_0 \\ & - \int_{S_{oc}^{(1)}} \mathbf{U}^{(1)*} \cdot \mathbf{\Pi}^{(1)} \cdot \mathbf{N}^{(1)} dS_0 \\ & + \int_{S_{oc}^{(1)}} \mathbf{U}^{(1)*} \left[ \langle \hat{\lambda}_N + \epsilon_N g \rangle \mathbf{v} - \left( 1 - \left\langle 1 - \frac{\mu(\hat{\lambda}_N + \epsilon_N g)}{\|\hat{\lambda}_T + \epsilon_T \mathbf{V}_T\|} \right\rangle \right) (\hat{\lambda}_T + \epsilon_T \mathbf{V}_T) \right] dS_0 \\ & - \int_{S_{oc}^{(2)}} \mathbf{U}^{(2)*} \cdot \mathbf{\Pi}^{(2)} \cdot \mathbf{N}^{(2)} dS_0 \\ & - \int_{S_{oc}^{(1)}} \mathbf{U}^{(2)*} \left[ \langle \hat{\lambda}_N + \epsilon_N g \rangle \mathbf{v} - \left( 1 - \left\langle 1 - \frac{\mu(\hat{\lambda}_N + \epsilon_N g)}{\|\hat{\lambda}_T + \epsilon_T \mathbf{V}_T\|} \right\rangle \right) (\hat{\lambda}_T + \epsilon_T \mathbf{V}_T) \right] dS_0 \\ & = \sum_{i=1}^{2} \int_{\Omega_0^{(i)}} \rho_0^{(i)} \mathbf{U}^{(i)} \mathbf{U}^{(i)*} d\Omega_0, \end{aligned} \right\} \quad (32) $$

where $\mathbf{N}^{(i)}$ is the unit outward normal at point $\mathbf{X}^{(i)} \in S_0^{(i)}$.

Note that the last integral of the left-hand side of (32) is taken over $S_{oc}^{(1)}$. It can be recast into an integral over $S_{oc}^{(2)}$ by using the chain of bijections $\xi \to \mathbf{X} \to \mathbf{Y} \to \eta$ which allows one to change the integration variables:

$$ \begin{aligned} dS_0^{(1)} &= \|\mathbf{X}_{,\xi^1} \times \mathbf{X}_{,\xi^2}\| d\xi^1 d\xi^2 = \|\mathbf{X}_{,\xi^1} \times \mathbf{X}_{,\xi^2}\| \frac{D\xi}{D\eta} d\eta^1 d\eta^2, \\ dS_0^{(2)} &= \|\mathbf{Y}_{,\eta^1} \times \mathbf{Y}_{,\eta^2}\| d\eta^1 d\eta^2. \end{aligned} \quad (33) $$

Hence

$$ dS_0^{(1)} = \frac{\|\mathbf{X}_{,\xi^1} \times \mathbf{X}_{,\xi^2}\|}{\|\mathbf{Y}_{,\eta^1} \times \mathbf{Y}_{,\eta^2}\|} \frac{D\xi}{D\eta} dS_0^{(2)}. \quad (34) $$

By substituting (34) into (32), one gets relations (2), (3b), (23) and

$$ \begin{aligned} \mathbf{T}^{(1)} &= \mathbf{\Pi}^{(1)}\mathbf{N}^{(1)} \\ &= \langle \hat{\lambda}_N + \epsilon_N g \rangle \mathbf{v} - \left( 1 - \left\langle 1 - \frac{\mu(\hat{\lambda}_N + \epsilon_N g)}{\|\hat{\lambda}_T + \epsilon_T \mathbf{V}_T\|} \right\rangle \right) (\hat{\lambda}_T + \epsilon_T \mathbf{V}_T) \mathrm{on} S_{oc}^{(1)}. \end{aligned} \quad (35) $$

Comparing (35) with the definition (17), $\mathbf{T}^{(1)} = T_N \mathbf{v} - \mathbf{T}_T$, and making use of (24) and (28) lead to equality (22) between the contact tractions and the multipliers. Next, the normal and tangential contact laws (19) and (20) are derived from (22), (27) and (31).

• It remains to consider the points on $S_{oc}^{(1)}$ where $\hat{\lambda}_T + \epsilon_T \mathbf{V}_T = \mathbf{0}$. Equality (28) is then replaced by $\hat{\lambda}_T = \hat{\lambda}_T + \epsilon_T \mathbf{V}_T$, which entails $\hat{\lambda}_T = \mathbf{V}_T = \mathbf{0}$. The subsequent reasoning is analogous yet simpler than the case $\hat{\lambda}_T + \epsilon_T \mathbf{V}_T \neq \mathbf{0}$. It can easily be checked that the tangential contact law (20) as well as (Relations 22) and (23) are satisfied at the points where $\hat{\lambda}_T + \epsilon_T \mathbf{V}_T = \mathbf{0}$. The whole Proposition 1 has been proved. ■

### 4. Semi-discrete equation system

In this section, the weighted residual relationship (21) will be discretized in space by means of the finite element method. The discretization process will be described in the 3D framework only but the results for 2D problems can be derived in a similar manner. The following shorthand notations will be used for brevity:

$$ \hat{\lambda}_N = \hat{\lambda}_N + \epsilon_N g, \quad \hat{\lambda}_T = \hat{\lambda}_T + \epsilon_T \mathbf{V}_T. \quad (36) $$

We shall confine ourselves to the terms pertaining to contact in equality (21), denoted $\mathcal{W}^{*-\mathrm{contact}}$, since the other terms can be discretized in a standard way. Let the contactor surface $S_{oc}^{(1)}$ be subdivided into a number of isoparametric finite elements $e^{(1)}$ (the superscript ‘(1)’ is used throughout to remind of the contactor surface $S_{oc}^{(1)}$), then expression $\mathcal{W}^{*-\mathrm{contact}}$ can be recast as the sum of the element contact terms $(\mathcal{W}^{*-\mathrm{contact}})^{e^{(1)}}$ related to elements $e^{(1)}$ and evaluated using the Gauss quadrature rule:

$$ \begin{aligned} (\mathcal{W}^{*-\mathrm{contact}})^{e^{(1)}} &= \sum_{\text{quadrature points}} \left\{ \left[ \langle \hat{\lambda}_N \rangle \mathbf{v} - \left( 1 - \left\langle 1 - \frac{\mu(\hat{\lambda}_N)}{\|\hat{\lambda}_T\|} \right\rangle \right) \hat{\lambda}_T \right] \right. \\ & \qquad \left. \cdot (\mathbf{U}^{(1)*}(\mathbf{X}) - \mathbf{U}^{(2)*}(\mathbf{Y}(\mathbf{X}))) + (\hat{\lambda}_N - \langle \hat{\lambda}_N \rangle) \frac{\hat{\lambda}_N^*}{\epsilon_N} \right. \\ & \qquad \left. + \left[ \hat{\lambda}_T - \left( 1 - \left\langle 1 - \frac{\mu(\hat{\lambda}_N)}{\|\hat{\lambda}_T\|} \right\rangle \right) \hat{\lambda}_T \right] \frac{\hat{\lambda}_T^*}{\epsilon_T} \right\} j(\xi), \end{aligned} \quad (37) $$

where all terms are computed at Gauss points and $j(\xi)$ is the area jacobian times the weight corresponding to the parent coordinates $\xi = (\xi_1, \xi_2) \in \mathbb{R}^2$ of quadrature point $\mathbf{X} \in e^{(1)}$, see Fig. 1:
$$
j(\xi) = \| \mathbf{X}_{\xi_1}(\xi) \times \mathbf{X}_{\xi_2}(\xi) \| \cdot \text{(weight at location } \xi). \quad (38)
$$

In the sequel, the following notational convention is adopted for
the matrix representations: curly brackets {} designate a column
vector, the angle brackets <> an arrow vector, and the square
brackets [ ] a general matrix. In any element e<sup>(1)</sup> ⊂ S<sub>oc</sub><sup>(1)</sup>, the refer-
ence coordinates X, the displacement U<sup>(1)</sup> and the multipliers λ
are interpolated by

$$
\begin{align}
\{\mathbf{X}\} &= [\mathcal{N}^{(1)}(\xi)]^{e^{(1)}} \{\mathbf{X}\}^{e^{(1)}}, \nonumber \\
\{\mathbf{U}^{(1)}\} &= [\mathcal{N}^{(1)}(\xi)]^{e^{(1)}} \{\mathbf{U}^{(1)}\}^{e^{(1)}}, \tag{39} \\
\{\boldsymbol{\lambda}\} &= [\mathcal{N}^{(1)}(\xi)]^{e^{(1)}} \{\boldsymbol{\lambda}\}^{e^{(1)}}. \nonumber
\end{align}
$$

In (39b) for instance, {$\mathbf{U}^{(1)}$} is the column vector with the three
components of the displacement $\mathbf{U}^{(1)}$ related to a fixed Cartesian
orthonormal basis ($\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3$). The column vector {$\mathbf{U}^{(1)}$}e^{(1)} contains
the 3nne^{(1)} nodal displacement components associated with ele-
ment e^{(1)}, where nne^{(1)} is the number of nodes of element e^{(1)}.
The matrix [$\mathcal{N}^{(1)}(\xi)$]e^{(1)} contains the shape functions of element
e^{(1)} and is of dimension 3 × 3nne^{(1)}.

The normal and tangential components of multiplier $\lambda$ can be
derived from (39c) via:

$$
\begin{align}
\lambda_N = \langle \mathbf{v} \rangle [\mathcal{N}^{(1)}(\xi)]^{e^{(1)}} \{\boldsymbol{\lambda}\}^{e^{(1)}}, \tag{40} \\
\{\boldsymbol{\lambda}_T\} = -[\mathbf{I} - \mathbf{v} \otimes \mathbf{v}] [\mathcal{N}^{(1)}(\xi)]^{e^{(1)}} \{\boldsymbol{\lambda}\}^{e^{(1)}}, \nonumber
\end{align}
$$

where {$\lambda_T$} is a 3-component column vector.

Like the contactor surface, the target surface $S_{oc}^{(2)}$ is subdivided
into a number of isoparametric finite elements. Given a quadrature
point $\mathbf{X}$ in element $e^{(1)}$, let $e^{(2)}$ (the or an) element of $S_{oc}^{(2)}$
containing the reference point $\mathbf{Y}$ corresponding to the projection
$\mathbf{y} = \text{proj}_{\xi^{(2)}} \mathbf{x} \in S_{oc}^{(2)}$. In element $e^{(2)}$, the reference coordinates $\mathbf{Y}$
and the displacement $\mathbf{U}^{(2)}$ are interpolated by

$$
\begin{align}
\{\mathbf{Y}\} &= [\mathcal{N}^{(2)}(\boldsymbol{\eta})]^{e^{(2)}} \{\mathbf{Y}\}^{e^{(2)}}, \tag{41} \\
\{\mathbf{U}^{(2)}\} &= [\mathcal{N}^{(2)}(\boldsymbol{\eta})]^{e^{(2)}} \{\mathbf{U}^{(2)}\}^{e^{(2)}}, \nonumber
\end{align}
$$

where $\boldsymbol{\eta} = (\eta_1, \eta_2) \in \mathbb{R}^2$ is the parent coordinate of $\mathbf{Y}$ and the super-
script (2) is used to remind of element $e^{(2)} \subset S_{oc}^{(2)}$. The column vector
$\{\mathbf{U}^{(2)}\}^{e^{(2)}}$ contains the 3nne$^{(2)}$ nodal displacement components associated
with element $e^{(2)}$, where nne$^{(2)}$ is the number of nodes of ele-
ment $e^{(2)}$. The 3 × 3nne$^{(2)}$ matrix [$\mathcal{N}^{(2)}(\boldsymbol{\eta})$]$^{e^{(2)}}$ contains the shape
functions of element $e^{(2)}$.

As in [3], given a quadrature point X in element e<sup>(1)</sup> &sub; S<sub>oc</sub><sup>(1)</sup>, there
is at least one element e<sup>(2)</sup> &sub; S<sub>oc</sub><sup>(2)</sup> containing the reference point Y of
y = proj<sub>&xi;<sup>(2)</sup></sub> x. Both element e<sup>(2)</sup> and the number of nodes nne<sup>(2)</sup>,
which depend on the projection y, vary as a function of the quad-
rate point X. Moreover, the whole set of quadrature points of ele-
ment e<sup>(1)</sup> generally involves several elements e<sup>(2)</sup> in S<sub>oc</sub><sup>(2)</sup>.

All the virtual quantities are interpolated in a similar way. Even-
tually, the element contact expression (W<sup>-</sup><sub>contact</sub>)<sup>e(1)</sup> in (37) on the
element e<sup>(1)</sup> level can be written in the discretized form:

$$
\begin{equation}
\begin{split}
(\mathcal{W}_{\text{contact}}^{*})^{e^{(1)}} = \sum_{\text{quadrature points}} & \left( \left\langle \langle \mathbf{U}^{(1)*} \rangle^{e^{(1)}}, \langle \mathbf{U}^{(2)*} \rangle^{e^{(2)}} \right\rangle \right) \{\mathbf{\Phi}_{\text{contact}}\}^{e} \\
& - (\mathbf{A}^{*})^{e^{(1)}} \{\mathbf{R}_{\mathbf{A}}\}^{e} j(\xi),
\end{split}
\tag{42}
\end{equation}
$$

where {$\mathbf{\Phi}_{\text{contact}}$}^e is the element contact force vector and {$\mathbf{R}_A$}^e the element residual vector for the multipliers, defined as

$$
\{\mathbf{\Phi}_{\text{contact}}\}^e = \left\{ \begin{array}{l} [\mathcal{N}^{(1)}(\xi)]^{e^{(1)T}} \left( \langle \hat{\lambda}_N \rangle \{\mathbf{v}\} - \left(1 - \left\langle 1 - \frac{\mu(\hat{\lambda}_N)}{\|\hat{\lambda}_T\|} \right\rangle\right) \{\hat{\lambda}_T\} \right) \\ - [\mathcal{N}^{(2)}(\boldsymbol{\eta})]^{e^{(2)T}} \left( \langle \hat{\lambda}_N \rangle \{\mathbf{v}\} - \left(1 - \left\langle 1 - \frac{\mu(\hat{\lambda}_N)}{\|\hat{\lambda}_T\|} \right\rangle\right) \{\hat{\lambda}_T\} \right) \end{array} \right\}, \quad (43)
$$

$$
\{\mathbf{R}_A\}^e = [\mathcal{N}^{(1)}(\xi)]^{e^{(1)T}} \left( \frac{1}{\epsilon_N} \left( \langle \hat{\lambda}_N \rangle - \hat{\lambda}_N \right) \{\mathbf{v}\} \right. \\
\left. \qquad + \frac{1}{\epsilon_T} \left\{ \hat{\lambda}_T - \left( 1 - \left\langle 1 - \frac{\mu(\hat{\lambda}_N)}{\|\hat{\lambda}_T\|} \right\rangle \right) \hat{\lambda}_T \right\} \right). \quad (44)
$$

Vector {$\mathbf{\Phi}_{\text{contact}}$}^e in (43) has nne^{(1)} + nne^{(2)} components, where
nnec^{(2)} varies as a function of the considered quadrature point; vector
$\{\mathbf{R}_A\}^e$ has nne^{(1)} components. Note that

$$
1 - \left\langle 1 - \frac{\mu \langle \hat{\lambda}_N \rangle}{\| \hat{\lambda}_T \|} \right\rangle =
\begin{cases}
0 & \text{if } \hat{\lambda}_N < 0 \\
& \text{(algorithmic gap),} \\
\frac{\mu \hat{\lambda}_N}{\| \hat{\lambda}_T \|} & \text{if } \hat{\lambda}_N \ge 0 \text{ and } \| \hat{\lambda}_T \| \ge \mu \hat{\lambda}_N \\
& \text{(algorithmic contact with slip),} \\
1 & \text{if } \hat{\lambda}_N \ge 0 \text{ and } \| \hat{\lambda}_T \| < \mu \hat{\lambda}_N \\
& \text{(algorithmic contact with stick).}
\end{cases}
\tag{45}
$$

In (45), use has been made of the adjective ‘algorithmic’ in order
to emphasize the difference with the ‘real’ situation. For instance,
the ‘algorithmic gap’ case means $\hat{\lambda}_N = \hat{\lambda}_N + \epsilon_N g < 0$, in contrast
with the real ‘gap’ case which merely means $g < 0$; and the ‘algo-
rithmic slip’ case means $\|\hat{\lambda}_T\| \ge \mu \hat{\lambda}_N$, as opposed to the ‘slip’ case
which means $\|\hat{\lambda}_T\| = \mu \hat{\lambda}_N$.

By adding the element contributions in (42), one obtains the fol-
lowing semi-discrete coupled equation system:

$$
\begin{equation}
\begin{aligned}
& [\mathbf{M}]\{\dot{\mathbf{U}}\} + \{\Psi(\mathbf{U})\} = \{\Phi(\mathbf{U})\} + \{\mathbf{\Phi}_{\text{contact}}(\mathbf{U}, \mathbf{\Lambda})\}, \\
& \{\mathbf{R}_A(\mathbf{U}, \mathbf{\Lambda})\} = \{\mathbf{0}\},
\end{aligned}
\tag{46}
\end{equation}
$$

where vectors {$\mathbf{\Phi}_{\text{contact}}(\mathbf{U}, \mathbf{\Lambda})$} and {$\mathbf{R}_A(\mathbf{U}, \mathbf{\Lambda})$} are the assembly of the element vectors {$\mathbf{\Phi}_{\text{contact}}$}e$^{j(\xi)}$ and {$\mathbf{R}_A$}e$^{j(\xi)}$ given in (43) and (44), respectively. The unknowns in (46) are the nodal displacement vector {$\mathbf{U}$} = {$\begin{Bmatrix} \{\mathbf{U}^{(1)}\} \\ \{\mathbf{U}^{(2)}\} \end{Bmatrix}$} in the two bodies 1 and 2, and the multiplier vector {$\mathbf{\Lambda}$} on the reference contact surface $S_{oc}^{(1)}$. Except for the above-mentioned vectors {$\mathbf{\Phi}_{\text{contact}}(\mathbf{U}, \mathbf{\Lambda})$} and {$\mathbf{R}_A(\mathbf{U}, \mathbf{\Lambda})$} which are specific to contact, the other vectors and matrices are standard and result from the discretization of terms other than $\mathcal{W}^{*-\text{contact}}$ in relationship (21): $[\mathbf{M}]$ is the mass matrix, {$\Psi(\mathbf{U})$} the internal force vector resulting from the strain work in the bodies, {$\mathbf{\Phi}(\mathbf{U})$} the external force vector (excluding contact force) which may depend on displacement {$\mathbf{U}$} if follower loads (e.g., a pressure) are applied on either body. The nonlinear semi-discrete Eq. (46) will be solved as a whole using the standard full Newton-Raphson method to simultaneously obtain the nodal displacements and the multipliers, as was done in [4, 7, 13, 8] (Curnier and co-workers even used a so-called generalized Newton method which is more sophisticated).

* The slip velocity $\mathbf{V}_T = \dot{\eta}^\gamma \mathbf{a}_\gamma$ defined in (15) and involved in (43) and (44) has to be computed here by means of an adequate temporally discrete expression. For this purpose, let us first note that Proposition 1 remains valid if $\mathbf{V}_T$ is replaced with $\mathbf{V}_T\Delta t$, where $\Delta t$ is an arbitrary finite time interval. The quantity $\mathbf{V}_T\Delta t$ having the dimension of a length will be referred to as the slip, its use instead of the slip velocity is more consistent in the numerical context because the penalty parameters $\epsilon_N$ and $\epsilon_T$ involved in the weak form then have the same dimension of a
force by unit volume, i.e. the same SI units of N/m³. Let us now look at a typical iteration of the current time step n and introduce some new notations as summarized in Fig. 2. For the sake of brevity, any quantity computed at the current step n will not be indexed by n, whereas any quantity related to the previous step will be distinguished by subscript n-1. Consider two times t and t<sub>n-1</sub>, which are respectively the current time corresponding to the current iteration of the current step n and the time corresponding to the previous step (n-1). Given a particle on the contactor surface S<sub>c</sub><sup>(1)</sup> with reference position X, its position at time t<sub>n-1</sub> is x<sub>n-1</sub> = φ<sup>(1)</sup>(X, t<sub>n-1</sub>) and its current position is x = φ<sup>(1)</sup>(X, t). We denote the following quantities related to time t<sub>n-1</sub>:

- **y**<sub>n-1</sub> the projection of **x**<sub>n-1</sub> on the target surface S<sub>c</sub><sup>(2)</sup>(t<sub>n-1</sub>) at time t<sub>n-1</sub>,

- **Y**<sub>n-1</sub> and η<sub>n-1</sub> the inverse images of **y**<sub>n-1</sub> under the chain of parametrizations (5).

It should be emphasized that (i) **y** and **y**<sub>n-1</sub> are not the respective positions at times t and t<sub>n-1</sub> of a same particle in body 2, and (ii) the current position of the particle **Y**<sub>n-1</sub> is φ<sup>(2)</sup>(**Y**<sub>n-1</sub>, t), in general distinct from **y** = φ<sup>(2)</sup>(**Y**<sub>n-1</sub>, t). With these notations in hand, using the chain of parametrizations (5), which gives **y** = ψ<sup>(2)</sup>(η, t) and ψ<sup>(2)</sup>(η<sub>n-1</sub>, t) = Φ<sup>(2)</sup>(**Y**<sub>n-1</sub>, t), and applying the backward Euler integration scheme with time step Δt = t - t<sub>n-1</sub> leads to the following expression for the slip:

$$ \mathbf{V}_T \Delta t = \mathbf{y} - \mathbf{\Phi}^{(2)}(\mathbf{Y}_{n-1}, t). \quad (47) $$

It should be noted that other expressions such as V<sub>T</sub>Δt = ψ<sup>(2)</sup>(η, t) - ψ<sup>(2)</sup>(η<sub>n-1</sub>, t) would not be suitable as surface S<sub>c</sub><sup>(2)</sup> is divided into finite elements and the parametrization (5) is in fact defined piecewisely. It may happen that the projection points **y**<sub>n-1</sub> and **y** do not belong to the same element on S<sub>c</sub><sup>(2)</sup> although they are close to each other; in this case their local parent coordinates η<sub>n-1</sub> and η may be quite distant even though the time step Δt tends to zero. The algebraic expression (47) does not involve the parent coordinates and thus avoids the aforementioned discontinuity problem. It is adopted to compute the slip in the contact force vector (43) and the residual (44), as well as the tangent stiffness in the next section.

* The important issues of convergence and stability of a mixed finite element formulation were investigated in the literature for continuum problems and later for contact problems. They are related to the inf-sup condition which is expressed in contact problems as a compatibility condition between the displacement and multiplier spaces on the contact surfaces.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Illustration of the slip V<sub>T</sub>Δt.</figcaption></figure>

Following Chapelle and Bathe's developments [23] on incompressible elasticity, El-Abbasi and Bathe [24] examined the small deformation frictionless contact problem and proposed an inf-sup test procedure where a series of refined meshes is used in order to evaluate the evolution of the numerical inf-sup value. A further mathematical study was carried out by Bathe and Brezzi in [25], giving rise to precise relations between the degrees of the polynomial interpolations for the displacements and the multipliers. In general, the mathematical study of the inf-sup condition is confined to the small deformations framework since it is difficult when the structure undergoes large deformations and the contact surface is unknown a priori. No attempt has been made in this work in order to prove-analytically or numerically-that the proposed method passes the inf-sup condition. This issue should be investigated in the future. In the meanwhile, it can be observed through the numerical examples presented in Section 6 that the convergence is achieved when the mesh is refined and there is neither locking nor spurious contact pressure.

## 5. Contact tangent stiffness matrix

This section is devoted to the derivation of the tangent stiffness matrix required for the Newton-Raphson solution of the coupled equation system (46). We will confine ourselves to that part of the tangent stiffness related to contact, since the remaining part related to non contact terms can be linearized in a standard way. Let δU and δλ be the independent variations of the displacement and the multiplier, respectively defined by

$$ \delta \mathbf{U} = \delta \mathbf{U}^i \mathbf{e}_i, \quad \delta \lambda = \delta \lambda_N \mathbf{v} - \delta \lambda_T, \quad \delta \lambda_T = \delta \lambda_{T2} \mathbf{a}^\alpha \perp \mathbf{v}, \quad (48) $$

where Latin indices take the values 1, 2, 3, Greek indices take the values 1, 2, repeated indices imply summation and the variations δU<sup>i</sup>, i ∈ {1, 2, 3}, δλ<sub>N</sub>, δλ<sub>T1</sub>, δλ<sub>T2</sub> are arbitrary. One has to compute the exact linearization of ψ<sup>*contact</sup> corresponding to these variations and extract the contact stiffness [K<sub>contact</sub>] from the relation:

$$ \forall \{\delta \mathbf{U}\}, \{\delta \mathbf{\Lambda}\}, \quad \delta \mathcal{W}^{*}_{\text{contact}} = - \left\{ \left\{ \begin{matrix} \mathbf{U} \\ \mathbf{\Lambda}^* \end{matrix} \right\} \right\}^T [\mathbf{K}_{\text{contact}}] \left\{ \begin{matrix} \delta \mathbf{U} \\ \delta \mathbf{\Lambda} \end{matrix} \right\}. \quad (49) $$

The variation (δψ<sup>*contact</sup>)<sup>e(1)</sup> over any element e<sup>(1)</sup> ⊂ S<sub>oc</sub><sup>(1)</sup> is obtained by taking the directional derivatives of (37) in the directions δU<sup>(i)</sup>, i ∈ {1, 2}, δλ<sub>N</sub> and δλ<sub>T</sub>. The variations of the basic kinematic variables are expressed in terms of the displacement variations via [2,3]:

$$ \begin{align} \delta g &= - (\delta \mathbf{U}^{(1)}(\mathbf{X}) - \delta \mathbf{U}^{(2)}(\mathbf{Y}(\mathbf{X}, t))) \cdot \mathbf{v}, \nonumber \\ \delta \eta^\alpha &= \bar{\mathbf{a}}^{\alpha\beta} \mathbf{a}_\beta (\delta \mathbf{U}^{(1)}(\mathbf{X}) - \delta \mathbf{U}^{(2)}(\mathbf{Y}(\mathbf{X}, t))) - g \bar{\mathbf{a}}^{\alpha\beta} \mathbf{v} \delta \mathbf{U}_\beta^{(2)}(\eta), \nonumber \\ \delta \mathbf{v} &= -(\kappa_{\alpha\gamma} \bar{\mathbf{a}}^{\gamma\beta} \mathbf{a}_\beta) (\delta \mathbf{U}^{(1)}(\mathbf{X}) - \delta \mathbf{U}^{(2)}(\mathbf{Y}(\mathbf{X}, t))) - (\bar{\mathbf{a}}^{\alpha\beta} \mathbf{a}_\beta \otimes \mathbf{v}) \delta \mathbf{U}_\alpha^{(2)}(\eta), \tag{50} \end{align} $$

where δU<sup>(2)</sup><sub>α</sub>(η) denotes the partial derivative with respect to η<sup>α</sup> of the composite function δU<sup>(2)</sup>(η) = δU<sup>(2)</sup>(Y(η)).

As for the variation of the slip, it is derived from (47):

$$ \delta(\mathbf{V}_T \Delta t) = \delta \mathbf{U}^{(2)}(\mathbf{Y}) - \delta \mathbf{U}^{(2)}(\mathbf{Y}_{n-1}) + \mathbf{W} \cdot (\delta \mathbf{U}^{(1)}(\mathbf{X}) - \delta \mathbf{U}^{(2)}(\mathbf{Y})) - \mathbf{Z}^\alpha \delta \mathbf{U}_{,\alpha}^{(2)}(\eta), \quad (51) $$

where W and Z<sup>α</sup>, α ∈ {1, 2}, are non symmetric second-order tensors defined as

$$ \mathbf{W} = \bar{\mathbf{a}}^{\alpha\beta} \mathbf{a}_\alpha \otimes \mathbf{a}_\beta, \quad \mathbf{Z}^\alpha = g \bar{\mathbf{a}}^{\alpha\beta} \mathbf{a}_\beta \otimes \mathbf{v}. \quad (52) $$
Like the displacement, the variations of the displacement and the multiplier are obtained by interpolating the nodal variations

$$
\begin{align}
\{\delta \mathbf{U}^{(1)}\} &= [\mathcal{N}^{(1)}(\xi)]^{e(1)} \{\delta \mathbf{U}^{(1)}\}^{e(1)}, \nonumber \\
\{\delta \mathbf{U}^{(2)}\} &= [\mathcal{N}^{(2)}(\eta)]^{e(2)} \{\delta \mathbf{U}^{(2)}\}^{e(2)}, \nonumber \\
\{\delta \boldsymbol{\lambda}\} &= [\mathcal{N}^{(1)}(\xi)]^{e(1)} \{\delta \boldsymbol{\lambda}\}^{e(1)}. \tag{53}
\end{align}
$$

In the above relationship, {$\delta\lambda$} for instance is the 3-component
column vector representing $\delta\lambda$ in the fixed basis ($\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3$) and the
column vector {$\delta\lambda$}^{e(1)} with $3nne^{(1)}$ components contains the trip-
lets of the nodal components of $\delta\lambda$ in the same basis. One also
needs to interpolate the variation {$\delta\mathbf{U}^{(2)}(\mathbf{Y}_{n-1})$} involved in (51):

$$
\{\delta \mathbf{U}^{(2)}(\mathbf{Y}_{n-1})\} = [\mathcal{N}^{(2)}(\boldsymbol{\eta}_{n-1})]^{e(2)} \{\delta \mathbf{U}^{(2)}\}^{e(2)}, \quad (54)
$$

where $e_{n-1}^{(2)}$ denotes the element in the reference target surface $S_{oc}^{(2)}$
containing point $\mathbf{Y}_{n-1}$, in general different from element $e^{(2)}$.

Eventually, relations (50)–(54) enable one to express the ele-
ment linearization $(\delta \mathcal{W}_{\text{contact}}^{*})^{e(1)}$ in the form:

$$
(\delta \mathcal{W}_{\text{contact}}^{*})^{e(1)} = - \sum_{\text{quadrature points}} \begin{Bmatrix} \{\mathbf{U}^{(1)*}\}^{e(1)} \\ \{\mathbf{U}^{(2)*}\}^{e(2)} \\ \{\mathbf{\Lambda}^{*}\}^{e(1)} \end{Bmatrix}^T \begin{Bmatrix} \{\delta \mathbf{U}^{(1)}\}^{e(1)} \\ \{\delta \mathbf{U}^{(2)}\}^{e(2)} \\ \{\delta \mathbf{\Lambda}\}^{e(1)} \end{Bmatrix}, \quad (55)
$$

where the element contact stiffness $[\mathbf{K}_{\text{contact}}]^{e}$ is a rectangular matrix of the form:

$$
[\mathbf{K}_{\text{contact}}]^{e} = \begin{bmatrix} [\mathbf{K}_{UU}]^{e} & [\mathbf{K}_{U\Lambda}]^{e} \\ [\mathbf{K}_{\Lambda U}]^{e} & [\mathbf{K}_{\Lambda\Lambda}]^{e} \end{bmatrix} \quad (56)
$$

The explicit expression of the above matrix is given in the Appendix. It should be emphasized that (i) despite the superscript *e* which stands for element, $[\mathbf{K}_{\text{contact}}]^{e}$ is actually defined for each quadrature point of element $e^{(1)} \subset S_{oc}^{(1)}$ and (ii) although the element matrices $[\mathbf{K}_{\text{contact}}]^{e}j(\xi)$ are rectangular, their assembly eventually results in the square contact stiffness $[\mathbf{K}_{\text{contact}}]$ of the form:

$$
[\mathbf{K}_{\text{contact}}] = \begin{bmatrix} [\mathbf{K}_{UU}] & [\mathbf{K}_{U\Lambda}] \\ [\mathbf{K}_{\Lambda U}] & [\mathbf{K}_{\Lambda\Lambda}] \end{bmatrix} \quad (57)
$$

**6. Numerical examples**

This section presents the numerical examples which are carried
out on the basis of the previous finite element formulation in order
to assess the efficiency of the proposed formulation. Although the
semi-discrete equation system (46) holds in dynamics, from now
on we will confine ourselves to those quasistatic cases that display
enough specific features of the contact phenomenon to consider
the formulation validated. The dynamic cases present additional
physical and numerical difficulties and are deliberately disre-
garded at this stage. All numerical results presented here have
been obtained by a home made finite element code.

6.1. Contact patch test

The first example we consider is related to the contact patch
test which was originally conceived in [26,27], with the aim of
checking whether the finite element contact formulation allows
transmission of a constant pressure through an arbitrary noncon-
forming contact surface.

In this work, we used the patch test proposed in [28,24] and
shown in Fig. 3a. An elastic body is pressed against another one
by a uniform distributed load p = 1000 N/m² applied on its top sur-
face. Assuming a frictionless interface (μ = 0) between the two
bodies, the exact solution of the patch test entails a constant con-
tact pressure equal to p along the contact line. The patch test is
passed if the numerical results yield such a uniform stress distribu-
tion, regardless of the mesh used. In practice, one has to verify the
test using some typical meshes.

Here the problem is solved assuming that the two bodies are
made of an isotropic hyperelastic Saint-Venant Kirchhoff material
characterized by the strain energy function in terms of the Green
strain **E**:

$$
w(\mathbf{E}) = \frac{1}{2} \lambda (\mathrm{tr} \, \mathbf{E})^2 + \bar{\mu} \mathrm{tr} (\mathbf{E}^2), \quad (58)
$$

where the Lamé constants $\lambda$, $\bar{\mu}$ are related to Young's modulus $E$ and
Poisson's ratio $\nu$ by the same relations as in linear elasticity:

$$
\lambda = \frac{E\nu}{1 + \nu} \frac{1 - \nu}{1 - 2\nu}, \quad \bar{\mu} = \frac{E}{2(1 + \nu)}. \quad (59)
$$

The problem is modeled in plane strain conditions, taking the
dimensions of the bodies $a = 5$ m and $b = 1$ m, and the material

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Contact patch test. (a) Problem definition. (b) Regular mesh. (c) Distorted mesh.</figcaption></figure>
properties $E = 10^6$ N/m² and $v = 0$. The structure is discretized using two meshes – the first one is regular and the second distorted – as shown in Fig. 3b–c, in both cases the mesh is non-conforming along the contact interface. The two bodies are discretized with 8-node quadrilaterals (Q8), so that the elements on the contact surface are 3-node line elements. The numerical quadrature is performed using $3 \times 3$ Gauss quadrature points in the Q8 elements and 7 Gauss points in line elements belonging to the contact surface. The penalty parameter is $\epsilon_N = 10^{12}$ N/m³.

Fig. 4 shows the deformed configurations and contact stress distributions along the contact line when the upper, respectively the lower, body is selected to be the contactor. The arrows represent the contact stresses exerted by the target body upon the contactor. In all cases, the contact stresses are found to be constant and equal to $p$ up to 5 or 6 significant digits.

All the nodes on the contact interface have the same vertical displacement $U_y = -0.0010015$ m. The stress field inside the two bodies is constant, with the only non-zero stress component $\Sigma_{yy} = -1001$ N/m²($\Sigma$ is the second Piola–Kirchhoff stress tensor). These values agree with the analytical solution of the nonlinear elasticity (nonlinear elasticity must be used for the comparison since the numerical results are obtained by a nonlinear finite element code).

## 6.2. Hertzian contact problem

In this example, we consider a Hertzian contact problem in order to investigate the accuracy of the proposed contact formulation. Referring to Fig. 5a, a 2D elastic cylinder with radius $R = 10$ m is pressed against a rigid foundation by a force $F$ which is perpendicular to the rigid plane and applied on the symmetry axis of the cylinder. The cylinder is made of the same isotropic hyperelastic Saint–Venant Kirchhoff material as in (58), here we take $E = 100$ N/m², $v = 0.3$. The computations are performed assuming plane strain conditions and either a frictionless ($\mu = 0$) or a frictional ($\mu = 0.3$) contact between the cylinder and the foundation.

In the frictionless case, the analytical solution assuming the linear small strain theory of elasticity is well known. The pressure distribution in the contact region is given by

$$p(x) = \frac{2F}{\pi b^2} \sqrt{b^2 - x^2}, \qquad (60)$$

where $b$ is the half contact region width:

$$b = 2\sqrt{\frac{FR(1-v^2)}{\pi E}}. \qquad (61)$$

<figure><img src="image_5.png" /><figcaption>Fig. 4. Deformed configuration and contact stress distribution in the patch test: (a–b) when the upper body is selected to be contactor, (c–d) when the lower body is selected to be contactor (the displacement field is magnified by 100).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. 2D Hertzian contact problem of an elastic cylinder on a rigid foundation. (a) Problem definition. (b) Finer mesh.</figcaption></figure>
The computation is carried out with two values of the total load: F = 1 N/m and F = 2 N/m. Due to symmetry, only half of the structure is considered. Two meshes are considered: a coarse one containing 198 elements and 455 nodes, and a finer one containing 345 elements and 773 nodes (the latter is shown in Fig. 5b). In both cases, the rigid foundation is chosen as the target and modeled by one single element with fixed nodes. Clearly, the meshes are nonconforming.

The numerical quadrature employs 7 points in the triangles and the elements belonging to the contact line. The penalty parameter is $\epsilon_N = 100$ N/m$^3$. The (relative) convergence tolerances are taken to be $10^{-6}$ for the force residual $\{\mathbf{R}(\mathbf{U}, \mathbf{\Lambda})\} = \{\mathbf{\Psi}(\mathbf{U})\} - \{\mathbf{\Phi}(\mathbf{U})\} - \{\mathbf{\Phi}_{\text{contact}}(\mathbf{U}, \mathbf{\Lambda})\}$ in (46a) and the multiplier residual $\{\mathbf{R}_A(\mathbf{U}, \mathbf{\Lambda})\}$ in (46b). Fig. 6 plotting the numerically computed contact pressure acting on the cylinder shows that the finite element solution converges as the mesh is refined.

It should be noted that for a given mesh and a given load, there is in general one element of the cylinder which overlaps the contact area and the free surface, as can be seen in Fig. 6. In such an element the contact stress vector-or the multiplier $\lambda$-differs from zero in the contact zone and vanishes outside, and the interpolation of $\lambda$ from the nodal multipliers, see Relation (39c), is not precise. Despite this difficulty, the numerical results obtained here are in good agreement with the theory. With the finer mesh and $F = 1$ N/m, the maximal pressure at $x = 0$ is 2% lower than the theoretical value. With the larger load $F = 2$ N/m, the numerical results deviate more from the theory, the maximal pressure at $x = 0$ is 2.8% lower than the theoretical value.

In actual fact, the contact stresses obtained here are the nominal Piola-Kirchhoff ones (i.e. the multipliers according to (22)) whereas the analytical pressure in (60) comes from the linear elasticity theory. Although the deformation is small, there still is a slight difference between these stresses. An elaborate calculation shows that the Cauchy stresses corresponding to the numerical Piola-Kirchhoff ones are twice closer to the theoretical values. Having said this, we decide to only present the Piola-Kirchhoff stresses throughout this work.

In the frictional case, the computations are made with the fine mesh only. The analytical solution was obtained by Spence [29,30], assuming that the influence of friction on the normal stress is negligible in order to uncouple the equations governing the normal and tangential stresses. The contact region is subdivided into a stick zone $|x| \le c$ and a slip zone $c \le |x| \le b$. The extent $b$ of the contact region is given by the same relation (61) as in the friction-

less case, i.e. $b = 0.3404$ m if $F = 1$ N/m. The extent $c$ of the stick region is obtained by solving the following relation giving $c/b$ as a function of the coefficient of friction $\mu$ and Poisson's ratio $\nu$

$$ \frac{K(c/b)}{K'(c/b)} = \frac{\tan^{-1} \mu}{\frac{1}{2} \frac{1+\beta}{1-\beta}}, \qquad (62) $$

where $\beta = \frac{1-2\nu}{2(1-\nu)}$, $K(k)$ denotes the complete elliptic integral of the first kind, $K(k) = \int_0^{\pi/2} \frac{d\phi}{\sqrt{1-k^2\sin^2\phi}}$, $0 \le k \le 1$, and $K'(k) = K(\sqrt{1-k^2})$. With $F = 1$ N/m, $\mu = 0.3$ and $\nu = 0.3$, it is found that $c/b = 0.7006$, so $c = 0.2385$ m.

The normal pressure is given by the same relation (60) as in the frictionless case while the tangential contact stresses is computed via:

$$ q(x) = \begin{cases} \frac{2\mu Fx}{\pi b^2} \left[ \frac{1}{K(c/b)} \int_{x/b}^{c/b} \frac{K(\zeta_b \sin^{-1} \frac{t}{c/b})}{t^2 \sqrt{1-t^2}} dt + \frac{b}{c} \sqrt{1-(\frac{\zeta_b}{c})^2} \right], & |x| \le c \\ \mu \cdot p(x), & c \le |x| \le b \end{cases} \qquad (63) $$

where $K(k, \psi)$ is the incomplete elliptic integral of the first kind, $K(k, \psi) = \int_0^\psi \frac{d\phi}{\sqrt{1-k^2\sin^2\phi}}$, $0 \le k \le 1$,

The total load $F = 1$ N/m is applied in 50 equal increments, the penalty parameters are $\epsilon_N = \epsilon_T = 100$ N/m$^3$. In each time step, the solution is obtained after 3 iterations in average, with a maximum of 7 iterations at the first time step. The shape of the stresses along the contact interface is displayed in Fig. 7a, which clearly shows the progressive transition from stick to slip over the contact line. The normal $T_y$ and tangential $T_x$ contact stresses are plotted in Fig. 7b, showing a good agreement with Spence's theoretical solution, except in the central stick region where the friction seems to have the effect of raising the normal contact pressure. The value at $x = 0$ of the normal stress exceeds the analytical one by about 4%, showing the same trend as in [12] where a similar problem of the contact between a rigid punch and an elastic halfspace was investigated.

In order to check whether the finite element procedure enables one to satisfy the local inequality constraints, the ratio $T_x/\mu T_y$ is computed at every node on the contact line for $F = 0.5$ and $1$ N/m and it is verified that the Coulomb inequality $|T_x| \le \mu T_y$ holds everywhere, see Fig. 7c. At the sliding nodes, the equality $|T_x| = \mu T_y$ is met within 0.1%.

<figure><img src="image_7.png" /><figcaption>Fig. 6. Contact stress distribution in the frictionless Hertzian contact problem. (a) Coarse mesh. (b) Fine mesh. ◦ or •: numerical results; solid lines: theoretical results.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Contact stress distribution in the frictional Hertzian contact problem ($F = 1$ N/m, $\mu = 0.3$): (a) Spatial distribution along the contact line, (b) ◦ or •: numerical results, solid lines: theoretical results. (c) Check of the Coulomb friction law ($F = 0.5$ and $1$ N/m).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Frictional Hertzian contact ($F = 1$ N/m, $\mu = 0.3$). Influence of the loading rate on the contact stresses.</figcaption></figure>

Eventually, two additional computations are carried out with only one and 10 increments instead of 50. As shown in Fig. 8, the contact width, the normal contact pressure as well as the frictional shear in the slip region can be recovered with good precision using one or 10 time steps, whereas the frictional shear distribution in the stick zone is rather sensitive to the load increment and requires up to 50 time steps to be accurate. This result is in full accord with Pietrzak's observation in [12].

### 6.3. Contact between two beams

In this example, we are concerned with the quasistatic contact between two beams, one of which being subjected to either a dead or a follower force, with a view to showing the ability of the proposed formulation to deal with large deformation contact. The two beams have identical reference geometries, with rectangular cross section, a length of $L = 1$ m along x, a height of 0.05m along y and a thickness of 0.1m along z, see Fig. 9. They are both made of an isotropic hyperelastic Saint-Venant Kirchhoff material with $E = 2.10^6$ N/m$^2$ and $v = 0.3$.

The beams are clamped at end $X = 0$. In the reference configuration, the lower face of the upper beam (contactor) coincide with the upper face of the lower beam (target).
<figure><img src="image_10.png" /><figcaption>Fig. 9. Reference configuration of the two beam problem: (a) 3D view of the mesh (b) x-y in-plane view.</figcaption></figure>

Each beam is modeled with 15 20-node brick (i.e. hexahedral) elements (H20) and the contact surface with 15 8-node quadrilaterals (Q8). The numerical quadrature employs 3 × 3 × 3 Gauss points in the H20's and 6 × 6 points in the Q8 elements belonging to the contact surface.

In the following, two types of loading, a dead and a follower force, are applied on the contactor. For each loading, the computations are carried out assuming either frictionless or contact, $\mu = 0.3$. In any case, the penalty parameters are $\epsilon_N = \epsilon_T = 10^5$ N/m$^3$.

### 6.3.1. Dead load

A dead surface load $-F_y$ is uniformly distributed over the cross section $X = 1$ m of the upper beam. Force $F$ is progressively increased from zero using an arc length method with the arc length regularly incremented by a factor 1.1. The (relative) convergence tolerances are taken to be $10^{-6}$ for the force residual ${R(U, \Lambda)}$ in (46), $10^{-6}$ for the arc length residual and $10^{-10}$ for the multiplier residual ${R_A(U, \Lambda)}$ in (46b).

In the frictionless case, the computation is carried out in 29 steps, the resultant force at the last step is $F = 62.62$ N and the convergence is achieved in 7-8 iterations in average. Fig. 10 shows the deformed configurations and the contact traction vectors acting on the upper beam. Only x-y in-plane views are shown since the complete 3D views are not legible. In the frictional case, the entire calculation requires 38 steps, the final resultant force is $F = 61.20$ N and the convergence achieved in 11-12 iterations in average. The corresponding deformed configurations and the contact tractions are shown in Fig. 11.

The contact immediately takes place after force $F$ is applied and persists over the entire lower surface of the upper beam until the last step. During the first steps, the contact tractions are higher near the end cross section where force $F$ is applied, whereas far from this section, the contact tractions are smaller and barely visible in Figs. 10 and 11. As the deformations grow up, the contact tractions become significant in the neighborhood of the built-in section. It should be mentioned that there are possible spatial oscillations of the contact stresses due to the use of quadratic shape functions which take negative values at certain points and due to coarse meshes. Thus, during the bending process, there may be some nodes on the contact surface where the contact law is not satisfied, namely $T_N < 0$ or $||T_T|| > \mu \cdot T_N$ in frictional case. Fortunately, these violations only occur at the nodes where $T_N$ and $||T_T||$ are very small, specifically smaller than about 1% of the maximal $T_N$ and $||T_T||$ values on the whole contact surface. Moreover, the negative contact pressures disappear if the mesh is fine enough. In Figs. 10 and 11, the contact tractions violating the contact law are hardly visible. In presence of the friction, the contact state is rather complex. Roughly speaking, it is found that there

<figure><img src="image_11.png" /><figcaption>Fig. 10. Contact between two beams. Dead load and frictionless case: (a) $F = 1.625$ N, (b) $F = 6.914$ N, (c) $F = 17.17$ N, (d) $F = 29.21$ N, (e) $F = 41.11$ N, (f) $F = 62.62$ N. The dead load $F$ is not drawn to scale.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Contact between two beams. Dead load and frictional case $\mu = 0.3$: (a) $F = 1.769$ N, (b) $F = 6.914$ N, (c) $F = 19.78$ N, (d) $F = 34.36$ N, (e) $F = 49.28$ N, and (f) $F = 61.20$ N. The dead load *F* is not drawn to scale.</figcaption></figure>

is slip at most nodes and every step, and stick at some nodes and some steps.

Fig. 12 shows that the deflection at point A (central point of the cross section $X = 1$ m of the upper beam, as shown in Fig. 9) versus force *F* is almost the same, with or without friction. If the friction prevented the slip to occur at some steps, force *F* necessary to reach the same deflection would be higher than in the frictionless case and the two plots would be distinct.

### 6.3.2. Follower force

The dead force considered above is now replaced with a uniform pressure distributed over the last Q8 element of the upper

face of the upper beam, next to the cross section $X = 1$ m. The pressure is a follower force in the sense that it remains continually perpendicular to the area on which it is applied, so that its direction changes as the beam deforms, contrary to the dead force which keeps a fixed direction. We choose the same automatic control of the pressure and convergence tolerances as in the dead load case.

In the frictionless case, the computation is carried out in 34 steps leading to the final pressure resultant $P = 73.84$ N, the convergence is achieved in 9–10 iterations in average and the deformed configurations together with the contact stresses are shown in Fig. 13. In the frictional case ($\mu = 0.3$), the computation is carried out in 43 steps, the final pressure resultant is $P = 74.83$ N, the convergence is achieved in 10–11 iterations in average and the results are shown in Fig. 14.

As shown in Fig. 13, at the beginning the contact occurs on the whole lower surface of the upper beam. At an advanced deformed state, the absence of friction leads to a release over a large portion of the former contact surface, as shown in Fig. 11e and f. On the other hand, in the presence of sufficiently high friction, the contact is maintained all over the lower surface of the upper beam, as shown in Fig. 14.

As in the dead load case, the contact state in the presence of friction is rather complex. It is found that at every step there is slip at most nodes, except the last step where there is more stick than slip.

Fig. 15 shows the deflection at point A of the upper beam versus pressure resultant *P*. Note that under the follower loading there is a limit point for the deflection, the maximum absolute value for the deflection being equal to about three-quarter of the beam length.

### 6.4. Contact between two membranes

This example relating to the contact between two pressurized membranes involves large displacement, large strain and large slip as well. The problem is modeled in plane strain with respect to the

<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Contact between two beams. Deflection at point A versus dead load *F* (*I* = the second moment of area of the upper beam). ○: no friction, ●: with friction.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Fig. 13. Contact between two beams. Follower load and frictionless case: (a) P = 1.607 N, (b) P = 6.757 N, (c) P = 18.88 N, (d) P = 33.35 N, (e) P = 46.44 N, and (f) P = 73.84 N. The follower load is not drawn to scale.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 14. Contact between two beams. Follower load and frictional case $\mu = 0.3$: (a) P = 1.757 N, (b) P = 6.743 N, (c) P = 20.38 N, (d) P = 36.97 N, (e) P = 53.79 N, and (f) P = 74.83 N. The follower load is not drawn to scale.</figcaption></figure>

x-y plane. In the reference configuration, both membranes are 1 mm thick and represented by line segments with a width of 1 m along x. The membranes are $\Delta y = 1$ m apart in the direction

y and their centers are shifted by a distance $\Delta x = 1$ m along x, as shown in Fig. 16. The edges of the membranes parallel to z-axis are fixed.
<figure><img src="image_16.png" /><figcaption>Fig. 15. Contact between two beams. Follower load: deflection at point A versus pressure resultant P, ◦: no friction, •: with friction.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 16. Reference configuration of the two membrane problem.</figcaption></figure>

The membranes are made of an isotropic hyperelastic material of the neo-Hookean type, defined by the strain energy function in terms of the invariants $I_1 = \text{tr}C$ and $I_3 = \det C$ of the right Cauchy-Green tensor $C$ [31]:

$$w(\mathbf{C}) = \frac{\bar{\mu}}{2}(I_1 - 3) - \frac{\bar{\mu}}{2}\ln I_3 + \frac{\lambda}{8}(\ln I_3)^2. \quad (64)$$

We take $E = 2.10^6$ N/m$^2$, $\nu = 0.3$, where Young's modulus $E$ and Poisson's ratio $\nu$ are related to constants $\lambda$ and $\bar{\mu}$ by the same relations (59). The second Piola-Kirchhoff stress tensor $\Sigma$ is

<figure><img src="image_18.png" /><figcaption>Fig. 17. Contact between two pressurized membranes. Frictionless case: (a) contact onset, $p = 2049$ N/m$^2$, (b) deformed configuration at pressure $p = 4023$ N/m$^2$.</figcaption></figure>

$$\Sigma = 2 \frac{\partial w}{\partial \mathbf{C}} = \mu (\mathbf{I} - \mathbf{C}^{-1}) + \frac{\lambda}{2} \ln (\det \mathbf{C}) \cdot \mathbf{C}^{-1}. \quad (65)$$

The membranes are subjected to the same pressure value $p$ which brings them into contact. The computations are carried out assuming either no or with friction, $\mu = 0.5$, a rather significant coefficient of friction in order to simulate the contact between rubber membranes.

Each membrane is divided into 20 3-node line elements of the membrane type (i.e. with no bending stiffness), in plane strain with respect to $x-y$ plane and plane stress with respect to the membrane thickness direction. Here, the mesh of the contactor (upper membrane) is also used as the mesh of the contact surface. In order to avoid the singularity of the stiffness matrix at the beginning of the computation-which is a common phenomenon in membrane problems-a small artificial prestress is applied at the first step and removed after.

The numerical quadrature is carried using 3 Gauss points in the line elements to compute the internal force vector, 6 Gauss points to compute the contact force vector. The penalty parameters are $\epsilon_N = \epsilon_T = 2.10^4$ N/m$^3$. The pressure $p$ is controlled by the arc length method with the arc length regularly incremented by 1.05. The convergence tolerances are $10^{-6}$ for the force residual $\{\mathbf{R}\}$, $10^{-6}$ for the arc length residual and $10^{-8}$ for the multiplier residual $\{\mathbf{R}_1\}$.

Without friction, the computation is made in 43 steps, the maximum pressure value is $p = 4023$ N/m$^2$ and the convergence achieved in 4 iterations in average. With friction, the computation is made in 124 steps, the maximum pressure value is $p = 5420$ N/m$^2$ and the convergence achieved in 4 iterations in average.

Fig. 17 shows the deformed configurations and the contact stresses exerted on the upper membrane in the frictionless case. Let us mark out two particular particles belonging respectively to the two membranes, identified by letters A and B in Fig. 17. At the contact onset ($p = 2049$ N/m$^2$), their locations $A_1$ and $B_1$ are almost the same, see Fig. 17a. When the pressure is $p = 4023$ N/m$^2$, their locations $A_2$ and $B_2$ are quite distinct as indicated in Fig. 17b. The distance $A_2B_2$ is a (poor) approximation of the slip $\mathbf{V}_T\Delta t$, computed at the particle marked by $A_1$ and between the two time instants corresponding to the two quoted pressures, see (47) and Fig. 2.

Fig. 18 shows the deformed configurations and the contact stresses exerted on the upper membrane in the frictional case. Note that the contact stresses in Figs. 17 and 18 are not drawn at the same scale in order to optimize the drawings.

While the shape of final deformed configuration is similar to that in the frictionless case, the displacements of the particles are quite different. To see this, consider the above-mentioned pair of particles again, with almost identical locations $A_1$ and $B_1$ at the beginning of the contact (pressure $p = 2049$ N/m$^2$). This time, due to significant friction, it is found that, up to pressure

<figure><img src="image_19.png" /></figure>
<figure><img src="image_20.png" /><figcaption>**Fig. 18.** Contact between two pressurized membranes. Frictional case: (a) contact onset $p = 2049$ N/m$^2$, (b) deformed configuration at pressure $p = 5420$ N/m$^2$.</figcaption></figure>

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

<figure><img src="image_22.png" /><figcaption>**Fig. 19.** Contact between two pressurized membranes. Deflection of point C versus the pressure.</figcaption></figure>

$p = 5420$ N/m$^2$ corresponding to an advanced deformed shape, these two particles almost stick together, as shown in Fig. 18b. Otherwise stated, the slip velocity at the particle marked by $A_1$ is virtually zero.

An in-depth analysis at different steps shows that there is stick in the neighborhood of particle $A_1$, and there may be slip at contact nodes outside this neighborhood.

Fig. 19 plots the deflection of the middle point of membrane 2 (point C in Fig. 16) versus pressure $p$. As predicted, the deflection is much smaller in the presence of friction. The resultant contact forces $F_x$ and $F_y$ between the membranes along x- and y-direction, respectively, together with the magnitude $|F| = \sqrt{F_x^2 + F_y^2}$, are shown in Fig. 20 versus pressure $p$. It can be seen that for a given pressure the magnitude $|F|$ is a little smaller in frictional than in frictionless case whereas the force direction, determined by $F_x$ and $F_y$, differs significantly.

### 6.5. Conical extrusion of a cylinder

This section deals with the quasistatic extrusion of a cylinder into a rigid conical die, Fig. 21. The geometry data are taken from [1,3]: the wall angle of the die is 5° and in the reference configuration, the cylinder is 25.4 cm long and 5.08 cm in radius.

In this example, use should be made of a suitable constitutive model in order to overcome computational difficulties due to extremely high compressive stresses in the bodies at some advanced deformed stage. In the classical elastoplasticity model, the elastic response is represented by the Saint-Venant Kirchhoff law which has an unstable behavior under compression, so that in the present example the cylinder is found unable to move forward when it is compressed enough. Modification of the plastic part does not help to improve the situation and this difficulty only disappears when the Saint-Venant Kirchhoff law is replaced with an appropriate one, e.g., by using the hyperelastic stored energy given in [3] p37. Here, we choose the so-called Mooney rate-independent elastoplasticity model instead, which is described in [12] and proved to

<figure><img src="image_23.png" /><figcaption>**Fig. 20.** Contact between two pressurized membranes. Contact forces $F_x, F_y$ and $|F|$ versus the pressure.</figcaption></figure>
behave correctly even for huge deformations both in compression and tension. This model is chosen merely because it is easier to be implemented in our finite element code. The question about which elastoplasticity model best predicts the extrusion response when compared to experimental results is crucial but out of the scope of this work. Here, the extrusion computations using the chosen

elastoplasticity model only aim to show the ability of the formula-
tion to successfully deal with contact in large deformation with a
nonlinear constitutive law of the elastoplastic type.

The cylinder is made of an elastoplastic material with a linear isotropic hardening. The elastoplastic formulation and the material data are the same as in [12]: Young's modulus $E = 100000$ MPa, Poisson's ratio $\nu = 0.3$, yield stress $\sigma_0 = 300$ MPa and hardening modulus $H = 3560$ MPa. Note that the geometry in [12] is not considered since it is not described with enough precision (the die is not conical and its geometry is rather complex).

The problem is modeled in axisymmetry. The cylinder (taken as the contactor) is discretized with 20 8-node quadrilaterals (Q8) interconnected by 74 nodes. The 10 elements on the contact surface are 3-node line elements. The rigid die (the target) is modeled with one single Q8 element. The numerical quadrature is performed using 3 × 3 Gauss quadrature points in the Q8 elements, 7 Gauss points in line elements belonging to the contact surface. The computation are made assuming either frictionless or frictional contact between the cylinder and the die wall ($\mu = 0.1$). The penalty parameters are $\epsilon_N = 10^{11}$ N/m$^3$, $\epsilon_T = 5.10^{10}$ N/m$^3$.

<figure><img src="image_24.png" /><figcaption>Fig. 21. Reference configuration of the conical extrusion problem.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Fig. 22. Deformed configurations and plastic contours for the conical extrusion problem. Frictionless case.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Fig. 23. Contact stresses on the extruded cylinder. Frictionless case.</figcaption></figure>
The numerical computation is carried out by regularly incre-
menting the prescribed axial displacement on the left end of the
cylinder by factor 1.1. The convergence tolerances are 10⁻⁶ for
the force residual {R} and the multiplier residual {R₄}.

In the frictionless case, the computation is achieved in 31 steps,
the axial displacement of the pushed end is 18.2 cm at the last step
and the convergence requires 4-5 iterations in average. Fig. 22 dis-
plays the deformed configurations together with the contours of
the effective plastic strain while Fig. 23 shows the contact stresses
exerted on the extruded cylinder. The plastic strain is concentrated
near the cylinder head (where the contact stresses are highest, too)
and reaches about 130% at some points.

With friction, the computation is achieved in 109 steps, the
axial displacement at the last step is 20.4 cm and the conver-
gence requires 4 iterations in average. Fig. 24 displays the de-
formed configurations together with the contours of the
effective plastic strain while Fig. 25 shows the contact stresses.

From the numerical results it is found that all the nodes on
the contact line are sliding at every loading step and the equality
$|T_T| = \mu \cdot T_N$ is accurately rendered to the last displayed digit
(6th digit).

Note that at the final step the contact stresses near the cylinder
head increase more than twice as compared with the frictionless
case, as shown in Fig. 23d and Fig. 25d. The plastic strain is concen-
trated near the cylinder head where the contact stresses are stron-
gest, with a maximal value of about 150%.

Fig. 26 plots the extrusion force *F* versus the displacement *U* of the left end of the cylinder, in frictional and frictionless cases. As can be predicted, the extrusion force increases as the cylinder is moving forward and the contact surface between the cylinder and the die is larger. Also, given an amount prescribed displacement *U*, *force F* is always higher in frictional than frictionless case. When the end displacement *U* is about 18 cm, *F* is about 4 times as big as in the frictionless case. The extrusion force *F* in the frictional

<figure><img src="image_27.png" /><figcaption>Fig. 24. Deformed configurations and plastic contours for the conical extrusion problem. Frictional case, $\mu = 0.1$.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Fig. 25. Contact stresses on the extruded cylinder. Frictional case, $\mu = 0.1$.</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>Fig. 26. Thrust force versus the axial displacement of the cylinder.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>Fig. 27. Volume change in the extruded cylinder.</figcaption></figure>

case plotted in Fig. 26 shows the same trend as that obtained in [3]
whereas the numerical values differ rather significantly since the
material model and the material data are not identical.

Fig. 27 shows that the volume of the cylinder diminishes grad-
ually during the extrusion process and the relative volume change
is maximal in absolute value at the last step, equal to −6.8% in the
frictionless case and −15.9% in the frictional case. The volume de-
crease is basically the same in both cases as long as the prescribed
displacement is less than about half the reference length of the cyl-
inder. When the prescribed displacement exceeds this value, the
volume decreases more rapidly in the presence of friction.

7. Conclusions

The main features of the present formulation can be summa-
rized as follows:

* A weighted residual relationship is proposed in (21) as an extension of the standard virtual work principle to cope with the large deformation contact problem with Coulomb friction. It is inspired by the augmented Lagrangian formulation of

Pietrzak and Curnier [13] and amounts to rewrite the aug-
mented Lagrangian formulation of Alart and Curnier [4] in a
generic weak form. Proposition 1 shows that the proposed
weak form is equivalent to the strong form of the contact
problem. It is noteworthy that equality (21) yields the inequal-
ity constraints of contact. Another interesting fact is that Rela-
tion (23) which expresses the equilibrium of the traction
vectors at the contact interface is found to be a consequence
of the weak form (21).Apart from the analogy with an aug-
mented Lagrangian principle, the weak form (21) is stated as
a genuine weighted residual relationship in that it involves
arbitrary (regular) virtual displacements and multipliers. As a
consequence, it can be discretized by means of the finite ele-
ment method in a simple way.

* The finite element discretization process holds for any finite element type on the contact surface. It is only required that the finite elements are isoparametric, so that the current and reference coordinates can be interpolated in the same way.

* Any discrete contact quantity – namely, the element contact force vector {$\Phi_{contact}$}^e in (43), the element residual vector for the multipliers {$R_A$}^e in (44) or the element contact stiffness matrice {$K_{contact}$}^e in (56) – is defined at a quadrature point $\mathbf{X}$ of an element $e^{(1)}$ belonging to the contactor surface $S_{oc}^{(1)}$. Its numerical computation involves the nodal quantities associated to two elements: element $e^{(1)}$ and the (or any) element $e^{(2)}$ in the target surface $S_{oc}^{(2)}$ containing the projection point $\mathbf{Y}$ of the quadrature point $\mathbf{X}$ in question. Element $e^{(2)}$, which depends on $\mathbf{Y}$, varies as a function of $\mathbf{X}$. This feature is common to the formulation given by Laursen [3].

* The temporally discrete relation (47) proposed for the slip $V_T \Delta t$ does not involve the parent coordinate $\eta$ of the projection point $\mathbf{Y}$ and thus avoids difficulties due to the discontinuity in variable $\eta$ related to the finite discretization of $S_{oc}^{(2)}$. Instead, it involves variables $\mathbf{y}$ and $\mathbf{Y}_{n-1}$ which have intrinsic kinematic meanings independent of the finite discretization of $S_{oc}^{(2)}$, as clearly shown in Fig. 2.

* The element contact stiffness matrices related to relative tangential displacement, [$k_{contact}^4$]$^e$ and [$k_{contact}^5$]$^e$ in (A6) and (A7), are rectangular since they also involve the projection point $\mathbf{Y}_{n-1}$ determined at the previous time step. Nevertheless, their assembly eventually yields the global contact stiffness matrix which is square as usual.

* As regards the smoothness of the kinematic functions, it should be noted that the present formulation involves differentiations of order 2 at most (see, for instance, the definition of the curvature $\kappa_{x\theta}$ in (14)), whereas other formulations may require higher order differentiations.

* As the weak form proposed in the present formulation already incorporates the contact laws, no integration of the frictional equations are required on the local level of the numerical solution procedure. This feature resembles, in some ways, an implicit property in Curnier and co-workers' formulations [4,13] and differs from most of other developments in the literature.

* The numerical examples presented in quasistatistics shows the ability of the proposed formulation to deal with large deformation contact. Investigations are in progress to deal with numerical examples in dynamic contact.

Acknowledgments

The authors would like to thank Prof. Nguyen Quoc Son, of Lab-
oratoire de Mecanique des Solides, Ecole Polytechnique in France,
for his instructive remarks which have enabled them to refine
Proposition 1.
# Appendix. The element contact stiffness matrix

This Appendix gives the explicit expressions for the four blocks of the element contact stiffness matrix (56), namely $[\mathbf{K}_{UU}]^e$, $[\mathbf{K}_{UA}]^e$, $[\mathbf{K}_{AU}]^e$ and $[\mathbf{K}_{AA}]^e$.

(1) The rectangular matrix $[\mathbf{K}_{UU}]^e$ in (56) pertains to the displacements $\{\mathbf{U}^{(1)}\}^{e(1)}$ on $e^{(1)}$, $\{\mathbf{U}^{(2)}\}^{e(2)}$ on $e^{(2)}$ and $\{\mathbf{U}^{(2)}\}^{e_{n-1}^{(2)}}$ on $e^{(2)}$. It is of dimension $(3nne^{(1)} + 3nne^{(2)}) \times (3nne^{(1)} + 3nne^{(2)} + 3nne^{(2)}_{n-1})$ where $nne_{n-1}^{(2)}$ is the number of nodes of element $e_{n-1}^{(2)}$, and is given by

$$[\mathbf{K}_{UU}]^e = H(\hat{\lambda}_N) \epsilon_N [\mathbf{k}_{\text{contact}}^1]^e + \langle \hat{\lambda}_N \rangle [\mathbf{k}_{\text{contact}}^2]^e - H\left(1 - \frac{\mu(\hat{\lambda}_N)}{\|\hat{\lambda}_T\|}\right) \mu H(\hat{\lambda}_N) \epsilon_N [\mathbf{k}_{\text{contact}}^3]^e + H\left(1 - \frac{\mu(\hat{\lambda}_N)}{\|\hat{\lambda}_T\|}\right) \mu \langle \hat{\lambda}_N \rangle \epsilon_T [\mathbf{k}_{\text{contact}}^4]^e + \left(1 - \langle 1 - \frac{\mu(\hat{\lambda}_N)}{\|\hat{\lambda}_T\|} \rangle\right) (\epsilon_T [\mathbf{k}_{\text{contact}}^5]^e + [\mathbf{k}_{\text{contact}}^6]^e), \quad (A1)$$

where symbol $H$ stands for the Heaviside function:

$$H\left(1 - \frac{\mu(\hat{\lambda}_N)}{\|\hat{\lambda}_T\|}\right) = \begin{cases} 1 & \text{if } \hat{\lambda}_N < 0 \quad (\text{algorithmic gap}) \\ 1 & \text{if } \hat{\lambda}_N \ge 0 \text{ and } \|\hat{\lambda}_T\| \ge \mu\hat{\lambda}_N \quad (\text{algorithmic contact with slip}) \\ 0 & \text{if } \hat{\lambda}_N \ge 0 \text{ and } \|\hat{\lambda}_T\| < \mu\hat{\lambda}_N \quad (\text{algorithmic contact with stick}) \end{cases} \quad (A2)$$

and the matrices $[\mathbf{k}_{\text{contact}}^1]^e$ to $[\mathbf{k}_{\text{contact}}^6]^e$ are listed in the sequel. The matrices $[\mathbf{k}_{\text{contact}}^1]^e$ and $[\mathbf{k}_{\text{contact}}^6]^e$ are square, of dimension $3nne^{(1)} + 3nne^{(2)}$ and symmetric:

$$[\mathbf{k}_{\text{contact}}^1]^e = \begin{bmatrix} [\mathcal{N}^{(1)}]^{e(1)T} \{\boldsymbol{\nu}\} < \boldsymbol{\nu} > [\mathcal{N}^{(1)}]^{e(1)} \\ -[\mathcal{N}^{(2)}]^{e(2)T} \{\boldsymbol{\nu}\} < \boldsymbol{\nu} > [\mathcal{N}^{(1)}]^{e(1)} \end{bmatrix} \begin{bmatrix} -[\mathcal{N}^{(1)}]^{e(1)T} \{\boldsymbol{\nu}\} < \boldsymbol{\nu} > [\mathcal{N}^{(2)}]^{e(2)} \\ [\mathcal{N}^{(2)}]^{e(2)T} \{\boldsymbol{\nu}\} < \boldsymbol{\nu} > [\mathcal{N}^{(2)}]^{e(2)} \end{bmatrix}, \quad (A3)$$

$$[\mathbf{k}_{\text{contact}}^2]^e = \begin{bmatrix} [\mathcal{N}^{(1)}]^{e(1)T} [\kappa_{\alpha\gamma} \bar{a}^{\gamma\beta} \mathbf{a}^{\alpha} \otimes \mathbf{a}_{\beta}] [\mathcal{N}^{(1)}]^{e(1)} & -[\mathcal{N}^{(1)}]^{e(1)T} [\kappa_{\alpha\gamma} \bar{a}^{\gamma\beta} \mathbf{a}^{\alpha} \otimes \mathbf{a}_{\beta}] [\mathcal{N}^{(2)}]^{e(2)} \\ -[\mathcal{N}^{(2)}]^{e(2)T} [\kappa_{\alpha\gamma} \bar{a}^{\gamma\beta} \mathbf{a}^{\alpha} \otimes \mathbf{a}_{\beta}] [\mathcal{N}^{(1)}]^{e(1)} & +[\mathcal{N}^{(1)}]^{e(1)T} [\bar{a}^{\alpha\beta} \mathbf{a}_{\beta} \otimes \boldsymbol{\nu}] [\mathcal{N}^{(2)}]^{e(2)} \\ +[\mathcal{N}^{(2)}, \alpha]^{e(2)T} [\bar{a}^{\alpha\beta} \boldsymbol{\nu} \otimes \mathbf{a}_{\beta}] [\mathcal{N}^{(1)}]^{e(1)T} & [\mathcal{N}^{(2)}]^{e(2)T} [\kappa_{\alpha\gamma} \bar{a}^{\gamma\beta} \mathbf{a}^{\alpha} \otimes \mathbf{a}_{\beta}] [\mathcal{N}^{(2)}]^{e(2)} \\ & -[\mathcal{N}^{(2)}]^{e(2)T} [\bar{a}^{\alpha\beta} \mathbf{a}_{\beta} \otimes \boldsymbol{\nu}] [\mathcal{N}^{(2)}, \alpha]^{e(2)} \\ & -[\mathcal{N}^{(2)}, \alpha]^{e(2)T} [\bar{a}^{\alpha\beta} \boldsymbol{\nu} \otimes \mathbf{a}_{\beta}] [\mathcal{N}^{(2)}]^{e(2)} \\ & -[\mathcal{N}^{(2)}, \alpha]^{e(2)T} [g\bar{a}^{\alpha\beta} \boldsymbol{\nu} \otimes \boldsymbol{\nu}] [\mathcal{N}^{(2)}, \beta]^{e(2)} \end{bmatrix}, \quad (A4)$$

In (A4), the notation $[\mathbf{T}]$ for any second order tensor $\mathbf{T}$ designates the 3 × 3 matrix representing $\mathbf{T}$ in the fixed basis $(\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3)$. For instance, if $\mathbf{T} = T_{\alpha}^{\beta} \mathbf{a}^{\alpha} \otimes \mathbf{a}_{\beta} = T^{\alpha\beta} \mathbf{a}_{\alpha} \otimes \mathbf{a}_{\beta}$, then the $(p,q)$-th element of matrix $[\mathbf{T}]$ is $T_{pq} = \mathbf{e}_p \cdot \mathbf{T} \cdot \mathbf{e}_q = T_{\alpha}^{\beta}(\mathbf{a}^{\alpha})_p(\mathbf{a}_{\beta})_q = T^{\alpha\beta}(\mathbf{a}_{\alpha})_p(\mathbf{a}_{\beta})_q$, where $(\mathbf{a}^{\alpha})_p$ is the p-th component of $\mathbf{a}^{\alpha}$ in the fixed basis. As $\kappa_{\alpha\gamma} \bar{a}^{\gamma\beta} \mathbf{a}^{\alpha} \otimes \mathbf{a}_{\beta} = \bar{a}^{\alpha\gamma} K_{\gamma\beta} \mathbf{a}_{\alpha} \otimes \mathbf{a}^{\beta}$ is a symmetric tensor, so is matrix $[\mathbf{k}_{\text{contact}}^2]^e$.

Matrix $[\mathbf{k}_{\text{contact}}^3]^e$ in (A1) is not symmetric and given by

$$[\mathbf{k}_{\text{contact}}^3]^e = \begin{bmatrix} [\mathcal{N}^{(1)}]^{e(1)T} \frac{\{\hat{\lambda}_T\}}{\|\hat{\lambda}_T\|} < \boldsymbol{\nu} > [\mathcal{N}^{(1)}]^{e(1)} & -[\mathcal{N}^{(1)}]^{e(1)T} \frac{\{\hat{\lambda}_T\}}{\|\hat{\lambda}_T\|} < \boldsymbol{\nu} > [\mathcal{N}^{(2)}]^{e(2)} \\ -[\mathcal{N}^{(2)}]^{e(2)T} \frac{\{\hat{\lambda}_T\}}{\|\hat{\lambda}_T\|} < \boldsymbol{\nu} > [\mathcal{N}^{(1)}]^{e(1)} & [\mathcal{N}^{(2)}]^{e(2)T} \frac{\{\hat{\lambda}_T\}}{\|\hat{\lambda}_T\|} < \boldsymbol{\nu} > [\mathcal{N}^{(2)}]^{e(2)} \end{bmatrix}. \quad (A5)$$
The contribution of point $\eta_{n-1}$ is found in the rectangular matrices $[\mathbf{k}_{\text{contact}}]^e$ and $[\mathbf{k}_{\text{contact}}^5]^e$ of dimension $(3nne^{(1)} + 3nne^{(2)}) \times (3nne^{(1)} + 3nne^{(2)} + 3nne_{n-1}^{(2)})$:

$$ [\mathbf{k}_{\text{contact}}^4]^e = \left[ \begin{array}{c|c} -[\mathcal{N}^{(1)}]^{e(1)T} \frac{\{\hat{\lambda}_T\} <\hat{\lambda}_T>}{\|\hat{\lambda}_T\|^3} [\mathbf{W}][\mathcal{N}^{(1)}]^{e(1)} & [\mathcal{N}^{(1)}]^{e(1)T} \frac{\{\hat{\lambda}_T\} <\hat{\lambda}_T>}{\|\hat{\lambda}_T\|^3} \frac{[\mathcal{N}^{(1)}]^{e(1)T} \{\hat{\lambda}_T\} <\hat{\lambda}_T>}{\|\hat{\lambda}_T\|^3} [\mathcal{N}^{(2)}(\eta_{n-1})]^{e(2)}_{n-1} \\ -[\mathcal{N}^{(2)}]^{e(2)T} \frac{\{\hat{\lambda}_T\} <\hat{\lambda}_T>}{\|\hat{\lambda}_T\|^3} [\mathbf{W}][\mathcal{N}^{(1)}]^{e(1)} & -[\mathcal{N}^{(2)}]^{e(2)T} \frac{\{\hat{\lambda}_T\} <\hat{\lambda}_T>}{\|\hat{\lambda}_T\|^3} \frac{[\mathcal{N}^{(2)}]^{e(2)T} \{\hat{\lambda}_T\} <\hat{\lambda}_T>}{\|\hat{\lambda}_T\|^3} [\mathcal{N}^{(2)}(\eta_{n-1})]^{e(2)}_{n-1} \end{array} \right]. \quad (A6) $$

$$ [\mathbf{k}_{\text{contact}}^5]^e = \left[ \begin{array}{c|c} [\mathcal{N}^{(1)}]^{e(1)T} [\mathbf{W}][\mathcal{N}^{(1)}]^{e(1)} & [\mathcal{N}^{(1)}]^{e(1)T} \frac{[\mathbf{I} - \mathbf{W}][\mathcal{N}^{(2)}]^{e(2)} - [\mathbf{Z}^{\alpha}][\mathcal{N}_{,\alpha}^{(2)}]^{e(2)}}{[\mathcal{N}^{(2)}]^{e(2)T} [\mathbf{W} - \mathbf{I}][\mathcal{N}^{(2)}]^{e(1)}} & -[\mathcal{N}^{(1)}]^{e(1)T} [\mathcal{N}^{(2)}(\eta_{n-1})]^{e(2)}_{n-1} \\ -[\mathcal{N}^{(2)}]^{e(2)T} [\mathbf{W}][\mathcal{N}^{(1)}]^{e(1)} & [\mathcal{N}^{(2)}]^{e(2)T} \frac{[\mathcal{N}^{(2)}]^{e(2)T} [\mathbf{I} - \mathbf{W}][\mathcal{N}^{(2)}]^{e(2)} + [\mathbf{Z}^{\alpha}][\mathcal{N}_{,\alpha}^{(2)}]^{e(2)}}{[\mathcal{N}^{(2)}]^{e(2)T} [\mathbf{W} - \mathbf{I}][\mathcal{N}^{(2)}]^{e(1)}} & [\mathcal{N}^{(2)}]^{e(2)T} [\mathcal{N}^{(2)}(\eta_{n-1})]^{e(2)}_{n-1} \end{array} \right]. \quad (A7) $$

In (A6) and (A7), $[\mathbf{W}]$, $[\mathbf{Z}^\alpha]$ are 3 × 3 matrices representing the tensors $\mathbf{W} \cdot \mathbf{Z}^\alpha$ defined by (52), in the fixed basis $(\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3) : W_{ij} = \mathbf{e}_i \cdot \mathbf{W} \cdot \mathbf{e}_j = \bar{a}^{\alpha\beta}(\mathbf{a}^\alpha)_i(\mathbf{a}^\beta)_j, Z_{ij}^\alpha = \mathbf{e}_i \cdot \mathbf{Z}^\alpha \cdot \mathbf{e}_j = g\bar{a}^{\alpha\beta}(\mathbf{a}^\mu)_i v_{ij}$. The last matrix in (A1), $[\mathbf{k}_{\text{contact}}^6]^e$, is square of the same dimension as $[\mathbf{k}_{\text{contact}}^1]^e$, yet not symmetric:

$$ [\mathbf{k}_{\text{contact}}^6]^e = \left[ \begin{array}{c|c} [\mathbf{0}] & [\mathbf{0}] \\ -[\mathcal{N}^{(2)}, \alpha]^{e(2)T} \{\hat{\lambda}_T\} <\bar{a}^{\alpha\beta} \mathbf{a}_\beta> [\mathcal{N}^{(1)}]^{e(1)} & [\mathcal{N}^{(2)}, \alpha]^{e(2)T} \{\hat{\lambda}_T\} \\ (<\bar{a}^{\alpha\beta} \mathbf{a}_\beta> [\mathcal{N}^{(2)}]^{e(2)} + g <\bar{a}^{\alpha\beta} \boldsymbol{\nu}> [\mathcal{N}_{,\beta}^{(2)}]^{e(2)}) & [\mathcal{N}^{(2)}]^{e(2)T} [\mathcal{N}^{(2)}(\eta_{n-1})]^{e(2)}_{n-1} \end{array} \right]. \quad (A8) $$

(2) The matrices $[\mathbf{K}_{\text{U}U}]^e$ and $[\mathbf{K}_{\text{U}U}]^e$ in (56) representing the coupling between the displacements and the multipliers are rectangular, of dimension $(3nne^{(1)} + 3nne^{(2)}) \times 3nne^{(1)}$ and $3nne^{(1)} \times (3nne^{(1)} + 3nne^{(2)})$, respectively:

$$ [\mathbf{K}_{\text{U}U}]^e = \left[ \begin{array}{c|c} [\mathcal{N}^{(1)}]^{e(1)T} & [-H(\hat{\lambda}_N)\{\mathbf{v}\}\langle\mathbf{v}\rangle + H\left(1 - \frac{\mu\langle\hat{\lambda}_N\rangle}{\|\hat{\lambda}_T\|}\right)\mu H(\hat{\lambda}_N) \frac{\{\hat{\lambda}_T\}}{\|\hat{\lambda}_T\|} \langle\mathbf{v}\rangle \\ -[\mathcal{N}^{(2)}]^{e(2)T} & +H\left(1 - \frac{\mu\langle\hat{\lambda}_N\rangle}{\|\hat{\lambda}_T\|}\right)\mu\langle\hat{\lambda}_N\rangle \frac{\{\hat{\lambda}_T\}}{\|\hat{\lambda}_T\|}^3 \langle\mathbf{v}\rangle - \left(1 - \left\langle 1 - \frac{\mu\langle\hat{\lambda}_N\rangle}{\|\hat{\lambda}_T\|}\right\rangle\right) [\mathbf{I} - \mathbf{v} \otimes \mathbf{v}] \right] [\mathcal{N}^{(1)}]^{e(1)} \quad (A9) $$

and

$$ [\mathbf{K}_{\text{U}U}]^e = [\mathcal{N}^{(1)}]^{e(1)T} \left[ -H(\hat{\lambda}_N) \{\mathbf{v}\} \langle\mathbf{v}\rangle \left\langle [\mathcal{N}^{(1)}]^{e(1)}, -[\mathcal{N}^{(2)}]^{e(2)} \right\rangle + H\left(1 - \frac{\mu(\hat{\lambda}_N)}{\|\hat{\lambda}_T\|}\right) \mu H(\hat{\lambda}_N) \frac{\epsilon_N}{\epsilon_T} \frac{\{\hat{\lambda}_T\}}{\|\hat{\lambda}_T\|} \langle\mathbf{v}\rangle \right. \\ \times \left\langle [\mathcal{N}^{(1)}]^{e(1)}, -[\mathcal{N}^{(2)}]^{e(2)} \right\rangle + H\left(1 - \frac{\mu(\hat{\lambda}_N)}{\|\hat{\lambda}_T\|}\right) \mu \langle\hat{\lambda}_N\rangle \frac{\{\hat{\lambda}_T\}}{\|\hat{\lambda}_T\|} \left( [\mathbf{W}] [\mathcal{N}^{(1)}]^{e(1)}, [\mathbf{I} - \mathbf{W}] [\mathcal{N}^{(2)}]^{e(2)} + [\mathbf{Z}^\alpha] [\mathcal{N}_{,\alpha}^{(2)}]^{e(2)} \right. \\ \left. - [\mathcal{N}^{(2)}(\eta_{n-1})]^{e(2)}_{n-1} - \left(1 - \frac{\mu(\hat{\lambda}_N)}{\|\hat{\lambda}_T\|}\right) [\mathbf{I} - \mathbf{v} \otimes \mathbf{v}] \left( [\mathbf{W}] [\mathcal{N}^{(1)}]^{e(1)}, [\mathbf{I} - \mathbf{W}] [\mathcal{N}^{(2)}]^{e(2)} + [\mathbf{Z}^\alpha] [\mathcal{N}_{,\alpha}^{(2)}]^{e(2)} \right. \right. \\ \left. - [\mathcal{N}^{(2)}(\eta_{n-1})]^{e(2)}_{n-1} \right] \right]. \quad (A10) $$
(3) Lastly, the matrix $[K_{AA}]^e$ in (56) is square of order $3nne^{(1)}$:

$$[K_{AA}]^e = [\mathcal{N}^{(1)}]^{e(1)T} \left[ \frac{1}{\epsilon_N} (H(\hat{\lambda}_N) - 1) \{\mathbf{v}\} \langle \mathbf{v} \rangle - H \left(1 - \frac{\mu(\hat{\lambda}_N)}{\|\hat{\lambda}_T\|} \right) \mu H(\hat{\lambda}_N) \frac{1}{\epsilon_T} \frac{\mu(\hat{\lambda}_N)}{\|\hat{\lambda}_T\|} \langle \mathbf{v} \rangle - \frac{\mu(\hat{\lambda}_N)}{\|\hat{\lambda}_T\|} \frac{\mu}{\epsilon_T} \langle \hat{\lambda}_N \rangle \frac{\{\hat{\lambda}_T\} \langle \hat{\lambda}_T \rangle}{\|\hat{\lambda}_T\|^3} \right. \\ \left. - \left\langle 1 - \frac{\mu(\hat{\lambda}_N)}{\|\hat{\lambda}_T\|} \right\rangle \frac{1}{\epsilon_T} [\mathbf{I} - \mathbf{v} \otimes \mathbf{v}] \right] [\mathcal{N}^{(1)}]^{e(1)} \quad (A11)$$

## References

[1] Kikuchi N, Oden JT. A treatise on variational inequalities and finite element methods for contact problems in elasticity. Philadelphia: SIAM Publications; 1988.

[2] Wriggers P. Computational contact mechanics. Wiley; 2002.

[3] Laursen TA. Computational contact and impact mechanics – fundamentals of modeling interfacial phenomena in nonlinear finite element analysis. Springer; 2003.

[4] Alart P, Curnier A. A mixed formulation for frictional contact problems prone to Newton like solution methods. Comput Meth Appl Mech Eng 1991;92:353-75.

[5] Zavarise G, Wriggers P, Schrefler BA. On augmented Lagrangian algorithms for thermomechanical contact problems with friction. Int J Numer Meth Eng 1995;38:2929-49.

[6] Glowinski R, Le Tallec P. Augmented lagrangian and operator-splitting methods in nonlinear mechanics. SIAM, Studies in Applied Mathematics; 1989.

[7] Heegaard J-H, Curnier A. An augmented Lagrangian method for discrete large-slip contact problems. Int J Numer Meth Eng 1993;36:569-93.

[8] Zavarise G, Wriggers P. A superlinear convergent augmented Lagrangian procedure for contact problems. Eng Comput 1999;16:88-119.

[9] Simo JC, Laursen TA. An augmented Lagrangian treatment of contact problems involving friction. Comput Struct 1992;42:97-116.

[10] Laursen TA, Simo JC. Algorithmic symmetrization of Coulomb frictional problems using augmented Lagrangians. Comput Meth Appl Mech Eng 1993;108:133-46.

[11] Wriggers P. Finite element methods for contact problems with friction. Tribol Int 1996;29:651-8.

[12] Pietrzak G. Continuum mechanics modelling and augmented Lagrangian formulation of large deformation frictional contact problems. Ph.D. Thesis. Ecole Polytechnique FTdTrale de Lausanne; 1997.

[13] Pietrzak G, Curnier A. Large deformation frictional contact mechanics: continuum formulation and augmented Lagrangian treatment. Comput Meth Appl Mech Eng 1999;177:351-81.

[14] Feng Z-Q, Peyraut F, Labed N. Solution of large deformation contact problems with friction between Blatz-Ko hyperelastic bodies. Int J Eng Sci 2003;41:2213-25.

[15] Mijar AR, Arora JS. An augmented Lagrangian optimization method for contact analysis problems, 1: formulation and algorithm. Struct Multidiscip Optim 2004;28:99-112.

[16] Mijar AR, Arora JS. An augmented Lagrangian optimization method for contact analysis problems, 2: numerical evaluation. Struct Multidiscip Optim 2004;28: 113-26.

[17] Giannakopoulos AE. The return mapping method for the integration of friction constitutive relations. Comput. Struct. 1989;32:157-67.

[18] Bathe K-J. Finite element procedures. Prentice Hall; 1996.

[19] Curnier A, He Qi-Ch, Klarbring A. Continuum mechanics modelling of large deformation contact with friction. In: Raous M, Jean M, Moreau JJ, editors. Contact mechanics. New York: Plenum; 1995. p. 145-58.

[20] Klarbring A. Large displacement frictional contact: a continuum framework for finite element discretisation. Eur J Mech A-Solid 1995;14:237-53.

[21] Le van A, Nguyen THT. Une formulation variationnelle du problème de contact avec frottement de Coulomb. Comptes Rendus Mécanique (Académie des Sciences, Paris) 2008;33:606-11.

[22] Moreau JJ. Unilateral contact and dry friction in finite freedom dynamics. In: Moreau JJ, Panagiotopoulos PD, editors. Nonsmooth mechanics and applications CISM courses and lectures. Springer Verlag; 1988. p. 1-82.

[23] Chapelle D, Bathe K-J. The inf-sup test. Comput Struct 1993;47:537-45.

[24] El Abbasi N, Bathe K-J. Stability and patch test performance of contact discretizations and a new solution algorithm. Comput Struct 2001;79: 1473-86.

[25] Bathe K-J, Brezzi F. Stability of finite element mixed interpolations for contact problems. Rendiconti Lincei – Matematica E Applicazioni 1991;12: 167-83.

[26] Simo JC, Wriggers P, Taylor RL. A perturbed Lagrangian formulation for the finite element solution of contact problems. Comput Meth Appl Mech Eng 1985;50:163-80.

[27] Taylor RL, Papadopoulos P. On a patch test for contact problems in two dimensions. In: Wriggers P, Wagner W, editors. Computational methods in nonlinear mechanics. Berlin: Springer; 1991. p. 690-702.

[28] Crisfield MA. Re-visiting the contact patchtest. Int J Numer Meth Eng 2000;48:435-49.

[29] Spence DA. An eigenvalue problem for elastic contact with finite friction. Proc Camb Philos Soc 1973;73:249-68.

[30] Johnson KL. Contact mechanics. Cambridge University Press; 1985.

[31] Bonet J, Wood RD. Nonlinear continuum mechanics for finite element analysis. Cambridge University Press; 1997.