# Variationally consistent discretization schemes and numerical algorithms for contact problems

Barbara Wohlmuth

► To cite this version:

Barbara Wohlmuth. Variationally consistent discretization schemes and numerical algorithms for contact problems. Acta Numerica, 2011, 20, pp.569 - 734. 10.1017/S0962492911000079. hal-01382364

HAL Id: hal-01382364

https://hal.science/hal-01382364

Submitted on 17 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
# Variationally consistent discretization schemes and numerical algorithms for contact problems

Barbara Wohlmuth

Technische Universität München
Fakultät für Mathematik M2,
Boltzmannstr. 3, 85748 Garching, Germany
E-mail: wohlmuth@ma.tum.de
URL: www-m2.ma.tum.de

We consider variationally consistent discretization schemes for mechanical contact problems. Most of the results can also be applied to other variational inequalities, such as those for phase transition problems in porous media, for plasticity or for option pricing applications from finance. The starting point is to weakly incorporate the constraint into the setting and to reformulate the inequality in the displacement in terms of a saddle-point problem. Here, the Lagrange multiplier represents the surface forces, and the constraints are re-stricted to the boundary of the simulation domain. Having a uniform inf-sup bound, one can then establish optimal low-order *a priori* convergence rates for the discretization error in the primal and dual variables. In addition to the abstract framework of linear saddle-point theory, complementarity terms have to be taken into account. The resulting inequality system is solved by rewrit-ing it equivalently by means of the non-linear complementarity function as a system of equations. Although it is not differentiable in the classical sense, semi-smooth Newton methods, yielding super-linear convergence rates, can be applied and easily implemented in terms of a primal-dual active set strat-egy. Quite often the solution of contact problems has a low regularity, and the efficiency of the approach can be improved by using adaptive refinement techniques. Different standard types, such as residual- and equilibrated-based *a posteriori* error estimators, can be designed based on the interpretation of the dual variable as Neumann boundary condition. For the fully dynamic setting it is of interest to apply energy-preserving time-integration schemes. However, the differential algebraic character of the system can result in high oscillations if standard methods are applied. A possible remedy is to modify the fully discretized system by a local redistribution of the mass. Numeri-cal results in two and three dimensions illustrate the wide range of possible applications and show the performance of the space discretization scheme, non-linear solver, adaptive refinement process and time integration.
# 1. Introduction

In many industrial applications or engineering problems, contact between deformable elastic bodies plays a crucial role. As examples we mention incremental forming processes, the simulation of rolling wheels, braking pads on tyres and roller bearings. Although early theoretical results go back to Hertz (1882), there are still many open problems, and the numerical simulation of dynamic contact problems remains challenging.

These problems are discussed in several monographs on contact mechanics such as Fischer-Cripps (2000), Johnson (1985) and Kikuchi and Oden (1988). More recent theoretical results on existence and uniqueness can be found in Eck, Jarušek and Krbec (2005) and Han and Sofonea (2000), and on mathematical models and numerical simulation techniques in Laursen (2002), Willner (2003), Wriggers (2006), Wriggers and Nackenhorst (2007) and the references therein. One of the main challenges relates to the fact that the actual contact zone is not known *a priori* and has to be identified by use of an iterative solver. Moreover, the transition between contact and non-contact is characterized by a change in the type of the boundary condition, and thus possibly results in a solution of reduced regularity. From the mathematical point of view, contact problems can be formulated as free boundary value problems and analysed within the abstract framework of variational inequalities (Facchinei and Pang 2003a, 2003b, Glowinski, Lions and Trémolières 1981, Harker and Pang 1990, Kinderlehrer and Stampacchia 2000).

This work is an overview of theoretical and numerical results obtained in recent years. Most of the numerical examples are therefore taken from the original papers which are cited in the reference list. The numerical implementation is based on different software codes. In particular, DUMUX (Flemisch, Fritz, Helmig, Niessner and Wohlmuth 2007), DUNE (Bastian, Blatt, Dedner, Engwer, Klöfkorn, Kornhuber, Ohlberger and Sander 2008), NETGEN (Schöberl 1997), PARAVIEW (Ahrens, Geveci and Law 2005),
PARDISO (Karypis and Kumar 1998, Schenk and Gärtner 2004, 2006), and UG (Bastian, Birken, Johannsen, Lang, Neuß, Rentz-Reichert and Wieners 1997) have been used.

The structure of this paper is as follows. In Section 2, the governing equations and corresponding inequality constraints for frictional contact problems are stated. Different equivalent formulations are discussed and a weak saddle-point formulation is presented. Section 3 is devoted to space discretization. Special emphasis is placed on uniformly inf-sup stable pairings and a suitable approximation of the dual cone. Optimal *a priori* estimates for the discretization error in the displacement and in the surface traction are given for low-order finite elements in Section 4. Here, we restrict ourselves to very simple contact settings with given friction bounds and do not take non-matching contact zones into account, but allow for non-matching meshes.

In Section 5, we survey different solver techniques for the non-linear inequality system. Of special interest are so-called semi-smooth Newton methods applied to an equivalent non-linear system of equations. We discuss in detail the structure of the systems to be solved after consistent linearization. In particular, in the case of no friction the Newton solver can easily be implemented as a standard primal-dual active set strategy, updating in each iteration step the type and the value of the boundary condition node-wise.

Section 6 is devoted to different aspects of adaptive refinement. Bearing in mind that the mechanical role of the discrete Lagrange multiplier is that of a surface traction, different error indicators can easily be designed. However, the analysis is quite challenging and only a few theoretical results exist, taking into account possibly non-matching meshes and the inequality character of the formulation. Here, we provide upper and lower bounds for a simplified setting and comment on possible generalizations.

Section 7 is devoted to aspects of time integration. For many applications structure-preserving time-integration schemes are of special interest. In this context, energy preservation is of crucial importance. Unfortunately most of the standard techniques result either in very high oscillations in the dual variable or in numerical dissipation. We apply a newly combined time and space integration scheme which is motivated by a reduction of the index of the differential algebraic system.

Finally, in Section 8 we illustrate the flexibility of the proposed approach by considering applications from different areas. In particular, an example from finance shows that the Lagrange multiplier approach based on a $(d-1)$-dimensional $H^{1/2}$-duality pairing can also be applied to obstacle-type inequalities reflecting a $d$-dimensional $H^1$-duality pairing. Of special interest are examples where $d$- and $(d-1)$-dimensional constraints are imposed, such as phase transition problems in heterogeneous porous media and elasto-plastic mechanical contact problems.
## 2. Problem setting for mechanical contact

In many applications involving several deformable bodies, frictional contact has to be considered in conjunction with inelastic material behaviour such as plasticity. A great deal of research has been done on both of these topics: see for example Boieri, Gastaldi and Kinderlehrer (1987), Eck *et al.* (2005), Johnson (1985), Laursen (2002), Willner (2003), Wriggers (2006) and the references therein for an overview of contact problems. Characteristically, this type of application leads to a constrained minimization problem or more generally to a variational inequality (Harker and Pang 1990, Haslinger, Hlaváček, Nečas and Lovíšek 1988, Kikuchi and Oden 1988, Kinderlehrer and Stampacchia 2000). Mathematical analyses of variational inequalities and constrained minimization problems can also be found in Facchinei and Pang (2003a, 2003b), Geiger and Kanzow (2002), Glowinski (1984), Glowinski, Lions and Trémolières (1981) and Haslinger, Hlaváček and Nečas (1996). We refer to the recent monographs by Han and Reddy (1999) and Han and Sofonea (2002) and the references therein for an overview of the mathematical theory and numerical analysis for inequality problems in continuum mechanics.

Our formulation will be based on a primal-dual pair of variables. In addition to the displacement which represents the primal variable, the surface traction on the possible contact zone is introduced as dual variable: see, e.g., Christensen, Klarbring, Pang and Strömberg (1998). This new pair of variables has to be admissible, i.e., satisfy the inequality constraints arising from the non-penetration condition and the friction law.

In this section, we provide the setting of a quasi-static frictional contact problem between elastic bodies. Figure 2.1 shows the stress components $\sigma_{xx}$ and $\sigma_{xy}$ for two different situations in the case of three elastic bodies in contact.

Here, the contact of a deformable body with a rigid obstacle has been taken into account as well as the contact between deformable bodies. A fully symmetric situation is shown in Figure 2.1(a,b), and no Dirichlet boundary condition is imposed. The rigid body motions are fixed by the non-penetration condition and a zero tangential displacement of the centre of the upper circle. In Figure 2.1(c,d), there is an additional rigid obstacle on the right of the three circles. Then all rigid body modes are automatically fixed by the contact conditions.

To simplify the notation, we restrict our attention to two bodies, linear elasticity in the compressible range and a given constant Coulomb friction coefficient. However, most of our algorithmic results can easily be extended to more complex situations. We refer to the early papers by Laursen and Simo (1993a), Oden, Becker, Lin and Demkowicz (1985), Puso and Laursen (2004a, 2004b), Puso, Laursen and Solberg (2008), Yang and
<figure><img src="image_2.png" /><figcaption>Figure 2.1. Stress component $\sigma_{xx}$ and $\sigma_{xy}$ for a symmetric (a,b) and a non-symmetric setting (c,d).</figcaption></figure>

Laursen (2008a, 2008b) and Yang, Laursen and Meng (2005) for large deformation contact discretizations on non-matching meshes and to the recent contributions on solvers by Gitterle, Popp, Gee and Wall (2010), Popp, Gitterle, Gee and Wall (2010), Popp, Gee and Wall (2009) and Krause and Mohr (2011). Numerical examples also illustrate the performance of these approaches in more general formulations, e.g., for nearly incompressible materials and the inclusion of thermal effects with a temperature-dependent friction coefficient.

The two bodies in the reference configuration are given by two open bounded domains $\Omega^s$ and $\Omega^m \subset \mathbb{R}^d$, $d = 2, 3$, with Lipschitz boundary $\partial\Omega^s$ and $\partial\Omega^m$, respectively. The notation is adapted to the standard mortar framework, i.e., the upper index s stands for the slave side, and the index m refers to the master side. The contact conditions will be imposed weakly in terms of Lagrange multipliers defined on the slave side. Thus the displacement on the slave side has to follow the displacement of the master side in the event that the constraints are active. This observation motivates the terminology. The boundary $\partial\Omega^k$ is partitioned into three open disjoint measurable parts $\Gamma_D^k, \Gamma_N^k$, and $\Gamma_C^k$, $k \in \{\text{m}, \text{s}\}$. Dirichlet conditions will be set on $\Gamma_D^k$ and Neumann data on $\Gamma_N^k$. For simplicity of notation, we assume firstly that meas($\Gamma_D^k$) > 0, $k \in \{\text{s}, \text{m}\}$, secondly that $\Gamma_D^s$ is compactly embedded in $\partial\Omega^s \setminus \bar{\Gamma}_C^s$, and thirdly that the actual contact zone $B_n \subset \Gamma_C^s$ is compactly embedded in $\Gamma_C^s$. The first assumption on the Dirichlet boundary part means that Korn's inequality holds on each body, and thus that we do not have to deal with extra rigid body motions. The second assumption guarantees that the trace space restricted to $\Gamma_C^s$ does not see any boundary condition originating from $\Gamma_D^s$. As we will see later, the third assumption on the actual contact zone guarantees that the support of the surface traction on $\Gamma_C^s$ is compactly embedded in $\Gamma_C^s$.

Throughout this paper, we use the standard notation for the Sobolev space $H^s(\omega)$, $s \ge 1$, where $\omega$ is a suitable subdomain of $\Omega^k$, $k \in \{\text{s}, \text{m}\}$, and denote the associated norm with $\| \cdot \|_{s;\omega}$. The broken $H^s$-norm on $\Omega$ is given by $\|v\|_{s;\Omega}^2 := \|v^s\|_{s;\Omega^s}^2 + \|v^m\|_{s;\Omega^m}^2$ for $v := (v^s, v^m) \in H^s(\Omega^s) \times H^s(\Omega^m)$.
On $(d-1)$-dimensional manifolds $\gamma$ such as $\Gamma_C^s$, we use the Sobolev space $H^s(\gamma)$, $s \ge 0$, and its dual space $H^{-s}(\gamma)$. We point out that in our notation $H^{-1/2}(\gamma)$ is not the dual space of $H_{00}^{1/2}(\gamma)$ but of $H^{1/2}(\gamma)$. The dual norm is defined in the standard way by

$$ \|\mu\|_{-s;\gamma} := \sup_{v \in H^s(\gamma)} \frac{\langle \mu, v \rangle_{s;\gamma}}{\|v\|_{s;\gamma}}, \quad s \ge 0, \qquad (2.1) $$

where $\langle \cdot, \cdot \rangle_{s;\gamma}$ denotes the duality pairing. We note that the second assumption on the Dirichlet boundary part allows one to work with $H^{1/2}(\Gamma_C^s)$. Otherwise, we would have to consider the more complex $H_{00}^{1/2}(\Gamma_C^s)$ space, as in the mortar framework with cross-points (Bernardi, Maday and Patera 1993, 1994). We refer to the recent monograph on the theory and implementation of mortar methods by Lacour and Ben Belgacem (2011).

In what follows, we shall frequently use the generic constants $0 < c, C < \infty$, which are independent of the mesh size but possibly depend on the regularity of the domain or the mesh. Vectorial quantities are written in bold, e.g., $\mathbf{x}$, $\mathbf{y}$, and for simplicity of notation $\mathbf{x}\mathbf{y}$ stands for the scalar product between $\mathbf{x}$ and $\mathbf{y}$. Tensorial quantities are represented by bold greek symbols.

## 2.1. Problem formulation in its strong form

For the moment we restrict ourselves on each body to a homogeneous isotropic linearized Saint Venant–Kirchhoff material and also to the small strain assumption. Then, the strain-displacement relation is defined by $\epsilon(\mathbf{v}) := 1/2(\nabla\mathbf{v} + (\nabla\mathbf{v})^\top)$ and the constitutive equation for the stress tensor is given in terms of the fourth-order Hooke tensor $C$ by

$$ \sigma(\mathbf{v}) := \lambda \operatorname{tr}(\epsilon(\mathbf{v}))\mathbf{Id} + 2\mu\epsilon(\mathbf{v}) := C\epsilon(\mathbf{v}). \qquad (2.2) $$

Here $\operatorname{tr}$ denotes the trace operator and $\mathbf{Id}$ the identity in $\mathbb{R}^{d \times d}$. The positive coefficient $\lambda$ and the shear modulus $\mu$ are the Lamé parameters, which are assumed to be constant in each subdomain $\Omega^k$, $k \in \{\text{s}, \text{m}\}$, but have possibly quite different values on the slave and master side. We note that the Lamé parameters can be easily calculated from the Poisson ratio and Young's modulus.

Then, the linearized elastic equilibrium condition for the displacement $\mathbf{u} := (\mathbf{u}^\text{m}, \mathbf{u}^\text{s})$ can be written as

$$ \begin{align} - \operatorname{div} \sigma(\mathbf{u}) &= \mathbf{f} \quad \text{in } \Omega, \nonumber \\ \mathbf{u} &= \mathbf{u}_\text{D} \quad \text{on } \Gamma_\text{D} := \Gamma_\text{D}^\text{m} \cup \Gamma_\text{D}^\text{s}, \tag{2.3} \\ \sigma(\mathbf{u})\mathbf{n} &= \mathbf{f}_\text{N} \quad \text{on } \Gamma_\text{N} := \Gamma_\text{N}^\text{m} \cup \Gamma_\text{N}^\text{s}, \nonumber \end{align} $$

where $\mathbf{n}$ stands for the outer unit normal vector, which is almost everywhere well-defined. Here, the volume force $\mathbf{f}$, the Neumann data $\mathbf{f}_\text{N}$, and
the Dirichlet condition $\mathbf{u}_D$ are assumed to be in $(L^2(\Omega))^d$, $(L^2(\Gamma_N))^d$ and $(C(\Gamma_D^s) \cap H^{1/2}(\Gamma_D^s))^d \times (C(\Gamma_D^m) \cap H^{1/2}(\Gamma_D^m))^d$, respectively. Moreover, we assume that $C_{\text{reg}} < \infty$ exists such that

$$ \sup_{\substack{\mathbf{v} \in (H^1(\Omega^s))^d \\ \mathbf{v}|_{\Gamma_D^s} = \mathbf{0}}} \frac{\int_{\Omega^s} \mathbf{f}\mathbf{v} \, dx + \int_{\Gamma_N^s} \mathbf{f}_N \mathbf{v} \, ds}{\|\mathbf{v}\|_{1;\Omega^s}} \le C_{\text{reg}} \sup_{\substack{\mathbf{v} \in (H^1(\Omega^s))^d \\ \mathbf{v}|_{\Gamma_D^s \cup \Gamma_C^s} = \mathbf{0}}} \frac{\int_{\Omega^s} \mathbf{f}\mathbf{v} \, dx + \int_{\Gamma_N^s} \mathbf{f}_N \mathbf{v} \, ds}{\|\mathbf{v}\|_{1;\Omega^s}}. \quad (2.4) $$

These regularity assumptions on the data can be considerably weakened but for most examples these hold.

In addition to (2.3), we have to satisfy the contact constraints on $\Gamma_C^s$: the linearized non-penetration condition in the normal direction and the friction law in the tangential direction. These constraints can be formulated by means of the displacement and the surface traction $\lambda := -\sigma(\mathbf{u}^s)\mathbf{n}^s$. The linearized non-penetration condition reads as

$$ [u_n] \le g, \quad \lambda_n \ge 0, \quad \lambda_n([u_n] - g) = 0, \quad (2.5) $$

where $g \in H^{1/2}(\Gamma_C^s)$ is the linearized gap function between the two deformable bodies. The linearized setting can be expressed in terms of the normal contributions with respect to the reference configuration. Here $\lambda_n := \lambda \mathbf{n}^s$ is the normal component of the boundary stress, and $[u_n] := (\mathbf{u}^s - \mathbf{u}^m \circ \chi)\mathbf{n}^s$ is the jump of the mapped boundary displacements, where $\chi(\cdot)$ denotes a suitable mapping from $\Gamma_C^s$ onto $\Gamma_C^m$.

In addition to (2.5), we have to satisfy the quasi-static Coulomb law

$$ \|\boldsymbol{\lambda}_t\| \le \nu\lambda_n, \quad [\dot{\mathbf{u}}_t]\boldsymbol{\lambda}_t - \nu\lambda_n \|[\dot{\mathbf{u}}_t]\| = 0, \quad (2.6) $$

where the tangential components are defined by $\boldsymbol{\lambda}_t := \boldsymbol{\lambda} - \lambda_n \mathbf{n}^s$ and $[\mathbf{u}_t] := [\mathbf{u}] - [u_n]\mathbf{n}^s$, $[\mathbf{u}] := \mathbf{u}^s - \mathbf{u}^m \circ \chi$, $\nu \ge 0$ is the friction coefficient, and $\| \cdot \|$ stands for the Euclidean norm. We note that for the dynamic case, inertia terms have to be included, and the volume mass density of the two bodies has to be taken into account: see Section 7. For the moment we focus on the static Coulomb law, i.e., we replace the tangential velocity by the tangential displacement in (2.6). This problem type has then to be solved in each time step if an implicit time-integration scheme is used.

Figure 2.2 illustrates the notation and the situation for finite deformations. In that case the Jacobian of the deformation mapping has to be taken into account and the non-penetration has to be formulated with respect to the actual configuration.

We recall that in the reference configuration $\Gamma_C^s$ and $\Gamma_C^m$ do not have to be matching, and thus the displacement from the master side has to be projected onto the slave side. Moreover, the contact surface tractions on the master and slave body have to be in equilibrium in the actual configuration. In the case of linear elasticity, reference and actual configuration can be
<figure><img src="image_3.png" /><figcaption>Figure 2.2. Illustration of the notation.</figcaption></figure>

identified, and thus the approach is simplified considerably compared to the
case of finite deformations.

We note that the Signorini problem with Coulomb friction was introduced
by Duvaut and Lions (1976); see also Fichera (1964) and Lions and Stam-
pacchia (1967). Although it is widely studied, not all aspects are yet fully
understood and some open questions remain. It is well known that frictional
contact problems do not necessarily have a unique weak solution. Many of
the early results are designed for the case of one elastic body in contact
with a rigid foundation. First existence results can be found in Demkowicz
and Oden (1982), Jarušek (1983) and Nečas, Jarušek and Haslinger (1980).
An alternative proof based on a penalization technique is given in Eck and
Jarušek (1998). Examples for non-uniqueness are given in Ballard (1999)
and Hild (2003, 2004), and uniqueness criteria are studied, e.g., in Ballard
and Basseville (2005), Hild and Renard (2006) and Renard (2006). We re-
fer to the recent monograph by Eck et al. (2005) for an excellent overview
of existence and uniqueness results for these contact problems. Roughly
speaking, for $\nu \ge 0$ small enough, the existence and uniqueness of a weak
solution is guaranteed.

Figure 2.3 shows the influence of the friction coefficient $\nu$ on the normal
and tangential component of the surface traction for the classical Hertz
problem (Hertz 1882). In between the maximal and minimal value of the
tangential contact traction, the two bodies stick together. On the rest of the
actual contact zone, a relative tangential displacement occurs. As can be
seen directly from (2.6) for $\nu = 0$, we have $\lambda_t = 0$, which is also obtained by
the numerical scheme. We observe that the size of the slippy contact zone
is largely influenced by $\nu$. The smaller $\nu$, the larger is the slippy contact
zone, whereas the contact radius is not very sensitive with respect to $\nu$.
<figure><img src="image_4.png" /><figcaption>Figure 2.3. Surface traction of a Hertz contact problem with Coulomb friction for different values of the friction coefficient $\nu \in \{0, 0.05, 0.15, 0.3, 0.5\}$.</figcaption></figure>

## 2.2. Formulation of the weak problem

We start with a detailed discussion of different but equivalent formulations of the contact problem with no friction, i.e., $\nu = 0$. Then the problem is equivalent to a standard variational inequality of the first kind or to an energy-minimization problem on a convex set. The frictionless contact problem with a linearized non-penetration condition can be stated as follows. Find $\mathbf{u} \in \mathbf{K}$ such that

$$J(\mathbf{u}) = \inf_{\mathbf{v} \in \mathbf{K}} J(\mathbf{v}), \qquad (2.7)$$

where the convex set $\mathbf{K}$ is given by all admissible solutions, and the energy is defined by $J(\mathbf{v}) := \frac{1}{2}a(\mathbf{v}, \mathbf{v}) - f(\mathbf{v})$. Here, the bilinear form $a(\cdot, \cdot)$ is given, for $\mathbf{v}, \mathbf{w} \in \mathbf{V} := \mathbf{V}^m \times \mathbf{V}^s := (H^1(\Omega^m))^d \times (H^1(\Omega^s))^d$, by

$$a(\mathbf{w}, \mathbf{v}) := a_m(\mathbf{w}^m, \mathbf{v}^m) + a_s(\mathbf{w}^s, \mathbf{v}^s), \\ a_k(\mathbf{w}^k, \mathbf{v}^k) := \int_{\Omega^k} \sigma(\mathbf{w}^k) : \epsilon(\mathbf{v}^k) \mathrm{d}\mathbf{x}, \quad k \in \{\mathrm{m}, \mathrm{s}\},$$

and the linear form $f(\cdot)$ is defined for all $\mathbf{v} \in \mathbf{V}$ in terms of

$$f(\mathbf{v}) := f_s(\mathbf{v}^s) + f_m(\mathbf{v}^m), \quad f_k(\mathbf{v}^k) := \int_{\Omega^k} \mathbf{f}\mathbf{v}^k \mathrm{d}\mathbf{x} + \int_{\Gamma_N^k} \mathbf{f}_N \mathbf{v}^k \mathrm{d}\mathbf{s}.$$

In addition to the Hilbert space $\mathbf{V}$, we introduce the subset $\mathbf{K} \subset \mathbf{V}$:

$$\mathbf{K} := \{\mathbf{v} \in \mathbf{V} \mid \mathbf{v} = \mathbf{u}_D \text{ on } \Gamma_D \text{ and } [v_n] \le g \text{ on } \Gamma_C^s\}.$$

By definition $\mathbf{K}$ is a closed convex non-empty set, $f(\cdot)$ is a continuous linear functional on $\mathbf{V}$, and the bilinear form $a(\cdot, \cdot)$ is continuous on $\mathbf{V} \times \mathbf{V}$ and elliptic with respect to the Hilbert space $\mathbf{V}_0 := \{\mathbf{v} \in \mathbf{V}, \mathbf{v} = \mathbf{0} \text{ on } \Gamma_D\}$.
Thus the constrained minimization problem (2.7) has a unique solution: see, e.g., Glowinski (1984). Due to the symmetry of the bilinear form $a(\cdot, \cdot)$, (2.7) can be equivalently written as a variational inequality of the first kind, i.e., find $\mathbf{u} \in \mathbf{K}$ such that

$$a(\mathbf{u}, \mathbf{v} - \mathbf{u}) \geq f(\mathbf{v} - \mathbf{u}), \quad \mathbf{v} \in \mathbf{K}, \qquad (2.8)$$

or as a variational inequality of the second kind, i.e., find $\mathbf{u} \in \mathbf{V}$ such that

$$a(\mathbf{u}, \mathbf{v} - \mathbf{u}) + \chi_{\mathbf{K}}(\mathbf{v}) - \chi_{\mathbf{K}}(\mathbf{u}) \geq f(\mathbf{v} - \mathbf{u}), \quad \mathbf{v} \in \mathbf{V}, \qquad (2.9)$$

where $\chi_{\mathbf{K}}$ is the indicator functional of $\mathbf{K}$, i.e., $\chi_{\mathbf{K}}(\mathbf{v}) := \infty$ if $\mathbf{v} \notin \mathbf{K}$ and zero otherwise. We refer to Brezis (1971), Duvaut and Lions (1976), Fichera (1964), Glowinski (1984), Glowinski et al. (1981), Kinderlehrer and Stampacchia (2000) and the references therein for an abstract mathematical framework on inequalities as well as for the so-called Signorini problem and its physical and mechanical interpretation. It is easy to see that in the special case of a variational inequality of the first kind with the convex set being a Hilbert space, existence and uniqueness of a solution follow directly from the Lax-Milgram theorem. Alternatively to the pure displacement-based formulation, the Signorini problem can be characterized in terms of the contact pressure as unknown variable; see, e.g., Demkowicz (1982).

We note that the convex set $\mathbf{K}$ can be characterized in terms of a dual cone. To do so, we introduce the dual space $\mathbf{M} := (M)^d := (W')^d := \mathbf{W}'$ of the trace space $\mathbf{W} := W^d := (H^{1/2}(\Gamma_{\mathbf{C}}^{\mathrm{s}}))^d$ and define the bilinear form $b(\cdot, \cdot)$ by

$$b(\boldsymbol{\mu}, \mathbf{v}) := \langle \boldsymbol{\mu}, [\mathbf{v}] \rangle_{\Gamma_{\mathbf{C}}^{\mathrm{s}}}, \quad \mathbf{v} \in \mathbf{V}, \boldsymbol{\mu} \in \mathbf{M},$$

where $\langle \cdot, \cdot \rangle_{\Gamma_{\mathbf{C}}^{\mathrm{s}}}$ stands for the $H^{1/2}$-duality pairing on $\Gamma_{\mathbf{C}}^{\mathrm{s}}$, and $[\mathbf{v}] := \mathbf{v}^{\mathrm{s}} - \mathbf{v}^{\mathrm{m}} \circ \chi$. It is assumed that $\chi$ is smooth enough such that for $\mathbf{v} \in \mathbf{V}$ we have $[\mathbf{v}] \in \mathbf{W}$. In terms of the bilinear form $b(\cdot, \cdot)$, the closed non-empty convex cone $\mathbf{M}^{+}$ is set to

$$\mathbf{M}^{+} := \{\boldsymbol{\mu} \in \mathbf{M} \mid \langle \boldsymbol{\mu}, \mathbf{w} \rangle_{\Gamma_{\mathbf{C}}^{\mathrm{s}}} \geq 0, \mathbf{w} \in \mathbf{W}_{n}^{+}\}, \qquad (2.10a)$$

$$\mathbf{W}_{n}^{+} := \{\mathbf{w} \in \mathbf{W} \mid w_n \in W^{+}\}, \quad W^{+} := \{w \in W \mid w \geq 0\}. \qquad (2.10b)$$

Here, we assume that $\mathbf{n}^{\mathrm{s}}$ is smooth enough such that for $\mathbf{w} \in \mathbf{W}$ and $w \in W$ we also have $w_n \in W$ and $w\mathbf{n}^{\mathrm{s}} \in \mathbf{W}$, respectively. Then the definitions of $\mu_n \in M$ for $\boldsymbol{\mu} \in \mathbf{M}$ and of $\mu_n \mathbf{n}^{\mathrm{s}} \in \mathbf{M}$ for $\mu_n \in M$ given by

$$\langle \mu_n, w \rangle_{\Gamma_{\mathbf{C}}^{\mathrm{s}}} := \langle \boldsymbol{\mu}, w\mathbf{n}^{\mathrm{s}} \rangle_{\Gamma_{\mathbf{C}}^{\mathrm{s}}}, \quad w \in W, \\ \langle \mu_n \mathbf{n}^{\mathrm{s}}, \mathbf{w} \rangle_{\Gamma_{\mathbf{C}}^{\mathrm{s}}} := \langle \mu_n, w_n \rangle_{\Gamma_{\mathbf{C}}^{\mathrm{s}}}, \quad \mathbf{w} \in \mathbf{W},$$

respectively, are consistent in the sense that $\langle \mu_n, w_n \rangle_{\Gamma_{\mathbf{C}}^{\mathrm{s}}} = \langle \mu_n \mathbf{n}^{\mathrm{s}}, w_n \mathbf{n}^{\mathrm{s}} \rangle_{\Gamma_{\mathbf{C}}^{\mathrm{s}}}$ for all $\boldsymbol{\mu} \in \mathbf{M}$ and $\mathbf{w} \in \mathbf{W}$. As a result the bilinear form $b(\cdot, \cdot)$ can be split
in a well-defined normal and tangential part,

$$b(\boldsymbol{\mu}, \boldsymbol{v}) = b_n(\boldsymbol{\mu}, \boldsymbol{v}) + b_t(\boldsymbol{\mu}, \boldsymbol{v}) \quad (2.11)$$

with $b_n(\boldsymbol{\mu}, \boldsymbol{v}) := b(\mu_n \mathbf{n}^s, \boldsymbol{v})$, $b_t(\boldsymbol{\mu}, \boldsymbol{v}) := b(\mu_t, \boldsymbol{v})$, $\mu_t := \boldsymbol{\mu} - \mu_n \mathbf{n}^s$.

The situation that $\mathbf{n}^s$ is only piecewise well-defined on $\Gamma_C^s$ can be handled by decomposing the contact zone into non-overlapping subparts $\gamma_j \subset \Gamma_C^s$, defining all quantities with respect to $\gamma_j$ and using product spaces and broken duality pairings.

Now, observing that $\mathbf{K}$ can be written as

$$\mathbf{K} = \{\boldsymbol{v} \in \mathbf{V} \mid \boldsymbol{v} = \mathbf{u}_D \text{ on } \Gamma_D \text{ and } b(\boldsymbol{\mu}, \boldsymbol{v}) \le g(\boldsymbol{\mu}), \boldsymbol{\mu} \in \mathbf{M}^+\},$$

where $g(\boldsymbol{\mu}) := \langle \mu_n, g \rangle_{\Gamma_C^s}$, $\boldsymbol{\mu} \in \mathbf{M}$, we obtain the saddle-point formulation of (2.8). Find $(\mathbf{u}, \boldsymbol{\lambda}) \in \mathbf{V}_D \times \mathbf{M}^+$ such that

$$\begin{aligned} a(\mathbf{u}, \boldsymbol{v}) + b(\boldsymbol{\lambda}, \boldsymbol{v}) &= f(\boldsymbol{v}), & \boldsymbol{v} &\in \mathbf{V}_0, \\ b(\boldsymbol{\mu} - \boldsymbol{\lambda}, \mathbf{u}) &\le g(\boldsymbol{\mu} - \boldsymbol{\lambda}), & \boldsymbol{\mu} &\in \mathbf{M}^+, \end{aligned} \quad (2.12)$$

with the convex set $\mathbf{V}_D := \{\boldsymbol{v} \in \mathbf{V}, \boldsymbol{v} = \mathbf{u}_D \text{ on } \Gamma_D\}$.

**Lemma 2.1.** The three inequality formulations (2.8), (2.9) and (2.12) are equivalent in the sense that if $(\mathbf{u}, \boldsymbol{\lambda})$ solves (2.12), then $\mathbf{u}$ is the solution of (2.8) and (2.9), and if $\mathbf{u}$ solves (2.8) or (2.9), then $(\mathbf{u}, \boldsymbol{\lambda})$, with $\boldsymbol{\lambda} \in \mathbf{M}$,

$$\langle \boldsymbol{\lambda}, \mathbf{w} \rangle_{\Gamma_C^s} := f(\mathcal{H}\mathbf{w}) - a(\mathbf{u}, \mathcal{H}\mathbf{w}), \quad \mathbf{w} \in \mathbf{W} \quad (2.13)$$

satisfies (2.12). Here we have used $\mathcal{H}\mathbf{w} := (\mathbf{0}, \mathcal{H}_s\mathbf{w})$, where $\mathcal{H}_s$ is the harmonic extension onto $\mathbf{V}_0^s := \{\boldsymbol{v} \in \mathbf{V}^s \mid \boldsymbol{v} = \mathbf{0} \text{ on } \Gamma_D^s\}$ with respect to the bilinear form $a_s(\cdot, \cdot)$.

*Proof.* For convenience of the reader we recall some of the basic steps and refer to the monographs by Glowinski (1984) and Glowinski et al. (1981) for further details. In particular, we comment on the formula (2.13) for the Lagrange multiplier. The equivalence between (2.8) and (2.9) is standard.

Let $(\mathbf{u}, \boldsymbol{\lambda})$ be a solution of (2.12); then for all $\boldsymbol{\mu} \in \mathbf{M}^+$ we have $\boldsymbol{\mu}+\boldsymbol{\lambda} \in \mathbf{M}^+$ and thus $\mathbf{u} \in \mathbf{K}$, and moreover $b(\boldsymbol{\lambda}, \mathbf{u}) = g(\boldsymbol{\lambda})$. For $\boldsymbol{v} \in \mathbf{K}$ we find $\boldsymbol{v}-\mathbf{u} \in \mathbf{V}_0$ and thus $a(\mathbf{u}, \boldsymbol{v}-\mathbf{u}) = f(\boldsymbol{v}-\mathbf{u}) - b(\boldsymbol{\lambda}, \boldsymbol{v}-\mathbf{u}) \ge f(\boldsymbol{v}-\mathbf{u}) + g(\boldsymbol{\lambda}) - g(\boldsymbol{\lambda})$.

Let $\mathbf{u}$ be the solution of (2.8). Then, for all $\mathbf{w} \in \mathbf{W}_n^+$, we have $\boldsymbol{v} := \mathbf{u} - \mathcal{H}\mathbf{w} \in \mathbf{K}$. Now the definition (2.13) of $\boldsymbol{\lambda}$ yields $\langle \boldsymbol{\lambda}, \mathbf{w} \rangle_{\Gamma_C^s} = a(\mathbf{u}, \boldsymbol{v} - \mathbf{u}) - f(\boldsymbol{v} - \mathbf{u}) \ge 0$ for all $\mathbf{w} \in \mathbf{W}_n^+$ and thus $\boldsymbol{\lambda} \in \mathbf{M}^+$. Moreover, observing that $\boldsymbol{v}_{\pm} := \mathbf{u} \pm \mathcal{H}((g - [u_n])\mathbf{n}^s)$ is in $\mathbf{K}$, we get $a(\mathbf{u}, \mathcal{H}((g - [u_n])\mathbf{n}^s)) = f(\mathcal{H}((g - [u_n])\mathbf{n}^s))$, from which we conclude that

$$0 = \langle \boldsymbol{\lambda}, (g - [u_n])\mathbf{n}^s \rangle_{\Gamma_C^s} = \langle \lambda_n, g \rangle_{\Gamma_C^s} - b_n(\boldsymbol{\lambda}, \mathbf{u}) = g(\boldsymbol{\lambda}) - b(\boldsymbol{\lambda}, \mathbf{u}).$$

Then the second line of (2.12) holds by the definition of $\mathbf{K}$. To see that (2.13) also satisfies the first line of (2.12), we set $\mathbf{w} := \mathbf{u} \pm (\boldsymbol{v} - \mathcal{H}[\boldsymbol{v}]) \in \mathbf{K}$
for $\mathbf{v} \in \mathbf{V}_0$ and use $\mathbf{w}$ as a test function in (2.8), resulting in

$$0 = a(\mathbf{u}, \mathbf{v} - \mathcal{H}[\mathbf{v}]) - f(\mathbf{v} - \mathcal{H}[\mathbf{v}]) = a(\mathbf{u}, \mathbf{v}) - f(\mathbf{v}) + \langle \boldsymbol{\lambda}, [\mathbf{v}] \rangle_{\Gamma_C^s}. \quad \square$$

**Remark 2.2.** We note that the saddle-point formulation (2.12) also has a unique solution. The uniqueness of the displacement is already established by Lemma 2.1. For the uniqueness of the surface traction a suitable inf-sup condition has to be satisfied. By definition, $\mathbf{M}$ is the dual space of $\mathbf{W}$, which is the trace space of $\mathbf{V}_0^s := \{\mathbf{v} \in \mathbf{V}^s | \mathbf{v} = \mathbf{0} \text{ on } \Gamma_D^s\}$, and the extension theorem yields that

$$\begin{aligned} \inf_{\boldsymbol{\mu} \in \mathbf{M}} \sup_{\mathbf{v} \in \mathbf{V}_0} \frac{b(\boldsymbol{\mu}, \mathbf{v})}{\|\boldsymbol{\mu}\|_{-\frac{1}{2}; \Gamma_C^s} \|\mathbf{v}\|_{1; \Omega}} &\ge \inf_{\boldsymbol{\mu} \in \mathbf{M}} \sup_{\mathbf{v} \in \mathbf{V}_0^s} \frac{b(\boldsymbol{\mu}, \mathbf{v})}{\|\boldsymbol{\mu}\|_{-\frac{1}{2}; \Gamma_C^s} \|\mathbf{v}\|_{1; \Omega^s}} \\ &\ge C_{\inf} \inf_{\boldsymbol{\mu} \in \mathbf{M}} \sup_{\mathbf{w} \in \mathbf{W}} \frac{\langle \boldsymbol{\mu}, \mathbf{w} \rangle_{\Gamma_C^s}}{\|\boldsymbol{\mu}\|_{-\frac{1}{2}; \Gamma_C^s} \|\mathbf{w}\|_{\frac{1}{2}; \Gamma_C^s}} = C_{\inf}. \end{aligned}$$

The case of Coulomb friction is more involved; we refer to Eck *et al.* (2005) for existence and regularity results and only mention that, for $\nu$ small enough, a unique solution exists. For contact problems in viscoelasticity we refer to Eck and Jarušek (2003) and Han and Sofonea (2002). In particular, the admissible solution space depends on the solution itself and cannot be characterized without knowledge of the contact pressure.

After these preliminary remarks, we can now easily extend our saddle-point formulation (2.12) for $\nu = 0$ to $\nu \ge 0$. We observe that $\mathbf{M}^+$ defined by (2.10a) can also be written as

$$\mathbf{M}^+ = \{\boldsymbol{\mu} \in \mathbf{M} | \boldsymbol{\mu}_n \in \mathbf{M}^+, \boldsymbol{\mu}_t = \mathbf{0}\},$$

with $\mathbf{M}^+ := \{\boldsymbol{\mu} \in \mathbf{M} | \langle \boldsymbol{\mu}, \mathbf{w} \rangle_{\Gamma_C^s} \ge 0, \mathbf{w} \in \mathbf{W}^+\}$. For $\nu > 0$ the tangential part of the surface traction, in general, does not vanish, and thus one has to work with a vectorial Lagrange multiplier, which is not necessary for a frictionless contact problem. Replacing the convex cone $\mathbf{M}^+$ in (2.12) by

$$\mathbf{M}(\lambda_n) := \{\boldsymbol{\mu} \in \mathbf{M} | \langle \boldsymbol{\mu}, \mathbf{v} \rangle_{\Gamma_C^s} \le \langle \nu \lambda_n, \|\mathbf{v}_t\| \rangle_{\Gamma_C^s}, \mathbf{v} \in \mathbf{W} \text{ with } -\nu_n \in \mathbf{W}^+\}, \quad (2.14)$$

we obtain the weak saddle-point formulation of a static Coulomb problem between two linearly elastic bodies as follows. Find $(\mathbf{u}, \boldsymbol{\lambda}) \in \mathbf{V}_D \times \mathbf{M}(\lambda_n)$ such that

$$\begin{aligned} a(\mathbf{u}, \mathbf{v}) + b(\boldsymbol{\lambda}, \mathbf{v}) &= f(\mathbf{v}), & \mathbf{v} \in \mathbf{V}_0, \\ b(\boldsymbol{\mu} - \boldsymbol{\lambda}, \mathbf{u}) &\le g(\boldsymbol{\mu} - \boldsymbol{\lambda}), & \boldsymbol{\mu} \in \mathbf{M}(\lambda_n). \end{aligned} \quad (2.15)$$

In the case of the quasi-static version, one has to replace $b(\boldsymbol{\mu} - \boldsymbol{\lambda}, \mathbf{u})$ in the second line of (2.15) by $b_n(\boldsymbol{\mu} - \boldsymbol{\lambda}, \mathbf{u}) + b_t(\boldsymbol{\mu} - \boldsymbol{\lambda}, \dot{\mathbf{u}})$. Comparing (2.12) and (2.15), we find that the only, but essential, difference is the solution cone for the Lagrange multiplier $\boldsymbol{\lambda}$. The key idea for the proof of existence is to
define a series of solutions $(\mathbf{u}^k, \lambda^k) \in \mathbf{V}_D \times \mathbf{M}(\lambda_n^{k-1})$ and apply Tikhonov's fixed-point theorem: see, e.g., Eck et al. (2005).

In the following we will frequently make use of the Karush-Kuhn-Tucker (KKT) conditions.

**Lemma 2.3.** Let $(\mathbf{u}, \lambda) \in \mathbf{V}_D \times \mathbf{M}(\lambda_n)$ be the solution of (2.15); then the non-penetration KKT condition

$$ \lambda_n \in M^+, \quad g - [u_n] \in W^+, \quad b_n(\lambda, \mathbf{u}) = \langle \lambda_n, g \rangle_{\Gamma_C^s} \quad (2.16) $$

holds. Moreover, under suitable regularity, the Coulomb law in its weak form

$$ \lambda \in \mathbf{M}(\lambda_n), \quad b_t(\lambda, \mathbf{u}) = \langle \nu\lambda_n, \|[\mathbf{u}_t]\| \rangle_{\Gamma_C^s} \quad (2.17) $$

is satisfied.

*Proof.* We observe that the constraint on $\lambda_n$ in (2.16) follows directly from the definition of $\mathbf{M}(\lambda_n)$. Using the additive splitting (2.11) and setting as test function $\mu := \lambda \pm \lambda_n \mathbf{n}^s$ in (2.15), we get the equality in (2.16). Observing that $W^+$, defined in (2.10b), can also be characterized by $W^+ = \{w \in W \mid \langle \mu, w \rangle_{\Gamma_C^s} \ge 0, \mu \in M^+\}$, it trivially holds that $g - [u_n] \in W^+$. For (2.17) we assume that there exists a $\chi$ such that $\|\chi\| \le 1$, $\chi[\mathbf{u}_t] = \|[\mathbf{u}_t]\|$, and we have $\chi\mathbf{v} \in W$ for $\mathbf{v} \in \mathbf{W}$. Then, we get $\mu := \lambda_n \mathbf{n}^s + \nu\lambda_n\chi \in \mathbf{M}(\lambda_n)$ and $b(\mu - \lambda, \mathbf{u}) = b(\nu\lambda_n\chi, \mathbf{u}) - b(\lambda_t, \mathbf{u}) = \langle \nu\lambda_n, \|[\mathbf{u}_t]\| \rangle_{\Gamma_C^s} - b_t(\lambda, \mathbf{u}) \le 0$, from which the equality in (2.17) follows from the definition (2.14). $\square$

**Remark 2.4.** The special case of a contact problem between one elastic body and a rigid obstacle can be obtained from the two-body situation. A rigid body can be regarded as an infinitively stiff elastic body, and thus the limit case $\lambda^m, \mu^m \to \infty$ results in a one-body case where formally $\mathbf{u}^m = \mathbf{0}$.

To conclude this section, we briefly comment on numerical stability issues in elasticity and on the extension to the case of a solution-dependent friction coefficient.

## 2.3. Nearly incompressible materials

In the nearly incompressible case, the Poisson ratio tends to 0.5 and thus the ratio between $\lambda$ and $\mu$ tends to infinity. The definition of the bilinear form $a(\cdot, \cdot)$ in terms of the linearized stress (2.2) shows that the continuity constant depends on $\max(\lambda, \mu)$, while the coercivity constant depends on $\mu$. As a consequence, a priori estimates for standard low-order finite elements involve large constants, and volumetric locking can be observed numerically. To handle such a case appropriately, special discretization techniques are required. Methods associated with the enrichment or enhancement of the strain or stress field by the addition of carefully chosen basis functions have proved to be highly effective and popular. The key work dealing with
<figure><img src="image_5.png" /><figcaption>Figure 2.4. Von Mises stress and deformed mesh of a contact between a soft nearly incompressible material (lower body) and two hard compressible ones (upper bodies): standard low-order (a,b) and Hu–Washizu-based (c,d) discretization.</figcaption></figure>

enhanced assumed strain formulations is Simo and Rifai (1990). Figure 2.4 shows the setting of a contact problem between a nearly incompressible soft rubber-like material and a compressible hard one. In Figure 2.4(a,b) standard conforming low-order finite elements are applied, whereas in (c,d) special low-order Hu–Washizu-based elements are used. Here the starting point is the Hu–Washizu formulation (Hu 1955, Washizu 1955), in which the unknown variables are displacement, strain, and stress. This formulation (see also Felippa (2000) for some historical comments) can serve as the point of departure for the development of enhanced strain formulations; see also, e.g., Braess, Carstensen and Reddy (2004), Kasper and Taylor (2000a, 2000b), Simo and Armero (1992), Simo, Armero and Taylor (1993) and Simo and Rifai (1990).

For both discretization schemes the numerically obtained von Mises stress and a zoom of the deformed meshes are depicted in Figure 2.4. As can be clearly observed in the case of the standard scheme, volumetric locking occurs, resulting in a very stiff response of the soft material; only a modified scheme can provide a good approximation. In the numerical experiment, we use a pure displacement-based formulation obtained from local static condensation of a three-field formulation: see Lamichhane, Reddy and Wohlmuth (2006) for details.

To get a better feeling for the influence of Poisson’s ratio on the quality of the discretization, we consider the classical Hertz contact problem between a circle and a half-plane, which is approximated by a rectangle. In that situation the maximum of the contact pressure as well as the contact radius can be computed analytically in terms of the material parameters (Johnson 1985)

$$p_{\max} = \frac{2f}{\pi r_{\text{cont}}}, \quad r_{\text{cont}} = 2\sqrt{\frac{fr(1 - \nu_{\text{Poisson}}^2)}{E\pi}},$$
<table><caption>Table 2.1. Comparison of standard scheme (Q1) with a displacement-based Hu–Washizu method (HW) for different values of the Poisson ratio.</caption><thead><tr><th rowspan="2">Poisson ratio<br>ν<sub>Poisson</sub></th><th colspan="3">Contact pressure: p<sub>max</sub></th><th colspan="3">Contact radius: r<sub>cont</sub></th></tr><tr><th>Q1</th><th>HW</th><th>'Exact'</th><th>Q1</th><th>HW</th><th>'Exact'</th></tr></thead><tbody><tr><td>0.1</td><td>16.7840</td><td>16.7708</td><td>17.1564</td><td>0.3437</td><td>0.3750</td><td>0.3709</td></tr><tr><td>0.45</td><td>19.0973</td><td>19.0460</td><td>18.9238</td><td>0.3125</td><td>0.3125</td><td>0.3363</td></tr><tr><td>0.49</td><td>19.9255</td><td>19.5690</td><td>19.3374</td><td>0.2812</td><td>0.3125</td><td>0.3291</td></tr><tr><td>0.499</td><td>21.5690</td><td>19.7007</td><td>19.4394</td><td>0.2500</td><td>0.3125</td><td>0.3274</td></tr><tr><td>0.4999</td><td>31.0458</td><td>19.7144</td><td>19.4497</td><td>0.1875</td><td>0.3125</td><td>0.3272</td></tr><tr><td>0.49999</td><td>52.2592</td><td>19.7158</td><td>19.4508</td><td>0.1250</td><td>0.3125</td><td>0.3272</td></tr><tr><td>0.499999</td><td>58.4622</td><td>19.7159</td><td>19.4509</td><td>0.0938</td><td>0.3125</td><td>0.3272</td></tr></tbody></table>

where *r* is the radius of the circle, *f* the applied point load, *ν*<sub>Poisson</sub> the Poisson number and *E* Young's modulus.

A quantitative comparison of the two discretization schemes is given in Table 2.1. In the compressible range, both schemes provide quite good and accurate numerical approximations even for coarse meshes. However, the situation is drastically changed if Poisson's ratio tends to 0.5. From 0.49 on, from row to row, λ is increased by a factor of 10. In the case of standard conforming low-order elements, the contact radius tends to zero, and thus the maximal contact pressure tends to infinity. In contrast, the analytical solution as well as the Hu–Washizu-based formulation yield convergence to a non-zero contact radius, and the maximum of the contact pressure remains finite. The limit of the standard scheme is an unphysical point contact, with the surface traction being a delta distribution.

From now on, we assume that we are in the compressible range and that we do not have to face numerical problems due to the material parameters.

## 2.4. Thermo-mechanical contact problem

Coupled contact problems where the coefficient of friction depends on the solution itself are quite difficult to analyse. Although non-trivial from the theoretical point of view, these generalized settings do fit perfectly well into the computational framework. Figure 2.5 shows the temperature distribution of a sliding body undergoing thermo-mechanical contact. The discretization in space is based on non-matching meshes and no re-meshing has to be done. For the time integration we apply a simple mid-point rule in combination with standard mass lumping techniques; see Hüeber and Wohlmuth (2009) for details and further numerical results for this example.

In addition to the displacement, the temperature *T* is a primal variable, and a bi-directionally coupled thermo-mechanical system has to be con-
<figure><img src="image_6.png" /><figcaption>Figure 2.5. Temperature distribution at time $t^{k-1/2}$ for $k = 14, 26, 38, 50$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 2.6. Friction coefficient and temperature at the contact nodes on the cutting line $y = 0$ for $k = 14, 26, 38, 50$.</figcaption></figure>

More precisely, the relative temperature enters in terms of the thermal expansion coefficient in the definition of the mechanical stress. In addition, we have to consider the first and second law of thermodynamics. The heating from the Joule effect adds the source term $\text{div } \dot{\mathbf{u}}$ to the heat equation, and thus a fully coupled system is obtained: see, e.g., Fluegge (1972) and Willner (2003). From the theoretical point of view, thermo-mechanical contact problems have been analysed in Eck (2002) and Eck and Jarušek (2001).

Moreover a friction coefficient $\nu(T) \ge 0$, which is monotone decreasing in $T$, modelling a thermal softening effect, has been applied. It tends to zero in the critical case that the temperature tends to the damage temperature: see, e.g., Laursen (2002).

Figure 2.6 shows the evolution of the temperature-dependent friction coefficient and the temperature at the nodes in contact for different time steps. The dashed horizontal line marks the static coefficient of friction, and the two vertical lines indicate the actual contact zone. Due to the heating of the two bodies, the temperature increases over time and thus $\nu(T)$ decreases.
### 3. Variationally consistent space discretization

In the past, penalty methods and simple node-to-node coupling concepts have been widely used and are quite often integrated in commercial software codes. The starting point of penalty techniques is the observation that contact without friction can easily be formulated as a minimization problem on a constrained space. This approach is closely related to (2.9) when replacing the proper convex and lower semi-continuous indicator functional by a regularized one with finite values. Working on the larger unconstrained space, but incorporating the restriction as an additional finite energy contribution, a non-linear variational system of equations is obtained. Although relatively easy to implement, oscillations possibly occur, and the numerical results are very sensitive to the choice of the penalty parameter. A too-small penalty parameter gives considerable penetration and a poor approximation of the contact forces. A too-large penalty parameter yields a badly conditioned system which has to be solved by a suitable non-linear scheme. Figure 3.1 illustrates the influence of the penalty parameter on the deformed geometry.

For a small penalty parameter (see Figure 3.1(a)), significant penetration of the rigid obstacle into the elastic material occurs, and the von Mises stress is highly underestimated. In Figure 3.1(b), a very large penalty parameter is applied which gives numerical results of high quality. However, the number of multigrid steps then required to solve the system is quite high compared to the case of a small penalty parameter.

Figure 3.2 illustrates the difference between a penalty and a Lagrange multiplier-based approach. Figures 3.2(a) and 3.2(d) show, respectively, the influence of the penalty parameter $\epsilon$ on the penetration and the contact stress. As can be seen, for small penalty parameters there is significant penetration and the approximation of the contact stress is rather poor. In the limit $\epsilon \to 0$, the solution tends to the unconstrained one and the

<figure><img src="image_8.png" /><figcaption>Figure 3.1. Comparison of the deformed geometry and the von Mises stress for a three-dimensional one-sided contact problem: small penalty parameter (a) and large penalty parameter (b).</figcaption></figure>
<figure><img src="image_9.png" /></figure>

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

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

<figure><img src="image_12.png" /><figcaption>Figure 3.2. Penalty versus Lagrange multiplier approach.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 3.3. Node-to-node interpolation versus variationally consistent coupling.</figcaption></figure>
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>Figure 3.4. Non-optimal error decay for a node-to node coupling: scalar case (a) and linear elasticity (b).</figcaption></figure>

contact stress vanishes. In the limit $\epsilon \to \infty$, the correct constrained solution is recovered but then the condition number of the system tends to be extremely large. Figures 3.2(b) and 3.2(c) compare the numerical results of a penalty approach with a large value for $\epsilon$ with a weakly consistent Lagrange multiplier-based formulation. Both results show the same good quality of approximation, but the latter approach has a much better condition number and thus is more suitable for fast iterative solvers such as multigrid or domain decomposition methods.

Nowadays penalty techniques and simple node-to-node coupling strategies are increasingly replaced by variationally consistent methods which pass suitable patch tests in the case of non-matching meshes. The admissibility of the discrete solution is then formulated in a weak variational framework. Displacement and surface traction form a primal-dual pair of unknown variables and have to be discretized.

Figure 3.3 illustrates the difference between a simple node-to-node coupling strategy and a variationally consistent approach. A constant force can only be mapped correctly from the slave to the master side if a weak coupling is applied. The simple node-to-node coupling yields poor numerical results if non-matching meshes are used, whereas a sliding of the mesh does not influence the approximation quality in the case of a variationally consistent scheme.

Figure 3.4 shows the quantitative error decay for a node-to-node coupling in the case of non-matching meshes. As a patch test, we select a linear solution which can be represented exactly by standard low-order finite elements. In Figure 3.4(a), the scalar-valued Laplace operator is considered, whereas in Figure 3.4(b) the results for the vector-valued system of linear elasticity are presented; see Dohrmann, Key and Heinstein (2000) for the parameter specifications. In both cases, the exact solution cannot be reproduced, and the error decay is sub-optimal. A first-order error decay can be only observed for the $L^2$-norm but not for the $H^1$-norm. However, in the
case of a variationally consistent weak formulation, the exact solution can
be recovered and the error is equal to zero on all meshes. From the theo-
retical point of view, a node-to-node coupling is associated with a discrete
Lagrange multiplier being represented by a linear combination of delta dis-
tributions, which is not compatible with the required $H^{1/2}$-duality pairing.
In the following, we restrict ourselves to discrete Lagrange multiplier spaces
being $L^2$-conforming.

To obtain a stable and well-posed discrete setting, a uniform inf-sup condi-
tion has to be satisfied. Roughly speaking this means that the trace space of
the discrete displacement has to be well balanced with the finite-dimensional
space for the surface traction, also called the Lagrange multiplier. A neces-
sary condition is that the dimension of the Lagrange multiplier space is less
than or equal to the dimension of the jump of the trace spaces. There exists
a large variety of different construction principles, all leading to optimal
*a priori* estimates in the case of standard variational equalities. Quite often
such a condition is numerically verified by the Bathe-Chapelle inf-sup test
(Chapelle and Bathe 1993). A mathematically rigorous analysis can be per-
formed within the abstract framework of mortar settings on non-matching
meshes; see, e.g., Ben Belgacem and Maday (1997) and Bernardi, Maday
and Patera (1993, 1994). Early theoretical results on uniform stable dis-
cretization schemes for contact problems without friction can be found in
Ben Belgacem (2000) and Ben Belgacem, Hild and Laborde (1997, 1999).
When a vector-valued Lagrange multiplier is used, there is no algebraic dif-
ference between a contact problem with Coulomb friction and one without.
Thus, quite often solvers and error estimators designed for contact problems
without friction naturally apply to contact problems with Coulomb friction.
However, we recall that from the theoretical point of view there is a possi-
bly considerable difference, as for existence and uniqueness results; see, e.g.,
Eck *et al.* (2005) and Kikuchi and Oden (1988).

In this section, we illustrate the fact that a weakly consistent discretiza-
tion based on a biorthogonal set of displacement traces and surface tractions
is well suited to the numerical simulation of contact problems. While simple
node-to-node coupling strategies are known to show locking effects as well
as unphysical oscillations when applied to non-matching meshes, variation-
ally consistent formulations based on uniform inf-sup stable pairings do not
exhibit this behaviour. Moreover, the biorthogonality of the basis functions
of such a pairing yields a stable node-to-segment coupling concept where
the simple interpolation is replaced by a quasi-projection. This is quite
attractive because of the locality of the coupling constraints.

The discretization of the system is based on the saddle-point formulation
(2.15). Both cases, the frictionless case and that with Coulomb friction,
can be handled within the same abstract framework. In the case $\nu = 0$,
we do not work with the primal variational inequality (2.8) but also use
the primal-dual variational inequality setting. A low-order pair of primal-dual variables for the displacement **u** and the surface traction **λ** on the contact zone will be applied. As usual in the mortar context, the Lagrange multiplier space is associated with the (d-1)-dimensional surface mesh on $\Gamma_C^s$ inherited from the volume mesh on the slave side. In addition, the degrees of freedom from the master side will not be required for the inf-sup condition to hold. Thus the inf-sup constant is independent of the ratio between the mesh sizes of the master and slave sides and also independent of the non-matching character of the meshes, which is quite attractive in dynamic situations when sliding geometries occur.

In the linear saddle-point theory, it is well established (Brezzi and Fortin 1991, Nicolaides 1982) that a priori estimates in terms of the best approximation error of the primal and dual variables can be obtained if stability and continuity of the relevant bilinear forms are given. The norm for the displacement is the product $H^1$-norm, and for the surface traction, the $H^{-1/2}$-norm defined as the dual norm of the $H^{1/2}$-norm on $\Gamma_C^s$. Thus, to obtain first-order estimates for the best approximation error, the natural choice for the displacement is the lowest-order conforming finite element space on each of the two subdomains, whereas for the Lagrange multiplier several interesting choices exist. Basically all existing possibilities from the mortar literature can be used, e.g., piecewise constants associated with the dual mesh or low-order conforming finite elements.

### 3.1. A pairing, not uniformly stable

Before going into the details of the discretization, we consider a counterexample. We note that element-wise constants for the Lagrange multiplier do not yield optimal estimates. Although the best approximation error of element-wise constants with respect to the $L^2$-norm is of order one, this combination of primal and dual variables is not uniformly inf-sup stable. A mesh-dependent inf-sup constant results in a reduced convergence rate. Figure 3.5 illustrates this non-uniformly stable pairing in a one-dimensional setting.

Let the unit interval $I := (0,1)$ be decomposed into $N_l := 2^l, l \in \mathbb{N}$, sub-intervals $I_i := (i-1,i)/N_l, i = 1, \dots, N_l$, of equal length and

$$W_l := \{v \in C(I); v|_{I_i} \in P_1(I_i), i = 1, \dots, N_l\},$$

$$M_l := \{v \in L^2(I); v|_{I_i} \in P_0(I_i), i = 1, \dots, N_l\}.$$

Then each $\mu_l \in M_l$ can be written as $\mu_l = \sum_{i=1}^{N_l} \alpha_i \psi_i$, with $\alpha_i \in \mathbb{R}$, and where $\psi_i$ stands for the characteristic function of the sub-interval $I_i$. Each $v_l \in W_l$ has the form $v_l = \sum_{i=0}^{N_l} b_i \phi_i$, with $b_i \in \mathbb{R}$, and where $\phi_i$ denotes the standard hat function associated with the node $x_i := i/N_l, i = 0, \dots, N_l$.
<figure><img src="image_16.png" /><figcaption>Figure 3.5. Nodal finite element $\phi_i$ and Lagrange multiplier basis function $\psi_i$ (a), special choice of $\mu_l$ and $v_l$ for $l = 4$ and $l = 5$ (b,c).</figcaption></figure>

**Lemma 3.1.** The pairing $(W_l, M_l)$ is not uniformly inf-sup stable with respect to the $H^{\frac{1}{2}}$-duality. Moreover there exists a $c_{\inf} > 0$ independent of the number of sub-intervals $N_l$ such that

$$ \inf_{\mu_l \in M_l} \sup_{v_l \in W_l} \frac{\int_0^1 \mu_l v_l \, ds}{\|\mu_l\|_{-\frac{1}{2};I} \|v_l\|_{\frac{1}{2};I}} \geq \frac{c_{\inf}}{N_l}, \quad (3.1) $$

and the estimate in (3.1) is sharp.

*Proof.* To see that the non-uniform inf-sup condition given by (3.1) holds, it is sufficient to define a Fortin operator $F_l : H^{\frac{1}{2}}(I) \to W_l$ such that $\int_I \mu_l F_l v \, ds = \int_I \mu_l v \, ds$ for all $\mu_l \in M_l$ and $\|F_l v\|_{\frac{1}{2};I} \leq c N_l \|v\|_{\frac{1}{2};I}$ for all $v \in H^{\frac{1}{2}}(I)$; see, e.g., Brezzi and Fortin (1991). Let $v \in H^{1/2}(I)$ be given. In a first step, we set $w_1 := \sum_{i=0}^{N_l} \gamma_i \phi_i$ with

$$ \gamma_i := \frac{2 \int_{I_i \cup I_{i+1}} v(2\phi_i - \phi_{i-1} - \phi_{i+1}) \, ds}{|I_i \cup I_{i+1}|}, \quad i = 0, \dots, N_l, \quad I_0 := I_{N_l+1} := \emptyset. $$

In a second step, we define recursively $\Delta\gamma_i := 2N_l \int_{I_i} (v - w_1) \, ds - \Delta\gamma_{i-1}$, $i = 1, \dots, N_l$ and $\Delta\gamma_0 := 0$, $w_2 := \sum_{i=0}^{N_l} \Delta\gamma_i \phi_i$. We note from the definition of $w_1$ that it is $H^{\frac{1}{2}}$-stable and has $L^2$-approximation properties. In terms of $w_1$ and $w_2$, we set $F_l v := w_1 + w_2$. Then by construction, $F_l v$ and $v$ have the same mean value on each sub-interval $I_i$, $i = 1, \dots, N_l$. Introducing the matrix $B_l \in \mathbb{R}^{N_l \times N_l}$ as follows, we get for its inverse

$$ B_l := \begin{pmatrix} 1 \\ 1 & 1 \\ 1 & 1 \\ 1 & 1 \\ 1 & 1 \\ \vdots & \vdots \end{pmatrix}, \quad B_l^{-1} = \begin{pmatrix} 1 \\ -1 & 1 \\ 1 & -1 & 1 \\ -1 & 1 & -1 & 1 \\ 1 & -1 & 1 & -1 & 1 \\ \vdots & \vdots & \vdots & \vdots & \vdots \end{pmatrix}. $$
Moreover, a standard inverse estimate for finite elements and the fact that the Euclidean norm of $B_l^{-1}$ is bounded in terms of $N_l$, i.e., $\|B_l^{-1}\| \le CN_l$, gives

$$ \begin{aligned} \|F_l v\|_{\frac{1}{2};I} &\le \|w_1\|_{\frac{1}{2};I} + \|w_2\|_{\frac{1}{2};I} \le C(\|v\|_{\frac{1}{2};I} + \sqrt{N_l}\|w_2\|_{0;I}) \\ &\le C(\|v\|_{\frac{1}{2};I} + \sqrt{N_l}\|B_l^{-1}\|\|\Pi_l(v - w_1)\|_{0;I}) \\ &\le C(\|v\|_{\frac{1}{2};I} + \sqrt{N_l}N_l\|v - w_1\|_{0;I}) \le CN_l\|v\|_{\frac{1}{2};I}, \end{aligned} $$

where $\Pi_l$ stands for the $L^2$-projection onto $M_l$.

To show that the estimate (3.1) is sharp, we have to specify a $\mu_l \in M_l$ such that no better bound can be obtained. Let us consider the choice $\alpha_i := (-1)^i(i-1)(N_l - i)$: see Figure 3.5(b,c). The definition of the dual norm (2.1) yields $\|\mu_l\|_{-\frac{1}{2};I} \ge c \sum_{i=1}^{N_l} |\alpha_i|/(N_l\sqrt{N_l}) \ge cN_l\sqrt{N_l}$. Then, for $l \ge 2$ a straightforward computation and a standard inverse estimate shows

$$ \begin{aligned} \int_I \mu_l v_l \, ds &= \sum_{i=1}^{N_l} (-1)^i (i-1)(N_l - i)(b_{i-1} + b_i)/(2N_l) \\ &= \sum_{i=1}^{N_l-1} (-1)^i (2i - N_l)b_i/(2N_l) \\ &= \frac{1}{2N_l} \sum_{i=1}^{N_l/4} (N_l + 1 - 4i)((b_{2i-1} - b_{2i}) + (b_{N_l-2i} - b_{N_l+1-2i})) \\ &\quad + \frac{1}{2N_l} \sum_{i=1}^{N_l/2-1} (b_i - b_{N_l/2+i}) \\ &\le C(|v_l|_{1;I} + \|v_l\|_{0;I}) \le C\sqrt{N_l}\|v_l\|_{\frac{1}{2};I} \le \frac{C}{N_l}\|v_l\|_{\frac{1}{2};I}\|\mu_l\|_{-\frac{1}{2};I}. \quad \square \end{aligned} $$

**Remark 3.2.** A possible remedy would be to use a coarser mesh for the Lagrange multiplier space. The pairing $(W_l, M_{l-1})$ is uniformly inf-sup stable. Alternatively, the space $W_l$ can be enriched by locally supported bubble functions, as is done in Brezzi and Marini (2001) and Hauret and Le Tallec (2007). Here we do not follow these possibilities, but use only Lagrange multiplier spaces being defined on the same mesh as the trace space of the slave side and having the same nodal degrees of freedom.

## 3.2. Stable low-order discretization

On each subdomain $\Omega^k$, $k \in \{\text{m}, \text{s}\}$, independent families of shape-regular triangulations $\mathcal{T}_l^k$, $l \in \mathbb{N}_0$, will be used, and we set $\mathcal{T}_l := \mathcal{T}_l^m \cup \mathcal{T}_l^s$ and $\bar{\Omega}_l^k := \cup_{T \in \mathcal{T}_l^k} \bar{T}$. The maximum element diameter of the triangulation $\mathcal{T}_l$ is
<figure><img src="image_17.png" /><figcaption>Figure 3.6. Discontinuous biorthogonal basis functions satisfying (3.2)–(3.5).</figcaption></figure>

denoted by $h_l$. The restriction of $\mathcal{T}_l^k$ to $\Gamma_{C;l}^k := \cup_{F \in \mathcal{F}_l^k} \bar{F}$, where $\mathcal{F}_l^k$ stands for the set of all contact faces on the side $k$, defines a $(d-1)$-dimensional surface mesh. Moreover, the surface mesh on $\Gamma_{C;l}^m$ will be mapped by $\chi_l^{-1}$ onto $\Gamma_{C;l}^s$, resulting in possibly non-matching meshes on the contact zone.

For the displacement, we use standard low-order conforming finite elements and for the surface traction dual finite elements which reproduce constants

$$ 
\begin{aligned}
V_l &:= V_l^m \times V_l^s, & V_l^k &:= (V_l^k)^d, & V_l^k &:= \operatorname{span}_{p \in \mathcal{P}_l^k} \{\phi_p\}, \\
M_l &:= M_l^s, & M_l^k &:= (M_l^k)^d, & M_l^k &:= \operatorname{span}_{p \in \mathcal{P}_{C;l}^k} \{\psi_p\},
\end{aligned}
$$

where $\mathcal{P}_l^k$ stands for all vertices of $\mathcal{T}_l^k$ not being on $\bar{\Gamma}_D^k$, and $\mathcal{P}_{C;l}^k$ is the set of all vertices on $\bar{\Gamma}_C^k$, $k \in \{\text{m}, \text{s}\}$. Moreover $\phi_p$ denotes the standard conforming nodal basis function associated with the vertex $p$. The basis functions $\psi_p \in M_l^k$ are required to have the following properties.

* Locality of the support:

$$ \operatorname{supp} \psi_p = \operatorname{supp} \phi_p|_{\Gamma_{C;l}^k}, \quad p \in \mathcal{P}_{C;l}^k, \qquad (3.2) $$

* Local biorthogonality relation:

$$ \int_F \psi_p \phi_q \, ds = \delta_{pq} \int_F \phi_q \, ds \ge 0, \quad p, q \in \mathcal{P}_{C;l}^k, F \in \mathcal{F}_l^k, \qquad (3.3) $$

* Best approximation property:

$$ \inf_{\mu_l \in M_l^k} \|\mu - \mu_l\|_{-\frac{1}{2}; \Gamma_{C;l}^k} \le C_{\text{best}} h_l |\mu|_{\frac{1}{2}; \Gamma_{C;l}^k}, \quad \mu \in H^{\frac{1}{2}}(\Gamma_{C;l}^k), \qquad (3.4) $$

where $C_{\text{best}} < \infty$ does not depend on the mesh size.

* Uniform inf-sup condition:

$$ \sup_{v_l \in W_l^k} \frac{\int_{\Gamma_{C;l}^k} \mu_l v_l \, ds}{\|v_l\|_{\frac{1}{2}; \Gamma_{C;l}^k}} \ge c_{\inf} \|\mu_l\|_{-\frac{1}{2}; \Gamma_{C;l}^k}, \quad \mu_l \in M_l^k, \qquad (3.5) $$
<figure><img src="image_18.png" /><figcaption>Figure 3.7. Continuous biorthogonal basis functions satisfying (3.2)–(3.5).</figcaption></figure>

where $c_{inf} > 0$ does not depend on the mesh size and the discrete trace space $W_l^k$ is given by $W_l^k := \text{span}_{p \in \mathcal{P}_{C;l}^k} \{\phi_p |_{\Gamma_{C;l}^k}\}$.

From the local relation (3.3) follows directly a global one, i.e.,

$$ \int_{\Gamma_{C;l}^k} \psi_p \phi_q \, ds = \delta_{pq} m_p := \delta_{pq} \int_{\Gamma_{C;l}^k} \phi_q \, ds \geq 0, \quad p, q \in \mathcal{P}_{C;l}^k, \quad (3.6) $$

with $m_p > 0$. We note that there exists no set of non-negative basis functions satisfying (3.6). Moreover (3.3) in combination with (3.4) automatically yields $\sum_{p \in \mathcal{P}_{C;l}^k} \psi_p = 1$ and

$$ \int_F \psi_p \, ds = \int_F \phi_p \, ds, \quad p \in \mathcal{P}_{C;l}^k, F \in \mathcal{F}_l^k. \quad (3.7) $$

The most popular choice of a dual Lagrange multiplier (Wohlmuth 2000, 2001) is obtained by an element-wise biorthogonalization process of the local nodal finite elements followed by a node-wise glueing step to reduce the number of degrees of freedom. As a result piecewise-linear but discontinuous basis functions are created for a one-dimensional contact zone. This technique works for the lowest-order finite elements on all types of surface meshes. In the case of higher-order elements, we have to use Gauss-Lobatto nodes. This approach is restricted to tensorial meshes where the element mapping is affine: see Lamichhane and Wohlmuth (2007). Alternatively, for low-order elements, as in the present discussion, we can apply piecewise-constant or -quadratic basis functions on a sub-mesh: see Figure 3.6(b,c), respectively.

In some applications it is of interest to work with continuous basis functions. Although the local construction of biorthogonal basis functions, node-wise defined, quite often results in a set of discontinuous basis functions, continuous ones do exist. Figure 3.7 shows basis functions which satisfy (3.2)–(3.5) and which in addition are continuous. In Figure 3.7(a), the basis function is cubic on each element. For the construction of the basis function shown in Figure 3.7(b,c), we use sub-elements and a piecewise quadratic and
linear approach, respectively. As will be seen, the cubic dual Lagrange multiplier goes hand in hand with the definition of $H(\text{div})$-conforming mixed finite elements for linear elasticity and fits well into the construction of a *posteriori* error estimators based on element-wise lifting techniques: see Section 6.

The following two remarks briefly comment on more general formulations. In particular, the construction of a dual Lagrange multiplier basis possibly depends on the geometry in the case of cylinder coordinates or of sliding meshes.

**Remark 3.3.** If a three-dimensional situation is reduced to a two-dimensional setting by introducing cylinder coordinates and exploiting symmetry arguments, then duality has to be formulated with respect to a weighted scalar product. The distance to the symmetry axis enters as weight in the local biorthogonality relation (3.3). As example, we consider the piecewise affine case. Let $r_1$ and $r_2$ be the distances of the two face nodes to the symmetry axis. Then the dual Lagrange multiplier restricted to a face can be written as

$$ \psi_1 = \frac{2r_1 + r_2}{r_1^2 + 4r_1r_2 + r_2^2} ((r_1 + 3r_2)\phi_1 - (r_1 + r_2)\phi_2), $$

$$ \psi_2 = \frac{r_1 + 2r_2}{r_1^2 + 4r_1r_2 + r_2^2} ((3r_1 + r_2)\phi_2 - (r_1 + r_2)\phi_1), $$

where $\phi_1$ and $\phi_2$ are the two nodal basis function associated with the face. We note that $\psi_1 + \psi_2 = 1$, and for $r_1 = r_2$ we fall back to the piecewise affine dual Lagrange multiplier depicted in Figure 3.6(a).

**Remark 3.4.** A typical benchmark problem for large deformation contact is a small cube sliding over a larger block. If the large block is defined as slave side, then the integral over the face $F \in \mathcal{F}_l^s$ in the local biorthogonality relation (3.3) has to be replaced by the integral over $F \cap \Gamma_{C;l}^m$. Let us consider the 2D reference case $F = (0,1)$ and $(s,1) \subset \Gamma_{C;l}^m$ with $s \in (0,1)$. Then a straightforward calculation shows

$$ \psi_1 = \frac{2(1 + s + s^2)}{1 - s} \phi_1 - (1 + 2s)\phi_2, $$

$$ \psi_2 = -\frac{(1 + 2s)(1 + s)}{1 - s} \phi_1 + 2(1 + s)\phi_2, $$

where $\phi_1$ is the nodal basis function associated with the endpoint $p = 0$ and $\phi_2$ is associated with $p = 1$.

The weak problem formulation will be based on suitable subsets of $V_l$ and $M_l$. For the displacement, we only impose the Dirichlet condition on the space and use no constraint related to the non-penetration condition.
The convex set $V_{l;D}$ is given by

$$V_{l;D} := V_{l;D}^{m} \times V_{l;D}^{s}, \quad V_{l;D}^{k} := \sum_{p \in \mathcal{P}_{D;l}^{k}} u_{D}(p)\phi_{p} + V_{l}^{k},$$

where $\mathcal{P}_{D;l}^{k}$ stands for the set of all vertices of the actual mesh $\mathcal{T}_{l}^{k}$ on $\bar{\Gamma}_{D}^{k}$. To handle the contact conditions (2.5) and (2.6) in a weakly consistent form, we have to impose constraints on the Lagrange multiplier space. The non-penetration condition restricts the normal part of $M_{l}$, and the Coulomb law requires a solution-dependent inequality bound for the tangential part.

Let the discrete solution $\lambda_{l} \in M_{l}$ be given by

$$\lambda_{l} = \sum_{p \in \mathcal{P}_{C;l}^{s}} \gamma_{p} \psi_{p}, \quad \gamma_{p} \in \mathbb{R}^{d},$$

and denote the discrete normal component by

$$\lambda_{l}^{n} := \sum_{p \in \mathcal{P}_{C;l}^{s}} \gamma_{p}^{n} \psi_{p}, \quad \gamma_{p}^{n} := \gamma_{p} n_{p}^{s}.$$

Here $n_{p}^{s}$ stands for a discrete normal vector associated with the node $p$. In the case of a non-planar contact surface, it can be obtained as a weighted combination of the adjacent element centre normals. We observe that $\lambda_{l}^{n} \in M_{l}^{s}$ but, in general, it is not equal to $\lambda_{l}n^{s}$. Based on the discrete normal surface traction $\lambda_{l}^{n}$, we then define the convex set

$$M_{l}(\lambda_{l}^{n}) := \left\{ \mu_{l} = \sum_{p \in \mathcal{P}_{C;l}^{s}} \beta_{p} \psi_{p}, \beta_{p} \in \mathbb{R}^{d}, \beta_{p}^{n} \geq 0, \|\beta_{p}^{t}\| \leq \nu \gamma_{p}^{n} \right\}, \quad (3.8)$$

as an approximation for the solution-dependent cone $M(\lambda^{n})$ defined in (2.14). Here $\beta_{p}^{n} := \beta_{p}n_{p}^{s}$, $\beta_{p}^{t} := \beta_{p} - \beta_{p}^{n}n_{p}^{s}$. We assume that $\Gamma_{C;l}^{s}$ is large enough such that $\gamma_{p} = 0$ for $p \in \partial\Gamma_{C;l}^{s}$.

**Remark 3.5.** If $\nu = 0$, the convex set $M_{l}(\lambda_{l}^{n})$ is solution-independent, and its definition reduces to

$$M_{l}^{+} := \left\{ \mu_{l} = \sum_{p \in \mathcal{P}_{C;l}^{s}} \beta_{p} \psi_{p}, \beta_{p} \in \mathbb{R}^{d}, \beta_{p}^{n} \geq 0, \beta_{p}^{t} = 0 \right\}.$$

In this case it is sufficient to work with a scalar-valued Lagrange multiplier space as is often done in the literature. Here, we use a vector-valued Lagrange multiplier to be in the same abstract framework for all $\nu \geq 0$.

If the standard nodal Lagrange multiplier basis is used to define $M_{l}$, there are then two natural but different ways to discretize $M^{+}$; see also, e.g., Hild
<figure><img src="image_19.png" /><figcaption>Figure 3.8. Two elements of $M_{l;2}^+ \setminus M_{l;1}^+$: local support (a) and global support (b).</figcaption></figure>

and Renard (2010). The first yields $M_{l;1}^+ \subset M^+$, whereas the second choice is based on the definition of $M^+$ as a dual cone:

$$M_{l;1}^+ := \left\{ \mu_l \in M_l \mid \mu_l = \sum_{p \in \mathcal{P}_{C;l}^s} \beta_p \psi_p, \beta_p \ge 0, p \in \mathcal{P}_{C;l}^s \right\},$$

$$M_{l;2}^+ := \left\{ \mu_l \in M_l \mid \int_{\Gamma_{C;l}^s} \mu_l \phi_p \, ds \ge 0, p \in \mathcal{P}_{C;l}^s \right\}.$$

We note that these two definitions yield two different spaces, with $M_{l;1}^+$ being a proper subspace of $M_{l;2}^+$. Figure 3.8 shows elements of $M_{l;2}^+$ which are clearly not in $M_{l;1}^+$. In Figure 3.8(a), the element $\mu_l \in M_l$ depends on the parameter value $a$ and, as a straightforward computation shows, $\mu_l \in M_{l;2}^+ \setminus M_{l;1}^+$ if and only if $a \in [-0.5; 0)$. The function in Figure 3.8(b) is obviously not in $M_{l;1}^+$. Testing it with all nodal basis functions $\phi_p$, $p \in \mathcal{P}_{C;l}^s$, yields the non-negative values (marked with bullets) and thus $\mu_l \in M_{l;2}^+$.

Both choices of $M_{l;i}^+$, $i = 1, 2$, can be applied in the discrete setting. The first one yields a conforming approach, whereas in the second one the non-conformity of $M_l^+$ has to be taken into account in the *a priori* estimates. The difference in the spaces stems from the fact that the matrix given by $\int_{\Gamma_{C;l}^s} \psi_q \phi_p \, ds$ is then the standard mass matrix, which is not an M-matrix. The situation is different if our Lagrange multiplier basis, satisfying (3.6), is applied. Then the mass matrix is diagonal and positive definite and both definitions yield the same space $M_l^+$. However, using a biorthogonal basis automatically results in a non-conforming approach, i.e., $M_l^+ \not\subset M^+$. To enforce conformity in that situation is not a good idea, since then the locality of the elements in $M_l^+ \cap M^+$ is lost; see also Figure 3.9.

Having the conforming finite element space $V_l \subset V_0$ and the non-conforming closed convex cone $M_l(\lambda_l^n)$, we can formulate the discrete weak version of (2.15) as follows. Find $\mathbf{u}_l := (\mathbf{u}_l^m, \mathbf{u}_l^s) \in V_{l;D}$, $\lambda_l \in M_l(\lambda_l^n)$ such that

$$\begin{aligned} a_l(\mathbf{u}_l, \mathbf{v}_l) + b_l(\boldsymbol{\lambda}_l, \mathbf{v}_l) &= f_l(\mathbf{v}_l), & \mathbf{v}_l &\in \mathbf{V}_l, \\ b_l(\boldsymbol{\mu}_l - \boldsymbol{\lambda}_l, \mathbf{u}_l) &\le g_l(\boldsymbol{\mu}_l - \boldsymbol{\lambda}_l), & \boldsymbol{\mu}_l &\in \mathbf{M}_l(\lambda_l^n). \end{aligned} \quad (3.9)$$
<figure><img src="image_20.png" /><figcaption>Figure 3.9. Elements in $M_l^+ \cap M^+$ for different refinement levels.</figcaption></figure>

The mesh-dependent bilinear and linear forms are obtained from the associated continuous ones in a natural way by replacing the volume and surface terms by the corresponding discrete analogue, i.e., summing over the volume and surface elements of the mesh.

### 3.3. Coupling in terms of the mortar projection

Although the given variational setting is a two-body formulation with possibly non-matching meshes, we can reformulate the contact conditions in a way similar to a one-body system. To do so, we introduce the mortar projection $\Pi_l := (\Pi_l)^d$ onto $W_l^s$ and the dual mortar projection $\Pi_l^* := (\Pi_l^*)^d$ onto $M_l^s$; see, e.g., Bernardi, Maday and Patera (1993, 1994). We recall that due to the assumption on the Dirichlet boundary part, no modification at the endpoints is required. For $w \in H^{1/2}(\Gamma_{C;l}^s)$ and $\mu \in H^{-1/2}(\Gamma_{C;l}^s)$ we set

$$ \int_{\Gamma_{C;l}^s} \Pi_l w \mu_l \, ds := \langle \mu_l, w \rangle_{\Gamma_{C;l}^s}, \quad \mu_l \in M_l^s, \qquad (3.10a) $$

$$ \int_{\Gamma_{C;l}^s} \Pi_l^* \mu w_l \, ds := \langle \mu, w_l \rangle_{\Gamma_{C;l}^s}, \quad w_l \in W_l^s. \qquad (3.10b) $$

We note that $\Pi_l$ and $\Pi_l^*$ restricted to $W_l^s$ and $M_l^s$ is the identity, respectively. In terms of $\Pi_l$, we can write $\Pi_l[u_l] = \sum_{p \in \mathcal{P}_{C;l}^s} \alpha_p \phi_p$. Moreover, we define $g_p := \int_{\Gamma_{C;l}^s} g_l \psi_p \, ds/m_p$ with a suitable approximation $g_l$ for the linearized gap, where $m_p$ is specified in (3.6). The following lemma shows that Lemma 2.3 has a node-wise discrete analogue.

**Lemma 3.6.** For each node $p \in \mathcal{P}_{C;l}^s$ the following discrete node-wise KKT conditions hold for non-penetration:

$$ 0 \le \gamma_p^n, \quad \alpha_p^n \le g_p, \quad \gamma_p^n(\alpha_p^n - g_p) = 0. \qquad (3.11) $$

Moreover, a discrete static Coulomb law holds for each node:

$$ \|\gamma_p^t\| \le \nu\gamma_p^n, \quad \alpha_p^t \gamma_p^t - \nu\gamma_p^n \|\alpha_p^t\| = 0. \qquad (3.12) $$
<figure><img src="image_21.png" /><figcaption>Figure 3.10. Discrete approximations $\Omega_l$ of the domain $\Omega$.</figcaption></figure>

*Proof.* We observe that the constraints on $\gamma_p^n$ in (3.11) and on $\gamma_p^t$ in (3.12) follow directly from the definition of $M_l(\lambda_l^n)$. The biorthogonality (3.6) plays an essential role in the proof. Using $\boldsymbol{\mu}_l = \boldsymbol{\lambda}_l \pm \gamma_p^n \mathbf{n}_p^s \psi_p \in M_l(\lambda_l^n)$ as test function in (3.9), we find $\pm \gamma_p^n b(\mathbf{n}_p^s \psi_p, \mathbf{u}_l) = \pm \gamma_p^n \alpha_p^n m_p \le g_p \gamma_p^n m_p$ and thus the complementarity condition in (3.11). (3.12) obviously holds for $\boldsymbol{\alpha}_p^t = \mathbf{0}$. For $\boldsymbol{\alpha}_p^t \ne \mathbf{0}$, we set $\boldsymbol{\mu}_l = \boldsymbol{\lambda}_l - \boldsymbol{\gamma}_p^t \psi_p + \nu \gamma_p^n \mathbf{e}_t \psi_p \in M_l(\lambda_l^n)$ with $\mathbf{e}_t := \boldsymbol{\alpha}_p^t / \|\boldsymbol{\alpha}_p^t\|$ as the test function in (3.9) and get (3.12). $\square$

The discrete contact conditions (3.11) and (3.12) only involve quantities associated with the slave nodes and thus have the same structure as a one-body system. However, to compute $\boldsymbol{\alpha}_p$, we do have to evaluate the mortar projection $\boldsymbol{\Pi}_l$ applied on $[\mathbf{u}_l]$. Its algebraic representation can be obtained from the entries of the mass matrix associated with $b_l(\cdot, \cdot)$. For its implementation, we not only have to map the mesh elements on the possible contact zone but also the basis functions.

In Figure 3.10, we show different possible matching and non-matching situations in the case of a non-planar contact surface. We note that in this situation standard triangulations do not resolve the domain exactly. In contrast to the continuous setting where $\Gamma_C^s = \Gamma_C^m$, we find, in the discrete setting, that the possible contact zones on the master and the slave sides are not the same, i.e., $\Gamma_{C;l}^s \ne \Gamma_{C;l}^m$.

Using hierarchical tree structures or front tracking techniques, the assembly of the surface-based coupling matrices between master and slave side can be realized quite efficiently and is of lower complexity, whereas in 3D a naive approach results in a higher complexity compared to the assembling process of the volume contributions. For the integration in 3D we use quadrature formulas on surface sub-triangles. Figure 3.11 illustrates different steps of the projection and partitioning procedure. This algorithm goes back to Puso (2004) (see also Puso et al. (2008)) and, alternatively, the recent papers by Dickopf and Krause (2009a, 2009b).

An analysis and a numerical study of the influence of curvilinear interfaces in the mortar situation can be found in Flemisch, Melenk and Wohlmuth (2005a). The abstract setting of blending elements (Gordon and Hall 1973a, 1973b) plays a key role in establishing optimal upper bounds
<figure><img src="image_22.png" /><figcaption>Figure 3.11. Element-wise mapping from the master to the slave side.</figcaption></figure>

of the consistency error. From the theoretical point of view, the mapping between $\Gamma_{C;l}^s$ and $\Gamma_{C;l}^m$ has to be globally smooth. From the computational point of view, an element-wise smooth but possibly discontinuous mapping is more attractive and works well in practice. Replacing the $H^{1/2}$-norm on the contact interface by a weighted mesh-dependent $L^2$-norm, this simplification can also be theoretically justified. A 3D analysis, in the case of planar interfaces and mesh-dependent norms, can be found in Braess and Dahmen (1998).

We note that for non-planar interfaces, a slave side associated with the coarser mesh and a vectorial partial differential equation, e.g., linear elasticity, possibly poor numerical results can be observed in the pre-asymptotic range, if the Lagrange multiplier with respect to its Cartesian coordinates is discretized by dual or piecewise constant basis functions: see, e.g., Flemisch, Puso and Wohlmuth (2005b). This effect can be explained by the observation that a constant normal surface and a zero tangential force cannot very well be approximated in terms of the Lagrange multiplier space. Firstly decomposing the Lagrange multiplier into its normal and tangential parts, and secondly discretizing yield much better results. However, in that case we have to be quite careful to handle rigid body motions correctly. These oscillations do not occur for the Laplace operator, with a finer slave mesh side, standard linear Lagrange multipliers, or quadratic dual Lagrange multipliers. For contact problems, non-penetration and the friction law are directly expressed in terms of the normal and tangential component of the surface traction. Thus a discretization of the traction in its locally rotated coordinate system seems to be quite attractive, in particular for finite deformations, and has already been used to define $M_l(\lambda_l^n)$. For small deformations and a constant contact normal, both approaches give the same results.

## 4. Optimal *a priori* error estimates

Abstract error estimates for variational inequalities can be found, e.g., in Brezzi, Hager and Raviart (1977), Falk (1974), Glowinski (1984) and Glowinski et al. (1981) and a priori bounds for the discretization error of unilateral contact problems are given, e.g., in Haslinger and Hlaváček (1981) and Haslinger et al. (1996). It is well known that the finite element solution of a variational inequality may have a reduced convergence order, compared to that of the best approximation. This holds true for higher-order finite
elements but may also be true for low-order approaches. A proof purely based on standard techniques will generally yield only $O(\sqrt{h_l})$ a priori bounds. We refer to the monograph by Han and Reddy (1999) for an introduction into this area for applications in plasticity.

Mortar techniques for contact problems without friction have been introduced in Ben Belgacem, Hild and Laborde (1998), Ben Belgacem (2000) and Lhalouani and Sassi (1999). We also refer to Ben Belgacem and Renard (2003), Coorevits, Hild, Lhalouani and Sassi (2001), Hild (2000) and Hild and Laborde (2002), where standard Lagrange multiplier spaces have been considered and analysed. In early papers on mortar, unilateral contact problems have quite often been considered, taking no friction into account, and using a scalar-valued Lagrange multiplier. The choice of the contact pressure as Lagrange multiplier is motivated by the fact that in that case the tangential component of the surface traction is zero. There is a series of papers on a priori estimates for two-body contact problems with no friction on non-matching meshes starting with order 1/4 bounds for the discretization error (Ben Belgacem et al. 1998). A priori error estimates for the displacements in the $H^1$-norm and for the Lagrange multiplier in the $H^{-1/2}$-norm of order 3/4 have been established; see, e.g., Ben Belgacem, Hild and Laborde (1999), Ben Belgacem and Renard (2003), Coorevits et al. (2001) and Lhalouani and Sassi (1999), under an $H^2$-regularity assumption. Using additional quite strong and restrictive regularity assumptions on the Lagrange multiplier, order one has been shown; see, e.g., Coorevits et al. (2001) and Hild (2000). These first a priori results have been considerably improved over the last decade. Under suitable assumptions on the actual contact zone and a $H^2$-regular solution quasi-optimal, i.e., $h_l \sqrt{|\log h_l|}$ and $h_l \sqrt[4]{|\log h_l|}$ a priori estimates can be found in Ben Belgacem (2000) and Ben Belgacem and Renard (2003). Most of the theoretical results are obtained for standard Lagrange multipliers, no friction and in the two-dimensional setting. Here, we apply these techniques to vector-valued dual Lagrange multiplier spaces and provide a priori error estimates for the displacement in the $H^1$-norm and for the surface traction in the $H^{-1/2}$-norm. In 3D, only sub-optimal bounds can be obtained for a problem with non-trivial friction. In 2D, we follow the lines of Hüeber, Matei and Wohlmuth (2005b) and Hüeber and Wohlmuth (2005a) and establish in a simplified problem setting optimal a priori bounds under some regularity assumption on the actual contact part and on the sticky zone.

We assume that no variational crimes are committed, i.e., the discrete bilinear and linear forms are exact. In particular, this implies that no quadrature error occurs and that $\Omega_l^k = \Omega^k$; we refer to Ciarlet (1991, 1998) for a rigorous mathematical analysis of the influence of quadrature formulas. Moreover, we assume a zero gap, i.e.,

$$g = 0, \Gamma_C := \Gamma_C^s = \Gamma_{C;l}^s \subset \Gamma_C^m = \Gamma_{C;l}^m$$
<figure><img src="image_23.png" /><figcaption>Figure 4.1. Decomposition into 4 and 25 subdomains (a,b) and error decay in the $H^1$-norm (c).</figcaption></figure>

and a constant unit vector $\mathbf{n} := \mathbf{n}^s$ on the possible contact zone. The a priori analysis of a curvilinear interface for the classical linear- and scalar-valued mortar case can be found in Flemisch et al. (2005a). Figure 4.1 shows that the number of subdomains in the linear mortar setting does not influence the constants in the a priori bound. The unit square is decomposed into $l^2$, $l = 2, 3, 4, 5$ subdomains; the interface is given by a sinus wave function. Due to the curvilinear character of the interface, a mesh-dependent mapping between the discrete master and slave interface is required.

In the case of a non-linear contact problem, the same approach can be applied. For simplicity of notation, we do not provide any technical detail here and restrict ourselves to simple geometrical settings such as, e.g., a square on a rectangle. Moreover, we will work with globally quasi-uniform meshes, such that there exists a regularity constant $c_{\text{reg}} > 0$ so that, for all nodes $p$ on $\Gamma_C$, we have $B_p(c_{\text{reg}}h_l) \subset \text{supp } \phi_p|_{\Gamma_C}$, where $B_p(c_{\text{reg}}h_l)$ is the $(d-1)$-dimensional ball with centre $p$ and radius $c_{\text{reg}}h_l$.

Most importantly, we replace for the rest of this section the Coulomb friction law (2.6) by the more simple Tresca law with a constant friction bound $\mathcal{F}$, i.e.,

$$ \|\lambda_t\| \le \mathcal{F}, \quad [\dot{u}_t]\lambda_t - \mathcal{F}\|[\dot{u}_t]\| = 0. \quad (4.1) $$

In the following, we will frequently make use of the Tresca version of Lemmas 2.3 and 3.6. Following the lines of the proof, we find for the normal components the complementarity conditions

$$ b_n(\lambda, \mathbf{u}) = 0 = b_n(\lambda_l, \mathbf{u}_l), \quad (4.2) $$

in addition to $\lambda_n \in M^+$, $\lambda_l^n \in M_l^+$ and $-[u_n], -\Pi_l[u_l^n] \in W^+$. We note that due to the possibly non-matching meshes, $-[u_l^n]$ is, in general, not in $W^+$. The friction law (4.1) guarantees that the tangential components satisfy

$$ b_t(\lambda, \mathbf{u}) = \int_{\Gamma_C} \mathcal{F}\|[\mathbf{u}_t]\| \, \mathrm{d}s, \quad b_t(\lambda_l, \mathbf{u}_l) = \int_{\Gamma_C} \mathcal{F}\|\Pi_l[\mathbf{u}_l^t]\| \, \mathrm{d}s, \quad (4.3) $$
<figure><img src="image_24.png" /><figcaption>Figure 4.2. Absolute value in 2D $|| \cdot ||$ (a), discrete absolute value $|| \cdot ||_l$ (b) and the difference $|| \cdot || - || \cdot ||_l$ (c).</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Figure 4.3. Numerical results for a contact problem with Tresca friction: normal displacement and contact pressure (a) and surface traction in normal and tangential direction (b).</figcaption></figure>

where the mesh-dependent Euclidean norm $|| \cdot ||_l$ is defined by

$$ ||\mathbf{\Pi}_l[\mathbf{u}_l^t]||_l := \sum_{p \in \mathcal{P}_{C;l}^s} ||\boldsymbol{\alpha}_p^t|| \phi_p \in W_l^+. $$

Figure 4.2 illustrates the difference between $|| \cdot ||$ and $|| \cdot ||_l$. The dashed line shows an element in $W_l$ if $\Gamma_C$ is a straight line. It is easy to verify that $||\mathbf{\Pi}_l[\mathbf{u}_l^t]|| - ||\mathbf{\Pi}_l[\mathbf{u}_l^t]||_l$ is equal to zero at all vertices of the slave-side mesh and that, in general, we do not have $||\mathbf{\Pi}_l[\mathbf{u}_l]|| \in W_l$. Moreover, we find that $||\mathbf{\Pi}_l[\mathbf{u}_l^t]|| - ||\mathbf{\Pi}_l[\mathbf{u}_l^t]||_l \le 0$.

Figure 4.3 shows a simple numerical example for a contact problem with a non-constant Tresca friction bound and a non-zero gap. Further, we observe that for such examples the discrete complementarity conditions hold true. Figure 4.3(a) shows that the coefficients $\gamma_p^n$ are non-zero only on the actual discrete contact zone and thus (4.2) is satisfied. In Figure 4.3(b), the tangential displacement vanishes if the tangential stress component is strictly below its given bound, and thus (4.3) holds.

In addition to the best approximation properties of the discrete Lagrange multiplier space (3.4), our proof relies on the properties of $\Pi_l$ and $\Pi_l^*$: see
Wohlmuth (2001). For $1 \le s \le 2$, we have

$$ \|\mu - \Pi_l^* \mu\|_{-\frac{1}{2}; \Gamma_C} \le C h_l^{s-1} |\mu|_{s-\frac{3}{2}; \Gamma_C}, \quad \mu \in H^{s-\frac{3}{2}}(\Gamma_C), \quad (4.4a) $$

$$ \|w - \Pi_l w\|_{\frac{1}{2}; \Gamma_C} \le C h_l^{s-1} |w|_{s-\frac{1}{2}; \Gamma_C}, \quad w \in H^{s-\frac{1}{2}}(\Gamma_C). \quad (4.4b) $$

**Remark 4.1.** We note that our regularity assumptions on the data and geometry guarantee that for $u \in (H^s(\Omega))^d$, $1 \le s \le 2$, we automatically have $\lambda = -\sigma(u^s)n \in (H^{s-3/2}(\Gamma_C))^d$, and moreover $|\lambda|_{s-3/2; \Gamma_C} \le C|u|_{s; \Omega}$. For $s = 2$ this is obvious; for $s = 1$ this does not hold for all $f \in V'_0$ and $f_N \in \{w = v|_{\Gamma_N}, v \in V_0\}$. However, by increasing the regularity of the data, as we do here, this holds true with a constant depending on $C_{\text{reg}}$ (see (2.4)), and then a standard interpolation argument yields the result for all $s \in [1, 2]$.

### 4.1. Upper bound for the discretization error

The starting point is the following abstract lemma. A similar lemma can be found in Hild and Laborde (2002) without friction and with quadratic finite elements associated with standard Lagrange multipliers. Here, we also have to consider the friction part and take into account the tangential component of the bilinear form $b(\cdot, \cdot)$. Introducing the error $E_l := (u - u_l, \lambda - \lambda_l)$ and its associated norm $\|E_l\|_{V \times M}^2 := \|u - u_l\|_{1; \Omega}^2 + \|\lambda - \lambda_l\|_{-1/2; \Gamma_C}^2$, the standard saddle-point theory and the complementarity conditions (4.2) provide a first a priori result.

**Lemma 4.2.** Let $(u, \lambda) \in V \times M(\mathcal{F})$ be the solution of (2.15) with $M(\lambda_n)$ replaced by $M(\mathcal{F})$ and let $(u_l, \lambda_l) \in V_l \times M_l(\mathcal{F})$ be the solution of the discrete formulation (3.9) with $M_l(\lambda_l^n)$ replaced by $M_l(\mathcal{F})$. Then, we have

$$ \|E_l\|_{V \times M} \le C \left\{ \begin{aligned} & \inf_{v_l \in V_l} \|u - v_l\|_{1; \Omega} + \inf_{\mu_l \in M_l} \|\lambda - \mu_l\|_{-\frac{1}{2}; \Gamma_C} \\ & + \max(b_n(\lambda_l, u), 0)^{\frac{1}{2}} + \max(b_n(\lambda, u_l), 0)^{\frac{1}{2}} \\ & + \max(b_t(\lambda_l - \lambda, u), 0)^{\frac{1}{2}} + \inf_{w_l \in W_l^s} \|[u^t] - w_l\|_{\frac{1}{2}; \Gamma_C} \end{aligned} \right\}. $$

*Proof.* Introducing $e_l := u - u_l$, we find for the error $e_l$ in the energy norm, and for $v_l \in V_l$

$$ \begin{aligned} a(\mathbf{e}_l, \mathbf{e}_l) &= a(\mathbf{e}_l, \mathbf{u} - \mathbf{v}_l) - b(\boldsymbol{\lambda}, \mathbf{v}_l - \mathbf{u}_l) + b(\boldsymbol{\lambda}_l, \mathbf{v}_l - \mathbf{u}_l) \\ &= a(\mathbf{e}_l, \mathbf{u} - \mathbf{v}_l) - b(\boldsymbol{\lambda} - \boldsymbol{\lambda}_l, \mathbf{v}_l - \mathbf{u}) - b(\boldsymbol{\lambda} - \boldsymbol{\lambda}_l, \mathbf{e}_l). \end{aligned} $$

Then Korn's inequality, which holds on both subdomains by assumption, and the continuity of the bilinear forms $a(\cdot, \cdot)$ and $b(\cdot, \cdot)$ yield an upper bound for the $H^1$-error of the displacement:

$$ \|\mathbf{e}_l\|_{1; \Omega}^2 \le C\left(\|\mathbf{e}_l\|_{1; \Omega} + \|\boldsymbol{\lambda} - \boldsymbol{\lambda}_l\|_{-\frac{1}{2}; \Gamma_C}\right) \|\mathbf{u} - \mathbf{v}_l\|_{1; \Omega} - b(\boldsymbol{\lambda} - \boldsymbol{\lambda}_l, \mathbf{e}_l). \quad (4.5) $$
Using standard techniques from the saddle-point framework and applying the discrete inf-sup condition (3.5), we get

$$
\begin{align*}
\|\boldsymbol{\mu}_l - \boldsymbol{\lambda}_l\|_{-\frac{1}{2};\Gamma_C} &\le C \sup_{\mathbf{w}_l \in \mathbf{V}_l} \frac{b(\boldsymbol{\mu}_l - \boldsymbol{\lambda}_l, \mathbf{w}_l)}{\|\mathbf{w}_l\|_{1;\Omega}} \\
&= C \sup_{\mathbf{w}_l \in \mathbf{V}_l} \frac{b(\boldsymbol{\mu}_l - \boldsymbol{\lambda}, \mathbf{w}_l) + a(\mathbf{u}_l - \mathbf{u}, \mathbf{w}_l)}{\|\mathbf{w}_l\|_{1;\Omega}} \\
&\le C (\|\boldsymbol{\mu}_l - \boldsymbol{\lambda}\|_{-\frac{1}{2};\Gamma_C} + \|\mathbf{u}_l - \mathbf{u}\|_{1;\Omega}).
\end{align*}
$$

Then the triangle inequality and Young's inequality applied on $\|\boldsymbol{\lambda}-\boldsymbol{\lambda}_l\|_{-\frac{1}{2};\Gamma_C}$ and (4.5), respectively, give

$$
\|\mathbf{E}_l\|_{\mathbf{V} \times \mathbf{M}}^2 \le C \left( \inf_{\mathbf{v}_l \in \mathbf{V}_l} \|\mathbf{u} - \mathbf{v}_l\|_{1;\Omega}^2 + \inf_{\boldsymbol{\mu}_l \in \mathbf{M}_l} \|\boldsymbol{\lambda} - \boldsymbol{\mu}_l\|_{-\frac{1}{2};\Gamma_C}^2 - b(\boldsymbol{\lambda} - \boldsymbol{\lambda}_l, \mathbf{e}_l) \right).
$$

We use the additive decomposition of $b(\cdot, \cdot)$ into $b_n(\cdot, \cdot) + b_t(\cdot, \cdot)$ and recall (4.2)

$$
b(\boldsymbol{\lambda}_l - \boldsymbol{\lambda}, \mathbf{e}_l) = b_n(\boldsymbol{\lambda}_l, \mathbf{u}) + b_n(\boldsymbol{\lambda}, \mathbf{u}_l) + b_t(\boldsymbol{\lambda}_l - \boldsymbol{\lambda}, \mathbf{u}) - b_t(\boldsymbol{\lambda}_l - \boldsymbol{\lambda}, \mathbf{u}_l).
$$

Each of the first three terms on the right can be bounded by the maximum
of zero and the term itself. To bound the last term, we use (4.3) and the
$H^{1/2}$-stability (see (4.4b)) of the mortar projection $\Pi_l$ defined by (3.10a).
For all $\boldsymbol{\mu}_l \in \mathbf{M}_l$ and $\mathbf{w}_l \in \mathbf{W}_l^s$, we have

$$
\begin{align*}
b_t(\boldsymbol{\lambda} - \boldsymbol{\lambda}_l, \mathbf{u}_l) &= \langle \boldsymbol{\lambda}^t - \boldsymbol{\lambda}_l^t, [\mathbf{u}_l^t] - \boldsymbol{\Pi}_l[\mathbf{u}_l^t] \rangle_{\Gamma_C} + \langle \boldsymbol{\lambda}^t - \boldsymbol{\lambda}_l^t, \boldsymbol{\Pi}_l[\mathbf{u}_l^t] \rangle_{\Gamma_C} \\
&\leq \langle \boldsymbol{\lambda}^t - \boldsymbol{\mu}_l, [\mathbf{u}_l^t] - \boldsymbol{\Pi}_l[\mathbf{u}_l^t] \rangle_{\Gamma_C} + \langle \mathcal{F}, \|\boldsymbol{\Pi}_l[\mathbf{u}_l^t]\| - \|\boldsymbol{\Pi}_l[\mathbf{u}_l^t]\|_l \rangle_{\Gamma_C} \\
&\leq \langle \boldsymbol{\lambda}^t - \boldsymbol{\mu}_l, [\mathbf{u}_l^t] - \boldsymbol{\Pi}_l[\mathbf{u}_l^t] \rangle_{\Gamma_C} \\
&\leq C \|\boldsymbol{\lambda} - \boldsymbol{\mu}_l\|_{-\frac{1}{2};\Gamma_C} (\|\mathbf{u} - \mathbf{u}_l\|_{1;\Omega} + \|[\mathbf{u}^t] - \mathbf{w}_l\|_{\frac{1}{2};\Gamma_C}).
\end{align*}
$$

Now Young's inequality gives the required bound.

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

The first two terms in the upper bound of Lemma 4.2 are the best ap-
proximation errors. They reflect the quality of the approximation of the
spaces **V**<sub>*l*</sub> and **M**<sub>*l*</sub>. The third, fourth and fifth term are consistency errors
of the approach. We remark that the term max(*b*<sub>*n*</sub>(*λ*, **u**<sub>*l*</sub>), 0) takes into ac-
count the discrete penetration of the two bodies on the actual contact set.
The term max(*b*<sub>*n*</sub>(*λ*<sub>*l*</sub>, **u**), 0) can be greater than zero if the discrete Lagrange
multiplier λ<sup>*n*</sup><sub>*l*</sub> is negative on a part of the actual contact set. We recall that
*M*<sub>*l*</sub><sup>+</sup> is not a subspace of *M*<sup>+</sup>, and thus λ<sup>*n*</sup><sub>*l*</sub> does not have to be non-negative.
To some extent this term measures the non-conformity of λ<sup>*n*</sup><sub>*l*</sub> with respect to
the physical requirement of a positive contact pressure. The fifth term sat-
isfies max(*b*<sub>*t*</sub>(*λ*<sub>*l*</sub> − *λ*, **u**), 0) ≤ max(∫<sub>Γ<sub>C</sub></sub>(‖*λ*<sub>*l*</sub><sup>*t*</sup>‖ − *F*)‖**u**<sub>*t*</sub>‖ ds, 0). Using that the
maximum of a nodal dual Lagrange multiplier is larger than one, ‖*γ*<sub>*p*</sub><sup>*t*</sup>‖ ≤ *F*
does not necessarily give ‖*λ*<sub>*l*</sub><sup>*t*</sup>‖ ≤ *F*, and thus this term is in general non-
zero and measures the violation of the friction law. Finally, the last term does not appear for contact problems without friction. A closer look into the proof reveals that the $H^{1/2}$-norm estimate on $\Gamma_C$ is too pessimistic for Coulomb problems, and it would be sufficient to consider the $H^{1/2}$-norm on the actual contact zone. Moreover, this term does not occur if we work with matching meshes.

## 4.2. Optimal a priori estimates

To prove optimal a priori error estimates under the $H^s$-regularity assumption for the displacements $\mathbf{u}$ with $1 \le s \le 2$, we have to consider in more detail the three terms in the upper bound of Lemma 4.2 which involve the bilinear form $b(\cdot, \cdot)$. We now give three lemmas providing upper bounds for these consistency errors.

**Lemma 4.3.** Let $(\mathbf{u}, \boldsymbol{\lambda}) \in \mathbf{V} \times \mathbf{M}(\mathcal{F})$ be the solution of (2.15) with $\mathbf{M}(\boldsymbol{\lambda}_n)$ replaced by $\mathbf{M}(\mathcal{F})$ and let $(\mathbf{u}_l, \boldsymbol{\lambda}_l) \in \mathbf{V}_l \times \mathbf{M}_l(\mathcal{F})$ be the solution of the discrete formulation (3.9) with $\mathbf{M}_l(\boldsymbol{\lambda}_l^n)$ replaced by $\mathbf{M}_l(\mathcal{F})$. Under the regularity assumption $\mathbf{u} \in (H^s(\Omega))^d$, $1 \le s \le 2$, we then have the a priori error estimate

$$b_n(\boldsymbol{\lambda}, \mathbf{u}_l) \le C(h_l^{2(s-1)}|\mathbf{u}|_{s;\Omega}^2 + h_l^{(s-1)}|\mathbf{u}|_{s;\Omega}\|\mathbf{u} - \mathbf{u}_l\|_{1;\Omega}).$$

*Proof.* For standard Lagrange multipliers, we refer to Hild and Laborde (2002). Although our dual basis functions of $\mathbf{M}_l$ are not positive, we can apply the same techniques. Using the discrete saddle-point formulation (3.9) and the definition of the mortar projection, we find, in terms of the approximation properties (4.4) and Remark 4.1, the upper bound

$$ \begin{aligned} b_n(\boldsymbol{\lambda}, \mathbf{u}_l) &= \langle \boldsymbol{\lambda}_n, [u_l^n] - \Pi_l[u_l^n] + \Pi_l[u_l^n] \rangle_{\Gamma_C} \le \langle \boldsymbol{\lambda}_n, [u_l^n] - \Pi_l[u_l^n] \rangle_{\Gamma_C} \\ &\le \langle \boldsymbol{\lambda}_n - \Pi_l^* \boldsymbol{\lambda}_n, [u_l^n] - \Pi_l[u_l^n] \rangle_{\Gamma_C} \\ &\le \|\boldsymbol{\lambda}_n - \Pi_l^* \boldsymbol{\lambda}_n\|_{-\frac{1}{2};\Gamma_C} \|[u_l^n] - \Pi_l[u_l^n]\|_{\frac{1}{2};\Gamma_C} \\ &\le C h_l^{s-1} |\boldsymbol{\lambda}_n|_{s-\frac{3}{2};\Gamma_C} (\|[u_l^n] - [u_n]\|_{\frac{1}{2};\Gamma_C} + \|[u_n] - \Pi_l[u_n]\|_{\frac{1}{2};\Gamma_C}) \\ &\le C(h_l^{s-1}|\mathbf{u}|_{s;\Omega}\|\mathbf{u} - \mathbf{u}_l\|_{1;\Omega} + h_l^{2(s-1)}|\mathbf{u}|_{s;\Omega}^2). \quad \square \end{aligned} $$

Before we focus on the terms $b_n(\boldsymbol{\lambda}_l, \mathbf{u})$ and $b_t(\boldsymbol{\lambda}_l - \boldsymbol{\lambda}, \mathbf{u})$, we consider a non-linear quasi-projection operator which preserves sign. This type of operator was originally introduced in Chen and Nochetto (2000). We also refer to Nochetto and Wahlbin (2002) for a negative result on the existence of higher-order sign-preserving operators and for a detailed discussion of the special role of extreme points. Let $\hat{S}_l: W \to W_l^s$ be a Clément-type operator which is defined node-wise by

$$\hat{S}_l w(p) := \frac{1}{|B_p(c_{\text{reg}}h_l)|} \int_{B_p(c_{\text{reg}}h_l)} w \, ds, \quad p \in \mathcal{P}_{C;l}^{\text{int}}.$$
Here $\mathcal{P}_{C;l}^{\text{int}} := \{p \in \mathcal{P}_{C;l}^{\text{s}}, p \notin \partial\Gamma_{\text{C}}\}$, and for all nodes $p \in \partial\Gamma_{\text{C}}$, we use a locally defined value depending only on the values of $w$ restricted to $B_p(c_{\text{reg}}h_l) \cap \Gamma_{\text{C}}$ such that $\hat{S}_l$ is $L^2$-stable and reproduces polynomials of degree one. Standard arguments show that for $1 \le s \le 2$ we get

$$ \|\hat{S}_l w - w\|_{\frac{1}{2};\Gamma_{\text{C}}} \le C h_l^{s-1} |w|_{s-\frac{1}{2};\Gamma_{\text{C}}}. $$

More importantly, $\hat{S}_l w$ preserves the sign of $w \in W^+$ in the mesh-dependent interior $\Gamma_{C;l}^{\text{int}} := \Gamma_{\text{C}} \setminus (\cup_{p \in \mathcal{P}_{C;l}^{\text{s}} \setminus \mathcal{P}_{C;l}^{\text{int}}} \text{supp } \phi_p)$, i.e., we have for $w \in W^+$ that $\hat{S}_l w(p) \ge 0$ for $p \in \mathcal{P}_{C;l}^{\text{int}}$.

In terms of the linear operator $\hat{S}_l$, we define the non-linear operator $S_l$:

$$ S_l w(p) := \begin{cases} \hat{S}_l w(p) & \text{supp } \phi_p \subset \text{supp } w \text{ or } p \in \partial\Gamma_{\text{C}}, \\ 0 & \text{otherwise,} \end{cases} \quad p \in \mathcal{P}_{C;l}^{\text{s}}, \quad (4.6) $$

the definition of which guarantees that $\text{supp } S_l w \cap \Gamma_{C;l}^{\text{int}} \subset \text{supp } w$. Moreover, for $w \in W^+$ we have $S_l w(p) \ge 0$ for $p \notin \partial\Gamma_{\text{C}}$.

**Assumption 4.4.** Let us define $\Sigma_l^n := \{x \in \Gamma_{\text{C}} | \text{dist}(x, \partial B_n) \le 2h_l\}$, where $B_n$ is the actual contact zone, i.e., $B_n := \text{supp } \lambda_n$. Then, we assume that $\Sigma_l^n$ and $B_n$ are compactly embedded in $\Gamma_{\text{C}}$ and $\Gamma_{C;l}^{\text{int}}$, respectively, and moreover that

$$ \|[u_n]\|_{0;\Sigma_l^n} \le C h_l^{s-\frac{1}{2}}|[u_n]|_{s-\frac{1}{2};\Gamma_{\text{C}}}. $$

Let us briefly comment on different aspects of this assumption. We note that for $h_l$ small enough, $\Sigma_l^n$ and $B_n$ are compactly embedded in $\Gamma_{\text{C}}$ and $\Gamma_{C;l}^{\text{int}}$, respectively, due to the assumption that $B_n$ is compactly embedded in $\Gamma_{\text{C}}$. This assumption can be weakened, but then the notation would become more technical. Setting $B_n^c := \Gamma_{\text{C}} \setminus B_n$ and defining $H_{00}^{s-1/2}(B_n^c) := \{w \in L^2(B_n^c) | w = v|_{B_n^c} \text{ for } v \in H^{s-1/2}(\Gamma_{\text{C}}) \text{ and } \text{supp } v \subset \overline{B_n^c}\}$, we get $[u_n] \in H_{00}^{s-1/2}(B_n^c)$ if $u \in (H^s(\Omega))^d$. Now, if $B_n^c$ is regular enough, the assumption is followed by a Poincaré–Friedrichs-type argument, together with suitable interpolation and a scaling. We refer to Li, Melenk, Wohlmuth and Zou (2010), where similar estimates for interfaces have been used, and to Melenk and Wohlmuth (2011), where these types of estimates are used to obtain quasi-optimal a priori $L^2$-norm estimates for the Lagrange multiplier in a linear mortar setting. In particular, an order-$h_l$ estimate is given for $H^1$-functions with vanishing trace. The assumption is naturally satisfied if $B_n$ is regular enough. If the boundary $\partial B_n$ is smooth enough, we can locally flatten $\partial B_n$, use the fact that Sobolev spaces are invariant under smooth changes of variables and apply the 1D Sobolev embedding result $\|v\|_{L^\infty(I)} \lesssim \|v\|_{H^s(I)}$, $s > \frac{1}{2}$ recursively, where $I$ is a fixed interval and $v \in H^s(I)$ (see, e.g., Adams (1975)).
<figure><img src="image_27.png" /><figcaption>Figure 4.4. Actual contact zone $B_n$ (a) and boundary strip $\Sigma_l^n$ (b).</figcaption></figure>

Figure 4.4 illustrates the definition of the strip $\Sigma_l^n$ which has a diameter of $4h_l$ perpendicular to $\partial B_n$.

In terms of these preliminary considerations, we can show the following bound for the consistency error of $\lambda_l^n$.

**Lemma 4.5.** Let $(\mathbf{u}, \boldsymbol{\lambda}) \in \mathbf{V} \times \mathbf{M}(\mathcal{F})$ be the solution of (2.15) with $\mathbf{M}(\lambda_n)$ replaced by $\mathbf{M}(\mathcal{F})$ and let $(\mathbf{u}_l, \lambda_l) \in \mathbf{V}_l \times \mathbf{M}_l(\mathcal{F})$ be the solution of the discrete formulation (3.9) with $\mathbf{M}_l(\lambda_l^n)$ replaced by $\mathbf{M}_l(\mathcal{F})$. Under Assumption 4.4 and the regularity assumption $\mathbf{u} \in (H^s(\Omega))^d$, $1 \le s \le 2$, we then have the *a priori* error estimate

$$b_n(\lambda_l, \mathbf{u}) \le C h_l^{s-1} |\mathbf{u}|_{s;\Omega} \|\boldsymbol{\lambda} - \lambda_l\|_{-\frac{1}{2};\Gamma_C}.$$

*Proof.* The operator $\hat{S}_l$ yields $\|\hat{S}_l[u_n] - [u_n]\|_{\frac{1}{2};\Gamma_C} \le C h_l^{s-1}|[u_n]|_{s-\frac{1}{2};\Gamma_C} \le C h_l^{s-1}|\mathbf{u}|_{s;\Omega}$ for $1 \le s \le 2$. Due to the fact that $\lambda_l^n \in M^+$ with $\gamma_p^n = 0$ for $p \in \partial\Gamma_C$, we find

$$\int_{\Gamma_C} \lambda_l^n S_l[u_n] \, ds = \sum_{p \in \mathcal{P}_{C;l}^{\text{int}}} \gamma_p^n S_l[u_n](p) m_p \le 0.$$

Then, (4.2), the construction of $S_l$ and the assumption $B_n \subset \Gamma_C^{\text{int}}$ yield

$$\begin{aligned} b_n(\lambda_l, \mathbf{u}) &= \int_{\Gamma_C} \lambda_l^n[u_n] \, ds \le \int_{\Gamma_C} \lambda_l^n([u_n] - S_l[u_n]) \, ds \\ &= \langle \lambda_l^n - \lambda_n, [u_n] - S_l[u_n] \rangle_{\Gamma_C} \le \|\lambda_l^n - \lambda_n\|_{-\frac{1}{2};\Gamma_C} \|[u_n] - S_l[u_n]\|_{\frac{1}{2};\Gamma_C} \\ &\le \|\lambda_l^n - \lambda_n\|_{-\frac{1}{2};\Gamma_C} (\|[u_n] - \hat{S}_l[u_n]\|_{\frac{1}{2};\Gamma_C} + \|S_l[u_n] - \hat{S}_l[u_n]\|_{\frac{1}{2};\Gamma_C}). \end{aligned}$$

As already noted the linear operator $\hat{S}_l$ has best approximation properties and thus it is sufficient to consider the second term on the right in more detail. The properties of the non-linear operator $S_l$ defined by (4.6) play a crucial role.

We observe that $S_l[u_n]$ and $\hat{S}_l[u_n]$ coincide in the two dark grey-shaded regions of Figure 4.4(a) and note that $S_l[u_n] - \hat{S}_l[u_n] = 0$ on $\Gamma_C \setminus \Sigma_l^n$. For the
<figure><img src="image_28.png" /><figcaption>Figure 4.5. Friction cone for $d = 2$ (a) and $d = 3$ (b).</figcaption></figure>

term $\|S_l[u_n] - \hat{S}_l[u_n]\|_{\frac{1}{2};\Gamma_C}$, we can now apply a standard inverse inequality and get

$$ \begin{aligned} \|S_l[u_n] - \hat{S}_l[u_n]\|_{\frac{1}{2};\Gamma_C}^2 &\le \frac{C}{h_l} \|S_l[u_n] - \hat{S}_l[u_n]\|_{0;\Gamma_C}^2 \\ &\le C \sum_{p \in \mathcal{P}_{\Sigma_l^n}} h_l^{d-2} (\hat{S}_l[u_n](p))^2 \le \frac{C}{h_l} \|[u_n]\|_{0;\Sigma_l^n}^2, \end{aligned} $$

where $\mathcal{P}_{\Sigma_l^n} := \{p \in \mathcal{P}_{\mathcal{C};l}^s \text{ such that } \text{supp } \phi \subset \Sigma_l^n\}$. Now Assumption 4.4 can be applied and yields the required bound. □

Now we combine the previous results and formulate a first optimal a priori error estimate for a two-body contact problem with no friction.

**Theorem 4.6.** Let $(u, \lambda) \in V \times M^+$ be the solution of (2.15) and let $(u_l, \lambda_l) \in V_l \times M_l^+$ be the solution of the discrete formulation (3.9) with $\nu = 0$. Under the Assumption 4.4 and the regularity assumption $u \in (H^s(\Omega))^d$, $1 \le s \le 2$, we then have the a priori error estimate

$$ \|u - u_l\|_{1;\Omega} + \|\lambda - \lambda_l\|_{-\frac{1}{2};\Gamma_C} \le Ch_l^{s-1} |u|_{s;\Omega}. $$

*Proof.* Using the well-known approximation property for the spaces $V_l$ and $M_l$, the proof is a direct consequence of Lemmas 4.2–4.5 by applying Young’s inequality and noting that $\lambda_t = \lambda_l^t = 0$. □

For a non-trivial given friction bound $\mathcal{F}$, the situation is more complex, and moreover there is a substantial difference between the two- and three-dimensional setting. In 2D, the tangential stress component can be identified with a scalar-valued functional, and we can follow the proof of Lemma 4.5. Figure 4.5 illustrates the difference. In 2D, the tangential surface traction $\lambda_t = \lambda_t t$ of a regular solution is either $\lambda_t = \mathcal{F}$ or $\lambda_t = -\mathcal{F}$ for sliding nodes. Moreover, the tangential displacement $[u_t] = [u_t]t$ can be separated into two parts $[u_t] = \max(0, [u_t]) + \min(0, [u_t]) =: u_t^+ - u_t^-$, $u_t^+$,
<figure><img src="image_29.png" /><figcaption>Figure 4.6. Tangential displacement for contact problems with friction: $d = 2$ (a) and $d = 3$ (b,c).</figcaption></figure>

$u_t^- \in W^+$: see Figure 4.6(a). In 3D, the tangential part $[u_t]$ is a vector-valued function, for which we cannot apply the friction law component-wise.

In Figure 4.6, we show typical numerical results for a contact problem with friction. We note that for these considerations there is no difference between a Tresca and a Coulomb problem. In Figure 4.6(a) a Hertz contact problem with a given but non-constant friction bound is simulated in 2D. In Figure 4.6(b), we show a visualization of the tangential displacement for $d = 3$ and Coulomb friction. The nodes in the centre marked with black bullets are sticky and do not carry a relative tangential displacement. For these nodes the tangential component of the stress satisfies the inequality of the Coulomb law strictly. The nodes situated in the outer ring do slide, and the sliding direction is possibly changing from node to node. As can be seen in the zoom in Figure 4.6(c), the sliding direction is, as required by the friction law, opposite to the tangential stress.

**Assumption 4.7. (2D setting)** Let us define

$$ \Sigma_l^{\pm;t} := \{x \in \Gamma_C \mid \mathrm{dist}(x, \partial B_t^{\pm}) \le 2h_l\}, $$

where $\partial B_t^+ \cup \partial B_t^-$ is the boundary of the actual sticky zone. More precisely, we set $\partial B_t^\pm := \partial \mathrm{supp} u_t^\pm \cap \Gamma_C$. Then, we assume that $\Sigma_l^{\pm;t}$ and $B_t^\pm$ are compactly embedded in $\Gamma_C$, and moreover that for $p \in \partial \Gamma_C \cap \mathrm{supp} u_t^\pm$ we have $\gamma_p^t = \pm \mathcal{F}$,

$$ \|u_t^\pm\|_{0;\Sigma_l^{\pm;t}} \le C h_l^{s-\frac{1}{2}} |[u_t]|_{s-\frac{1}{2};\Gamma_C}. $$

We note that for a Coulomb problem, we have $B_t \subset B_n$, and thus $B_t$ is automatically compactly embedded $\Gamma_C$ if $B_n$ is so. This is not necessarily the case for a Tresca friction problem. We point out that this assumption rules out the case $s = 2$ and $\lambda_t = \lambda_t t$ with $\lambda_t$ being a function which jumps from plus to minus of the friction bound. In that case $u_t^+$ and $u_t^-$ have a lower regularity than $u_t$. To be more precise $u_t \in H^{3/2}(\Gamma_C)$, whereas
$u_t^\pm \in H^{3/2-\epsilon}(\Gamma_C)$ for all $\epsilon > 0$, and thus only an order $h_l^{3/2-\epsilon}$ can be expected to hold true. Figure 4.7 shows the numerical solution for a Coulomb problem in 2D with two different friction coefficients. In Figure 4.7(a) the case $\nu = 0.8$ is shown whereas in Figure 4.7(b,c) the case $\nu = 0.3$ is presented. The close-up in Figure 4.7(c) reveals that $u_t$ is zero not only in the centre point but also in a non-trivial sub-interval, and thus Assumption 4.7 is satisfied for both cases.

The following lemma is the counterpart of Lemma 4.5 for the tangential component. Due to the partition $u_t = u_t^+ - u_t^-$ it only holds true for $d = 2$.

**Lemma 4.8.** Let $(\mathbf{u}, \boldsymbol{\lambda}) \in \mathbf{V} \times \mathbf{M}(\mathcal{F})$ be the solution of (2.15) with $\mathbf{M}(\lambda_n)$ replaced by $\mathbf{M}(\mathcal{F})$ and let $(\mathbf{u}_l, \boldsymbol{\lambda}_l) \in \mathbf{V}_l \times \mathbf{M}_l(\mathcal{F})$ be the solution of the discrete formulation (3.9) with $\mathbf{M}_l(\lambda_l^n)$ replaced by $\mathbf{M}_l(\mathcal{F})$. Under Assumption 4.7 and the regularity assumption $\mathbf{u} \in (H^s(\Omega))^2$, $1 \le s \le 2$, we then have the a priori error estimate for $d = 2$:

$$b_t(\boldsymbol{\lambda}_l - \boldsymbol{\lambda}, \mathbf{u}) \le C h_l^{s-1} |\mathbf{u}|_{s;\Omega} \|\boldsymbol{\lambda} - \boldsymbol{\lambda}_l\|_{-\frac{1}{2};\Gamma_C}.$$

*Proof.* The proof follows the lines of the proof of Lemma 4.5. We start with the observation that $b_t(\boldsymbol{\lambda}_l - \boldsymbol{\lambda}, \mathbf{u}) = \int_{\Gamma_C} (\lambda_l^t - \mathcal{F})u_t^+ ds - \int_{\Gamma_C} (\lambda_l^t + \mathcal{F})u_t^- ds$, where $\boldsymbol{\lambda}_l^t := \lambda_l^t \mathbf{t}$. Now we apply the operator $S_l$ to $u_l^+$ and $u_l^-$. Then, under Assumption 4.7 we get

$$\|u_t^+ - S_l u_t^+\|_{\frac{1}{2};\Gamma_C} + \|u_t^- - S_l u_t^-\|_{\frac{1}{2};\Gamma_C} \le C h_l^{s-1} \|[u_n]\|_{s-\frac{1}{2};\Gamma_C}.$$

The construction of $S_l$ yields that

$$
\begin{align*}
b_t(\boldsymbol{\lambda}_l - \boldsymbol{\lambda}, \mathbf{u}) &= \langle \lambda_l^t - \lambda_t, u_t^+ - S_l u_t^+ \rangle + \langle \lambda_l^t - \lambda_t, S_l u_t^- - u_t^- \rangle \\
&\quad + \langle \lambda_l^t - \lambda_t, S_l u_t^+ \rangle - \langle \lambda_l^t - \lambda_t, S_l u_t^- \rangle \\
&\le C h_l^{s-1} \|[u_n]\|_{s-\frac{1}{2};\Gamma_C} \|\lambda_l^t - \lambda_t\|_{-\frac{1}{2};\Gamma_C} \\
&\quad + \langle \lambda_l^t - \mathcal{F}, S_l u_t^+ \rangle - \langle \lambda_l^t + \mathcal{F}, S_l u_t^- \rangle \\
&= C h_l^{s-1} \|[u_n]\|_{s-\frac{1}{2};\Gamma_C} \|\lambda_l^t - \lambda_t\|_{-\frac{1}{2};\Gamma_C} \\
&\quad + \sum_{p \in \mathcal{P}_{C;l}^s} (\gamma_p^t - \mathcal{F}) S_l u_t^+(p) m_p \\
&\quad - \sum_{p \in \mathcal{P}_{C;l}^s} (\gamma_p^t + \mathcal{F}) S_l u_t^-(p) m_p.
\end{align*}
$$

Moreover, by definition of $M_l(\mathcal{F})$, it is easy to see that $\gamma_p^t - \mathcal{F} \le 0$ and
<figure><img src="image_30.png" /><figcaption>Figure 4.7. Tangential displacement for different friction coefficients.</figcaption></figure>

$\gamma_p^t + \mathcal{F} \geq 0$, and thus we get

$$ \sum_{p \in \mathcal{P}_{C;l}^{\text{int}}} [(\gamma_p^t - \mathcal{F})S_l u_t^+(p) - (\gamma_p^t + \mathcal{F})S_l u_t^-(p)] m_p \leq 0. $$

For $p \in \mathcal{P}_{C;l}^{\text{s}} \setminus \mathcal{P}_{C;l}^{\text{int}}$, we have that $(\gamma_p^t - \mathcal{F})S_l u_t^+(p) = 0 = (\gamma_p^t + \mathcal{F})S_l u_t^-(p)$. We note that these last arguments are typical for $d = 2$ but cannot be applied in 3D. $\square$

We now combine the previous results and formulate an optimal a priori error estimate for a two-body contact problem with Tresca friction in 2D.

**Theorem 4.9.** Let $(\mathbf{u}, \boldsymbol{\lambda}) \in \mathbf{V} \times \mathbf{M}(\mathcal{F})$ be the solution of (2.15), with $\mathbf{M}(\lambda_n)$ replaced by $\mathbf{M}(\mathcal{F})$, and let $(\mathbf{u}_l, \boldsymbol{\lambda}_l) \in \mathbf{V}_l \times \mathbf{M}_l(\mathcal{F})$ be the solution of the discrete formulation (3.9), with $\mathbf{M}_l(\lambda_l^n)$ replaced by $\mathbf{M}_l(\mathcal{F})$. Under Assumptions 4.4, 4.7 and the regularity assumption $\mathbf{u} \in (H^s(\Omega))^2$, $1 \leq s \leq 2$, we then have the a priori error estimate for $d = 2$:

$$ \|\mathbf{u} - \mathbf{u}_l\|_{1;\Omega} + \|\boldsymbol{\lambda} - \boldsymbol{\lambda}_l\|_{-\frac{1}{2};\Gamma_C} \leq C h_l^{s-1} |\mathbf{u}|_{s;\Omega}. $$

*Proof.* Using the well-known approximation property for the spaces $\mathbf{V}_l$, $\mathbf{M}_l$ and $\mathbf{W}_l^s$, the proof is a direct consequence of Lemmas 4.2–4.8 by applying Young’s inequality. $\square$

Let us briefly comment on the three-dimensional case. Lemma 4.8 is the only one where we have explicitly used a 2D construction. All the other results hold true for $d = 2$ and $d = 3$. In 3D, we do still get a priori bounds for the discretization error, although the optimal order is lost if $\mathbf{u} \in (H^s(\Omega))^d$, $3/2 < s \leq 2$. We introduce the operators $\mathbf{Z}_l^* : \mathbf{M} \to \mathbf{W}_l^s$ and $\mathbf{Z}_l : \mathbf{W} \to \mathbf{W}_l^s$ by $\mathbf{Z}_l^* := (Z_l^*)^d$ and $\mathbf{Z}_l := (Z_l)^d$, $Z_l^* : \mathbf{M} \to \mathbf{W}_l^s$ and $\mathbf{Z}_l : \mathbf{W} \to \mathbf{W}_l^s$:

$$ Z_l^* \mu := \sum_{p \in \mathcal{P}_{C;l}^{\text{s}}} \frac{\langle \mu, \phi_p \rangle_{\Gamma_C}}{m_p} \phi_p, \quad Z_l w := \sum_{p \in \mathcal{P}_{C;l}^{\text{s}}} \frac{\langle \phi_p, w \rangle_{\Gamma_C}}{m_p} \phi_p. \quad (4.7) $$
Obviously, both $Z_l$ and $Z_l^*$ are locally defined and reproduce constants. Moreover, $Z_l$ is $L^2$- and $H^1$-stable. From the $H^{1/2}$-stability we get the $H^{-1/2}$-stability of $Z_l^*$ with the same stability constant

$$ \|Z_l^* \mu\|_{-\frac{1}{2};\Gamma_C} = \sup_{w \in W} \frac{\langle Z_l^* \mu, w \rangle_{\Gamma_C}}{\|w\|_{\frac{1}{2};\Gamma_C}} = \sup_{w \in W} \frac{\langle \mu, Z_l w \rangle_{\Gamma_C}}{\|w\|_{\frac{1}{2};\Gamma_C}}. $$

In terms of these operators, we obtain the following non-optimal a priori bound for a Tresca friction problem in 3D with non-trivial but constant friction coefficient.

**Theorem 4.10.** Let $(\mathbf{u}, \boldsymbol{\lambda}) \in \mathbf{V} \times \mathbf{M}(\mathcal{F})$ be the solution of (2.15), with $\mathbf{M}(\lambda_n)$ replaced by $\mathbf{M}(\mathcal{F})$, and let $(\mathbf{u}_l, \boldsymbol{\lambda}_l) \in \mathbf{V}_l \times \mathbf{M}_l(\mathcal{F})$ be the solution of the discrete formulation (3.9), with $\mathbf{M}_l(\lambda_l^n)$ replaced by $\mathbf{M}_l(\mathcal{F})$. If $\mathbf{u} \in (H^s(\Omega))^d$, $1 \le s \le 3/2$, we then have the a priori error estimate

$$ \|\mathbf{u} - \mathbf{u}_l\|_{1;\Omega} + \|\boldsymbol{\lambda} - \boldsymbol{\lambda}_l\|_{-\frac{1}{2};\Gamma_C} \le Ch_l^{s-1}|\mathbf{u}|_{s;\Omega}. $$

*Proof.* We have to re-examine the two terms $b_n(\boldsymbol{\lambda}_l, \mathbf{u})$ and $b_t(\boldsymbol{\lambda}_l - \boldsymbol{\lambda}, \mathbf{u})$ in the upper bound of Lemma 4.2. To do so, we apply the operators $Z_l^*$ and $\mathbf{Z}_l^*$ and remark that $Z_l^* \lambda_l^n = \sum_{p \in \mathcal{P}_{C;l}^s} \gamma_p^n \phi_p \in W^+$ and

$$ \|\mathbf{Z}_l^* \boldsymbol{\lambda}_l^t\| - \mathcal{F} = \left\| \sum_{p \in \mathcal{P}_{C;l}^s} \gamma_p^t \phi_p \right\| - \mathcal{F} \le \sum_{p \in \mathcal{P}_{C;l}^s} \|\gamma_p^t\| \phi_p - \mathcal{F} \le \sum_{p \in \mathcal{P}_{C;l}^s} \mathcal{F} \phi_p - \mathcal{F} = 0. $$

These preliminary observations in combination with (3.7) yield for the normal part

$$ \begin{aligned} b_n(\boldsymbol{\lambda}_l, \mathbf{u}) &= \langle \lambda_l^n - Z_l^* \lambda_l^n, [u_n] \rangle_{\Gamma_C} + \langle Z_l^* \lambda_l^n, [u_n] \rangle_{\Gamma_C} \\ &\le \langle \lambda_l^n - Z_l^* \lambda_l^n, [u_n] \rangle_{\Gamma_C} = \langle \lambda_l^n - Z_l^* \lambda_l^n, [u_n] - \Pi_{0;l}[u_n] \rangle_{\Gamma_C} \\ &\le Ch_l^{s-\frac{1}{2}} \|\lambda_l^n - Z_l^* \lambda_l^n\|_{0;\Gamma_C} |[u_n]|_{s-\frac{1}{2};\Gamma_C} \\ &\le Ch_l^{s-1} \|\lambda_l^n - Z_l^* \lambda_l^n\|_{-\frac{1}{2};\Gamma_C} |\mathbf{u}|_{s;\Omega}, \end{aligned} $$

where $\Pi_{0;l}$ is the $L^2$-projection onto element-wise constants. Here we have also used additionally the inverse estimate for $\|\lambda_l^n - Z_l^* \lambda_l^n\|_{0;\Gamma_C}$, which results from standard inverse estimates for finite elements and the fact that $Z_l^* \lambda_l^n \ge 0$.

Keeping in mind that $b_t(\boldsymbol{\lambda}, \mathbf{u}) = \int_{\Gamma_C} \mathcal{F} \|[\mathbf{u}]_t\| \mathrm{d}s$, the tangential part can be estimated in the same way:

$$ \begin{aligned} b_t(\boldsymbol{\lambda}_l - \boldsymbol{\lambda}, \mathbf{u}) &= b_t(\boldsymbol{\lambda}_l - \mathbf{Z}_l^* \boldsymbol{\lambda}_l, \mathbf{u}) + b_t(\mathbf{Z}_l^* \boldsymbol{\lambda}_l - \boldsymbol{\lambda}, \mathbf{u}) \\ &\le b_t(\boldsymbol{\lambda}_l - \mathbf{Z}_l^* \boldsymbol{\lambda}_l, \mathbf{u}) \le Ch_l^{s-1} \|\boldsymbol{\lambda}_l^t - \mathbf{Z}_l^* \boldsymbol{\lambda}_l^t\|_{-\frac{1}{2};\Gamma_C} |\mathbf{u}|_{s;\Omega}. \end{aligned} $$

In a last step, we have to consider $\|\boldsymbol{\lambda}_l - \mathbf{Z}_l^* \boldsymbol{\lambda}_l\|_{-\frac{1}{2};\Gamma_C}$ in more detail and
bound it. The stability of $Z_l^*$ in the $H^{-1/2}$-norm gives

$$ \begin{aligned} \|\lambda_l - Z_l^* \lambda_l\|_{-\frac{1}{2}; \Gamma_C} &\le C \|\lambda_l - \lambda\|_{-\frac{1}{2}; \Gamma_C} + \|\lambda - Z_l^* \lambda\|_{-\frac{1}{2}; \Gamma_C} \\ &\le C (\|\lambda_l - \lambda\|_{-\frac{1}{2}; \Gamma_C} + h_l^{s-1} |\mathbf{u}|_{s; \Omega}). \end{aligned} $$

Now, Lemma 4.3 in combination with Young's inequality yields the required a priori bound. □

**Remark 4.11.** We note that in contrast to Theorem 4.9 no additional assumptions on the actual contact zones are made in Theorem 4.10.

The main advantage of the dual Lagrange multiplier space is the possibility of computing $\lambda_l$ by a local post-process from the discrete displacement. Taking the local residual and using a simple scaling directly yield the coefficients. However, for the visualization in general $Z_l^* \lambda_l$ is plotted and not $\lambda_l$. However, both quantities have the same order of convergence. For $1 \le s \le 2$ and $\mathbf{u} \in (H^s(\Omega))^d$ we obtain in terms of the $H^{-1/2}$-stability of $Z_l^*$ and its approximation property

$$ \begin{aligned} \|\lambda - Z_l^* \lambda_l\|_{-\frac{1}{2}; \Gamma_C} &\le \|\lambda - Z_l^* \lambda\|_{-\frac{1}{2}; \Gamma_C} + \|Z_l^*(\lambda - \lambda_l)\|_{-\frac{1}{2}; \Gamma_C} \\ &\le C (\|\lambda - \lambda_l\|_{-\frac{1}{2}; \Gamma_C} + h_l^{s-1} |\mathbf{u}|_{s; \Omega}). \end{aligned} $$

**Remark 4.12.** Quite often in the context of mortar methods (see, e.g., Braess and Dahmen (2002)), one prefers to work with a norm that is easier to handle than the $H^{-1/2}$-norm. Thus this is replaced by a weighted mesh-dependent $L^2$-norm,

$$ \|\mu\|_{M_l}^2 := \sum_{F \in \mathcal{F}_{C;l}^s} h_F \|\mu\|_{0;F}^2 \quad \mu \in L^2(\Gamma_C^s) \text{ or } \mu \in (L^2(\Gamma_C^s))^d. \quad (4.8) $$

Then all our theoretical results also cover $\|\lambda - \lambda_l\|_{M_l}$ if the regularity of the solution is good enough. The proof follows exactly the same lines as for the $H^{-1/2}$-norm and uses an inverse estimate. It is known from the linear mortar setting that in the case of weighted $L^2$-norms, a uniform inf-sup condition also holds. Moreover, one can replace, in the proof of the best approximation properties of the constrained space, the discrete harmonic extension by a discrete zero extension to the interior nodes.

## 4.3. Numerical results

We note that in all our numerical results the mesh-dependent norm (4.8) has been used to measure the discretization error in the Lagrange multiplier and $Z_l^* \lambda_l$ has been used to make the plots. To illustrate the convergence rates of low-order finite elements numerically, we consider two simple two-dimensional test settings. As a reference solution, we use the numerical solution on $\mathcal{T}_{l_{\max}+2}$.
<figure><img src="image_31.png" /></figure>

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

<figure><img src="image_33.png" /><figcaption>Figure 4.8. Problem setting (a), tangential Lagrange multiplier (b) and tangential displacement (c) on $\Gamma_{C;l}^s$.</figcaption></figure>

Our first example is a scalar-valued model problem for anti-plane friction; we refer to Hüeber et al. (2005b) for details regarding the problem formulation. The friction bound is set to $F = 0.6$. Figure 4.8 shows the geometry as well as the tangential Lagrange multiplier and the displacement on master and slave sides. In Figure 4.8(b), we can clearly observe the discrete complementarity (3.12) with $\nu\gamma_p^n$ being replaced by the given bound $\mathcal{F}$. In Figure 4.8(c) the tangential displacements are shown. For $\|\gamma_p^t\|$ strictly smaller than $\mathcal{F}$, the two bodies have to stick together. Sliding is possible only on a part of the contact zone where $\|\gamma_p^t\| = \mathcal{F}$.

Table 4.1 shows the convergence rates for this simplified frictional contact problem. We show the $L^2$-norm and the $H^1$-norm of the displacement error and the mesh-dependent norm of the error in the Lagrange multiplier. Our numerical results confirm the theoretical ones, and the computed rates are fairly close to the optimal order of convergence. Although we do not have a theoretical result for the $L^2$-norm, the numerical results also show a significantly better rate compared to the $H^1$-norm. Standard duality techniques such as the Aubin–Nitsche approach are tricky to apply in the setting of variational inequalities since they depend on regularity assumptions of the dual problem.

Our second test example is the classical Hertz problem (Hertz 1882, Johnson 1985, Kikuchi and Oden 1988) with Coulomb friction. Although our theoretical results do not cover the case of a Coulomb contact problem, Table 4.2 shows that also for this case we obtain almost optimal convergence rates. In that case, the convergence order for the error in the Lagrange multiplier in the weighted $L^2$-norm is close to the best approximation order of 3/2.

**Remark 4.13.** This effect is also numerically well observed for mortar problems in the linear setting (Wohlmuth 2001). A theoretical analysis can be found in the recent contribution by Melenk and Wohlmuth (2011), where it is shown that quasi-optimal $L^2$-norm estimates for the Lagrange multiplier hold in the linear mortar setting under suitable regularity assumptions.
<table><caption>Table 4.1. Convergence rates for a contact problem with given friction bound.</caption><thead><tr><th>Level</th><th>||<b>u</b><sub>l</sub> - <b>u</b><sub>ref</sub>||<sub>0;Ω</sub> / ||<b>u</b><sub>ref</sub>||<sub>0;Ω</sub></th><th>||<b>u</b><sub>l</sub> - <b>u</b><sub>ref</sub>||<sub>1;Ω</sub> / ||<b>u</b><sub>ref</sub>||<sub>1;Ω</sub></th><th>||<b>λ</b><sub>l</sub> - <b>λ</b><sub>ref</sub>||<sub>*</sub></th></tr></thead><tbody><tr><td>0</td><td>8.9343e-02</td><td>-</td><td>3.2317e-01</td><td>-</td></tr><tr><td>1</td><td>3.0960e-02</td><td>1.53</td><td>1.8933e-01</td><td>0.77</td></tr><tr><td>2</td><td>8.6519e-03</td><td>1.84</td><td>1.0343e-01</td><td>0.87</td></tr><tr><td>3</td><td>2.4272e-03</td><td>1.83</td><td>5.6431e-02</td><td>0.87</td></tr><tr><td>4</td><td>6.6657e-04</td><td>1.86</td><td>3.0414e-02</td><td>0.89</td></tr><tr><td>5</td><td>1.7918e-04</td><td>1.90</td><td>1.6132e-02</td><td>0.91</td></tr><tr><td>6</td><td>4.6522e-05</td><td>1.95</td><td>8.3113e-03</td><td>0.96</td></tr></tbody></table>

<table><caption>Table 4.2. Convergence rates for a Hertz contact problem with Coulomb friction, $\nu = 0.5$.</caption><thead><tr><th>Level</th><th>||<b>u</b><sub>l</sub> - <b>u</b><sub>ref</sub>||<sub>1;Ω</sub> / ||<b>u</b><sub>ref</sub>||<sub>1;Ω</sub></th><th>||<b>λ</b><sub>l</sub> - <b>λ</b><sub>ref</sub>||<sub>*</sub></th></tr></thead><tbody><tr><td>1</td><td>4.465867e-01</td><td>-</td></tr><tr><td>2</td><td>3.056095e-01</td><td>0.55</td></tr><tr><td>3</td><td>1.693210e-01</td><td>0.85</td></tr><tr><td>4</td><td>9.155008e-02</td><td>0.89</td></tr><tr><td>5</td><td>4.857727e-02</td><td>0.91</td></tr><tr><td>6</td><td>2.450933e-02</td><td>0.99</td></tr></tbody></table>

Finally, we briefly comment on higher-order elements; see also Belhachmi and Ben Belgacem (2000) for an analysis and Fischer and Wriggers (2006) and Puso et al. (2008) for simulation results in applications, and the influence of the choice of the Lagrange multiplier space. Recently, hp-techniques have also been applied for contact problems in combination with boundary elements (Chernov, Maischak and Stephan 2008). The p-version in the boundary element method for contact is discussed in Gwinner (2009). From the algorithmic point of view, higher-order elements can be easily applied. As we will see, quadratic elements do yield higher-order *a priori* estimates, but optimal quadratic order cannot be achieved. This results from the fact that there exists no monotonicity-preserving operator of higher order; see also Nochetto and Wahlbin (2002). In addition, the solution of a contact problem is, in general, not in $H^3(\Omega)$, and thus also, from the point of view of best approximation, no second-order error decay can be expected. Nevertheless, a higher-order *a priori* estimate can be obtained by replacing $V_l$ by quadratic finite elements. Quadratic finite elements and linear dual
<table><caption>Table 4.3. Relative error for the displacement in the $H^1$-norm for linear and quadratic finite elements.</caption><thead><tr><th>Level</th><th colspan="2">(i, j) = (1, 1)</th><th colspan="2">(i, j) = (2, 1)</th><th colspan="2">(i, j) = (2, 2)</th></tr></thead><tbody><tr><td>0</td><td>4.663632e-01</td><td>-</td><td>3.159307e-01</td><td>-</td><td>3.903263e-01</td><td>-</td></tr><tr><td>1</td><td>3.214737e-01</td><td>0.54</td><td>1.592747e-01</td><td>0.99</td><td>1.376072e-01</td><td>1.50</td></tr><tr><td>2</td><td>1.807130e-01</td><td>0.83</td><td>6.777325e-02</td><td>1.23</td><td>5.656398e-02</td><td>1.28</td></tr><tr><td>3</td><td>9.735853e-02</td><td>0.89</td><td>2.992646e-02</td><td>1.18</td><td>2.422295e-02</td><td>1.22</td></tr><tr><td>4</td><td>5.111965e-02</td><td>0.93</td><td>1.340727e-02</td><td>1.16</td><td>1.028243e-02</td><td>1.23</td></tr><tr><td>5</td><td>2.584391e-02</td><td>0.98</td><td>-</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

<table><caption>Table 4.4. Mesh-dependent $L^2$-error for the Lagrange multiplier in the linear and quadratic approach.</caption><thead><tr><th>Level</th><th colspan="2">(i, j) = (1, 1)</th><th colspan="2">(i, j) = (2, 1)</th><th colspan="2">(i, j) = (2, 2)</th></tr></thead><tbody><tr><td>0</td><td>5.845412e+01</td><td>-</td><td>5.849757e+01</td><td>-</td><td>1.323412e+02</td><td>-</td></tr><tr><td>1</td><td>4.999477e+01</td><td>0.23</td><td>4.129640e+01</td><td>0.50</td><td>3.621992e+01</td><td>1.87</td></tr><tr><td>2</td><td>2.121223e+01</td><td>1.24</td><td>1.814467e+01</td><td>1.19</td><td>1.389391e+01</td><td>1.38</td></tr><tr><td>3</td><td>8.378905e+00</td><td>1.34</td><td>7.316218e+00</td><td>1.31</td><td>5.230080e+00</td><td>1.41</td></tr><tr><td>4</td><td>3.269796e+00</td><td>1.36</td><td>2.813967e+00</td><td>1.38</td><td>2.015976e+00</td><td>1.38</td></tr><tr><td>5</td><td>1.168347e+00</td><td>1.48</td><td>-</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

Lagrange multipliers yield an order $h^{s-1}$, $1 \le s < \frac{5}{2}$, upper bound for the discretization error if the solution is $H^s$-regular. Replacing the linear Lagrange multiplier space by quadratic Lagrange multipliers does not give a higher order: see, e.g., Hild and Laborde (2002). Revising the proof of Theorem 4.9 shows that the crucial steps are Lemma 4.5 and Lemma 4.8. These parts do not yield estimates of order two even if the spaces used have higher-order best approximation properties. For a proof and more detailed numerical results, comparing quadratic finite elements with linear Lagrange multipliers and with quadratic Lagrange multipliers, we refer to Hüeber, Mair and Wohlmuth (2005a). As a test example, we choose the simple Hertz contact problem without friction.

Table 4.3 shows the convergence rates for the relative $H^1$-norm of the error in the displacement, whereas Table 4.4 refers to the error of the Lagrange multiplier in the mesh-dependent $L^2$-norm. Here, we illustrate the influence of linear and quadratic finite elements. The indices $i = 1$ and $i = 2$ stand for the use of standard conforming linear and quadratic finite elements for the displacement, respectively. The indices $j = 1$ and $j = 2$ indicate the
use of biorthogonal basis function of lowest and second order, respectively. We note that the pairing $(i, j) = (1, 2)$ is not uniformly inf-sup stable with respect to the slave side and thus is not considered. Moreover $j = 1$ already gives a best approximation property of the Lagrange multiplier space of order $3/2$, and thus we do not expect a qualitative increase in the case $j = 2$ compared to $j = 1$ for the Lagrange multiplier. Although the convergence order in the $H^1$-norm for $i = 2$ is not equal to two, it is much higher than compared to $i = 1$. A more efficient strategy, however, is to combine higher-order elements in the interior with adaptive refinement techniques on the contact part. We refer to the recent hp-strategy for a simplified Tresca problem in 2D (Dörsek and Melenk 2010).

The last test shows that the numerical solution is quite insensitive to the choice of the dual Lagrange multiplier basis. We test the discontinuous piecewise constant and linear one (see Figure 3.6(a,b)) and the continuous piecewise cubic one (see Figure 3.7(a)). Table 4.5 shows a comparison of the maximum contact pressure for different Lagrange multipliers. From the very first levels, the maximal value for all the three tested Lagrange multipliers is in very good agreement. Thus the choice of the Lagrange multiplier basis is not relevant as long as the conditions (3.2)–(3.5) are satisfied.

<table><caption>Table 4.5. Maximum contact pressure for different low-order dual Lagrange multipliers.</caption><thead><tr><th>Level</th><th>Linear</th><th>Constant</th><th>Cubic</th></tr></thead><tbody><tr><td>1</td><td>382.057</td><td>382.057</td><td>382.057</td></tr><tr><td>2</td><td>514.166</td><td>514.172</td><td>514.172</td></tr><tr><td>3</td><td>504.190</td><td>504.229</td><td>504.229</td></tr><tr><td>4</td><td>496.765</td><td>496.755</td><td>496.755</td></tr><tr><td>5</td><td>494.805</td><td>494.809</td><td>494.809</td></tr><tr><td>6</td><td>494.264</td><td>494.266</td><td>494.266</td></tr><tr><td>7</td><td>494.174</td><td>494.175</td><td>494.175</td></tr><tr><td>8</td><td>494.202</td><td>494.202</td><td>494.202</td></tr></tbody></table>

## 5. Semi-smooth Newton solver in space

Early numerical approaches for two-body contact problems on non-matching meshes and for contact problems with Coulomb friction often tried to weaken the non-linearity by suitable fixed-point or decoupling strategies. Following the proof of existence, a Coulomb friction problem can be reduced to a sequence of simplified problems with given bound for the tangential traction,
(a)

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

(b)

<figure><img src="image_35.png" /><figcaption>Figure 5.1. Section view of the problem geometry (a) and convergence rate (b).</figcaption></figure>

and thus fixed-point strategies naturally apply, where in each step a contact problem with given friction bound has to be solved numerically. Figure 5.1 shows the convergence rates of two different solvers for a Coulomb contact problem with a complex geometry in 3D. As expected, the simple fixed-point approach has a linear convergence rate, whereas an alternative solver exhibits a super-linear rate.

Using the concept of domain decomposition, a two-sided contact problem can be rewritten as a one-sided contact problem in addition to the equilibrium of contact forces. Thus Dirichlet–Neumann-type algorithms are suitable; see, e.g., Bayada, Sabil and Sassi (2002), Chernov, Geyn, Maischak and Stephan (2006), Chernov et al. (2008) and Krause and Wohlmuth (2002). More precisely, in each iteration step we solve on the master side a linear elasticity problem with given surface traction, and on the slave side we consider numerically a non-linear contact problem where the displacement of the master side acts as a rigid obstacle. The update of the interface data is realized globally after each cycle. From a theoretical point of view convergence can only be proved for a sufficiently small damping parameter (see Bayada, Sabil and Sassi (2008) and Eck and Wohlmuth (2003)), and in practice these methods require sophisticated damping strategies and are barely competitive. Figure 5.2 shows the undamped version applied to a long hexahedral bar between two cylinders. The length of the bar highly influences the convergence behaviour. If it is small, then the undamped algorithm shows fast convergence, whereas if the bar is long, no convergence at all is obtained unless a suitable damping parameter is applied.

Although simple to apply, these coupling strategies give rise to inner and outer iteration schemes and are therefore quite expensive. Thus there is a strong need for efficient solvers which focus on all non-linearities at the same time and tackle the fully coupled system.
<figure><img src="image_36.png" /><figcaption>Figure 5.2. Oscillation of an undamped Dirichlet–Neumann contact solver.</figcaption></figure>

Monotone multigrid methods have been shown to be very attractive due to guaranteed convergence, if the underlying system is equivalent to a constrained minimization problem. Early references on multigrid methods for variational inequalities or free boundary problems such as the obstacle problem are given by Brandt and Cryer (1983), Hackbusch and Mitelmann (1983), Hoppe (1987), Hoppe and Kornhuber (1994), Kornhuber (1994, 1996), Kornhuber and Krause (2001) and the monograph by Kornhuber (1997). Nowadays these techniques have been applied very successfully to more challenging contact problems including two bodies with non-matching meshes, finite deformations and complex geometries in 3D (Dickopf and Krause 2009a, Krause 2008, 2009, Krause and Mohr 2011, Wohlmuth and Krause 2003).

Domain decomposition-based solvers such as FETI techniques are also widely applied. An excellent overview of these techniques applied to variational inequalities can be found in the recent monograph by Dostál (2009); see also the original research papers of Dostál, Friedlander and Santos (1998), Dostál, Gomes Neto and Santos (2000), Dostál and Horák (2003), Dostál, Horák, Kučera, Vondrák, Haslinger, Dobiaš and Pták (2005), Dostál, Horák and Stefanica (2007, 2009) and Schöberl (1998). Different alternatives exist, e.g., interior point methods (Wright 1997), SQP algorithms (Pang and Gabriel 1993), the radial return mapping or the catching-up algorithm (Moreau 1977, Simo and Hughes 1998), as well as penalty or augmented Lagrangian approaches (Glowinski and Le Tallec 1989, Laursen 2002, Simo and Laursen 1992). We refer to the recent monograph by Ito and Kunisch (2008a) for an overview of Lagrange multiplier-based methods for variational problems.

Here we choose an abstract and very flexible framework within which many different applications can be handled. The starting point is the observation that most inequality constraints can be equivalently stated in terms of a non-linear system. This holds true not only for contact problems but also
for other problems involving variational inequalities. The weak form of the
underlying partial differential equation and the non-linear complementarity
(NCP) function then form a coupled non-linear system on which a Newton
scheme can be applied and easily combined with fast iterative solvers, such
as multigrid (Hackbusch 1985) or domain decomposition techniques (Toselli
and Widlund 2005), for the consistent linearized system. Due to the char-
acteristic lack of classical differentiability of the NCP function, the assump-
tions for standard Newton methods (Deuflhard 2004) are not satisfied, but
the so-called semi-smooth Newton methods (Facchinei and Pang 2003a, Ito
and Kunisch 2003, Hintermüller, Kovtunenko and Kunisch 2004) can be ap-
plied; see also Pang (1990) and Pang and Qi (1993). Early applications of
this type of method can be found in the engineering literature. The classical
radial return mapping in plasticity (see Moreau (1977) for an early variant
of it) can be handled within this abstract framework, but it has also been
successfully applied to contact problems for roughly two decades (Alart and
Curnier 1991, Christensen 2002a, 2002b, Christensen and Pang 1999, De
Saxcé and Feng 1991, Simo and Laursen 1992). It is well established that
the semi-smooth Newton method converges locally super-linearly: see, e.g.,
the monograph by Facchinei and Pang (2003b). Global convergence can
be shown for some special cases, e.g., the Laplace operator-based obstacle
problem: see Ito and Kunisch (2008a). A simplified Signorini problem has
been analysed in Ito and Kunisch (2008b). For contact problems no global
convergence holds, but the pre-asymptotic robustness can, in particular in
3D, be widely improved by a suitable local rescaling of an NCP function and
a local node-wise regularization of the Jacobian. In each Newton step, the
contact condition and its boundary type have to be updated locally, e.g., a
Robin-type condition applies in the case of a sliding node. As a consequence,
the semi-smooth Newton method can be implemented as a primal-dual ac-
tive set strategy (Hintermüller, Ito and Kunisch 2002, Hüeber, Stadler and
Wohlmuth 2008, Ito and Kunisch 2004). The use of active sets allows for
local static condensation of either the dual variable or the corresponding
primal degrees of freedom, such that only a system of the size of the dis-
placement has to be solved in each Newton step. One of the attractive
features of this class of algorithms is that it can be easily combined with
other types of non-linearities, such as non-linear material laws, for example.
No inner and outer iteration loop is required even in the presence of the
different types of non-linearities.

5.1. Equivalent formulation as a non-linear equation system

In a first step, we rewrite the inequality constraints associated with the
discrete static Coulomb friction problem as a non-linear system. For sim-
plicity of notation, we present the algebraic form only for homogeneous
Dirichlet boundary conditions. Thus, after discretization the weak formula-
tion (3.9) has the following algebraic structure. Find $(\tilde{\mathbf{u}}_l, \boldsymbol{\lambda}_l) \in \mathbb{R}^{N_l^V} \times \mathbb{R}^{N_l^M}$,
$N_l^V := \dim \mathbf{V}_l$, $N_l^M := \dim \mathbf{M}_l$, such that

$$
\begin{align}
\tilde{A}_l \tilde{\mathbf{u}}_l + \tilde{B}_l \boldsymbol{\lambda}_l &= \tilde{\mathbf{f}}_l, \tag{5.1} \\
\tilde{C}_l(\boldsymbol{\lambda}_l, \tilde{\mathbf{u}}_l) &= \mathbf{0}. \nonumber
\end{align}
$$

Here we use the same symbol for $\boldsymbol{\lambda}_l \in \mathbf{M}_l$ and its vector representation
$\boldsymbol{\lambda}_l \in \mathbb{R}^{N_l^M}$. The matrices $\tilde{A}_l \in \mathbb{R}^{N_l^V} \times N_l^V$, $\tilde{B}_l \in \mathbb{R}^{N_l^V} \times N_l^M$ and the right-
hand side $\tilde{\mathbf{f}}_l \in \mathbb{R}^{N_l^V}$ result from the bilinear forms $a_l(\cdot, \cdot)$, $b_l(\cdot, \cdot)$ and the
linear form $f_l(\cdot)$, respectively, and are assembled with respect to the nodal
basis functions $\phi_p, \psi_p$. We note that $d\#\mathcal{P}_{C;l}^s := N_l^s = N_l^M$. The NCP
function $\tilde{C}_l(\cdot, \cdot) \in \mathbb{R}^{N_l^M}$ reflects the non-penetration condition (3.11) and
the static Coulomb law (3.12). It has a node-wise form and can be written
as $(\tilde{C}_l(\boldsymbol{\lambda}_l, \tilde{\mathbf{u}}_l))_p = \tilde{C}_p(\gamma_p, \tilde{\mathbf{u}}_l) \in \mathbb{R}^d, p \in \mathcal{P}_{C;l}^s$.

There exist many different choices for NCP functions in the literature. Quite often generalizations of the Fischer-Burmeister approach (Fischer 1992) are used. We refer to Chen, Chen and Kanzow (2000) for a penalized version and to Chen (2007), Hu, Huang and Chen (2009), Kanzow, Yamashita and Fukushima (1997) and Sun and Qi (1999) for the introduction and analysis of a family of NCP functions. An excellent overview can be found in the monograph by Facchinei and Pang (2003a). Here we use a different type of NCP function, which is based not on the root function but on the max function: see Alart and Curnier (1991). The main advantage of this type of NCP function is that the generalized derivatives are extremely easy to compute. This is quite important if not only the simple non-penetration law is considered but also more complex situations such as Coulomb friction or finite deformations. We refer to the series of recent papers by Gitterle et al. (2010) and Popp et al. (2009, 2010), where the concept of dual La-grange multipliers and semi-smooth Newton schemes have been applied to finite deformation problems. Moreover, it can be implemented easily in terms of an active set strategy which switches off and on different types of non-linear boundary conditions on the possible contact zone.

As is standard in the case of radial return mappings, we introduce trial test vectors. Here we need two of them, one in the normal and one in the tangential direction, that is,

$$
\gamma_{p;\text{tr}}^n := \gamma_p^n + \frac{2\mu c_n}{m_p} (\alpha_p^n - g_p), \quad \gamma_{p;\text{tr}}^t := \gamma_p^t + \frac{2\mu c_t}{m_p} \alpha_p^t, \qquad (5.2)
$$

where $c_n$ and $c_t$ are two positive mesh-independent constants. Keeping
in mind that $\boldsymbol{\alpha}_p, p \in \mathcal{P}_{C;l}^s$, are the coefficients with respect to the nodal
basis functions on the slave side of $\mathbf{\Pi}_l[\mathbf{u}_l]$, we find that $\boldsymbol{\alpha}_p$ is a function
of $\tilde{u}_l$. In terms of these trial vectors, which depend on $\gamma_p$ and $\tilde{u}_l$, one can easily reformulate the inequality constraints (3.11) and (3.12) as equality conditions.

**Lemma 5.1.** The inequality constraints (3.11) and (3.12) for each node $p \in \mathcal{P}_{C;l}^s$ are equivalent to $\tilde{C}_p(\gamma_p, \tilde{u}_l) := C_p(\gamma_p, \alpha_p(\tilde{u}_l)) = \mathbf{0}$, where the normal component $C_p^n(\gamma_p^n, \alpha_p^n)$ of $C_p(\gamma_p, \alpha_p)$ is given by

$$C_p^n(\gamma_p^n, \alpha_p^n) := \gamma_p^n - \max(0, \gamma_{p;\text{tr}}^n), \qquad (5.3)$$

and the tangential component of $C_p(\gamma_p, \alpha_p)$ is defined by $C_p^t(\gamma_p, \alpha_p) :=$

$$\left\{ \begin{array}{ll} \gamma_p^t & \mathcal{F}_{p;\text{tr}} = 0, \\ \max(\|\gamma_{p;\text{tr}}^t\|, \mathcal{F}_{p;\text{tr}})^s \left( \gamma_p^t - \min \left( 1, \frac{\mathcal{F}_{p;\text{tr}}}{\|\gamma_{p;\text{tr}}^t\|} \right) \gamma_{p;\text{tr}}^t \right) & \text{otherwise,} \end{array} \right. \qquad (5.4)$$

where $s \ge 0$ stands for a scaling parameter and $\mathcal{F}_{p;\text{tr}} := \nu \max(0, \gamma_{p;\text{tr}}^n)$.

**Proof.** The proof is a rather straightforward calculation, but for convenience of the reader it is given. If $C_p^n(\gamma_p^n, \alpha_p^n) = 0$, then obviously $\gamma_p^n \ge 0$. If $\gamma_p^n = 0$, we get $\alpha_p^n - g_p \le 0$, and for $\gamma_p^n > 0$, we have $\alpha_p^n = g_p$, and thus (3.11) is satisfied. Let $C_p^n(\gamma_p^n, \alpha_p^n) = 0$ and $C_p^t(\gamma_p, \alpha_p) = \mathbf{0}$; then $\mathcal{F}_{p;\text{tr}} = \nu\gamma_p^n$. In the case $\mathcal{F}_{p;\text{tr}} = 0$, (3.12) trivially holds. For $\mathcal{F}_{p;\text{tr}} > 0$, the scaling factor $\max(\|\gamma_{p;\text{tr}}^t\|, \mathcal{F}_{p;\text{tr}})^s$ is non-zero, and thus from $C_p^t(\gamma_p, \alpha_p) = \mathbf{0}$ it follows that

$$\gamma_p^t = \min \left( 1, \frac{\mathcal{F}_{p;\text{tr}}}{\|\gamma_{p;\text{tr}}^t\|} \right) \gamma_{p;\text{tr}}^t.$$

If $\|\gamma_{p;\text{tr}}^t\| \le \nu\gamma_p^n$, then $\gamma_p^t - \gamma_{p;\text{tr}}^t = \mathbf{0}$, and thus $\alpha_p^t = \mathbf{0}$ and $\|\gamma_p^t\| \le \nu\gamma_p^n$. For $\|\gamma_{p;\text{tr}}^t\| > \nu\gamma_p^n$, we get $\|\gamma_p^t\| = \nu\gamma_p^n$. Moreover, $\gamma_{p;\text{tr}}^t = (1 + \beta)\gamma_p^t$ with some positive value for $\beta$. Using the definition (5.2) of $\gamma_{p;\text{tr}}^t$, we find that $\alpha_p^t$ points in the same direction as $\gamma_p^t$ and thus (3.12) is satisfied.

Let (3.11) be true; then a straightforward computation shows that we have $C_p^n(\gamma_p^n, \alpha_p^n) = 0$. The situation is more complex for the discrete Coulomb law. If (3.11) and (3.12) hold, then either $\alpha_p^t = \mathbf{0}$ or it points in the same direction as $\gamma_p^t$ or $\gamma_p^n = 0$. In the first two cases, we can thus write $\gamma_{p;\text{tr}}^t = (1 + \beta)\gamma_p^t$ with some non-negative $\beta$. For $\beta = 0$, we trivially find that $\gamma_{p;\text{tr}}^t = \gamma_p^t$ and thus $C_p^t(\gamma_p, \alpha_p) = \mathbf{0}$. For $\beta > 0$, we have $\|\gamma_p^t\| = \nu\gamma_p^n$ and moreover $\|\gamma_{p;\text{tr}}^t\| > \nu\gamma_p^n$, yielding

$$\gamma_p^t - \min \left( 1, \frac{\mathcal{F}_{p;\text{tr}}}{\|\gamma_{p;\text{tr}}^t\|} \right) \gamma_{p;\text{tr}}^t = \gamma_p^t - \frac{1}{1 + \beta} \gamma_{p;\text{tr}}^t = \gamma_p^t - \frac{1 + \beta}{1 + \beta} \gamma_p^t = \mathbf{0}.$$

The case that $\gamma_p^n = 0$ yields $\mathcal{F}_{p;\text{tr}} = 0$ and $\gamma_p^t = \mathbf{0}$ and thus $C_p^t(\gamma_p, \alpha_p) = \mathbf{0}$.

<figure><img src="image_37.png" /></figure>
**Remark 5.2.** We note that (3.11) and (5.3) are equivalent but (3.12) and (5.4) are not. In particular, the trial vector in the normal direction enters into the definition of $C_p^t(\gamma_p, \alpha_p)$.

## 5.2. Basis transformation: from nodal to constrained

As already mentioned, from the algebraic point of view each two-body problem can be rewritten formally as a one-body problem by introducing a new basis. The weak inequalities then result in node-wise inequalities for the slave side, even in the case of non-matching meshes. One of the main advantages of the choice of a dual Lagrange multiplier space is that the basis transformation is a local operator with a sparse matrix representation. This is not the case for standard Lagrange multipliers.

Although $\tilde{C}_p(\gamma_p, \tilde{u}_l)$ only depends on the nodal coefficient of $\lambda_l$ at $p$, it does not have this simple structure with respect to $\tilde{u}_l$. A more local and thus implementationally attractive representation can be obtained by a suitable basis transformation. This transformation was introduced in Wohlmuth and Krause (2001) to construct a multigrid scheme for mortar finite element discretizations. The V-cycle analysis of a level-independent convergence rate can be found in Wohlmuth (2005).

Introducing the nodal block structure for the displacement

$$ \tilde{u}_l = (\tilde{u}_l^i, \tilde{u}_l^m, \tilde{u}_l^s)^\top, \quad \tilde{u}_l^i \in \mathbb{R}^{N_l^i}, \quad \tilde{u}_l^m \in \mathbb{R}^{N_l^m}, \quad \tilde{u}_l^s \in \mathbb{R}^{N_l^s}, $$

where $N_l^m := d\#\mathcal{P}_{C;l}^m$ and $N_l^i := N_l^V - N_l^m - N_l^s$ (see also Figure 5.3), we obtain the following structure for $\tilde{A}_l$ and $\tilde{B}_l$:

$$ \tilde{A}_l = \begin{pmatrix} \tilde{A}_l^{i,i} & \tilde{A}_l^{i,m} & \tilde{A}_l^{i,s} \\ \tilde{A}_l^{m,i} & \tilde{A}_l^{m,m} & 0 \\ \tilde{A}_l^{s,i} & 0 & \tilde{A}_l^{s,s} \end{pmatrix}, \quad \tilde{B}_l = \begin{pmatrix} 0 \\ -M_l^\top \\ D_l \end{pmatrix}, $$

where the entries of the coupling matrices $D_l$ and $M_l$ are defined by

$$ d_{pp} := m_p \mathrm{Id}_{d \times d}, \quad m_{pq} := \int_{\Gamma_{C;l}^s} \psi_p \phi_q \circ \chi_l \, ds \, \mathrm{Id}_{d \times d}, \quad p \in \mathcal{P}_{C;l}^s, q \in \mathcal{P}_{C;l}^m. $$

The diagonal structure of $D_l$ is a consequence of (3.6). In the case of standard Lagrange multipliers $D_l$ has the band-structure of a $(d-1)$-dimensional mass matrix, and thus $D_l^{-1}$ is dense.

We recall that, with respect to the nodal basis functions, the coefficient $\alpha_p$ in (3.11) and (3.12) does depend on the coefficients of the master and the slave nodes. The definition (3.10) of the mortar projection yields that

$$ (\alpha_p)_{p \in \mathcal{P}_{C;l}^s} := \alpha = \tilde{u}_l^s - D_l^{-1} M_l \tilde{u}_l^m. $$

To eliminate the dependence of the coefficient vector on the master side, we
<figure><img src="image_38.png" /><figcaption>Figure 5.3. Partitioning of the nodes into the three blocks.</figcaption></figure>

use a basis transformation. The nodal basis function $\phi_p$, $p \in \mathcal{P}_{C;l}^m$, will be replaced by a constrained basis function

$$ \tilde{\phi}_p := \phi_p + \sum_{q \in \mathcal{P}_{C;l}^s} r_{pq} \phi_q, \qquad (5.5) $$

while all other basis functions remain unchanged, i.e., $\tilde{\phi}_p := \phi_p$, $p \in \mathcal{P}_l \setminus \mathcal{P}_{C;l}^m$.

**Lemma 5.3.** Setting $r_{pq} := ((D_l^{-1}M_l)_{qp})_{11}$, $q \in \mathcal{P}_{C;l}^s$, $p \in \mathcal{P}_{C;l}^m$ in (5.5) gives

$$ \Pi_l[\tilde{\phi}_p] = 0, \quad p \in \mathcal{P}_{C;l}^m. $$

*Proof.* The definition of the mortar projection (3.10a) shows that $\Pi_l$, restricted to $W_l^s$, is the identity and that

$$ \Pi_l(\phi_p \circ \chi_l) = \sum_{q \in \mathcal{P}_{C;l}^s} r_{pq} \phi_q, \quad p \in \mathcal{P}_{C;l}^m $$

and thus

$$ \Pi_l[\tilde{\phi}_p] = \sum_{q \in \mathcal{P}_{C;l}^s} r_{pq} \phi_q - \Pi_l(\phi_p \circ \chi_l) = 0. \quad \square $$

Algebraically this basis transformation can be realized very efficiently in a local pre-process. The coefficients with respect to the constrained basis multiplied by the matrix

$$ Q_l = \begin{pmatrix} \mathrm{Id} & 0 & 0 \\ 0 & \mathrm{Id} & 0 \\ 0 & D_l^{-1} M_l & \mathrm{Id} \end{pmatrix} \qquad (5.6) $$

yield the coefficients with respect to the nodal basis. If $v_l = \sum_{p \in \mathcal{P}_l} \tilde{\beta}_p \phi_p$, then $v_l = \sum_{p \in \mathcal{P}_l} \beta_p \tilde{\phi}_p$ with $(\tilde{\beta}_p)_{p \in \mathcal{P}_l} = Q_l(\beta_p)_{p \in \mathcal{P}_l}$. To obtain the stiffness matrices $A_l$, $B_l$ and the right-hand side $f_l$ with respect to the new constrained basis, we have to apply $Q_l$ in a suitable way on $\tilde{A}_l$, $\tilde{B}_l$ and $\tilde{f}_l$:

$$ A_l = Q_l^\top \tilde{A}_l Q_l, \quad B_l = Q_l^\top \tilde{B}_l, \quad f_l = Q_l^\top \tilde{f}_l. $$
Due to the definition (5.6) of $Q_l$, it is easy to verify that the block structure of the matrix $B_l$, with respect to this new basis, has the form $B_l = (0, 0, D_l)^\top$. The non-linear system to be solved is equivalent to (5.1), and can be written as

$$ 
\begin{aligned}
A_l \mathbf{u}_l + B_l \boldsymbol{\lambda}_l &= \mathbf{f}_l, \\
C_p(\gamma_p, \mathbf{u}_p) &= \mathbf{0}, \quad p \in \mathcal{P}_{C;l}^s.
\end{aligned} 
\qquad (5.7) $$

Here we have used that, with respect to the new basis, $\alpha_p = \mathbf{u}_p$, $p \in \mathcal{P}_{C;l}^s$. We point out that the same notation for $\mathbf{u}_l \in \mathbf{V}_l$ is applied as for its coefficient vector $\mathbf{u}_l \in \mathbb{R}^{N_l^V}$ with respect to the new constrained basis. Finally, we mention that, with respect to this new basis, the $p$-component of the NCP function only depends on the $p$-component of $\boldsymbol{\lambda}_l$ and $\mathbf{u}_l$. The coefficients in the new basis with respect to the slave nodes no longer specify the total nodal displacement but now describe the movement of the underlying finite element node relative to the master side. Having $\mathbf{u}_l$, the Lagrange multiplier can easily be obtained, from the node-wise residual, by a diagonal scaling

$$ \boldsymbol{\lambda}_l = D_l^{-1}(\mathbf{f}_l^s - A_l^{\text{si}} \mathbf{u}_l^{\text{si}} - A_l^{\text{sm}} \mathbf{u}_l^{\text{m}} - A_l^{\text{ss}} \mathbf{u}_l^{\text{s}}). \qquad (5.8) $$

**Remark 5.4.** Such a basis transformation can be carried out for all types of Lagrange multiplier spaces as long as a discrete inf-sup condition holds. In most cases, it will result in a dense block corresponding to the possibly global character of the mortar projection. However, in our situation the basis transformation is a local operator for all space dimensions and thus inexpensive. Moreover, for uniformly stable pairs $(\mathbf{V}_l, \mathbf{M}_l)$, the basis transformation does not influence the order of the condition number of $A_l$, which is then comparable to that of $\tilde{A}_l$.

### 5.3. Semi-smooth Newton solver

As can be easily seen from (5.3) and (5.4), the NCP function $C_p(\gamma_p, \mathbf{u}_p)$ is not globally but only piecewise smooth. Thus, to solve (5.7), a Newton-type solver can be implemented in terms of an active set strategy. The active and inactive sets are defined by the different cases of $C_p(\cdot, \cdot)$ and can be selected node-wise.

Let $(\boldsymbol{\lambda}_l^{k-1}, \mathbf{u}_l^{k-1})$ be the previous iterate. We then obtain the new iterate

$$ (\boldsymbol{\lambda}_l^k, \mathbf{u}_l^k) = (\boldsymbol{\lambda}_l^{k-1}, \mathbf{u}_l^{k-1}) + (\delta\boldsymbol{\lambda}_l^{k-1}, \delta\mathbf{u}_l^{k-1}) $$

of the semi-smooth Newton step by solving a linear system for the update $(\delta\boldsymbol{\lambda}_l^{k-1}, \delta\mathbf{u}_l^{k-1})$. For each node $p \in \mathcal{P}_{C;l}^s$, the local system

$$ DC_p(\gamma_p^{k-1}, \mathbf{u}_p^{k-1}) (\delta\gamma_p^{k-1}, \delta\mathbf{u}_p^{k-1})^\top = -C_p(\gamma_p^{k-1}, \mathbf{u}_p^{k-1}) \qquad (5.9) $$

has to be satisfied. Here, $DC_p$ denotes the Jacobian of the local NCP function $C_p(\cdot, \cdot)$. We note that due to our basis transformation, $DC_p$ can
be regarded as a $d \times 2d$ matrix. Defining the trial vectors of the increment $(\delta\gamma_p, \delta\mathbf{u}_p)$, in a way similar to (5.2), by

$$ \delta\gamma_{p;\text{tr}}^n := \delta\gamma_p^n + \frac{2\mu c_n}{m_p} \delta u_p^n, \quad \delta\gamma_{p;\text{tr}}^t := \delta\gamma_p^t + \frac{2\mu c_t}{m_p} \delta \mathbf{u}_p^t, $$

we introduce three sets in terms of $S_n(\gamma_p, \mathbf{u}_p) := S_n := \max(0, \gamma_{p;\text{tr}}^n)$:

$$ \begin{aligned} \mathcal{I}_n &:= \mathcal{I}_n(\gamma_p, \mathbf{u}_p) := \{p \in \mathcal{P}_{C;l}^s : S_n = 0\}, \\ \mathcal{I}_t &:= \mathcal{I}_t(\gamma_p, \mathbf{u}_p) := \{p \in \mathcal{P}_{C;l}^s : \|\gamma_{p;\text{tr}}^t\| < \nu S_n\}, \\ \mathcal{A} &:= \mathcal{A}(\gamma_p, \mathbf{u}_p) := \{p \in \mathcal{P}_{C;l}^s : \|\gamma_{p;\text{tr}}^t\| \ge \nu S_n, S_n > 0\}. \end{aligned} \quad (5.10) $$

Now, we can easily define the generalized derivative $DC_p(\cdot, \cdot)$ for each $p \in \mathcal{P}_{C;l}^s$ in the direction of the update $(\delta\gamma_p, \delta\mathbf{u}_p)$. We observe that $\mathcal{I}_n, \mathcal{I}_t$ and $\mathcal{A}$ form a non-overlapping partition of $\mathcal{P}_{C;l}^s$. On $\mathcal{I}_n$ the nodes are free and not in contact and no boundary forces apply. On $\mathcal{I}_t \cup \mathcal{A}$ the nodes are actually in contact and sliding on $\mathcal{A}$.

Before we discuss the general case, we focus on $\nu = 0$. Here $\mathcal{I}_t = \emptyset$, and the Newton algorithm simplifies considerably. Observing that the tangential component of $C_p(\cdot, \cdot)$ is then linear, we get

$$ DC_p(\gamma_p, \mathbf{u}_p)(\delta\gamma_p, \delta\mathbf{u}_p) = \begin{cases} (\delta\gamma_p^n, \delta\gamma_p^t)^\top & p \in \mathcal{I}_n, \\ (\delta\gamma_p^n - \delta\gamma_{p;\text{tr}}^n, \delta\gamma_p^t)^\top & p \in \mathcal{A}. \end{cases} $$

The Newton update (5.9) then gives $(\gamma_p)^k = \mathbf{0}$ for $p \in \mathcal{I}_n^k := \mathcal{I}_n(\gamma_p^{k-1}, \mathbf{u}_p^{k-1})$ and $(\gamma_p^t)^k = \mathbf{0}$, $(u_p^n)^k = g_p$ for $p \in \mathcal{A}^k := \mathcal{A}(\gamma_p^{k-1}, \mathbf{u}_p^{k-1})$. For all possible contact nodes, we have a homogeneous Neumann boundary condition in the tangential direction, whereas in the normal direction we have a Dirichlet or Neumann condition. Thus, the implementation can be easily realized as a primal-dual active set strategy.

The situation is more involved for $\nu > 0$. For a Coulomb problem with a non-trivial friction coefficient, we obtain three different situations for $C_p := C_p(\gamma_p, \mathbf{u}_p)$ and $DC_p := DC_p(\gamma_p, \mathbf{u}_p)(\delta\gamma_p, \delta\mathbf{u}_p)$. Using the differential

$$ \partial \frac{\mathbf{x}}{\|\mathbf{x}\|} = \frac{1}{\|x\|} \left( \mathrm{Id} - \frac{\mathbf{x} \otimes \mathbf{x}}{\|\mathbf{x}\|^2} \right), \quad \mathbf{x} \neq \mathbf{0}, $$

a straightforward computation shows the following.

* If $p \in \mathcal{I}_n$, then $C_p = (\gamma_p^n, \gamma_p^t)^\top$ and

$$ DC_p = \begin{pmatrix} \delta\gamma_p^n \\ \delta\gamma_p^t \end{pmatrix}. $$
* If $p \in \mathcal{I}_t$, then $C_p = (\gamma_p^n - \gamma_{p;tr}^n, (\nu\gamma_{p;tr}^n)^s(\gamma_p^t - \gamma_{p;tr}^t))^\top$,

$$DC_p = \left( \frac{\delta\gamma_p^n - \delta\gamma_{p;tr}^n}{(\nu\gamma_{p;tr}^n)^s(\delta\gamma_p^t - \delta\gamma_{p;tr}^t + \frac{s}{\gamma_{p;tr}^n}\delta\gamma_{p;tr}^n(\gamma_p^t - \gamma_{p;tr}^t))} \right).$$

* If $p \in \mathcal{A}$, then $C_p = (\gamma_p^n - \gamma_{p;tr}^n, \|\gamma_{p;tr}^t\|^s(\gamma_p^t + \alpha_2\gamma_{p;tr}^t))^\top$,

$$DC_p = \left( \frac{\delta\gamma_p^n - \delta\gamma_{p;tr}^n}{\|\gamma_{p;tr}^t\|^s(\delta\gamma_p^t + \alpha_1\gamma_p^t + \alpha_2\delta\gamma_{p;tr}^t + \alpha_3\gamma_{p;tr}^t)} \right),$$

where the factors $\alpha_1$, $\alpha_2$ and $\alpha_3$ are given by

$$\alpha_1 := s \frac{\gamma_{p;tr}^t \delta \gamma_{p;tr}^t}{\|\gamma_{p;tr}^t\|^2}, \quad \alpha_2 := - \frac{\nu \gamma_{p;tr}^n}{\|\gamma_{p;tr}^t\|},$$

$$\alpha_3 := (1-s) \frac{\nu \gamma_{p;tr}^n}{\|\gamma_{p;tr}^t\|^3} (\delta \gamma_{p;tr}^t \gamma_{p;tr}^t) - \frac{\nu \delta \gamma_{p;tr}^n}{\|\gamma_{p;tr}^t\|}.$$

For all cases, the normal component of the NCP function is linear; this is not the case for the tangential component. Only for $p \in \mathcal{I}_n$ do we have a linear tangential component. Here, we restrict ourselves to $s \in [0, 1]$ and note that for $s = 0$ or $s = 1$ some of the terms cancel.

**Lemma 5.5.** The semi-smooth Newton solver applied to (5.7) can be implemented as a primal-dual active set strategy, where in each Newton step we have to update for each node $p \in \mathcal{P}_{C;l}^s$ the type and the value of the boundary condition. Moreover, for $\nu > 0$ we have the following.

* Homogeneous Neumann conditions for $p \in \mathcal{I}_n^k := \mathcal{I}_n(\gamma_p^{k-1}, \mathbf{u}_p^{k-1})$:

$$(\gamma_p)^k = \mathbf{0} \qquad (5.11)$$

* Inhomogeneous Dirichlet conditions for $p \in \mathcal{I}_t^k := \mathcal{I}_t(\gamma_p^{k-1}, \mathbf{u}_p^{k-1})$:

$$(u_p^n)^k = g_p \quad \text{and} \quad (\mathbf{u}_p^t)^k = \frac{-s(\delta\gamma_{p;tr}^n)^{k-1}}{(\gamma_{p;tr}^n)^{k-1}} (\mathbf{u}_p^t)^{k-1}, \qquad (5.12)$$

where the condition in the tangential direction also depends on the update of the normal surface traction.

* Dirichlet conditions in the normal and Robin conditions in the tangential direction for $p \in \mathcal{A}^k := \mathcal{A}(\gamma_p^{k-1}, \mathbf{u}_p^{k-1})$:

$$(u_p^n)^k = g_p \quad \text{and} \qquad (5.13a)$$

$$(\mathrm{Id} + L_{p;s}^{k-1})(\gamma_p^t)^k + \hat{c}_t L_{p;s}^{k-1}(\mathbf{u}_p^t)^k - \frac{\nu(\gamma_p^n)^k(\gamma_{p;tr}^t)^{k-1}}{\|(\gamma_{p;tr}^t)^{k-1}\|} = \mathbf{g}_p^{k-1}, \quad (5.13b)$$

where $\mathbf{g}_p^{k-1} := L_{p;s}^{k-1}(\gamma_{p;tr}^t)^{k-1}$, the mesh-dependent scaling factor $\hat{c}_t$ is
given by $\hat{c}_t := \frac{2\mu c_t}{m_p}$ and the matrix $L_{p;s}^{k-1}$ by

$$L_{p;s}^{k-1} := \beta_p^{k-1} \left( \alpha_2^{k-1} \mathrm{Id} + \frac{sL_1^{k-1} + (s-1)\alpha_2^{k-1}L_2^{k-1}}{\|(\gamma_{p;\mathrm{tr}}^t)^{k-1}\|^2} \right),$$

$$\alpha_2^{k-1} := - \frac{\nu(\gamma_{p;\mathrm{tr}}^n)^{k-1}}{\|(\gamma_{p;\mathrm{tr}}^t)^{k-1}\|}, \quad \beta_p^{k-1} := 1, \quad \omega_p^{k-1} := 1,$$

$$L_1^{k-1} := \omega_p^{k-1} (\gamma_p^t)^{k-1} \otimes (\gamma_{p;\mathrm{tr}}^t)^{k-1}, \quad L_2^{k-1} := (\gamma_{p;\mathrm{tr}}^t)^{k-1} \otimes (\gamma_{p;\mathrm{tr}}^t)^{k-1}.$$

*Proof.* The proof is based on (5.9) and the partitioning (5.10). For each node $p \in \mathcal{I}_n^k$, we get

$$\binom{(\delta\gamma_p^n)^{k-1}}{(\delta\gamma_p^t)^{k-1}} = - \binom{(\gamma_p^n)^{k-1}}{(\gamma_p^t)^{k-1}}$$

and thus (5.11). From now on, to simplify the notation we suppress the upper index $k-1$ of the Newton iteration.

If $p \in \mathcal{I}_t^k$, then

$$\left( \left( \delta\gamma_p^n - \delta\gamma_{p;\mathrm{tr}}^n \right) + \frac{s}{\gamma_{p;\mathrm{tr}}^n} \delta\gamma_{p;\mathrm{tr}}^n (\gamma_p^t - \gamma_{p;\mathrm{tr}}^t) \right) = - \left( \begin{pmatrix} \gamma_p^n - \gamma_{p;\mathrm{tr}}^n \\ (\gamma_p^t - \gamma_{p;\mathrm{tr}}^t) \end{pmatrix} \right).$$

Observing that

$$\gamma_p^n - \gamma_{p;\mathrm{tr}}^n = \frac{2\mu c_n}{m_p} (g_p - u_p^n) \quad \text{and} \quad \delta\gamma_p^n - \delta\gamma_{p;\mathrm{tr}}^n = -\frac{2\mu c_n}{m_p} \delta u_p^n,$$

we get the normal part of (5.12). For the tangential part, we use

$$\gamma_p^t - \gamma_{p;\mathrm{tr}}^t = -\frac{2\mu c_n}{m_p} \mathbf{u}_p^t \quad \text{and} \quad \delta\gamma_p^t - \delta\gamma_{p;\mathrm{tr}}^t = -\frac{2\mu c_n}{m_p} \delta\mathbf{u}_p^t.$$

For $p \in \mathcal{A}^k$, we find the same condition as in (5.12) for the normal part and thus (5.13a). Using the tangential component of the system (5.9), we find

$$\partial\gamma_p^t + \alpha_1\gamma_p^t + \alpha_2\delta\gamma_{p;\mathrm{tr}}^t + \alpha_3\gamma_{p;\mathrm{tr}}^t = -\gamma_p^t - \alpha_2\gamma_{p;\mathrm{tr}}^t. \quad (5.14)$$

Unfortunately, the new unknown updates are hidden in the coefficients $\alpha_1$ and $\alpha_3$. Using the definition of the matrix $L_{p;s}$ and applying the formula $(\mathbf{x} \otimes \mathbf{y})\mathbf{z} = (\mathbf{y}\mathbf{z})\mathbf{x}$, we get

$$\begin{align*}
\alpha_1 \gamma_p^t + \alpha_3 \gamma_{p;\mathrm{tr}}^t &= (L_{p;s} - \alpha_2 \mathrm{Id}) \delta \gamma_{p;\mathrm{tr}}^t - \frac{\nu \delta \gamma_{p;\mathrm{tr}}^n}{\|\gamma_{p;\mathrm{tr}}^t\|} \gamma_{p;\mathrm{tr}}^t \\
&= (L_{p;s} - \alpha_2 \mathrm{Id}) \delta \gamma_{p;\mathrm{tr}}^t - \frac{\nu (\gamma_{p;\mathrm{tr}}^n)^k}{\|\gamma_{p;\mathrm{tr}}^t\|} \gamma_{p;\mathrm{tr}}^t - \alpha_2 \gamma_{p;\mathrm{tr}}^t.
\end{align*}$$

Using (5.13a), which guarantees that $(\gamma_{p;\mathrm{tr}}^n)^k = (\gamma_p^n)^k$, and inserting the equality for $\alpha_1\gamma_p^t + \alpha_3\gamma_{p;\mathrm{tr}}^t$ in (5.14), we get (5.13b). $\square$
**Remark 5.6.** If a non-local friction law or a simple combination of Coulomb friction and Tresca friction is applied, possibly more cases have to be considered, we refer to Hager and Wohlmuth (2009b) for an abstract framework. If the friction bound depends not only on the contact pressure, then nodes which are not in contact are not automatically surface traction-free nodes.

Table 5.1 shows the convergence rate of the semi-smooth Newton scheme with $s = 1$ applied to the 2D contact problem with given friction bound depicted in Figure 4.3. Here we have used a discrete mesh-dependent norm for the Lagrange multiplier. The third column shows the ratio of the error in step $k$ and $k-1$. As expected, from the theoretical point of view, this ratio tends to zero.

<table><caption>Table 5.1. Super-linear convergence for the semi-smooth Newton scheme applied to a contact problem with given friction bound.</caption><thead><tr><th>k</th><th>e<sup>k</sup> := ||λ<sub>l</sub><sup>k</sup> - λ<sub>l</sub>||</th><th>e<sup>k</sup>/e<sup>k-1</sup></th></tr></thead><tbody><tr><td>1</td><td>5.622e+02</td><td>—</td></tr><tr><td>2</td><td>2.553e+02</td><td>4.541e-01</td></tr><tr><td>3</td><td>1.087e+02</td><td>4.257e-01</td></tr><tr><td>4</td><td>3.515e+01</td><td>3.233e-01</td></tr><tr><td>5</td><td>5.761e+00</td><td>1.638e-01</td></tr><tr><td>6</td><td>2.362e-01</td><td>4.100e-02</td></tr><tr><td>7</td><td>3.360e-03</td><td>1.422e-02</td></tr><tr><td>8</td><td>1.760e-07</td><td>5.239e-05</td></tr><tr><td>9</td><td>1.004e-12</td><td>5.707e-06</td></tr></tbody></table>

For now, we do not comment on the solvability of the global system. However, we note that if convergence can be achieved, then the Newton update will converge to zero and thus (5.12) yields that the physical condition of a sticky node is satisfied, i.e., the relative tangential displacement is equal to zero. However, in contrast to the physical non-penetration condition, it is not directly imposed as one might expect from a primal-dual active set strategy. A value for $s$ strictly smaller than one provides an additional damping. Moreover, (5.12) and (5.13b) show that the Newton solver couples the tangential and normal parts in the boundary conditions.

Figure 5.4 shows the influence of the friction coefficient $\nu$ on the convergence rate of the semi-smooth Newton scheme applied to a 3D contact problem. A comparison between a fixed-point approach and the Newton solver is given. For all three cases, super-linear rates are obtained for the
<figure><img src="image_39.png" /></figure>

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

<figure><img src="image_41.png" /><figcaption>Figure 5.4. (a) $v = 0.2$, (b) $v = 0.4$, (c) $v = 0.6$.</figcaption></figure>

semi-smooth Newton algorithm. The number of Newton steps required is quite insensitive to the friction coefficient, whereas for the simple fixed-point approach more iteration steps are required for a larger $v$.

### 5.4. Influence of the scaling parameter

It is well known, e.g., in plasticity, that the classical radial return algorithm, which is equivalent to setting $s = 0$, does not converge for large load steps. A similar observation can be made for contact problems with Coulomb or Tresca friction if $s = 0$. To get a better feeling for the scaling parameter $s$, we consider a simplified one-dimensional model $h(r, u) := \frac{5}{3} \max(|r+100u|, 100)^s (r - 100(r+100u)/\max(100, |r+100u|))$. Figure 5.5 shows the NCP function $h(\cdot, \cdot)$ for $s = 0$ and $s = 1$. A cut through $u = \text{const.}$ shows that for $s = 0$ the function is almost constant and thus a Newton linearization involves a badly conditioned Jacobian. The situation is improved for $s = 1$. This observation motivates the use of the scaling parameter $s$. The value $s = 1$ is already introduced in Hüeber et al. (2008) and also successfully applied in Koziara and Bicanic (2008).

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

<figure><img src="image_43.png" /><figcaption>Figure 5.5. NCP function for $s = 0$ (a) and $s = 1$ (b).</figcaption></figure>
<figure><img src="image_44.png" /><figcaption>Figure 5.6. Geometry (a), error decay of the semi-smooth Newton iteration for different values of $s$ (b) and total number of iterations for $t_1$ to $t_4$ with respect to $s$ (c).</figcaption></figure>

Figure 5.6(a) illustrates the geometry of a dynamic two-body contact problem in 3D with Coulomb friction. In Figure 5.6(b) we show the convergence of the semi-smooth Newton method for different scaling parameters. As can be clearly seen for small $s$, we do not have convergence. Numerically optimal convergence can be observed for $s = 0.75$. The convergence is also more robust for $s = 1$ than for smaller values of $s$. Figure 5.6(c) shows that for $s \in [0.75, 1]$, the global number of iterations required is almost constant. For $s > 1$, we observe a linear increase in the number of iteration steps, whereas for $s \le 0.5$ we see divergence. From now on we will restrict ourselves to the case $s \in [0, 1]$.

## 5.5. Stabilization in the pre-asymptotic range

Numerical experience shows that the use of the Robin boundary condition in the form of (5.13b) does not necessarily yield a robust and stable algorithm. Similar observations have been made in totally different fields of applications. We refer to Chan, Golub and Mulet (1999) and Hintermüller and Stadler (2006) for the use of primal-dual active set strategies in the field of image restoration and to the early work by Andersen, Christiansen, Conn and Overton (2000).

To get a better feeling, we consider the definition of $L_{p;s}^{k-1}$ in more detail. Let us assume that the Newton iterates converge towards a solution satisfying (3.11) and (3.12). Then, for a sliding node, we find that $L_1^{k-1}$ and $-\alpha_2^{k-1}L_2^{k-1}$ tend to the same limit and thus $L_{p;s}^{k-1}$ tends to the limit $L_p^\infty$ independently of $s$:

$$L_p^\infty := \frac{1}{\|\gamma_{p;\text{tr}}^t\|\|\gamma_p^t\|} \gamma_p^t \otimes \gamma_p^t - \frac{\|\gamma_p^t\|}{\|\gamma_{p;\text{tr}}^t\|} \text{Id.}$$

Noting that a rank-one matrix of the form $\mathbf{x} \otimes \mathbf{x}$ has exactly one non-trivial eigenvalue, namely $\|\mathbf{x}\|^2$, we find that $\text{Id} + L_{p;s}^{k-1}$ tends to a symmetric matrix
with eigenvalues 1 and $1 - \|\gamma_p^t\|/\|\gamma_{p;tr}^t\|$, which are strictly larger than zero for a physically correct sliding node. If no sliding occurs then $\|\gamma_p^t\| = \|\gamma_{p;tr}^t\|$, and the matrix is singular. In that situation, we face a standard Dirichlet boundary condition for $u_p^t$. Recalling that $\gamma_p^t$ stands for the negative surface traction in the tangential direction acting on node $p$, we find in the limit of a non-degenerate Robin boundary condition that

$$ -(\gamma_p^t)^k + L_{p;R}(u_p^t)^k = g_R^k, \quad L_{p;R} := \frac{2\mu c_t}{m_p} (-L_p^\infty)^{\frac{1}{2}} (\mathrm{Id} + L_p^\infty)^{-1} (-L_p^\infty)^{\frac{1}{2}}. $$

Here the Robin boundary data vector $g_R^k$ itself depends on the solution $(\gamma_p^n)^k$ and the previous iterate. The matrix $L_{p;R}$ is symmetric and positive definite, and thus yields a well-posed Robin condition for elliptic systems, and unique solvability for the linearized system is established. Unfortunately, these observations do not necessarily hold true in the pre-asymptotic range. Then it might occur that $\mathrm{Id} + L_{p;s}^{k-1}$ is non-singular but that $-(\mathrm{Id} + L_{p;s}^{k-1})^{-1} L_{p;s}^{k-1}$ is not positive semi-definite. This may result in no convergence.

To stabilize our approach, we introduce two modifications such that in each Newton step we obtain a well-defined system. This will be done by modifying $\omega_p^{k-1}$ and $\beta_p^{k-1}$, introduced in Lemma 5.5, so that

$$ \omega_p^{k-1} := \frac{\nu(\gamma_{p;tr}^n)^{k-1}}{\max(\nu(\gamma_{p;tr}^n)^{k-1}, \|(\gamma_p^t)^{k-1}\|)}, \quad (5.15) $$

$$ \beta_p^{k-1} := \min\left(1, \frac{1}{s(1 - \chi_p^{k-1} \zeta_p^{k-1})}\right), \quad (5.16) $$

where $\chi_p^{k-1}$ and $\zeta_p^{k-1}$ are defined by

$$ \chi_p^{k-1} := \frac{(\gamma_p^t)^{k-1} (\gamma_{p;tr}^t)^{k-1}}{\|(\gamma_p^t)^{k-1}\| \|(\gamma_{p;tr}^t)^{k-1}\|}, \quad \zeta_p^{k-1} := \min\left(\frac{\|(\gamma_p^t)^{k-1}\|}{\nu(\gamma_{p;tr}^n)^{k-1}}, 1\right). $$

Let us briefly comment on the definition of the two parameters $\omega_p^{k-1}$ and $\beta_p^{k-1}$. As we will see by our numerical results, both parameters play important roles in the pre-asymptotic range but do tend to 1 within the first few iterates if the algorithm converges. Thus asymptotically the exact Newton method is recovered. The damping parameter $\omega_p^{k-1}$ can be regarded as a penalty term and only differs from 1 if the friction bound is violated. It is obvious that $\zeta_p^{k-1} \in [0, 1]$ and that $\chi_p^{k-1}$, being the cosine of the angle between $(\gamma_p^t)^{k-1}$ and $(\gamma_{p;tr}^t)^{k-1}$, is in $[-1, 1]$.

It is now easy to see that $\beta_p^{k-1}$ is equal to 1 if $s \in [0, 0.5]$ or if $\chi_p^{k-1} \ge 0$. The last condition is equivalent to the fact that the angle between the actual and the trial tangential stress is bounded by $\pi/2$. Figure 5.7 shows $\beta_p^{k-1}$ as a function of $\chi_p^{k-1} \in [-1, 1]$ and $\zeta_p^{k-1} \in [0, 1]$ for $s \in \{0.5, 0.75, 1\}$. For
<figure><img src="image_45.png" /><figcaption>Figure 5.7. The damping term $\beta_p^{k-1}$: (a) $s = 0.5$, (b) $s = 0.75$, (c) $s = 1$.</figcaption></figure>

the physically correct solution, the angle between the actual and the trial tangential stress is zero. As our numerical results show, $\beta_p^{k-1}$ is mostly equal to one, reflecting the fact that within the first few iterations the orientation of the actual and trial tangential stress is adapted in the physically correct direction.

**Lemma 5.7.** Let $\omega_p^{k-1}$ and $\beta_p^{k-1}$ be defined by (5.15) and (5.16), respectively. Then the boundary condition (5.13) is well-posed.

*Proof.* In the case that $\mathrm{Id} + L_{p;s}^{k-1}$ is singular, the Robin condition reduces to a Dirichlet condition for the component of $(\mathbf{u}_p^t)^k$ in the direction of the kernel $K_1$ of $\mathrm{Id} + L_{p;s}^{k-1}$. Moreover, if $L_{p;s}^{k-1}$ is singular, then the Robin condition reduces to a Neumann condition for the component of $(\gamma_p^t)^k$ in the direction of the kernel $K_2$ of $L_{p;s}^{k-1}$. Since $K_1 \cap K_2 = \emptyset$, the Robin condition degenerates, in the direction of an element of the kernel $K_1$ or $K_2$, to either a well-defined Dirichlet or Neumann condition.

We have now to consider the non-degenerate case in more detail. Recalling that $L_{p;s}^{k-1} = \beta_p^{k-1}(\alpha_2^{k-1}\mathrm{Id} + \mathbf{x} \otimes \mathbf{y})$, we find for a non-singular $\mathrm{Id} + L_{p;s}^{k-1}$ and a non-singular $L_{p;s}^{k-1}$,

$$ (\mathrm{Id} + L_{p;s}^{k-1})^{-1} L_{p;s}^{k-1} = \mathrm{Id} - (\mathrm{Id} + L_{p;s}^{k-1})^{-1} = \frac{\beta_p^{k-1}}{1 + \beta_p^{k-1} \alpha_2^{k-1}} \left( \alpha_2^{k-1} \mathrm{Id} + \frac{\mathbf{x} \otimes \mathbf{y}}{1 + \beta_p^{k-1} \alpha_2^{k-1} + \beta_p^{k-1} \mathbf{x} \mathbf{y}} \right) =: aZ. $$

We note that $0 < \beta_p^{k-1} \le 1$ and $\alpha_2^{k-1} \ge -1$ and thus the factor $a$ in front of the matrix $Z$ is positive. It is easy to verify that we obtain the two eigenvalues of the matrix $Z$ by

$$ \alpha_2^{k-1} < 0, \quad \frac{(1 + \beta_p^{k-1} \alpha_2^{k-1})(\alpha_2^{k-1} + \mathbf{x}\mathbf{y})}{1 + \beta_p^{k-1}(\alpha_2^{k-1} + \mathbf{x}\mathbf{y})}. $$
<figure><img src="image_46.png" /><figcaption>Figure 5.8. Scaling factors $\omega_p^{k-1}$ and $\beta_p^{k-1}$ at iteration step $k$ and global convergence rate for $s = 1$.</figcaption></figure>

If $-1 < e_Z := \beta_p^{k-1}(\alpha_2^{k-1} + \mathbf{xy}) < 0$, then $Z$ is negative definite. The definition of $L_{p;s}^{k-1}$ in Lemma 5.5 shows that we can set

$$ \mathbf{x} := s\omega_p^{k-1}(\gamma_p)^{k-1} + (s-1)\alpha_2^{k-1}(\gamma_{p;\mathrm{tr}})^{k-1}, \quad \mathbf{y} := \frac{(\gamma_{p;\mathrm{tr}})^{k-1}}{\|(\gamma_{p;\mathrm{tr}})^{k-1}\|^2}. $$

Then, the value $e_Z$ can be rewritten as

$$ \begin{aligned} e_Z &= \beta_p^{k-1} \left( \alpha_2^{k-1} + s\omega_p^{k-1}\chi_p^{k-1} \frac{\|(\gamma_p)^{k-1}\|}{\|(\gamma_{p;\mathrm{tr}})^{k-1}\|} + (s-1)\alpha_2^{k-1} \right) \\ &= s\beta_p^{k-1} \left( \omega_p^{k-1}\chi_p^{k-1} \frac{\|(\gamma_p)^{k-1}\|}{\|(\gamma_{p;\mathrm{tr}})^{k-1}\|} + \alpha_2^{k-1} \right) \\ &= s\beta_p^{k-1}\alpha_2^{k-1} \left( 1 - \omega_p^{k-1}\chi_p^{k-1} \frac{\|(\gamma_p)^{k-1}\|}{\nu(\gamma_{p;\mathrm{tr}})^{k-1}} \right) = s\beta_p^{k-1}\alpha_2^{k-1}(1 - \chi_p^{k-1}\zeta_p^{k-1}) \\ &\ge \alpha_2^{k-1} \ge -1. \end{aligned} $$

The case $e_Z = -1$ is ruled out by the assumption that $\mathrm{Id} + L_{p;s}^{k-1}$ is non-singular, since $e_Z + 1$ is an eigenvalue of this matrix.

These considerations show that the Robin boundary condition, given by (5.13b), reduces to a Neumann or Dirichlet condition in the direction of the kernel of $L_{p;s}^{k-1}$ or $\mathrm{Id} + L_{p;s}^{k-1}$, respectively. For all other cases it forms a well-defined Robin condition. $\square$

**Remark 5.8.** We note that the proposed modifications are just one possibility. An alternative choice is to symmetrize and rescale $L_{p;s}^{k-1}$. All these modifications work well as long as they correctly normalize and, in the limit, tend to the original version; see also Hüeber et al. (2008) and Hüeber (2008).

To illustrate the effect of the modification, we consider different test examples in 3D. Figure 5.8 shows the case of a Tresca friction problem. We select a representative node and depict the value for $\omega_p^{k-1}$ and $\beta_p^{k-1}$. As can be seen, the modification actually applies only in the first few iteration steps.
<figure><img src="image_47.png" /><figcaption>Figure 5.9. Scaling factors $\omega_p^{k-1}$, $\beta_p^{k-1}$ at iteration step $k$ (a,b) and convergence rates for different strategies (c).</figcaption></figure>

In Figure 5.9, we apply the modifications to a Coulomb problem and consider the influence of the modification on the local and global convergence. To do so, we select the initial guess in two different ways. In the first case, it is chosen randomly on a fine mesh, whereas in the second situation it is interpolated from the solution of the previous coarser level and thus can be expected to be good. For this nested approach the modification does not affect the iteration scheme and thus the same rates of convergence can be observed. This is in good agreement with the observation that asymptotically $\omega_p^{k-1}$ and $\beta_p^{k-1}$ tend to 1 in the case of convergence. The situation is considerably different if we apply a bad initial guess. Here, without the modification, no convergence at all can be observed, whereas the modified version still shows a reasonable rate. These two test examples show that the effect of the modification can be neglected if the start iteration is already good enough. However in the pre-asymptotic range, it is of great significance and highly enlarges the domain of convergence.

Our last test combines the scaling $s \in [0, 1]$ with the modification for the nodes $p \in A^k$. Here, we consider a dynamic Coulomb problem in 3D on non-matching meshes where the actual contact zone is not simply connected. Figure 5.10 shows the influence of the convergence rate as a function of $s$. As can be clearly seen in Figure 5.10(a), the number of Newton iterations required is quite insensitive provided that $s$ is large enough. For small $s$, no convergence at all can be obtained. Therefore, at time $t_2$, only the results for $s = 0.5$, $s = 0.75$ and $s = 1$ are plotted. The choices $s = 0$ and $s = 0.25$ do not yield a convergent scheme.

Finally Figure 5.11 shows that for a selected node both parameters $\omega_p^{k-1}$ and $\beta_p^{k-1}$ tend to 1 very fast in the case of $s > 0.5$. Moreover, the cosine of the angle between the tangential stress and the trial tangential stress tends to 1. This indicates that our algorithm is able to adjust the correct sliding direction within the first few iterates. The situation is changed drastically
<figure><img src="image_48.png" /></figure>

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

<figure><img src="image_50.png" /><figcaption>Figure 5.10. Number of required iterations (a) and convergence rates for two different time steps (b,c).</figcaption></figure>

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

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

<figure><img src="image_53.png" /><figcaption>Figure 5.11. Scaling factors $\chi_p^{k-1}$, $\beta_p^{k-1}$ and $\omega_p^{k-1}$ at iteration step $k$.</figcaption></figure>

if $s \le 0.5$. Then the correct sliding direction cannot be identified and $\omega_p^{k-1}$ and $\chi_p^{k-1}$ do not tend to 1, reflecting the fact that no convergence is obtained.

## 5.6. Mesh-dependent convergence rates

Although the abstract framework of semi-smooth Newton methods is very flexible and quite attractive for a large class of problems, there is one bottleneck. The convergence rate is, in general, mesh-dependent. Numerical results show that the number of Newton steps increases linearly with the refinement level. Several strategies exist to overcome this problem, depending on the type of application. One of the most efficient ones is to embed the solver in a nested iteration. This is extremely easy to realize with time-dependent problems or with adaptive refinement techniques. In both situations the initial guess can be interpolated from the previous time step or mesh. This simple pre-processing is quite often sufficient to obtain a level-independent number of non-linear solver steps. Alternatively, or additionally, we can combine the Newton approach with an inexact solution strategy. Quite often the arising linear system is not solved by a fast direct solver but by preconditioned Krylov or subspace correction methods such as multigrid or FETI techniques. Then it is only natural to make the Newton update after a small number of steps of the linear solver. To avoid
over-solving, during the first non-linear iteration steps, the number of inner iterations should be set dynamically.

We consider the same example as shown in Figure 4.3 to illustrate the effect of the mesh-dependent convergence rates. We use a hierarchy of uniformly refined meshes and compare three different strategies. The ‘exact’ one starts with a randomly chosen initial guess on each level and solves the resulting linearized system in each Newton step by a multigrid method with fixed small tolerance. By $K_l$, we denote the number of Newton steps required. As can be seen in Table 5.2, $K_l$ depends linearly on $l$, and thus the total number $N_l^{\text{MG}}$ of multigrid steps on level $l$ increases linearly with $l$. The situation is different if we apply an inexact strategy, using the same bad initial guess, but do a Newton update after each multigrid step. This results in an ‘inexact’ strategy where the correct Jacobian of the NCP function is used but the system stiffness matrix is replaced by its multigrid approximation. Here, $M_l$ stands for the number of iterations required to identify the correct active sets. As before, this number increases linearly, but the total number of multigrid steps is significantly reduced. Finally we combine this inexact strategy with a good initial guess obtained by interpolation from the solution on the previous level. We call this the ‘nested’ approach. Then the total number of multigrid steps is bounded independently of the refinement level and is comparable to the number required to solve one linearized problem.

<table><caption>Table 5.2. Number of total multigrid steps for different strategies.</caption><thead><tr><th colspan="2">Strategy</th><th colspan="2">Exact</th><th colspan="2">Inexact</th><th colspan="2">Nested</th></tr><tr><th>Level l</th><th>DOF</th><th>K<sub>l</sub></th><th>MG</th><th>M<sub>l</sub></th><th>N<sub>l</sub><sup>MG</sup></th><th>M<sub>l</sub></th><th>N<sub>l</sub><sup>MG</sup></th></tr></thead><tbody><tr><td>1</td><td>27</td><td>3</td><td>46</td><td>3</td><td>12</td><td>3</td><td>12</td></tr><tr><td>2</td><td>125</td><td>3</td><td>62</td><td>3</td><td>16</td><td>3</td><td>16</td></tr><tr><td>3</td><td>729</td><td>4</td><td>72</td><td>4</td><td>17</td><td>4</td><td>15</td></tr><tr><td>4</td><td>4913</td><td>6</td><td>86</td><td>6</td><td>17</td><td>6</td><td>14</td></tr><tr><td>5</td><td>35937</td><td>7</td><td>106</td><td>9</td><td>19</td><td>6</td><td>16</td></tr></tbody></table>

In more complex 3D situations, or for non-linear material laws, more than one multigrid step is required before an update can be performed. In a last test, we combine semi-smooth Newton techniques and inexact solvers for the linearized system with an overlapping two-scale domain decomposition method in 3D. We refer to Brunnen, Hager, Wohlmuth and Schmid (2008) and Brunnen and Wohlmuth (2009) for details of the model and the specification of the data. In addition to the contact formulation, plasticity effects are taken into account. The approach is motivated by possible
<figure><img src="image_54.png" /></figure>

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

<figure><img src="image_56.png" /><figcaption>Figure 5.12. Surface mesh of the two-scale domain decomposition approach (a), convergence rates for the non-linear Newton solver (b,c).</figcaption></figure>

applications to incremental metal cold forming processes. Here, the forming zone is small, but very mobile, and the work tool is contacting almost every point of the workpiece at some time in the process. To avoid expensive re-meshing and to reduce the complexity of the elasto-plastic constitutive equations, an operator splitting technique in space can be introduced. A small local but mobile subdomain, with a fine mesh and the fully non-linear contact and plasticity model, will interact with the global coarse mesh associated with a simplified model. To solve the fully coupled non-linear system, semi-smooth Newton techniques in combination with a block Gauss-Seidel solver are quite efficient. This is in particular true if inexact strategies are applied. To avoid over-solving during the first non-linear iteration steps, the number of inner iteration should be set dynamically. Here, the stopping criterion is based on Dembo, Eisenstat and Steinhaug (1982) and Eisenstat and Walker (1996).

Figure 5.12 shows the increase in efficiency of the non-linear solver if, within each Newton step, the linearized coupled domain decomposition system is not solved exactly but by a few block Gauss-Seidel steps. The number of inner iteration steps is set dynamically depending on the non-linear residual. As can be seen from Figure 5.12(c), during the first few Newton iterations, there is no need to solve the linearized system up to very high accuracy. To obtain the full efficiency of the approach, the accuracy of the linear solver has to be gradually improved during the non-linear solution process.

## 6. A *posteriori* error estimates and adaptivity

Adaptive techniques based on *a posteriori* error estimators play an important role in enhancing the performance of the numerical simulation algorithm and are well established for finite element methods: see the monographs by Ainsworth and Oden (2000), Babuška and Strouboulis (2001), Han (2005), Repin (2008), Verfürth (1996) and the references therein. For abstract variational inequalities we refer to Ainsworth, Oden and Lee (1993),
Bostan, Han and Reddy (2005), Erdmann, Frei, Hoppe, Kornhuber and Wiest (1993), Fuchs and Repin (2010), Liu and Yan (2000), Nochetto, von Petersdorff and Zhang (2010), Moon, Nochetto, von Petersdorff and Zhang (2007) and Suttmeier (2005), whereas obstacle-type problems are considered in Bildhauer, Fuchs and Repin (2008), Braess (2005), Chen and Nochetto (2000), French, Larsson and Nochetto (2001), Hoppe and Kornhuber (1994), Johnson (1992), Kornhuber and Zou (2011), Nochetto, Siebert and Veeser (2003, 2005) and Veeser (2001), and early approaches for contact problems can be found in Blum and Suttmeier (2000), Buscaglia, Duran, Fancello, Feijoo and Padra (2001), Carstensen, Scherf and Wriggers (1999), Lee and Oden (1994) and Wriggers and Scherf (1998). A residual-type error estimator is introduced and analysed in Hild and Nicaise (2005, 2007) and in Bostan and Han (2006) and Hild and Lleras (2009) for a one-sided contact problem without friction and with Coulomb friction, respectively. In addition to standard face and volume residual terms, extra terms reflecting the non-conformity of the approach are taken into account. For boundary element discretizations, we refer to Eck and Wendland (2003) and Maischak and Stephan (2005, 2007). Although the error estimator in Hild and Lleras (2009) provides a mesh-independent upper bound for the discretization error, not all terms result in optimal lower bounds. Thus the efficiency of the error estimator cannot be guaranteed from the theoretical point of view. Early results on hp-techniques for frictional contact problems can be found in Lee and Oden (1994), whereas in the recent contribution of Dörsek and Melenk (2010) a simplified Tresca problem with a given surface normal traction equal to zero has been studied.

As it turns out, the saddle-point approach (2.15) and its discrete version (3.9) provide an excellent starting point for the construction of an error indicator. In Wohlmuth (2007) an estimator was introduced for the case of no friction and non-matching meshes. A theoretical analysis shows that a constant-free global upper bound and local lower bounds for the error can be established. However, in contrast to the standard linear conforming setting, the additional higher-order term is solution-dependent and cannot be controlled within the adaptive refinement process.

These first results can be improved considerably by following a more general construction principle. Firstly, we use the equilibrium of the saddle-point approach to construct our indicator. Secondly, we consider the influence of the discretization (3.8) of (2.14). For the special case of a one-sided contact problem without friction, it is shown in Weiss and Wohlmuth (2009) that an indicator, constructed in this way, provides upper and local lower bounds for the discretization error and is thus an efficient error estimator. Moreover, the higher-order terms are standard data oscillation terms and can be controlled within the adaptive refinement strategy. The results are shown for a flux-based approach but can easily be generalized to cover the
case of a residual-based error estimator: see Hüeber and Wohlmuth (2010) for numerical results. Here, we follow these lines and extend the approach and the analysis to contact problems with friction.

In this section, we discuss different element-oriented error indicators, i.e., $η^2 := \sum_{T \in \mathcal{T}_l} η_T^2$, and analyse the influence of the variational inequality. We focus both on the theoretical bounds and on computational aspects. Due to the variational inequality character of the given problem, we have to include a term which measures the non-conformity of the Lagrange multiplier. This contribution can be decomposed into a contact and friction term. For ease of presentation and analysis, we restrict ourselves to the two-dimensional setting, polygonal domains, a zero gap and a constant contact normal. Moreover, we do not analyse the influence of non-matching meshes but do provide the construction of the error indicator.

As is standard, we use a data oscillation term,

$$ \xi_i^2 := \sum_{T \in \mathcal{T}_l} \frac{h_T^2}{2\mu} \|\mathbf{f} - \mathbf{\Pi}_i \mathbf{f}\|_{0;T}^2, \quad i = 0, 1. \qquad (6.1) $$

To keep the notation simple, we further assume that we are working with simplicial meshes and that the given boundary data is compatible with the discretization, i.e., $\mathbf{u}_l|_{\Gamma_D} = \mathbf{u}_D$ and $\mathbf{f}_N$ is piecewise cubic and continuous on each straight segment $γ_i$ of $\Gamma_N$. Moreover, we assume that $\mathbf{n}_j\mathbf{f}_N|_{γ_i}(p) = \mathbf{n}_i\mathbf{f}_N|_{γ_j}(p)$ with $\partial\gamma_i \cap \partial\gamma_j = p$ and that $\mathbf{n}\mathbf{f}_N|_{\Gamma_N}(p) = 0$ where $p = \partial\Gamma_C \cap \partial\Gamma_N$, and $\mathbf{n}$ is the outer unit normal on $\Gamma_C$. Otherwise, as is standard, additional boundary face terms have to be included, measuring the weighted $L^2$-norms of the boundary error: see, e.g., Repin, Sauter and Smolianski (2003). The weight $\sqrt{h_f/(2\mu)}$ for the Neumann term is the inverse of the weight for the Dirichlet term, reflecting the $H^{1/2}$-duality between displacement trace and surface tractions. Here $h_f$ stands for the diameter of the boundary face $f$. In the following, we use the piecewise cubic biorthogonal basis function shown in Figure 3.7(a).

As a first preliminary step, we reformulate the coupled problem (3.9) by introducing a weakly consistent Neumann force in $M_l^m$ on the master side. To do so, we introduce $\mathbf{\Pi}_l^{*;m} := (\mathbf{\Pi}_l^{*;m})^d$ as the dual mortar projection $\mathbf{\Pi}_l^{*;m}: M \to M_l^m$ with respect to the master side,

$$ \int_{\Gamma_C} \mathbf{\Pi}_l^{*;m} \mu v_l \, ds = \langle \mu, v_l \rangle_{\Gamma_C}, \quad v_l \in W_l^m. \qquad (6.2) $$

Keeping in mind that $M_l^m$ and $W_l^m$ reproduce constants and have a locally defined basis, a straightforward Bramble-Hilbert argument implies the approximation properties

$$ \|\mu - \mathbf{\Pi}_l^{*;m} \mu\|_{-\frac{1}{2};\Gamma_C}^2 \le \sum_{f \in \mathcal{F}_l^m} h_f \|\mu - \mathbf{\Pi}_l^{*;m} \mu\|_{0;f}^2, \quad \mu \in L^2(\Gamma_C), \qquad (6.3) $$