On Nitsche's method for elastic contact problems

TL;DR

This paper analyzes Nitsche's method for elastic contact problems, establishing quasi-optimality and error estimates, and demonstrating its robustness and efficiency through numerical experiments.

Contribution

It introduces three Nitsche's mortaring techniques for contact boundary stabilization with minimal regularity requirements and no saturation assumption.

Findings

Nitsche's method is robust for elastic contact problems.

A posteriori error estimators are effective.

Numerical experiments confirm theoretical results.

Abstract

We show quasi-optimality and a posteriori error estimates for the frictionless contact problem between two elastic bodies with a zero-gap function. The analysis is based on interpreting Nitsche's method as a stabilised finite element method for which the error estimates can be obtained with minimal regularity assumptions and without the saturation assumption. We present three different Nitsche's mortaring techniques for the contact boundary each corresponding to a different stabilising term. Our numerical experiments show the robustness of Nitsche's method and corroborates the efficiency of the a posteriori error estimators.

Figures

Equations

$$\bm{\varepsilon}(\bm{u}_{i})=\frac{1}{2}\Big{(}\nabla\bm{u}_{i}+(\nabla\bm{u}_{i})^{T}\Big{)}.$$

$${\mathbf{\sigma }}_i({\mathbf{u}}_i)=2{\mu }_i\mathbf{\epsilon }({\mathbf{u}}_i)+{\lambda }_i\mathrm{tr}\mathbf{\epsilon }({\mathbf{u}}_i)\mathbf{I},$$

$${\mathbf{\sigma }}_i({\mathbf{u}}_i){\mathbf{n}}_i={\mathbf{\sigma }}_{i,n}({\mathbf{u}}_i)+{\mathbf{\sigma }}_{i,t}({\mathbf{u}}_i).$$

$${\sigma }_{1,n}({\mathbf{u}}_1)={\mathbf{\sigma }}_{1,n}({\mathbf{u}}_1)\cdot {\mathbf{n}}_1,$$

$${\sigma }_{2,n}({\mathbf{u}}_2)=-{\mathbf{\sigma }}_{2,n}({\mathbf{u}}_2)\cdot {\mathbf{n}}_2,$$

$${\sigma }_{1,n}({\mathbf{u}}_1)={\sigma }_{2,n}({\mathbf{u}}_2),{\sigma }_{i,n}({\mathbf{u}}_i)\le 0,i=1,2.$$

$${\mathbf{u}}_1\cdot {\mathbf{n}}_1+{\mathbf{u}}_2\cdot {\mathbf{n}}_2\le 0,$$

$$u_n=-({\mathbf{u}}_1\cdot {\mathbf{n}}_1+{\mathbf{u}}_2\cdot {\mathbf{n}}_2)$$

$$\left\llbracket u_{n}\right\rrbracket\geq 0,$$

$$-\mathbf{div}{\mathbf{\sigma }}_i({\mathbf{u}}_i)$$

$${\mathbf{u}}_i$$

$${\mathbf{\sigma }}_i({\mathbf{u}}_i){\mathbf{n}}_i$$

$${\mathbf{\sigma }}_{i,t}({\mathbf{u}}_i)$$

$${\sigma }_{1,n}({\mathbf{u}}_1)-{\sigma }_{2,n}({\mathbf{u}}_2)$$

$$\displaystyle\left\llbracket u_{n}\right\rrbracket$$

$${\sigma }_{i,n}({\mathbf{u}}_i)$$

$$\displaystyle\left\llbracket u_{n}\right\rrbracket\sigma_{i,n}(\bm{u}_{i})$$

$$-\mathbf{div}{\mathbf{\sigma }}_i({\mathbf{u}}_i)$$

$${\mathbf{u}}_i$$

$${\mathbf{\sigma }}_i({\mathbf{u}}_i){\mathbf{n}}_i$$

$${\mathbf{\sigma }}_{i,t}({\mathbf{u}}_i)$$

$$\lambda +{\sigma }_{1,n}({\mathbf{u}}_1)$$

$$\lambda +{\sigma }_{2,n}({\mathbf{u}}_2)$$

$$\displaystyle\left\llbracket u_{n}\right\rrbracket$$

$$\lambda$$

$$\displaystyle\left\llbracket u_{n}\right\rrbracket\lambda$$

$${\mathbf{V}}_i=\{{\mathbf{w}}_i\in [H^1({\Omega }_i){]}^d:{\mathbf{w}}_i{|}_{{\Gamma }_{D,i}}=0\},$$

$$\Vert w{\Vert }_{\frac{1}{2},\Gamma }^2=\Vert w{\Vert }_{0,\Gamma }^2+{\int }_{\Gamma }{\int }_{\Gamma }\frac{|w(x)-w(y){|}^2}{|x-y{|}^d}\mathrm{d}x\mathrm{d}y.$$

$$\Vert \xi {\Vert }_{-\frac{1}{2},\Gamma }=\underset{w\in W}{\sup }\frac{⟨w,\xi ⟩}{\Vert w{\Vert }_{\frac{1}{2},\Gamma }}.$$

$$\Lambda =\{\xi \in H^{-\frac{1}{2}}(\Gamma ):⟨w,\xi ⟩\ge 0\forall w\in H^{\frac{1}{2}}(\Gamma ),w\ge 0\text{a.e. on Γ}\}$$

$$\mathcal{B}(\bm{w},\xi;\bm{v},\eta)=\sum_{i=1}^{2}(\bm{\sigma}_{i}(\bm{w}_{i}),\bm{\varepsilon}(\bm{v}_{i}))_{\Omega_{i}}-\left\langle\left\llbracket v_{n}\right\rrbracket,\xi\right\rangle-\left\langle\left\llbracket w_{n}\right\rrbracket,\eta\right\rangle,$$

$$\mathcal{L}(\mathbf{v})=\sum _{i=1}^2({\mathbf{f}}_i,{\mathbf{v}}_i{)}_{{\Omega }_i}.$$

$$\mathcal{B}(\mathbf{u},\lambda ;\mathbf{v},\eta -\lambda )\le \mathcal{L}(\mathbf{v})\forall (\mathbf{v},\eta )\in \mathbf{V}\times \Lambda .$$

$${\mathbf{V}}_{i,h}$$

$${\mathbf{V}}_h$$

$$Q_h$$

$${\Lambda }_h=\{{\eta }_h\in Q_h:{\eta }_h\ge 0\}.$$

$${\mathcal{B}}_h({\mathbf{w}}_h,{\xi }_h;{\mathbf{v}}_h,{\eta }_h)=\mathcal{B}({\mathbf{w}}_h,{\xi }_h;{\mathbf{v}}_h,{\eta }_h)-\alpha {\mathcal{S}}_h({\mathbf{w}}_h,{\xi }_h;{\mathbf{v}}_h,{\eta }_h),$$

$$\mathcal{S}_{h}(\bm{w}_{h},\xi_{h};\bm{v}_{h},\eta_{h})=\sum_{i=1}^{2}\sum_{E\in\mathcal{G}_{h}^{i}}\frac{h_{E}}{\mu_{i}}\Big{(}\xi_{h}+\sigma_{i,n}(\bm{w}_{i,h}),\eta_{h}+\sigma_{i,n}(\bm{v}_{i,h})\Big{)}_{E},$$

Full text

The financial support from the Portuguese government through FCT (Fundação para a Ciência e a Tecnologia), I.P., under the projects PTDC/MAT-PUR/28686/2017 and UTAP-EXPL/MAT/0017/2017, is gratefully acknowledged.

1. introduction

In this paper, we analyse the Nitsche method for elastic contact problems. Over the last decade, this method has been studied by a number of authors, see, e.g., [9, 6, 7, 10], and shown to be a robust and efficient method. The advantages are an easy implementation based on the displacement variables only and, when compared to mixed methods with Lagrange multipliers, the absence of an ”inf-sup” stability condition which renders a symmetric positive definite system instead of one with a saddle point structure.

From a theoretical point of view, the previously mentioned works suffer from two shortcomings. First, for the problem posed in $H^{1}$, the solution is typically assumed to be in $H^{s}$, with $s>3/2$. Second, the a posteriori error analyses are often based on a non-rigorous saturation assumption.

We have addressed these issues in our recent articles, cf. [12, 13]. Our approach dates back to [23] where different ways to enforce weakly the Dirichlet boundary conditions were discussed in the context of the so called stabilised mixed methods [2, 3] wherein the bilinear form of the original mixed finite element method is augmented with a properly weighted residual term to ensure stability. In [23], it was shown that the local elimination of the Lagrange multiplier leads essentially to a method introduced by Nitsche in the early age of the finite element analysis [22]. Since Nitsche’s method is straightforward both to analyse (under the additional smoothness assumption) and to implement, we started to advocate it, in particular for contact problems, cf. [24, 4].

What we have realised recently is that one should take full advantage of the relation between Nitsche’s and stabilised method when analysing the former. In fact, we were able to get rid of both the smoothness and the saturation assumption for the membrane obstacle problem in [12]. In this paper, we will continue on this path and perform an error analysis, both quasi-optimality and a posteriori, for a simplified two-body contact problem without friction. Besides the theoretical improvements, we present three versions of the Nitsche’s method where the changes in the material parameters between the bodies are taken into account. The simplest is a typical ”master-slave” approach where the contact surface of the stiffer body is chosen as the master part and the slave surface is then mortared by the Nitsche’s technique. In the two other variants, the material parameters appear as weights in the Nitsche formulation so that the methods decide by themselves which part is the master and which is the slave. In order to simplify the notation, analysis and implementation of the adaptive methods, we assume that the elastic bodies are initially in full contact, see, e.g., [17], and leave the case with a non-vanishing initial gap between the elastic bodies for a future work.

Although our analysis is built upon our earlier works, cf. [12, 14], we will present proofs of all the main theorems. We also note that the elastic contact problem literature is vast and therefore we only refer to the review paper [27], and to all the references therein, for the analysis and application of finite element methods arising from mixed formulations and to [20, 8], and to all the references therein, for the a posteriori error analyses of contact problems. We end the paper by presenting results of our computational experiments.

2. The contact problem

Let $\Omega_{i}\subset\mathbb{R}^{d}$, $i=1,2$, $d\in\{2,3\}$, denote two elastic bodies in their reference configuration and assume that the bodies are initially in contact. Moreover, assume that $\Omega_{i}$ are polygonal (polyhedral) domains and denote by $\Gamma=\partial\Omega_{1}\cap\partial\Omega_{2}$ their common boundary. The boundary $\partial\Omega_{i}$ is split into three disjoint sets $\Gamma_{D,i},\Gamma_{N,i}$ and $\Gamma_{C,i}$, with $\Gamma_{D,i}$ denoting the part where homogeneous Dirichlet data is given, $\Gamma_{N,i}$ the part with a Neumann boundary condition and $\Gamma_{C,i}$ the part where contact can occur, see Figure 1.

Letting $\bm{u}_{i}:\Omega_{i}\rightarrow\mathbb{R}^{d}$, $i=1,2$, be the displacement of the body $\Omega_{i}$, the infinitesimal strain tensor is defined as

[TABLE]

We assume homogenous isotropic bodies and a plain strain problem in the two dimensional case. The stress tensor is thus given by

[TABLE]

where $\mu_{i}>0$ is the shear modulus and $\lambda_{i}$ the second Lamé parameter of the body $\Omega_{i}$ and $\bm{I}$ denotes the $d$-dimensional identity tensor. We will exclude the possibility that the materials are nearly incompressible and hence it holds $\lambda_{i}\lesssim\mu_{i}$. (For nearly incompressible materials the standard approach of reformulating the problem in mixed form [5] should be used.)

By $\bm{n}_{i}\in\mathbb{R}^{d}$ we denote the outward unit normal to $\partial\Omega_{i}$, and define $\bm{n}=\bm{n}_{1}=-\bm{n}_{2}$. In what follows, $\bm{t}$ denotes any unit vector that satisfies $\bm{n}\cdot\bm{t}=0$.

We decompose the traction vector on $\partial\Omega_{i}$, $\bm{\sigma}_{i}(\bm{u}_{i})\bm{n}_{i}$, into its normal and tangential parts, viz.

[TABLE]

For the scalar normal tractions we use the sign convention

[TABLE]

and

[TABLE]

and note that on $\Gamma$ these tractions are either both zero or continuous and compressive, i.e. it holds that

[TABLE]

The physical non-penetration constraint on $\Gamma$ reads as

[TABLE]

which, defining

[TABLE]

can be written as

[TABLE]

where \left\llbracket\cdot\right\rrbracket denotes the jump over $\Gamma$.

We thus have the following problem.

Problem 1** (Strong formulation).**

Find $\bm{u}_{i}:\Omega_{i}\rightarrow\mathbb{R}^{d}$, $i=1,2$, $d\in\{2,3\}$, such that

[TABLE]

where $\bm{f}_{i}\in[L^{2}(\Omega_{i})]^{d}$ denotes the volume force on $\Omega_{i}$.

Letting $\lambda=-{\sigma_{1,n}(\bm{u}_{1})}=-{\sigma_{2,n}(\bm{u}_{2})}$ denote a Lagrange multiplier associated with the contact constraint, we obtain an equivalent mixed formulation in which the normal traction on the contact surface is an independent unknown.

Problem 2** (Mixed formulation).**

Find $\bm{u}_{i}:\Omega_{i}\rightarrow\mathbb{R}^{d}$, $i=1,2$, $d\in\{2,3\}$, and $\lambda:\Gamma\rightarrow\mathbb{R}$, such that

[TABLE]

To present a variational formulation for Problem 2.11, we introduce function spaces for the displacements

[TABLE]

and equip them with the usual norms $\|\cdot\|_{1,\Omega_{i}}$. Moreover, we write $\bm{V}=\bm{V}_{1}\times\bm{V}_{2}$ and assume that $\Gamma$ is a compact subset of $\partial\Omega_{i}\setminus\Gamma_{D,i}$ for $i=1,2$. Thus the normal components of the displacement traces on the contact zone are in $H^{\frac{1}{2}}(\Gamma)$ with the intrinsic norm in $H^{\frac{1}{2}}(\Gamma)$ defined by (cf., e.g., [25])

[TABLE]

The inequality constraint on $\Gamma$ is imposed by the Lagrange multiplier which belongs to $H^{-\frac{1}{2}}(\Gamma)$, the topological dual of $H^{\frac{1}{2}}(\Gamma)$, i.e. $H^{-\frac{1}{2}}(\Gamma)=H^{\frac{1}{2}}(\Gamma)^{\prime}$. The duality pairing is denoted by $\left\langle\cdot,\cdot\right\rangle:H^{\frac{1}{2}}(\Gamma)\times H^{-\frac{1}{2}}(\Gamma)\rightarrow\mathbb{R}$, and the norm is then

[TABLE]

Moreover, we define the positive part of $H^{-\frac{1}{2}}(\Gamma)$ as

[TABLE]

and introduce the bilinear and linear forms

[TABLE]

and

[TABLE]

The variational problem now reads as follows:

Problem 3** (Weak formulation).**

Find $(\bm{u},\lambda)\in\bm{V}\times\varLambda$ such that

[TABLE]

We refer to [16, 15] for the derivation of weak formulation from Problem 2.11 and for the proof of existence and uniqueness of solutions to problem (2.18).

3. Finite element method

Let the bodies $\Omega_{i}\subset\mathbb{R}^{d}$ be separately divided into sets of non-overlapping simplices $\mathcal{C}_{h}^{i}$, $i=1,2$. The $d-1$ dimensional facets of the elements in $\mathcal{C}_{h}^{i}$ are further divided into the set of interior facets $\mathcal{E}_{h}^{i}$, the set of facets on the contact boundary $\mathcal{G}_{h}^{i}$, and the set of facets on the Neumann boundary $\mathcal{N}_{h}^{i}$. We denote by $\mathcal{G}_{h}^{12}$ the boundary mesh on $\Gamma$ which is obtained by intersecting the facets of $\mathcal{G}_{h}^{1}$ and $\mathcal{G}_{h}^{2}$. In particular, each $E\in\mathcal{G}_{h}^{12}$ corresponds to a pair $(E_{1},E_{2})\in\mathcal{G}_{h}^{1}\times\mathcal{G}_{h}^{2}$ such that $E=E_{1}\cap E_{2}$. The finite element subspaces are

[TABLE]

where $P_{p}(K)$ denotes the polynomials of degree $p$ on $K$. Moreover, we introduce a subset of $\varLambda$, denoted by $\varLambda_{h}$, as the positive part of $Q_{h}$, i.e.

[TABLE]

Now, defining a stabilised bilinear form $\mathcal{B}_{h}$ through

[TABLE]

where $\alpha>0$ is a stabilisation parameter and

[TABLE]

we arrive at the following finite element formulation which is an extension of the mortar method introduced in [19, 14].

Problem 4** (Stabilised discrete formulation).**

Find $(\bm{u}_{h},\lambda_{h})\in\bm{V}_{h}\times\varLambda_{h}$ such that

[TABLE]

We will now derive an equivalent formulation wherein the Lagrange multiplier is not explicitly present. To this end, we start by defining $L^{2}(\Gamma)$-functions $\mathpzc{h}_{i}$ through

[TABLE]

and introduce the notation

[TABLE]

i.e. a convex combination of the discrete normal tractions. Furthermore, we let

[TABLE]

where

[TABLE]

Next, we will show that the discrete Lagrange multiplier $\lambda_{h}$ can be eliminated locally (i.e. element by element). This leads to a Nitsche formulation with the displacements as sole unknowns. Choosing $\bm{v}_{h}=\bm{0}$ in the variational inequality (3.7), gives

[TABLE]

which, in view of the notation defined above, can be written as

[TABLE]

Let then $E\in\mathcal{G}_{h}^{12}$ be an element on which $\lambda_{h}|_{E}>0$ and denote by $\phi_{E}$ one of the basis functions of $Q_{h}|_{E}$. Moreover, choose a test function $\eta_{h}$ in (3.13) in such a way that it vanishes at $\Gamma\setminus E$ and $\eta_{h}|_{E}=\lambda_{h}\pm\epsilon\phi_{E}$, with $\epsilon>0$ chosen small enough so that $\eta_{h}|_{E}>0$. It follows that

[TABLE]

and, since

[TABLE]

we conclude that

[TABLE]

This shows that

[TABLE]

where $(a)_{+}=\max(0,a)$ denotes the positive part of $a$.

The discrete contact region, defined as

[TABLE]

can now, in view of (3.17), be written as

[TABLE]

On the other hand, testing with $\bm{v}_{h}$ in (3.7) and using (3.17) yields

[TABLE]

It follows from (3.10) that

[TABLE]

and on $\Gamma_{c}(\bm{u}_{h})$ it holds that

[TABLE]

Therefore, defining the jump

[TABLE]

and the $L^{2}(\Gamma)$-function

[TABLE]

and substituting the above five expressions into (3.20), we obtain after rearranging terms the following Nitsche’s formulation for Problem 3.7 with $\bm{u}_{h}$ as the sole unknown.

Nitsche formulation 1**.**

Find $\bm{u}_{h}\in\bm{V}_{h}$ such that

[TABLE]

Remark 3.1**.**

Since $\sigma_{n}(\bm{u}_{i})$ vanishes on $\Gamma\setminus\Gamma_{c}(\bm{u}_{h})$, this set can be reinterpreted as being part of $\Gamma_{N,i}$, $i=1,2$. Consequently, the term

[TABLE]

can be dropped.

Next we present two other variants of Nitsche’s method. The first is the so called ”master-slave” formulation.

Assume that the material parameters satisfy $\mu_{1}\geq\mu_{2}$. The body $\Omega_{1}$ is the master part, $\Omega_{2}$ the slave, and the mortaring at the contact surface is only done for the latter, less rigid body, i.e. the stabilising term is now

[TABLE]

Repeating the steps above, we obtain $\lambda_{h}=(l_{h}(\bm{u}_{h}))_{+}\,,$ with

[TABLE]

The contact region $\Gamma_{c}(\bm{u}_{h})$ is given by (3.19), with $l_{h}(\bm{u}_{h})$ taken from (3.28), and we have the following method.

Nitsche formulation 2**.**

Find $\bm{u}_{h}\in\bm{V}_{h}$ such that

[TABLE]

Again, the term

[TABLE]

can be dropped, see Remark 3.1.

In the third alternative, we follow [18] and define the stabilising term through

[TABLE]

Repeating once more the above computations, we arrive at the following method.

Nitsche formulation 3**.**

Find $\bm{u}_{h}\in\bm{V}_{h}$ such that

[TABLE]

with $\Gamma_{c}(\bm{u}_{h})$ given by (3.19) $($and $l_{h}(\bm{u}_{h})$ as in (3.17)$)$.

Also here the term

[TABLE]

can be dropped.

4. Error analysis

The energy norm for the problem is

[TABLE]

Since we exclude nearly incompressible materials, it holds $\lambda_{i}\lesssim\mu_{i}$, and hence with our choice of boundary conditions the Korn inequality is valid in both regions, and we have the norm equivalence

[TABLE]

The error estimate will be given in the continuous norm

[TABLE]

but in the analysis we will also use the following mesh dependent norm

[TABLE]

Theorem 4.1** (Continuous stability).**

For every $(\bm{w},\xi)\in\bm{V}\times Q$ there exists $\bm{v}\in\bm{V}$ such that

[TABLE]

and

[TABLE]

Proof.

It is well-known that the inf-sup condition

[TABLE]

holds in both subdomains $\Omega_{i}$ (cf. [1]). Therefore

[TABLE]

Assume then that $(\bm{w},\xi)\in\bm{V}\times Q$ is given and let $\bm{v}_{i}=\bm{w}_{i}-\bm{q}_{i}$ where $\bm{q}_{i}\in\bm{V}_{i}$ solves the problem

[TABLE]

Choosing $\bm{z}_{i}=\bm{q}_{i}$ above, we obtain after summing

[TABLE]

Moreover, from (4.7), it follows that

[TABLE]

and thus

[TABLE]

Now, it is easy to see that

[TABLE]

and that $\|\bm{v}\|_{V}=\|\bm{w}-\bm{q}\|_{V}\lesssim{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\bm{w},\xi)\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}.$ ∎

Above and in the following we write $a\gtrsim b$ (or $a\lesssim b$) when $a\geq Cb$ (or $a\leq Cb$) for some positive constant $C$ independent of the finite element mesh.

To derive the discrete stability estimate, we need the following discrete trace inequality, easily shown by a scaling argument.

Lemma 4.1** (Discrete trace estimate).**

There exists $C_{I}>0$, independent of the mesh parameter $h$, such that

[TABLE]

Theorem 4.2** (Discrete stability).**

Suppose that $0<\alpha<C_{I}$. Then, for every $(\bm{w}_{h},\xi_{h})\in\bm{V}_{h}\times Q_{h}$, there exists $\bm{v}_{h}\in\bm{V}_{h}$ such that

[TABLE]

and

[TABLE]

Proof.

From the discrete trace estimate it follows that

[TABLE]

which proves the result in the mesh-dependent norm of $\xi_{h}$ for $0<\alpha<C_{I}$.

On the other hand, the continuous inf-sup condition (4.8) implies that for any $\xi_{h}\in Q_{h}$ there exists $\bm{v}\in\bm{V}$ such that

[TABLE]

This means that (cf. the proof of Lemma 3.2 in [12])

[TABLE]

where $C_{2},C_{3},C_{4}$ are positive constants and $I_{h}\bm{v}\in\bm{V}_{h}$ is the Clément interpolant of $\bm{v}$. Using again the discrete trace estimate and inequalities (4.12) and (4.13), we then obtain

[TABLE]

Now, it is straightforward to show (cf. [14]) that there exists $\delta>0$ such that

[TABLE]

and that $\|\bm{w}_{h}-\delta I_{h}\bm{v}\|_{V}\lesssim{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|(\bm{w}_{h},\xi_{h})\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{h}$. ∎

In our improved error analysis, we use techniques from the a posteriori error analysis. Let $\bm{f}_{i,h}\in\bm{V}_{i,h}$ be the $[L^{2}(\Omega_{i})]^{d}$ projection of $\bm{f}_{i}$, define on any $K\in\mathcal{C}_{h}^{i}$ the oscillation of $\bm{f}_{i}$ by

[TABLE]

and, for each $E\in\mathcal{G}_{h}^{i}$, let $K(E)\in\mathcal{G}_{h}^{i}$ denote the element such that $\partial K(E)\cap E=E$.

Lemma 4.2**.**

For any $(\bm{v}_{h},\eta_{h})\in\bm{V}_{h}\times Q_{h}$, it holds that

[TABLE]

Proof.

We follow the reasoning presented for the mortar method in [14]. It is clearly enough to prove the result in $\Omega_{1}$. Thus, let $b_{E}\in P_{d}(E)$, $E\in\mathcal{G}_{h}^{1}$, be the usual edge/facet bubble function and define $\tau_{E}$ on $K(E)\in\mathcal{C}_{h}^{1}$ through

[TABLE]

where $K(E)$ is such that $\overline{K(E)}\cap E=E$. It follows that

[TABLE]

Next, defining $\bm{\tau}\in\bm{V}_{1,h}$ in such a way that $\tau_{n}:=-\bm{\tau}\cdot\bm{n}=\sum_{E\in\mathcal{G}_{h}^{1}}\tau_{E}$ and testing problem (2.18) with $(\bm{v}_{1},\bm{v}_{2},\eta)=(-\bm{\tau},0,\lambda)$, we obtain

[TABLE]

Summing (4.15) over the edges in $\mathcal{G}_{h}^{1}$, gives then

[TABLE]

Inverse estimates imply that

[TABLE]

Now, one readily sees, using trace inequalities and the norm equivalence (4.2), that

[TABLE]

from which, using the standard estimates for interior residuals (cf. [26]) and the inverse estimate (4.16) to bound the last term, it follows that

[TABLE]

∎

We can now establish the quasi-optimality of the method.

Theorem 4.3**.**

For $0<\alpha<C_{I}$ it holds that

[TABLE]

Proof.

On account of the discrete stability estimate, there exists $\bm{w}_{h}\in\bm{V}_{h}$ such that

[TABLE]

and

[TABLE]

Using the bilinearity and (3.7), we obtain

[TABLE]

The terms above can be estimated as follows. First, continuity of the bilinear form $\mathcal{B}$ and inequality (4.18) yield

[TABLE]

Next, using the weak formulation (2.18) and the fact that \left\llbracket u_{n}\right\rrbracket\geq 0 and $\lambda_{h}\geq 0$, we obtain

[TABLE]

Finally, from the discrete trace estimate (4.9) it follows that

[TABLE]

Using Lemma 4.14, and collecting the above estimates, we arrive at the asserted error estimate. ∎

Remark 4.1**.**

We refrain from giving an a priori error estimate assuming a regular solution. The reasons are twofold. Firstly, contact singularities are inevitable and essential in contact problems. Secondly, to derive an a priori bound, one would need to estimate the term \sqrt{\langle\left\llbracket u_{n}\right\rrbracket,\eta_{h}\rangle}, with $\eta_{h}$ being the interpolant to $\lambda$. Besides, and perhaps most importantly, one of the main results of this paper is the fact that we do not need to assume that the solution belongs to $H^{s}$, with $s>3/2.$

For the a posteriori error analysis, we define the local estimators

[TABLE]

with $i=1,2$. The corresponding global estimator $\eta$ is then defined as

[TABLE]

In addition, we need an estimator $S$ defined only globally as

[TABLE]

Theorem 4.4** (A posteriori error estimate).**

It holds that

[TABLE]

Proof.

In view of the continuous stability estimate, there exists $\bm{v}\in\bm{V}$, with

[TABLE]

and

[TABLE]

Let $\tilde{\bm{v}}\in\bm{V}_{h}$ be the Clément interpolant of $\bm{v}$. From (3.7), it follows that

[TABLE]

Using the weak formulation (2.18), this gives

[TABLE]

Integrating by parts, we obtain for the first two terms above

[TABLE]

Moreover, using an inverse inequality for the $H^{1/2}(\Gamma)$-norm (cf. [11]) we get

[TABLE]

Finally, using the discrete trace estimate (4.9) and the standard bounds for the Clément interpolant, and recalling (4.31), we obtain for the stabilising term

[TABLE]

Estimate (4.38) follows from collecting the above bounds. ∎

The estimator $\eta$ bounds the error from below. For the proof of the following theorem we refer to [12].

Theorem 4.5** (A posteriori estimate – efficiency).**

It holds that

[TABLE]

The analysis of Methods 2 and 3 is analogous. In the a posteriori estimates the term

[TABLE]

is replaced by

[TABLE]

and

[TABLE]

for Method 2 and 3, respectively.

5. Computational experiments

All computations presented in this section were obtained using the Nitsche formulation 3 with the term (3.32) dropped. Had we considered other formulations, the results would have been practically identical. We also note that since the stabilized/Nitsche’s method is variationally conforming (as a mortaring method) it passes the patch test of [21], p. 425. This was confirmed numerically up to machine accuracy.

We consider the geometry given by

[TABLE]

and define the boundary conditions on the following subsets:

[TABLE]

Thus, the geometry is the one given in Figure 1. A nonmatching discretisation of the geometry is depicted in Figure 2. Initially, the material parameters are $E_{1}=E_{2}=1$ and $\nu_{1}=\nu_{2}=0.3$ and the loading is

[TABLE]

For this loading, the displacement is constrained on $\Gamma_{D,i}$, $i=1,2$, only in the horizontal direction which minimizes the effect of the singularities – other than the ones related to the contact boundary – on the rates of convergence. We consider both linear and quadratic elements, with $\alpha=10^{-2}$ and $\alpha=10^{-3}$, respectively.

The adaptively refined meshes are shown in Figure 3(a) and 3(b), and the global error estimator $\eta+S$ is plotted as a function of the number of degrees-of-freedom $N$ in Figure 3(c). Since $\eta+S$ is an upper bound for the total error, the results suggest that the total error of the quadratic solution is limited to $\mathcal{O}(N^{-0.5})$ when using uniform refinements and that adaptivity successfully improves the order of the discretisation error to $\mathcal{O}(N^{-1})$.

Next we fix also the vertical displacement on $\Gamma_{D,i}$, $i=1,2$, and consider the loading

[TABLE]

which causes the left block to bend slightly downwards and, as a consequence, the active contact region is a nontrivial subset of $\Gamma$. The active contact region is found via an iterative solution of the linearised problem, cf. [12]. See Figure 4(a) and 4(b) for the final meshes and contact stresses, and Figure 4(c) for the convergence rates. We observe that the singularity at the upper corner of the contact region is properly resolved by the adaptive meshing strategy and that the convergence is similar albeit less idealised as in the first example.

In Figure 5, we demostrate how the improved convergence rates can be obtained for $P_{2}$ elements even if the value of the Young’s modulus changes significantly over the contact boundary. In Figure 6, we demonstrate that the effect of the stabilisation parameter is small in the asymptotic limit. Finally, in Figure 7, we consider the loading

[TABLE]

which results in an active contact boundary consisting of two disjoint parts and a perfectly symmetric contact stress.