# Numerical implementation of two nonconforming finite element methods for unilateral contact

Patrick Hild

► To cite this version:

Patrick Hild. Numerical implementation of two nonconforming finite element methods for unilateral contact. Computer Methods in Applied Mechanics and Engineering, 2000, 184 (1), pp.99 - 123. 10.1016/S0045-7825(99)00096-1. hal-01390457

HAL Id: hal-01390457

https://hal.science/hal-01390457

Submitted on 2 Nov 2016

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Numerical implementation of two nonconforming finite
element methods for unilateral contact

Patrick Hild

Mathématiques pour l'Industrie et la Physique, Unité Mixte de Recherches CNRS UPS INSAT (U.M.R. 5640), Université Paul Sabatier, 118 route de Narbonne, 31062 Toulouse Cedex 4, France

We consider the finite element approximation of the unilateral contact problem between elastic bodies. We are interested in a practical problem which often occurs in finite element computations concerning two independently discretized bodies in unilateral contact. It follows that the nodes of both bodies located on the contact surface do not fit together. We present two different approaches in order to define unilateral contact on nonmatching meshes. The first is an extension of the mortar finite element method to variational inequalities that defines the contact in a global way. On the contrary, the second one expresses local node on segment contact conditions. In both cases, the theoretical approximation properties are given. Then, we implement and compare the two methods. © 2000 Elsevier Science S.A. All rights reserved.

Keywords: Unilateral contact; Nonmatching meshes; Mortar finite element method; Global contact condition; Local contact condition

# 1. Introduction and notations

In finite element procedures solving unilateral contact problems between deformable bodies, each solid is often discretized independently of the other. So, a finite element mesh does not coincide with the other one on the contact zone. For instance, constitutive or geometrical nonlinearities lead to the use of an incremental scheme and, at each step, the updated finite element meshes cannot fit together on the contact zone. In other respects, mesh adaptivity procedures used in a contact context generally lead to nonmatching meshes. So the question is to define a convenient discrete contact condition for nonmatching meshes.

On the one hand, the mathematical framework associated with the finite element approximation for contact problems can be found in [16]. In this reference, Haslinger et al. consider the case of matching meshes on the contact zone.

On the other hand, the mortar element domain decomposition method, introduced by Bernardi et al. [8], allows the handling of nonmatching meshes. This technique has been studied for many problems governed by variational equalities.

The first extension of the mortar finite element method to a variational inequality has been made for the unilateral contact problem by Ben Belgacem et al. [5]. In this reference, the authors also extend the finite element analysis of contact problems to nonmatching meshes.
The main aim of this paper is to carry out the first numerical experiments associated with the theoretical results already obtained. The paper is organized as follows. First, we introduce the model describing the unilateral contact without friction between two deformable elastic bodies. The associated weak formulation is exhibited.

In Section 3 we consider two finite element methods in order to solve the problem with independent meshes. The first approach is of global type and corresponds to an extension of the mortar domain decomposition method. The second approximation is of local type and uses classical node-on-segment contact conditions. We recall the most significant approximation properties obtained in previous papers and we establish new results proving the optimal convergence of the global approach and the convergence of the local approach.

In Section 4 we obtain the corresponding matrix formulations and we mention the algorithm contained in the finite element code CASTEM 2000. Section 5 is devoted to the studies in which the contact conditions are compared from a numerical point of view.

## 2. Setting of the problem

We consider two elastic bodies occupying in the initial configuration two subsets $\bar{\Omega}^\ell$ of $\mathbb{R}^2$, $\ell = 1, 2$. The boundary $\partial\Omega^\ell$ of the domain $\Omega^\ell$ is assumed to be “smooth” and consists of three nonoverlapping parts $\Gamma_u^\ell$, $\Gamma_g^\ell$ and $\Gamma_c^\ell$. The unit outward normal on $\partial\Omega^\ell$ is denoted $\mathbf{n}^\ell$. The body is submitted to volume forces $\mathbf{f}^\ell$ on $\Omega^\ell$ and to surface forces $\mathbf{g}^\ell$ on $\Gamma_g^\ell$. On $\Gamma_u^\ell$, the displacements $\mathbf{U}^\ell$ are prescribed. In the initial configuration, both bodies have a common portion $\Gamma_c = \Gamma_c^1 = \Gamma_c^2$ which will be considered as the candidate contact surface for the sake of simplicity.

The unilateral contact problem consists of finding the displacement field $\mathbf{u} = (\mathbf{u}^1, \mathbf{u}^2)$, where $\mathbf{u}^\ell = \mathbf{u}|_{\Omega^\ell}$ and the stress tensor field $\sigma = (\sigma^1, \sigma^2)$, $(\sigma^\ell = \sigma|_{\Omega^\ell})$ satisfying the following conditions (2.1) (2.5) for $\ell = 1, 2$:

$$ 
\begin{array}{ll}
\mathbf{div} \, \sigma^\ell(\mathbf{u}^\ell) + \mathbf{f}^\ell = 0 & \text{in } \Omega^\ell, \\
\sigma^\ell(\mathbf{u}^\ell)\mathbf{n}^\ell = \mathbf{g}^\ell & \text{on } \Gamma_g^\ell, \\
\mathbf{u}^\ell = \mathbf{U}^\ell & \text{on } \Gamma_u^\ell.
\end{array}
\qquad (2.1) $$

The symbol **div** denotes the divergence operator defined by $\mathbf{div} \, \sigma = (\partial\sigma_{ij}/\partial x_j)_i$, where the summation convention of repeated indices is adopted.

We will consider small strains hypothesis so that the strain tensor $\varepsilon(\mathbf{v})$ induced by a displacement field $\mathbf{v}$ is $\varepsilon(\mathbf{v}) = (\nabla\mathbf{v} + \mathbf{v}^T \nabla\mathbf{v})/2$. The stress tensor field $\sigma^\ell$ is linked to the displacement field $\mathbf{u}^\ell$ by the constitutive law of linear elasticity

$$ \sigma^\ell(\mathbf{u}^\ell) = C^\ell \varepsilon(\mathbf{u}^\ell), \qquad (2.2) $$

where $C^\ell = (c_{ij,kh}^\ell)_{1 \le i,j,k,h \le 2}$ is a fourth order tensor satisfying $c_{ij,kh}^\ell = c_{ji,kh}^\ell = c_{kh,ij}^\ell$. We assume that there exists constants $\alpha^\ell > 0$ verifying

$$ c_{ij,kh}^\ell \varepsilon_{ij} \varepsilon_{kh} \ge \alpha^\ell \varepsilon_{ij} \varepsilon_{ij} \quad \forall \varepsilon_{ij} = \varepsilon_{ji}. $$

The conditions on the contact zone $\Gamma_c$ are as follows:

$$ (\sigma^1(\mathbf{u}^1)\mathbf{n}^1) \cdot \mathbf{n}^1 = (\sigma^2(\mathbf{u}^2)\mathbf{n}^2) \cdot \mathbf{n}^2 = \sigma_n(\mathbf{u}), \qquad (2.3) $$

$$ \mathbf{u}^1 \cdot \mathbf{n}^1 + \mathbf{u}^2 \cdot \mathbf{n}^2 \le 0, \quad \sigma_n(\mathbf{u}) \le 0, \quad \sigma_n(\mathbf{u})(\mathbf{u}^1 \cdot \mathbf{n}^1 + \mathbf{u}^2 \cdot \mathbf{n}^2) = 0, \qquad (2.4) $$

$$ \sigma_T^1(\mathbf{u}^1) = \sigma_T^2(\mathbf{u}^2) = 0, \qquad (2.5) $$

where

$$ \sigma_T^\ell(\mathbf{u}^\ell) = \sigma^\ell(\mathbf{u}^\ell)\mathbf{n}^\ell - \sigma_n(\mathbf{u})\mathbf{n}^\ell, \quad 1 \le \ell \le 2. $$

Relations (2.3) represent the action and the reaction principle. Conditions (2.4) express unilateral contact between the two bodies and finally (2.5) states a contact without friction.
In order to obtain the variational formulation of the problem, we introduce the spaces $V(\Omega^\ell)$ ($\ell = 1, 2$)

$$V(\Omega^\ell) = \{ \mathbf{v} \in (H^1(\Omega^\ell))^2, \mathbf{v} = \mathbf{U}^\ell \text{ on } \Gamma_\mathbf{u}^\ell \},$$

where $H^1(\Omega^\ell)$ is the classical Sobolev space (see [1]). A vector field $\mathbf{v} \in V(\Omega^1) \times V(\Omega^2)$ is denoted
$\mathbf{v} = (\mathbf{v}^1, \mathbf{v}^2)$. Endowed with the standard inner product

$$(\mathbf{u}, \mathbf{v}) = (\mathbf{u}^1, \mathbf{v}^1)_{(H^1(\Omega^1))^2} + (\mathbf{u}^2, \mathbf{v}^2)_{(H^1(\Omega^2))^2},$$

$V(\Omega^1) \times V(\Omega^2)$ is a Hilbert space and the corresponding energy norm is denoted $\| \cdot \|$.

We define the bilinear form

$$a(\mathbf{u}, \mathbf{v}) = \sum_{\ell=1}^{2} \int_{\Omega^{\ell}} C^{\ell} \varepsilon(\mathbf{u}^{\ell}) \cdot \varepsilon(\mathbf{v}^{\ell}) \, d\Omega^{\ell}$$

for all $\mathbf{u}, \mathbf{v} \in V(\Omega^1) \times V(\Omega^2)$. Next, we denote $L(\cdot)$ the linear form which corresponds to the external loads:

$$L(\mathbf{v}) = \sum_{\ell=1}^{2} \left( \int_{\Omega^{\ell}} \mathbf{f}^{\ell} \cdot \mathbf{v}^{\ell} \, d\Omega^{\ell} + \int_{\Gamma_{\mathbf{g}}^{\ell}} \mathbf{g}^{\ell} \cdot \mathbf{v}^{\ell} \, d\Gamma^{\ell} \right).$$

The closed convex set $\mathbf{K}$ of admissible displacements is the subset of $V(\Omega^1) \times V(\Omega^2)$ which contains the
displacement fields satisfying the nonpenetration condition:

$$\mathbf{K} = \mathbf{v} = (\mathbf{v}^1, \mathbf{v}^2) \in V(\Omega^1) \times V(\Omega^2), \mathbf{v}^1 \cdot \mathbf{n}^1 + \mathbf{v}^2 \cdot \mathbf{n}^2 \le 0 \text{ on } \Gamma_c.$$

The variational inequality associated with the unilateral contact problem (2.1) (2.5) consists of finding $\mathbf{u}$ such that (see [13,16,21]):

$$\mathbf{u} \in \mathbf{K}, \quad a(\mathbf{u}, \mathbf{v} - \mathbf{u}) \ge L(\mathbf{v} - \mathbf{u}) \quad \forall \mathbf{v} \in \mathbf{K}. \qquad (2.6)$$

Using Stampacchia's theorem, we conclude that problem (2.6) has only one solution when (for example) $\Gamma_\mathbf{u}^\ell$, $\ell = 1, 2$, is of positive measure. Other conditions leading to existence and uniqueness results can be found in [16].

**Remark 2.1.** The more general and technical study involving an initial gap and the corresponding results are given in [16], Chapter 3, Sections 5 and 6.

### 3. The finite element approximation for nonmatching meshes

#### 3.1. The global and the local contact conditions

We suppose that $\Omega^1$ and $\Omega^2$ are domains with polygonal boundaries and we assume that $\Gamma_c$ is a straight line segment to simplify. Let the approximation parameter $h = (h_1, h_2)$ be a given pair of real positive numbers that will decay to 0. With each subdomain $\Omega^\ell$, we then associate a family of triangulations $T_h^\ell$, made of triangles denoted $\kappa$, the diameter of which does not exceed $h_\ell$. Therefore, we can write

$$\bar{\Omega}^{\ell} = \bigcup_{\kappa \in T_h^{\ell}} \bar{\kappa}.$$

The extreme points $c_1$ and $c_2$ of the contact part $\Gamma_c$ are supposed to belong to both meshes associated with $T_h^1$ and $T_h^2$. The contact zone $\Gamma_c$ inherits two independent families of discretizations arising from $T_h^1$ and $T_h^2$. The mesh $\mathcal{T}_{c,h}^\ell$ on $\Gamma_c$ is defined as the set of all the edges of $\kappa \in T_h^\ell$ on the contact zone. The set of the nodes associated with $\mathcal{T}_{c,h}^\ell$ is denoted $\xi_h^\ell$. In general $\xi_h^1$ and $\xi_h^2$ are not identical on account of the nonmatching meshes.
The space of the polynomials on $\kappa$ whose global degree is lower or equal to $q$, ($q$ nonnegative integer) is denoted $\mathbf{P}_q(\kappa)$. The finite element space used in $\Omega^\ell$ is then defined as (see [11]):

$$ \mathbf{V}_h(\Omega^\ell) = \{ \mathbf{v}_h^\ell \in (\mathscr{C}(\bar{\Omega}^\ell))^2 \ \forall \kappa \in T_h^\ell, \mathbf{v}_h^\ell|_{\kappa} \in (\mathbf{P}_1(\kappa))^2, \mathbf{v}_h^\ell|_{T_h^\ell} = \mathbf{U}_h^\ell \}, $$

where $\mathscr{C}(\bar{\Omega}^\ell)$ denotes the space of continuous functions on $\bar{\Omega}^\ell$ and $\mathbf{U}_h^\ell$ is a finite element approximation of $\mathbf{U}^\ell$.

In order to express the contact constraints (2.4) on $\Gamma_c$, we need to introduce some functional spaces. Let $W_h^\ell(\Gamma_c)$ be the range of $\mathbf{V}_h(\Omega^\ell)$ by the normal trace operator on $\Gamma_c$:

$$ W_h^\ell(\Gamma_c) = \{ \varphi_h = \mathbf{v}_h^\ell |_{\Gamma_c} \cdot \mathbf{n}^\ell, \mathbf{v}_h^\ell \in \mathbf{V}_h(\Omega^\ell) \}. \quad (3.1) $$

Next, we introduce the space of the Lagrange multipliers that will be useful to define a projection operator:

$$ M_h^\ell(\Gamma_c) = \{ \psi_h \in W_h^\ell(\Gamma_c), \psi_h|_{T \in \mathbf{P}_0(T)} \forall T \in \mathcal{T}_{c,h}^\ell, \text{ such that } \mathbf{c}_1 \text{ or } \mathbf{c}_2 \in T \}. $$

The notation $\pi_h^\ell$ stands for the projection operator on $W_h^\ell(\Gamma_c)$ defined for any function $\varphi \in \mathscr{C}(\bar{\Gamma}_c)$ as

$$ \pi_h^\ell \varphi \in W_h^\ell(\Gamma_c), \quad (\pi_h^\ell \varphi)(\mathbf{c}_i) = \varphi(\mathbf{c}_i) \quad \text{for } i = 1 \text{ and } 2, \quad \int_{\Gamma_c} (\varphi - \pi_h^\ell \varphi) \psi_h \, d\Gamma = 0 \quad \forall \psi_h \in M_h^\ell(\Gamma_c). \quad (3.2) $$

The condition $(\pi_h^\ell \varphi)(\mathbf{c}_i) = \varphi(\mathbf{c}_i)$, $i = 1, 2$, has been introduced in order to handle more general problems in a domain decomposition context (see [8]). It is easy to check that the classical $L^2$ projection operator on $W_h^\ell(\Gamma_c)$ does not satisfy such a condition. The approximation properties of $\pi_h^\ell$ are enumerated in Ben Belgacem, [4].

By using this projection operator (of global character), we are in a position to define the discrete admissibility convex cone $\mathbf{K}_h^{\text{glo}}$:

$$ \mathbf{K}_h^{\text{glo}} = \{ \mathbf{v}_h = (\mathbf{v}_h^1, \mathbf{v}_h^2) \in \mathbf{V}_h(\Omega^1) \times \mathbf{V}_h(\Omega^2), \mathbf{v}_h^1 \cdot \mathbf{n}^1 + \pi_h^1(\mathbf{v}_h^2 \cdot \mathbf{n}^2) \le 0 \text{ on } \Gamma_c \}. \quad (3.3) $$

Let us notice that the condition incorporated in $\mathbf{K}_h^{\text{glo}}$ is expressed in the space $W_h^1(\Gamma_c)$. Following the terminology of Bernardi et al. [8], $W_h^2(\Gamma_c)$ stands for the mortar space.

**Remark 3.1.** Of course, it is possible to give a symmetrical definition of the convex by choosing as mortar space $W_h^1(\Gamma_c)$ and using the projection operator $\pi_h^2$.

Let $\mathcal{I}_h^\ell$ denote the Lagrange interpolation operator ranging in $W_h^\ell(\Gamma_c)$. Then, we define the admissibility convex cone $\mathbf{K}_h^{\text{loc}}$ by using the interpolation operator of local character:

$$ \mathbf{K}_h^{\text{loc}} = \{ \mathbf{v}_h = (\mathbf{v}_h^1, \mathbf{v}_h^2) \in \mathbf{V}_h(\Omega^1) \times \mathbf{V}_h(\Omega^2), \mathbf{v}_h^1 \cdot \mathbf{n}^1 + \mathcal{I}_h^1(\mathbf{v}_h^2 \cdot \mathbf{n}^2) \le 0 \text{ on } \Gamma_c \}. \quad (3.4) $$

The discrete local contact conditions inserted in the definition of $\mathbf{K}_h^{\text{loc}}$ are similar to the classical node-on-segment conditions.

**Remark 3.2.** There are other possibilities of defining unilateral contact with nonmatching meshes (see [7,12,17]).

In addition, it is straightforward to check that $\mathbf{K}_h^{\text{glo}} \not\subset \mathbf{K}$ and $\mathbf{K}_h^{\text{loc}} \not\subset \mathbf{K}$. Therefore, both approximations are not “Hodge” conforming (see [11]). When matching meshes are used, the discrete unilateral constraints can be expressed in both cases merely by the natural node-on-node condition $\mathbf{v}_h^1 \cdot \mathbf{n}^1 + \mathbf{v}_h^2 \cdot \mathbf{n}^2 \le 0$ and the approximation becomes conforming ($\mathbf{K}_h^{\text{glo}} = \mathbf{K}_h^{\text{loc}} \subset \mathbf{K}$). This situation was extensively studied by Haslinger and Hlaváček [15] and Haslinger et al. [16].

The finite element problem issued from (2.6) is the following variational inequality: find $\mathbf{u}_h$ such that

$$ \mathbf{u}_h \in \mathbf{K}_h, \quad a(\mathbf{u}_h, \mathbf{v}_h - \mathbf{u}_h) \ge L(\mathbf{v}_h - \mathbf{u}_h), \quad \forall \mathbf{v}_h \in \mathbf{K}_h, \quad (3.5) $$

where $\mathbf{K}_h = \mathbf{K}_h^{\text{glo}}$ or $\mathbf{K}_h = \mathbf{K}_h^{\text{loc}}$.
Using again Stamppacchia's theorem, we conclude that problem (3.5) admits a unique solution under the assumptions mentioned in the previous section.

**Remark 3.3.** The finite element approximation in the general case of an initial gap between the bodies (with matching meshes) and the associated error estimations can be found in [16], chapter 3, Section 8.

## 3.2. Error estimation

We intend to give in the present part an estimate of the error committed on the exact solution by the global and the local finite element approximations.

In the next theorems, we adopt regularity assumptions which have been introduced by Brezzi et al. [9] for a Signorini problem, and used by Haslinger et al. [16] for the unilateral contact problem with matching meshes on the contact zone.

For technical reasons, we assume that the family of triangulations $T_h^\ell$ is regular (see [11]) and that $h_1/h_2$ is bounded. Moreover, we suppose that the measure of $\Gamma_u^\ell$ does not vanish and that $\mathbf{U}^\ell = 0$, $\ell = 1, 2$. We will make use of the standard Lebesgue and Sobolev spaces $L^\infty$, $W^{1,\infty}$, $(H^\ell)_{\tau \in \mathbb{R}_+}$; the detailed presentation of these spaces can be found in [1].

The approximation result associated with the global contact case is given in the following theorem.

**Theorem 3.1.** Suppose that the solution $\mathbf{u}$ of the continuous problem (2.6) is such that $\mathbf{u}^1 \in (H^2(\Omega^1))^2$, $\mathbf{u}^2 \in (H^2(\Omega^2))^2$, $\mathbf{u}^1 \cdot \mathbf{n}^1 \in W^{1,\infty}(\Gamma_c)$, $\mathbf{u}^2 \cdot \mathbf{n}^2 \in W^{1,\infty}(\Gamma_c)$ and $\sigma_n(\mathbf{u}) \in L^\infty(\Gamma_c)$. Suppose that the set of points of $\Gamma_c$ in which the change from $\mathbf{u}^1 \cdot \mathbf{n}^1 + \mathbf{u}^2 \cdot \mathbf{n}^2 < 0$ to $\mathbf{u}^1 \cdot \mathbf{n}^1 + \mathbf{u}^2 \cdot \mathbf{n}^2 = 0$ occurs is finite. Let $\mathbf{u}_h$ be the solution of the problem (3.5) with $\mathbf{K}_h = \mathbf{K}_h^{\text{glo}}$. Then

$$ \|\mathbf{u} - \mathbf{u}_h\| \le C(\mathbf{u})(h_1 + h_2), $$

where $C(\mathbf{u})$ is independent of $h$.

**Proof.** The starting point of the proof consists of a result in [6] which is the following:

$$ \|\mathbf{u} - \mathbf{u}_h\|^2 \le C^2(\mathbf{u})(h_1^2 + h_2^2) + C \int_{\Gamma_c} \sigma_n(\mathbf{u})(\mathcal{I}_h^1[\mathbf{u} \cdot \mathbf{n}] - [\mathbf{u} \cdot \mathbf{n}]) \, d\Gamma + C(\mathbf{u})h_1 (\|\mathbf{u} - \mathbf{u}_h\| + C(\mathbf{u})h_2), $$

where $[\mathbf{u} \cdot \mathbf{n}] = \mathbf{u}^1 \cdot \mathbf{n}^1 + \mathbf{u}^2 \cdot \mathbf{n}^2$. Then, by writing $2h_1h_2 \le h_1^2 + h_2^2$ and

$$ 2C(\mathbf{u})h_1 \|\mathbf{u} - \mathbf{u}_h\| \le \beta \|\mathbf{u} - \mathbf{u}_h\|^2 + \frac{1}{\beta} C^2(\mathbf{u}) h_1^2 $$

for any positive $\beta$, it comes out that

$$ \|\mathbf{u} - \mathbf{u}_h\|^2 \le C^2(\mathbf{u})(h_1^2 + h_2^2) + C \int_{\Gamma_c} \sigma_n(\mathbf{u})(\mathcal{I}_h^1[\mathbf{u} \cdot \mathbf{n}] - [\mathbf{u} \cdot \mathbf{n}]) \, d\Gamma \quad (3.6) $$

if $\beta$ is chosen small enough.

Using the condition $\sigma_n(\mathbf{u})(\mathbf{u}^1 \cdot \mathbf{n}^1 + \mathbf{u}^2 \cdot \mathbf{n}^2) = 0$ on $\Gamma_c$, and writing the integral term as a sum of integrals on the segments $t_h^1$ defined by the mesh of $\Omega^1$, we obtain

$$ \int_{\Gamma_c} \sigma_n(\mathbf{u})(\mathcal{I}_h^1[\mathbf{u} \cdot \mathbf{n}] - [\mathbf{u} \cdot \mathbf{n}]) \, d\Gamma = \int_{\Gamma_c} \sigma_n(\mathbf{u})(\mathcal{I}_h^1[\mathbf{u} \cdot \mathbf{n}]) \, d\Gamma = \sum_{t_h^1 \in \mathcal{F}_{c,h}} \int_{t_h^1} \sigma_n(\mathbf{u})(\mathcal{I}_h^1[\mathbf{u} \cdot \mathbf{n}]) \, d\Gamma. $$

It is obvious that the integrals which do not involve points of $\Gamma_c$ in which the change from $\mathbf{u}^1 \cdot \mathbf{n}^1 + \mathbf{u}^2 \cdot \mathbf{n}^2 < 0$ to $\mathbf{u}^1 \cdot \mathbf{n}^1 + \mathbf{u}^2 \cdot \mathbf{n}^2 = 0$ occurs are equal to zero. Therefore, it remains a finite number (independent of $h$) of integral terms which are bounded by using the regularity assumptions on the exact solution. This yields
$$
\begin{align*}
\int_{\Gamma_c} \sigma_n(\mathbf{u})(\mathcal{J}_h^1[\mathbf{u} \cdot \mathbf{n}] - [\mathbf{u} \cdot \mathbf{n}]) \, d\Gamma &= \sum_{\text{finite}} \int_{\Gamma_h^l} \sigma_n(\mathbf{u})(\mathcal{J}_h^1[\mathbf{u} \cdot \mathbf{n}]) \, d\Gamma \\
&\le Ch_1 \|\sigma_n(\mathbf{u})\|_{L^\infty(\Gamma_c)} \|\mathcal{J}_h^1[\mathbf{u} \cdot \mathbf{n}]\|_{L^\infty(\Gamma_c)} \\
&\le Ch_1^2 \|\sigma_n(\mathbf{u})\|_{L^\infty(\Gamma_c)} (\|\mathbf{u}^1 \cdot \mathbf{n}^1\|_{W^{1,\infty}(\Gamma_c)} + \|\mathbf{u}^2 \cdot \mathbf{n}^2\|_{W^{1,\infty}(\Gamma_c)}).
\end{align*}
$$

That concludes the proof. □

This theorem extends the result by Haslinger et al. established for matching meshes (see [16], Theorem 8.1).

**Remark 3.4.** The smoothness conditions $\mathbf{u}^l \cdot \mathbf{n}^l \in W^{1,\infty}(\Gamma_c)$, $l = 1, 2$, $\sigma_n(\mathbf{u}) \in L^\infty(\Gamma_c)$ as well as the condition on the finite number of points can be avoided. Indeed, under $H^2 \times H^2$ assumptions on the displacements, the convergence rate of the method is of the order $h_1^{3/4} + h_2$ (see [5]) as in the matching case (see [16]). The latter regularity assumptions can be again weakened, and under $H^v \times H^v$ ($3/2 < v \le 2$) assumptions, we obtain a convergence rate of order $h_1^{(v/2)-(1/4)} + h_2^{v-1}$ (see [6]).

**Remark 3.5.** It can be proved with a counterexample that the integral term of (3.6) cannot be bounded below $h_1^{3/2}$ under $H^2 \times H^2$ regularity assumptions (see [19]).

In the local contact case, we can only obtain the following convergence rate, which is suboptimal in the finite element sense.

**Theorem 3.2.** Let the assumptions of the previous theorem on $\mathbf{u}$ be fulfilled. Let $\mathbf{u}_h$ be the solution of problem (3.5) with $\mathbf{K}_h = \mathbf{K}_h^{\text{loc}}$. Then

$$ \|\mathbf{u} - \mathbf{u}_h\| \le C(\mathbf{u}) (\sqrt{h_1 + h_2}), $$

where $C(\mathbf{u})$ is independent of $h$.

**Proof.** By using an analogous estimate with that established in [6], we write

$$ \|\mathbf{u} - \mathbf{u}_h\|^2 \le C^2(\mathbf{u})(h_1^2 + h_2^2) + C \int_{\Gamma_c} \sigma_n(\mathbf{u})(\mathcal{J}_h^1[\mathbf{u} \cdot \mathbf{n}] - [\mathbf{u} \cdot \mathbf{n}]) \, d\Gamma + C(\mathbf{u}) \sqrt{h_1}(\|\mathbf{u} - \mathbf{u}_h\| + C(\mathbf{u})h_2). $$

Using the same arguments as in the previous theorem yields the result. □

**Remark 3.6.** The convergence rate of order $\sqrt{h_1 + h_2}$ comes from the poor approximation properties of the Lagrange interpolation operator in dual Sobolev spaces (see the counterexample in [19]) and therefore it has been proved that the estimates obtained in the analysis are optimal. Of course, one could dream that the analysis is inappropriate. The latter question seems to be open.

## 4. Matrix formulations

When solving the discretized unilateral contact problem, we use the finite element code CASTEM 2000 and a saddle-point formulation in which the multipliers are continuous functions on the contact zone and piecewise linear on the mesh of $\Omega^1$. A minimization type formulation for frictional contact problems can be found e.g. in [23]; for an augmented Lagrangian approach with nonmatching grids in linear elasticity, see [22]. A saddle-point formulation in which the multipliers are piecewise constant functions on the contact zone can be found in [16,25].
## 4.1. Preliminaries

At first, we intend to define the closed convex cone $\mathcal{M}_h$ of the discrete Lagrange multipliers. We set

$$ \mathcal{M}_h = \left\{ \lambda_h \in W_h^1(\Gamma_c), \int_{\Gamma_c} \lambda_h \varphi_h \, d\Gamma \le 0 \quad \forall \varphi_h \in W_h^1(\Gamma_c), \varphi_h \ge 0 \right\}, \quad (4.1) $$

where $W_h^1(\Gamma_c)$ is a space of continuous and piecewise linear functions defined in (3.1).

**Remark 4.1.** A function belonging to $\mathcal{M}_h$ is not necessarily nonpositive on $\Gamma_c$.

The following lemma is the tool used in order to bear out the choice of the Lagrange multipliers convex cone $\mathcal{M}_h$.

**Lemma 4.1.** Let $\varphi_h \in W_h^1(\Gamma_c)$. Then

$$ \varphi_h \le 0 \quad \text{if and only if} \quad \int_{\Gamma_c} \varphi_h \psi_h \, d\Gamma \ge 0 \quad \forall \psi_h \in \mathcal{M}_h. \quad (4.2) $$

**Proof.** Let us notice that the bilinear form $A$ defined on $W_h^1(\Gamma_c)$ by

$$ A(\theta_h, \rho_h) = \int_{\Gamma_c} \theta_h \rho_h \, d\Gamma $$

is an inner product on $W_h^1(\Gamma_c)$. Set

$$ \mathcal{N}_h = \{ \varphi_h \in W_h^1(\Gamma_c), \varphi_h \le 0 \}, $$

which is a closed convex cone. Then, we consider the polar cone of $\mathcal{N}_h$, denoted $\mathcal{N}_h^0$, and defined as follows (see [20]):

$$ \mathcal{N}_h^0 = \left\{ \psi_h \in W_h^1(\Gamma_c), \int_{\Gamma_c} \psi_h \varphi_h \, d\Gamma \le 0 \quad \forall \varphi_h \in \mathcal{N}_h \right\}. $$

The definition of $\mathcal{M}_h$ in (4.1) yields $\mathcal{N}_h^0 = -\mathcal{M}_h$. Using the property that the bipolar cone of a closed convex cone is the same convex cone, we deduce

$$ \begin{aligned} \mathcal{N}_h = (\mathcal{N}_h^0)^0 &= \left\{ \varphi_h \in W_h^1(\Gamma_c), \int_{\Gamma_c} \varphi_h \psi_h \, d\Gamma \le 0 \quad \forall \varphi_h \in \mathcal{N}_h^0 \right\}, \\ &= \left\{ \varphi_h \in W_h^1(\Gamma_c), \int_{\Gamma_c} \varphi_h \psi_h \, d\Gamma \ge 0 \quad \forall \varphi_h \in \mathcal{M}_h \right\}. \end{aligned} $$

Hence the lemma. □

## 4.2. The global contact case ($\mathbf{K}_h = \mathbf{K}_h^{\text{glo}}$)

Setting $\mathbf{V}_h = \mathbf{V}_h(\Omega^1) \times \mathbf{V}_h(\Omega^2)$, we consider the saddle-point problem on $\mathbf{V}_h \times \mathcal{M}_h$ associated with the following Lagrangian $\mathcal{L}^{\text{glo}}$:

$$ \mathcal{L}^{\text{glo}}(\mathbf{v}_h, \mu_h) = \frac{1}{2} a(\mathbf{v}_h, \mathbf{v}_h) - L(\mathbf{v}_h) - \int_{\Gamma_c} \mu_h (\mathbf{v}_h^1 \cdot \mathbf{n}^1 + \pi_h^1(\mathbf{v}_h^2 \cdot \mathbf{n}^2)) \, d\Gamma. \quad (4.3) $$
It is easy to verify that there exists a unique saddle-point $(\mathbf{u}_h, \lambda_h)$ on $\mathbf{V}_h \times \mathcal{M}_h$ and, using (4.2), it comes out that $\mathbf{u}_h$ is the unique solution of the variational problem (3.5) with $\mathbf{K}_h = \mathbf{K}_h^{\text{glo}}$. Moreover, it can be also proved that the multiplier $\lambda_h$ tends towards $\sigma_n(\mathbf{u})$ if $h = (h_1, h_2)$ decays to zero; this convergence result is beyond the scope of this paper and it will be established in a forthcoming study.

Denoting by $\mathbf{V}$ and $\mathbf{U}$ the vectors corresponding to the nodal values of $\mathbf{v}_h$ and $\mathbf{u}_h$, respectively, and by $\mathcal{M}$ and $\Lambda$ the vectors corresponding to the nodal values of $\mu_h$ and $\lambda_h$, respectively, the saddle-point problem (4.3) consists of finding $(\mathbf{U}, \Lambda)$ solution to

$$ \max_{A^1 M \le 0} \left( \min_{\mathbf{V}} \frac{1}{2}'\mathbf{V}\mathbf{K}\mathbf{V} - {}'\mathbf{V}\mathbf{F} - {}'(\mathbf{B}\mathbf{V})A^1 M \right), \quad (4.4) $$

where $K$ is the stiffness matrix, $F$ is a generalized load vector and $A^1$ is the mass matrix associated with the mesh of $\Omega^1$ on $\Gamma_c$. This means that $A^1$ is a $m$-by-$m$ matrix, where $m$ is the number of nodes of the mesh of $\Omega^1$ on $\Gamma_c$, and satisfying

$$ (A^1)_{i,j} = \int_{\Gamma_c} \psi_i \psi_j \, d\Gamma, \quad 1 \le i, j \le m, \quad (4.5) $$

where $\psi_i \in W_h^1(\Gamma_c)$ is equal to one on node number $i$ and to zero on the other nodes.

The matrix $B$ expresses the contact condition and requires the calculation of the projection operator $\pi_h^1$ mapping $W_h^2(\Gamma_c)$ into $W_h^1(\Gamma_c)$ (see (3.2)). In order to determine $B$, we need to give the matrix formulation of the condition $\mathbf{v}_h^1 \cdot \mathbf{n}^1 + \pi_h^1(\mathbf{v}_h^2 \cdot \mathbf{n}^2)$ incorporated in (3.3) and (4.3). Denoting by $n$ the number of nodes of $\Omega^2$ on $\Gamma_c$ and by $I_m$ the $m$-by-$m$ identity matrix, we have to find the $m$-by-$(m+n)$ matrix:

$$ (I_m | \Pi_h^1), \quad (4.6) $$

where $\Pi_h^1$ is the $m$-by-$n$ projection operator matrix. We denote by $\varphi_j \in W_h^2(\Gamma_c)$, $1 \le j \le n$, the function equal to one on node number $j$ and to zero on the other nodes. Then, we define $\theta_k \in W_h^1(\Gamma_c)$, $2 \le j \le m-1$, as follows:

$$ \theta_2 = \psi_1 + \psi_2, \quad \theta_k = \psi_k, \quad 3 \le k \le m-2, \quad \theta_{m-1} = \psi_{m-1} + \psi_m. $$

Using (3.2), we deduce $\Pi_h^1 = C^{-1}D$ where $C$ is the following $m$-by-$m$ matrix:

$$ \begin{aligned} C_{1,1} &= 1, \\ C_{1,j} &= 0, \quad 2 \le j \le m, \\ C_{i,j} &= \int_{\Gamma_c} \theta_i \psi_j \, d\Gamma, \quad 2 \le i \le m-1, \quad 1 \le j \le m, \\ C_{m,j} &= 0, \quad 1 \le j \le m-1, \\ C_{m,m} &= 1, \end{aligned} $$

and $D$ is the $m$-by-$n$ matrix verifying

$$ \begin{aligned} D_{1,1} &= 1, \\ D_{1,j} &= 0, \quad 2 \le j \le n, \\ D_{i,j} &= \int_{\Gamma_c} \theta_i \varphi_j \, d\Gamma, \quad 2 \le i \le m-1, \quad 1 \le j \le n, \\ D_{m,j} &= 0, \quad 1 \le j \le n-1, \\ D_{m,n} &= 1. \end{aligned} $$

In order to compute $D_{i,j}$, $2 \le i \le m-1$, $1 \le j \le n$, we denote by $\xi_h$ the set of the nodes located on $\Gamma_c$ such that $\xi_h = (\xi_h^1 \cup \xi_h^2) \setminus (\xi_h^1 \cap \xi_h^2)$. Denoting by $p$ the number of nodes in $\xi_h$ (one has $\max(m,n) \le p \le m+n-2$), we introduce the functions $(\chi_k)_{1 \le k \le p}$. The function $\chi_k$ is continuous on $\Gamma_c$, piecewise linear on the mesh defined by $\xi_h$, equal to one on node number $k$ and to zero on the other nodes of $\xi_h$. As a result, we write
$$D_{i,j} = \int_{\Gamma_c} \theta_i \varphi_j \, d\Gamma = \int_{\Gamma_c} \left( \sum_{k=1}^{p} (\alpha_i)_k \chi_k \right) \left( \sum_{k=1}^{p} (\beta_j)_k \chi_k \right) \, d\Gamma = \sum_{k=1}^{p} \sum_{k'=1}^{p} (\alpha_i)_k (\beta_j)_{k'} \int_{\Gamma_c} \chi_k \chi_{k'} \, d\Gamma.$$

The determination of the $(\alpha_i)_k$ and the $(\beta_j)_{k'}$ is done by taking the values of $\theta_i$ and $\varphi_j$ at the nodes of $\xi_h$. The $m$-by-$n$ projection matrix $C^{-1}D$ is then computed once for all. In the current bidimensional context, with the examples we consider (see the numerical studies), the computation of $C^{-1}D$ is not expensive. Nevertheless, if we want to adapt and extend to the three-dimensional case this global contact procedure, it will certainly be necessary to avoid the complete construction of $C^{-1}D$.

The solution $(\mathbf{U}, \mathcal{A})$ of (4.4) satisfies the relation $K\mathbf{U} - {}^t\mathcal{B}\mathbf{A}^1\mathcal{A} = \mathbf{F}$. So, setting $\Phi = \mathbf{A}^1\mathbf{M}$, the saddle-point problem (4.4) can be rewritten as a minimization problem of a quadratic functional with linear inequality constraints:

$$\min_{\Phi \le 0} \left( \frac{1}{2} {}^t \Phi B K^{-1} {}^t B \Phi + {}^t \Phi B K^{-1} F + \frac{1}{2} {}^t F K^{-1} F \right). \quad (4.7)$$

Since $m$ is the rank of $B$ and $K$ is symmetric and positive definite, it comes out that the matrix $BK^{-1}B$ is symmetric and positive definite. If $\Phi_0$ is the solution of the minimization problem (4.7), then $\mathcal{A} = (\mathbf{A}^1)^{-1}\Phi_0$ and the calculation of $\mathbf{U} = K^{-1}(\mathbf{F} + {}^t\mathcal{B}\Phi_0)$ is straightforward.

As already noticed in Remark 4.1, the components of the vector $\mathcal{A}$ (representing $\sigma_n(\mathbf{u})$) are not necessarily nonpositive. In a posteriori error estimation (see [12]), this “lack of nonpositiveness” must be added to the error.

The finite element code CASTEM 2000 solves the minimization problem (4.7) by using the iterative Frank and Wolfe algorithm (see [14]) which we recall hereafter.

Consider the problem of minimizing the functional $J: \mathbb{R}^n \to \mathbb{R}$ under linear constraints

$$\min_{\Phi \le 0} J(\Phi). \quad (4.8)$$

The method of Frank and Wolfe is iterative and generates a sequence of points $\Phi^0, \Phi^1, \dots, \Phi^k$, where $\forall k$, $\Phi^{k+1}$ is defined by using $\Phi^k$ as follows: solve the linear programming problem

$$(LP(\Phi^k)) \quad \min_{\Phi \le 0} {}^t(\nabla J(\Phi^k)) \cdot \Phi. \quad (4.9)$$

Let $y^k$ be an extremal point of $X = \{\Phi \in \mathbb{R}^n, C \le \Phi \le 0\}$ (where $|C|$ is chosen large enough such that $X$ contains a solution of (4.8)) and optimal solution to $(LP(\Phi^k))$. Then $\Phi^{k+1}$ is given by

$$J(\Phi^{k+1}) = \min_{\Phi \in [\Phi^k, \Phi^k]} J(\Phi).$$

If $J$ is continuously differentiable and if $J(\Phi) \to \infty$ as $\|\Phi\|_{\mathbb{R}^n} \to \infty$, then for every $\Phi^0 \in X$, the method converges towards a local minimum of $J(\Phi)$, $\Phi \in X$ (see [24]). Notice that for convex problems, the linear subproblem (4.9) provides a lower bound on the optimal objective value. The upper bound on the objective value is updated at each step and the algorithm is terminated when the relative difference between the bounds is smaller than a priori set parameter. The theoretical convergence rate of the algorithm is arithmetic (see [26,10]).

### 4.3. The local contact case ($\mathbf{K}_h = \mathbf{K}_h^{\text{loc}}$)

In this case, we consider the saddle-point problem on $\mathbf{V}_h \times \mathcal{M}_h$ associated with the Lagrangian $\mathcal{L}^{\text{loc}}$:

$$\mathcal{L}^{\text{loc}}(\mathbf{v}_h, \mu_h) = \frac{1}{2} a(\mathbf{v}_h, \mathbf{v}_h) - L(\mathbf{v}_h) - \int_{\Gamma_c} \mu_h (\mathbf{v}_h^1 \cdot \mathbf{n}^1 + \mathcal{I}_h^1 (\mathbf{v}_h^2 \cdot \mathbf{n}^2)) \, d\Gamma. \quad (4.10)$$

As in the global contact case, there exists a unique saddle-point $(\mathbf{u}_h, \lambda_h)$ on $\mathbf{V}_h \times \mathcal{M}_h$ and, using (4.2), we conclude that $\mathbf{u}_h$ is the unique solution of (3.5) with $\mathbf{K}_h = \mathbf{K}_h^{\text{loc}}$. The study of the convergence of $\lambda_h$ towards $\sigma_n(\mathbf{u})$ will be proposed in a following study.
The problem (4.10) consists then of solving the minimization problem

$$ \min_{\Phi \le 0} \left( \frac{1}{2} \Phi B' K^{-1} B' \Phi + \Phi B' K^{-1} F + \frac{1}{2} F K^{-1} F \right), \qquad (4.11) $$

where the matrix $B'$ involves now the local contact conditions. In fact, we need only to determine the matrix formulation of the Lagrange interpolation operator $\mathcal{I}_h^1$ mapping $W_h^2(\Gamma_c)$ into $W_h^1(\Gamma_c)$. This is done with the following $m$-by-$n$ matrix denoted $I_h^1$:

$$ (I_h^1)_{i,j} = (\mathcal{I}_h^1 \varphi_j)(a_i), \quad 1 \le i \le m, \ 1 \le j \le n, \qquad (4.12) $$

where $\varphi_j \in W_h^2(\Gamma_c)$ is equal to one on the node number $j$ and to zero on the others, and $a_i$ denotes node number $i$ on the mesh of $\Omega^1$ on $\Gamma_c$.

Problem (4.11) is solved in the finite element code like problem (4.7).

**Remark 4.2.** In the case where the bodies are supposed to come into contact after deformation, we must take into account of the initial gap. So, we replace in (4.7) and (4.11), $BV$ and $B'V$ by $BV - G$ and $B'V - G$, respectively, where $G$ is the $m$-vector whose components are the distances between the nodes of $\Omega^1$ on the candidate contact zone and the boundary of $\Omega^2$.

**Remark 4.3.** In order to show concretely the differences between the global contact conditions and the local node-on-segment conditions, we propose to illustrate their characteristics with the simple example depicted in Fig. 1. There are 7 equidistant nodes of $\Omega^1$ and 5 equidistant nodes of $\Omega^2$ on the contact zone. It is easy to see that the interpolation matrix $I_h^1$ defined in (4.12) is given by

$$ I_h^1 = \begin{pmatrix} 1.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 \\ 0.3333 & 0.6666 & 0.0000 & 0.0000 & 0.0000 \\ 0.0000 & 0.6666 & 0.3333 & 0.0000 & 0.0000 \\ 0.0000 & 0.0000 & 1.0000 & 0.0000 & 0.0000 \\ 0.0000 & 0.0000 & 0.3333 & 0.6666 & 0.0000 \\ 0.0000 & 0.0000 & 0.0000 & 0.6666 & 0.3333 \\ 0.0000 & 0.0000 & 0.0000 & 0.0000 & 1.0000 \end{pmatrix}, $$

whereas the projection matrix $\Pi_h^1$ introduced in (4.6) is as follows:

$$ \Pi_h^1 = \begin{pmatrix} 1.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 \\ 0.2947 & 0.7440 & -0.0379 & -0.0016 & 0.0008 \\ -0.0566 & 0.7799 & 0.2727 & 0.0080 & -0.0040 \\ 0.0152 & -0.0303 & 1.0303 & -0.0303 & 0.0152 \\ -0.0040 & 0.0080 & 0.2727 & 0.7799 & -0.0566 \\ 0.0008 & -0.0016 & -0.0379 & 0.7740 & 0.2947 \\ 0.0000 & 0.0000 & 0.0000 & 0.0000 & 1.0000 \end{pmatrix}. $$

Of course, the terms of the matrices $I_h^1$ and $\Pi_h^1$ are rounded numbers. The local character of the interpolation operator is given by the numerous terms (outside of the “diagonal”) of the matrix $I_h^1$ which are

<figure><img src="image_2.png" /><figcaption>Fig. 1. A simple example of nonmatching meshes.</figcaption></figure>
equal to zero. That merely means that the nodes which are distant do not interact in the definition of the node-on-segment conditions. On the contrary, the global character of the projection operator is shown by the nonzero terms of the matrix $\Pi_h^1$. The terms which are near to 1 represent strong interaction between close nodes of both bodies whereas the terms near to zero correspond to little interaction between distant nodes.

# 5. Numerical studies

In this section, we report numerical studies on several problems dealing with global and local contact conditions. The numerical experiments have been made at the Laboratoire de Mathématiques pour l'Industrie et la Physique of the Université de Toulouse, at the Laboratoire de Mécanique et Technologie of the Ecole Normale Supérieure de Cachan and more recently at the Laboratoire de Mathématiques of the Université de Savoie. In all three cases, we used the finite element code CASTEM 2000.

In accordance with the previous notations of Sections 3 and 4, we adopt this convention: the upper body always stands for $\Omega^1$.

## 5.1. Test 1: comparison of convergence rates for global, local and node-on-node conditions

In a previous section, we have considered the theoretical convergence rates of the discretized solutions towards the solution **u** of the continuous problem (2.6). This test consists of comparing these convergence rates in a numerical context.

Contact problems generally do not admit an analytical solution. Therefore, we must have a solution with finely discretized bodies at our disposal, which is a reference solution (denoted $\mathbf{u}_{h,\text{ref}}$) for error estimates. In order to obtain the convergence curve of the error, we build a family of nested meshes. This family is obtained by an algorithm which divides the triangles: we begin with a very coarse mesh and the following mesh is obtained by the natural subdivision of each triangle in four triangles. We then compute the finite element solution $\mathbf{u}_h$ on each mesh. As previously noticed, the error $\|\mathbf{u} - \mathbf{u}_h\|$ is approximated by $\|\mathbf{u}_{h,\text{ref}} - \mathbf{u}_h\|$. The latter expression is estimated by $\|\sigma_h((\mathcal{I}_h \mathbf{u}_{h,\text{ref}}) - \mathbf{u}_h)\|_*$ where $\mathcal{I}_h$ denotes the Lagrange interpolation operator, $\sigma_h((\mathcal{I}_h \mathbf{u}_{h,\text{ref}}) - \mathbf{u}_h)$ is the stress tensor field associated with $(\mathcal{I}_h \mathbf{u}_{h,\text{ref}}) - \mathbf{u}_h$ and $\|\cdot\|_*$ stands for the standard $L^2(\Omega^1 \cup \Omega^2)$-norm defined on the space of tensor fields. The most refined mesh is the reference mesh. For obvious reasons, the error can not be estimated by taking $\mathbf{u}_h = \mathbf{u}_{h,\text{ref}}$ and we choose the most refined mesh for error computations such that $h = 4h_{\text{ref}}$. We are interested in the rate of convergence denoted $\alpha$ such that $\|\mathbf{u}_{h,\text{ref}} - \mathbf{u}_h\|/\|\mathbf{u}_{h,\text{ref}}\| = Ch_1^\alpha$, where the notation $h_1$ represents the discretization parameter of the upper body.

* Comparison between matching and nonmatching meshes with global contact: We consider the contact problem of Fig. 2. In order to avoid singularities of Dirichlet Neumann type, we adopt symmetry conditions.

The length of the edges of the bodies in Fig. 2 is 1 mm and plane strain conditions are assumed. We choose a Poisson's ratio of 0.2 for both solids and Young's modulus $E_1 = 13\,000$ Mpa and $E_2 = 30\,000$ Mpa for the upper and the lower body, respectively. The applied loads on the two parts of the boundary of the upper body are 100 daN/mm².

In the matching case with node-on-node contact conditions, the reference solution is obtained with meshes corresponding to 66 564 d.f. and 65 536 triangular elements. The contact zone comprises 128 matching meshes and we use the node-on-node contact condition. The relative normal displacement on the contact zone for the reference problem is represented in Fig. 3. The error is estimated by using 5 nested meshes and the most refined comprises 32 matching meshes on the contact zone.

In the nonmatching case with global contact conditions, the reference solution is obtained with meshes corresponding to 28 804 d.f. and 28 160 triangular elements. The reference mesh contains 65 nodes of the upper body and 97 nodes of the lower body on the contact zone. The curve of the error is obtained with 4 nested meshes, the most refined having on the contact zone 17 nodes of the upper mesh and 25 of the lower mesh.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Reference problem.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. The relative normal displacement on the contact zone for the reference problem.</figcaption></figure>

Fig. 4 depicts the convergence rate of the relative error (in the energy norm) as a function of the discretization parameter $h_1$ of the upper body. The mean value of the convergence rate is $\alpha = 1.21$ in the matching case and $\alpha = 1.26$ in the nonmatching case. Let us notice that the study of the convergence in the $L^2$-norm of the error instead of the $H^1$-norm, yields the following mean values of the convergence rates: $\alpha = 1.75$ in the matching case and $\alpha = 1.69$ in the nonmatching case. The two dotted curves of Fig. 5 show the relative error in the energy norm as a function of the number of d.f. On this example, the convergence rates are not weakened when nonmatching meshes and global contact conditions are used, as already proved in a theoretical context.

* *Comparison between global and local contact conditions:* Now, we intend to compare the numerical convergence rates corresponding to global and local contact conditions.

Let us consider the problem of Fig. 2. Henceforward, we lay emphasis on the case of nonmatching meshes. We choose a nonmatching reference mesh corresponding to 64388 d.f. and 63488 triangular elements. On the contact zone, there are 65 nodes of the upper body and 161 nodes of the lower body.

For both contact conditions, the computation of the error is done with 4 nested meshes and the most refined mesh comprises 17 nodes of the upper solid and 41 nodes of the lower solid on the
<figure><img src="image_5.png" /><figcaption>Fig. 4. The convergence curves for matching meshes with node on node contact condition (o) and for nonmatching meshes with global contact condition (x).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. The convergence curves for node on node (+), global (o) and local (x) contact conditions (the two dotted lines correspond to the curves of Fig. 4).</figcaption></figure>

contact zone. Notice that the calculation using the two contact conditions is achieved with the same meshes.

The mean value of the convergence rate of the error (in the energy norm) as a function of the discretization parameter $h_1$ of the upper body is $\alpha = 1.26$ in the global case and $\alpha = 0.98$ in the local case. The convergence rates in the $L^2$-norm are as follows: $\alpha = 1.69$ with global contact and $\alpha = 1.53$ with local contact.
The two full curves in Fig. 5 represent the convergence rates of the error (in the energy norm) as a function of the number of d.f. associated with the nested meshes.

Owing to the test, it seems that the behaviour of the error is better in the global case than in the local case.

## 5.2. Test 2: a qualitative comparison between global and local contacts

The geometry of the problem and the finite element meshes are shown in Fig. 6. The geometrical and material characteristics of both bodies are the same as in the first test and we apply a uniform load of 100 daN/mm² on the top of the upper body.

We consider the discretization with nonmatching meshes depicted in Fig. 6 and we compare the two different contact conditions (global and local) with this configuration.

The difference between the solutions of global and local type is considerable. Using the same number of inequalities (13) corresponding to the number of nodes of the upper body on the contact zone, the global contact approach yields a very satisfactory solution in Fig. 7 (with an undetectable interpenetration and a straight contact zone), whereas the local contact approach shows a quite unacceptable solution (with an important interpenetration of the bodies) in Fig. 8.

Then, we consider the stress fields $\sigma^{\text{glo}}$ and $\sigma^{\text{loc}}$ (obtained from the displacement fields with the constitutive law). We intend to compare the component $\sigma_{yy}^{\text{glo}}$ of $\sigma^{\text{glo}}$ with the component $\sigma_{yy}^{\text{loc}}$ of $\sigma^{\text{loc}}$, where y denotes the vertical. The exact solution for this problem is a uniform $\sigma_{yy}$ field of value -100.

In the case of global contact, the obtained $\sigma_{yy}^{\text{glo}}$ field is quasi-uniform (minimum value = -100.00008, maximum value = -99.99990) as shown in Fig. 9. On the contrary, the local contact approach yields a field $\sigma_{yy}^{\text{loc}}$ which is not at all uniform (minimum value = -233, maximum value = -5), particularly near the contact zone (see Fig. 10).

Finally, Fig. 11 shows the Lagrange multipliers equal to $(A^1)^{-1}\Phi_0$. The mass matrix $A^1$ has been introduced in (4.5) and $\Phi_0$ is the solution of the minimization problem (4.7) in the global case and of (4.11) in the local case. These multipliers, defined on the contact zone, express the normal stresses (exact value =

<figure><img src="image_7.png" /><figcaption>Fig. 6. The reference problem and the meshes.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 7. Deformed configuration with global contact (amplification: 100).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Deformed configuration with local contact (amplification: 100).</figcaption></figure>

−100). Once again, we notice the quite good value of the global contact multiplier and the important irregularities of the multiplier obtained when using local contact conditions.

If we choose the symmetrical definition of the contact condition (3.3) (see Remark 3.1) and the symmetrical definition of the contact condition (3.4), we obtain the deformed configurations of Figs. 12 and 13. In this case, the difference between the global and the local contact conditions is obviously less significant than in the previous comparison but the global solution remains a bit better.

As in the symmetrical case, we consider the stress fields $\sigma^{glo}$ and $\sigma^{loc}$. In the case of global contact, the obtained $\sigma_{yy}^{glo}$ field is still quasi-uniform (minimum value = −100.0032, maximum value = −99.9965) as shown in Fig. 14. Notice that the local approach gives a suitable field $\sigma_{yy}^{loc}$ (minimum value = −101.86, maximum value = −95.935, see Fig. 15). Concerning the Lagrange multipliers depicted in Fig. 16, we can still notice the very good results given by the global approach and the satisfactory solution yielded by the local conditions.
<figure><img src="image_10.png" /></figure>

GIBI FECIT

<caption>Fig. 9. The $\sigma_{yy}^{glo}$ field obtained with global contact (minimum value 100.00008, maximum value 99.99990).</caption>

VAL = 120

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

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

GIBI FECIT

<caption>Fig. 10. The $\sigma_{yy}^{loc}$ field obtained with local contact (minimum value 233.77, maximum value 5.42).</caption>

VAL = 120

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

From this test, it becomes manifest that the local contact approach must be avoided, especially when defining the constraints on the coarser grid. Let us remark that the meshes have been chosen precisely to show a great difference between the two results. Talking of that, we notice that 13 inequalities of local type (or of node-on-segment type) describe in a very poor way the contact, whereas the 13 inequalities of global type are quite representative. This also explains the spectacular superiority of the global contact approach.

When defining the contact conditions on the finer grid, the difference between the two approaches is less significant but the global technique still leads to better results.
<figure><img src="image_14.png" /><figcaption>Fig. 11. Lagrange multipliers corresponding to global (o) and local (x) conditions.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 12. Deformed configuration with the symmetrical definition of global contact (amplification: 100).</figcaption></figure>

5.3. Test 3: *A case with an initial gap*

Fig. 17 shows the contact problem between an elastic half-disc and an elastic support. The aim of this example is to try to adapt the global contact procedure to a more general context than the previous ones.

In such a configuration, we have to consider nonmatching meshes, on account of the geometries of the bodies. Moreover, there is an initial gap and consequently, there are points of the boundaries initially not in contact which will come into contact after deformation. So, we define the contact by introducing an extended global condition which takes into account of the initial gap (see Remark 4.2). We choose end points
<figure><img src="image_16.png" /><figcaption>Fig. 13. Deformed configuration with the symmetrical definition of local contact (amplification: 100).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 14. The $\sigma_{yy}^{glo}$ field obtained with the symmetrical definition of global contact (minimum value 100.0032, maximum value 99.9965).</figcaption></figure>

$c_1 = (x_{c1}, y_{c1})$, $c_2 = (x_{c2}, y_{c2})$ on $\partial\Omega^1$ and $d_1 = (x_{d1}, y_{d1})$, $d_2 = (x_{d2}, y_{d2})$ on $\partial\Omega^2$ satisfying $x_{c1} = x_{d1}$ and $x_{c2} = x_{d2}$. The latter construction is achieved in order to have a common interface $\Gamma_c = [d_1, d_2]$, where we can project (in the vertical direction) the nodes of the arc of the circle $(c_1, c_2)$ and define the global contact conditions.

The half-disc is 20 mm in diameter and the length of an edge of the support is 40 mm. A Poisson's ratio of 0.2 for both solids, Young's modulus $E_1 = 25\,000$ Mpa for the upper body and $E_2 = 15\,000$ Mpa for the lower body are assumed. The applied loads on the top are 500 daN/mm$^2$.
<figure><img src="image_18.png" /><figcaption>Fig. 15. The $\sigma_{yy}^{loc}$ field obtained with the symmetrical definition of local contact (minimum value 101.86, maximum value 95.935).</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 16. Lagrange multipliers corresponding to symmetrical global (o) and symmetrical local (x) conditions.</figcaption></figure>

The initial and the deformed configurations are depicted in Fig. 18. Fig. 19 represents the initial and the deformed meshes near the contact zone, and we observe a deformed configuration which seems quite satisfactory, particularly on the contact part.
<figure><img src="image_20.png" /><figcaption>Fig. 17. Reference problem.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 18. The initial and the deformed configurations.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Fig. 19. The initial and the deformed meshes near the contact zone.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 20. The contact pressure on the arc of the circle ($c_1, c_2$) on $\partial\Omega^1$.</figcaption></figure>

The contact pressure (given by the Lagrange multiplier) on the arc of the circle ($c_1, c_2$) on $\partial\Omega^1$ is depicted in Fig. 20. By using the generalized load vectors on the nodes of $[d_1, d_2]$ on $\partial\Omega^2$, it becomes possible to obtain the contact pressure on $[d_1, d_2]$ as shown in Fig. 21.
<figure><img src="image_24.png" /><figcaption>Fig. 21. The contact pressure on the segment $[d_1, d_2]$ on $\partial\Omega^2$.</figcaption></figure>

### 5.4. Test 4: Taking into account quasi-matching meshes and strong variations of the contact pressure

The purpose of this last example is to show how the global contact conditions can take into account quasi-matching meshes and strong variations of the contact pressure.

Let us consider the contact problem of Fig. 22. The dimensions of $\Omega^1$ and $\Omega^2$ are 1 × 0.05 mm. A Poisson's ratio of 0.2 for both solids, Young's modulus $E_1 = 25\,000$ Mpa for the upper body and $E_2 = 15\,000$ Mpa for the lower body are assumed. The applied loads are 100 daN/mm². The mesh of $\Omega^1$ divides $\Gamma_c$ into 120 identical segments and the mesh of $\Omega^2$ divides $\Gamma_c$ into 121 identical segments as suggested in Fig. 23. The deformed meshes are shown in Fig. 24 and the generalized loads at the nodes of $\Omega^1$ on $\Gamma_c$ are depicted in Fig. 25. Finally, the multiplier, representing the contact pressure is obtained by using the latter loads (see Fig. 26). As already noticed in Remark 4.1, the multiplier is not always nonpositive.

As a result, this example shows that the global contact procedure takes into account strong variations of the contact pressure (and meshes which seem difficult to handle) in a satisfactory way.

<figure><img src="image_25.png" /><figcaption>Fig. 22. Setting of the problem.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>Fig. 23. The meshes (left part of the structure).</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>Fig. 24. The deformed configuration.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Fig. 25. The generalized loads at the nodes of Ω¹ on Γ_c.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Fig. 26. The multiplier representing the contact pressure.</figcaption></figure>
## 6. Conclusion and perspectives

In order to solve the unilateral contact problem between two elastic bodies, we have considered two finite element approximations (one global, one local) of order one using nonmatching meshes on the contact zone.

We have proved that the global approximation extends in an optimal way the results of Haslinger et al. [16], established in the case of matching meshes and we have obtained a convergence result in the case of the local approach. The methods have been compared and we come to the conclusion that the global approach could be a promising technique.

For frictional contact, the global contact condition can also be used, and the first theoretical studies can be found in [3,18]. In other respects, it would be interesting to use other minimization algorithms, in particular Newton like minimization techniques (see [2]). Finally, the extension of such a global contact technique to three dimensional problems by using the results of [4] is the following study which should be investigated.

## Acknowledgements

The author is deeply grateful to P. Laborde and F. Ben Belgacem from MIP-Toulouse for their useful advice and contribution. Many thanks to J.-P. Pelle and P. Coorevits from LMT-Cachan for their help in the first numerical experiments. I am also very thankful to the anonymous referee for the numerous remarks and suggestions which have led to significant improvements of this paper.

## References

[1] R.A. Adams, Sobolev Spaces, Academic Press, New York, 1975.

[2] P. Alart, Méthode de Newton généralisée en mécanique du contact, J. Math. Pures Appl. 76 (1997) 83-108.

[3] G. Bayada, M. Chambat, K. Lhalouani, T. Sassi, Éléments finis avec joints pour des problèmes de contact avec frottement de Coulomb non local, C. R. Acad. Sci. Paris 325 (Série I (1997) 1323-1328.

[4] F. Ben Belgacem, Discrétisations 3 D non conformes par la méthode de décomposition de domaine des éléments avec joints: Analyse mathématique et mise en œuvre pour le problème de Poisson, Thèse de l'Université Pierre et Marie Curie, Paris 6, 1993.

[5] F. Ben Belgacem, P. Hild, P. Laborde, Approximation of the unilateral contact problem by the mortar finite element method, C. R. Acad. Sci. Paris 324 (Série I (1997) 123-127.

[6] F. Ben Belgacem, P. Hild, P. Laborde, Extension of the mortar finite element method to a variational inequality modeling unilateral contact, Internal report of MIP, IR96.18, Math. Mod. and Meth. in the Appl. Sci. 9 (2), (1999) 287-303.

[7] F. Ben Belgacem, P. Hild, P. Laborde, The mortar finite element method for contact problems, Math. Comput. Model. 28 (1998) 263-271.

[8] C. Bernardi, Y. Maday, A.T. Patera, A new nonconforming approach to domain decomposition: The mortar element method, in: H. Brezis, J. L. Lions, Pitman (Eds.), Collège de France Seminar, 1994, pp. 13-51.

[9] F. Brezzi, W.W. Hager, P.A. Raviart, Error estimates for the finite element solution of variational inequalities, Numer. Math. 28 (1997) 431-443.

[10] M.D. Canon, C.D. Cullum, A tight upper bound on the rate of convergence of the Frank Wolfe algorithm, SIAM J. on Control 6 (1968) 509-516.

[11] P. G. Ciarlet, The Finite Element Method for Elliptic Problems, North Holland, Amsterdam, 1978.

[12] P. Coorevits, P. Hild, J. P. Pelle, Contrôle des calculs par éléments finis pour un problème de contact unilatéral, Internal report of LMT Cachan, IR203, submitted to Revue Européenne des Éléments Finis, 1998.

[13] G. Duvaut, J. L. Lions, Les inéquations en mécanique et en physique, Dunod, Paris, 1972.

[14] M. Frank, P. Wolfe, An algorithm for quadratic programming, Naval Research Logist Quarterly 3 (1956) 95-110.

[15] J. Haslinger, I. Hlaváček, Contact between elastic bodies, 2. Finite Element Analysis, Aplikace Matematiky 26 (1981) 263-290.

[16] J. Haslinger, I. Hlaváček, J. Nečas, Numerical methods for unilateral problems in solid mechanics, in: P.G. Ciarlet, J. L. Lions (Eds.), Handbook of Numerical Analysis, vol. IV, Part 2, North Holland, Amsterdam, 1996, pp. 313-485.

[17] P. Hild, Problèmes de contact unilatéral et maillages éléments finis incompatibles, Thèse de l'Université Paul Sabatier, Toulouse 3, 1998.

[18] P. Hild, Éléments finis non conformes pour un problème de contact unilatéral avec frottement, C. R. Acad. Sci. Paris 324 (Série I (1997) 707-710.

[19] P. Hild, A propos d'approximation par éléments finis optimale pour les problèmes de contact unilatéral, C. R. Acad. Sci. Paris 326 (Série I (1998) 1233-1236.
[20] J. B. Hiriart Urruty, C. Lemaréchal, Convex Analysis and Minimisation Algorithms, vol. I and II, in: Grundlehren der Mathematischen Wissenschaften (305,306), Springer, Berlin, 1993.

[21] N. Kikuchi, J.T. Oden, Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods, SIAM, Philadelphia, 1988.

[22] P. Le Tallec, T. Sassi, Domain decomposition with nonmatching grids: Augmented Lagrangian approach, Math. of Comp. 64 (1995) 1367–1396.

[23] C. Licht, E. Pratt, M. Raous, Remarks on a numerical method for unilateral contact including friction, in: International Series of Numerical Mathematics, vol. 101, Birkhauser, 1991, pp. 129–144.

[24] M. Minoux, Programmation mathématique, théorie et algorithmes, Tome 1, Dunod, 1983.

[25] K. Lhalouani, T. Sassi, Nonconforming mixed variational formulation and domain decomposition for unilateral problems, Internal report of Equipe d'Analyse Numérique Lyon/Saint Etienne, IR286, 1998.

[26] P. Wolfe, Convergence theory in nonlinear programming, in: J. Abadie (Ed.), Integer and Nonlinear Programming, North Holland, Amsterdam, 1970, pp. 1–36.