Weak imposition of Signorini boundary conditions on the boundary element method

TL;DR

This paper develops a boundary element method for inequality boundary conditions, specifically Signorini contact conditions, using a variational approach and provides a comprehensive error analysis with numerical validation.

Contribution

It introduces a novel weak imposition technique for Signorini boundary conditions within the boundary element framework, combining Calderón projector and augmented Lagrangian methods.

Findings

Complete a priori error analysis conducted.

Numerical examples validate the theoretical results.

Method effectively handles inequality boundary conditions.

Abstract

We derive and analyse a boundary element formulation for boundary conditions involving inequalities. In particular, we focus on Signorini contact conditions. The Calder\'on projector is used for the system matrix and boundary conditions are weakly imposed using a particular variational boundary operator designed using techniques from augmented Lagrangian methods. We present a complete numerical a priori error analysis and present some numerical examples to illustrate the theory.

Equations

$$-\Delta u$$

$$u$$

$$u⩽g_{\text{C}}\text{and}\frac{\partial u}{\partial \mathbf{\nu }}$$

$$\left(\frac{\partial u}{\partial \mathbf{\nu }}-{\psi }_{\text{C}}\left)\right(u-g_{\text{C}}\right)$$

$$G(\mathbf{x},\mathbf{y})=\frac{1}{4\text{textpi}|\mathbf{x}-\mathbf{y}|}.$$

$$(\mathcal{V}\mu )(\mathbf{x})$$

$$(\mathcal{K}v)(\mathbf{x})$$

$${\gamma }_{\text{D}}f(\mathbf{x})$$

$${\gamma }_{\text{N}}f(\mathbf{x})$$

$$u=-\mathcal{K}({\gamma }_{\text{D}}u)+\mathcal{V}({\gamma }_{\text{N}}u).$$

$$u_{\mu }^{\mathcal{V}}:=\mathcal{V}\mu$$

$$\Vert u_{\mu }^{\mathcal{V}}{\Vert }_{H^1(\Omega )}⩽c\Vert \mu {\Vert }_{H^{-1/2}(\Gamma )}.$$

$$u_v^{\mathcal{K}}:=\mathcal{K}v$$

$$\Vert u_v^{\mathcal{K}}{\Vert }_{H^1(\Omega )}⩽c\Vert v{\Vert }_{H^{1/2}(\Gamma )}.$$

$$(\mathsf{K}v)(\mathbf{x})$$

$$(\mathsf{W}v)(\mathbf{x})$$

$$\mathsf{C}:=\left(\begin{array}{cc}(1-\sigma )\mathsf{Id}-\mathsf{K}& \mathsf{V}\\ \mathsf{W}& \sigma \mathsf{Id}+{\mathsf{K}}^{\prime }\end{array}\right),$$

$$\sigma (\mathbf{x})=\underset{ϵ\to 0}{\lim }\frac{1}{4\text{textpi}{ϵ}^2}{\int }_{\mathbf{y}\in \Omega :|\mathbf{y}-\mathbf{x}|=ϵ}\mathrm{d}\mathbf{y}.$$

$$\mathsf{C}\left(\begin{array}{c}{\gamma }_{\text{D}}u\\ {\gamma }_{\text{N}}u\end{array}\right)=\left(\begin{array}{c}{\gamma }_{\text{D}}u\\ {\gamma }_{\text{N}}u\end{array}\right).$$

$${⟨{\gamma }_{\text{D}}u,\mu ⟩}_{\Gamma }$$

$${⟨{\gamma }_{\text{N}}u,v⟩}_{\Gamma }$$

$$\begin{gathered} \mathcal{C}[(u,\lambda ),(v,\mu )]:={⟨(\frac{1}{2}\mathsf{Id}-\mathsf{K})u,\mu ⟩}_{\Gamma }+{⟨\mathsf{V}\lambda ,\mu ⟩}_{\Gamma } \\ +{⟨(\frac{1}{2}\mathsf{Id}+{\mathsf{K}}^{\prime })\lambda ,v⟩}_{\Gamma }+{⟨\mathsf{W}u,v⟩}_{\Gamma }. \end{gathered}$$

$$\mathcal{C}[(u,\lambda ),(v,\mu )]={⟨u,\mu ⟩}_{\Gamma }+{⟨\lambda ,v⟩}_{\Gamma }.$$

$$\mathcal{A}[(u,\lambda ),(v,\mu )]:=-{⟨\mathsf{K}u,\mu ⟩}_{\Gamma }+{⟨\mathsf{V}\lambda ,\mu ⟩}_{\Gamma }+{⟨{\mathsf{K}}^{\prime }\lambda ,v⟩}_{\Gamma }+{⟨\mathsf{W}u,v⟩}_{\Gamma }.$$

$$\mathcal{A}[(u,\lambda ),(v,\mu )]=\frac{1}{2}{⟨u,\mu ⟩}_{\Gamma }+\frac{1}{2}{⟨\lambda ,v⟩}_{\Gamma }.$$

$$\mathbb{V}:=H^{1/2}(\Gamma )\times H^{-1/2}(\Gamma )$$

$$\Vert (v,\mu ){\Vert }_{\mathbb{V}}:=\Vert v{\Vert }_{H^{1/2}(\Gamma )}+\Vert \mu {\Vert }_{H^{-1/2}(\Gamma )}.$$

$$\left|\mathcal{A}[(w,\eta ),(v,\mu )]\right|$$

$$\alpha \left(|v{|}_{H_{\ast }^{1/2}(\Gamma )}^2+\Vert \mu {\Vert }_{H^{-1/2}(\Gamma )}^2\right)$$

$${\text{P}}_h^k(\Gamma )$$

$${\text{DP}}_h^l(\Gamma )$$

$${\widetilde{\text{DP}}}_h^l(\Gamma )$$

$$\underset{{\eta }_h\in {\text{DUAL}}_h^0(\Gamma )}{\inf }\Vert \mu -{\eta }_h{\Vert }_{H^{-1/2}(\Gamma )}$$

$${\mathbb{V}}_h:={\text{P}}_h^k(\Gamma )\times {\Lambda }_h^l,$$

$$R_{{\Gamma }_{\text{D}}}(u_h,{\lambda }_h):={\beta }_{\text{D}}^{1/2}(g_{\text{D}}-u_h).$$

Full text

{AMS}

65N38, 65R20, 74M15

1 Introduction

The application of Nitsche techniques to deal with variational inequalities has received increasing interest recently, starting from a series of works by Chouly, Hild and Renard for elasticity problems with contact [7]. Their approach goes back to an augmented Lagrangian formulation, that has first been introduced by Alart & Curnier [1].

In a previous paper [2], we have shown how Nitsche techniques can be used to impose Dirichlet, Neumann, mixed Dirichlet–Neumann or Robin conditions weakly within boundary element methods. By using the Calderón projector, we were able to derive a unified framework that can be used for different boundary conditions.

The purpose of this article is to extend these techniques to boundary conditions involving inequalities, such as Signorini contact conditions. In particular, we consider the Laplace equation with mixed Dirichlet and Signorini boundary conditions: Find $u$ such that

[TABLE]

Here $\Omega\subset\mathbb{R}^{3}$ denotes a polyhedral domain with outward pointing normal $\bm{\nu}$ and boundary $\Gamma:=\Gamma_{\textup{D}}\cup\Gamma_{\textup{C}}$. We assume for simplicity that the boundary between $\Gamma_{\textup{D}}$ and $\Gamma_{\textup{C}}$ coincides with edges between the faces of $\Gamma$. Whenever it is ambiguous, we will write $\bm{\nu}_{\bm{x}}$ for the outward pointing normal at the point ${\bm{x}}$. We assume that $g=\begin{cases}g_{\textup{D}}&\text{in }\Gamma_{\textup{D}}\\ g_{\textup{C}}&\text{in }\Gamma_{\textup{C}}\end{cases}\in L^{2}(\Gamma)$ and $\psi_{\textup{C}}\in H^{1/2}(\Gamma_{\textup{C}})$.

Observe that when $\Gamma_{\textup{C}}=\varnothing$, there exists a unique solution to 1 by the Lax–Milgram lemma. In the case that meas($\Gamma_{\textup{C}})>0$, the theory of Lions and Stampacchia [12] for variational inequalities yields existence and uniqueness of solutions. We assume that $u\in H^{3/2+\epsilon}(\Omega)$, for some $\epsilon>0$.

Boundary element methods for Signorini problems were first studied by Han [11]. A variational formulation involving the Calderón projector was presented in [10]. An alternative formulation is based on Steklov-Poincaré operators [20, 22]. The numerical approaches to solve such formulations include a penalty formulation [15], operator splitting techniques [17, 23] or semi-smooth Newton methods [20, 22]. The latter reference includes besides the usual energy norm estimates an $L^{2}(\Gamma)$-error estimate based on a duality argument. Maischak & Stephan [13] presented a posteriori error estimates and an $hp$-adaptive algorithm for the Signorini problem. A priori error estimates for a penalty-based $hp$ algorithm were shown by Chernov, Maischak & Stephan [6]. Recently, an augmented Lagrangian approach has been presented in combination with a semi-smooth Newton method [22], and variational inequalities have been successfully used for time-dependent contact problems [9].

We will consider an approach where the full Calderón projector is used and the boundary conditions are included by adding properly scaled penalty terms to the two equations. This results in formulations similar to the ones obtained for weak imposition of boundary conditions using Nitsche’s method [14]. The proposed framework is flexible and allows for the design of a range of different methods depending on the choice of weights and residuals.

An outline of the paper is as follows. In Section 2, we introduce the basic boundary operators that will be needed and review some of their properties. Then, in Section 3, we introduce the variational framework and review the results from [2] for the pure Dirichlet problem. In Section 4, we show how the framework can be applied to Signorini boundary conditions and the mixed problem 1. The method is analysed in Section 5. We conclude by showing some numerical experiments in Section 6.

2 Boundary operators

We define the Green’s function for the Laplace operator in $\mathbb{R}^{3}$ by

[TABLE]

In this paper, we focus on the problem in $\mathbb{R}^{3}$. Similar analysis can be used for problems in $\mathbb{R}^{2}$, in which case this definition should be replaced by $G({\bm{x}},{\bm{y}})=-\log|{\bm{x}}-{\bm{y}}|/2\text{\textpi}$.

In the standard fashion (see e.g. [19, Chapter 6]), we define the single layer potential operator, $\mathcal{V}:H^{-1/2}(\Gamma)\to H^{1}(\Omega)$, and the double layer potential, $\mathcal{K}:H^{1/2}(\Gamma)\to H^{1}(\Omega)$, for $v\in H^{1/2}(\Gamma)$, $\mu\in H^{-1/2}(\Gamma)$, and ${\bm{x}}\in\Omega\setminus\Gamma$ by

[TABLE]

We define the space $H^{1}(\Delta,\Omega):=\{v\in H^{1}(\Omega):\Delta v\in L^{2}(\Omega)\}$, and the Dirichlet and Neumann traces, $\gamma_{\textup{D}}:H^{1}(\Omega)\to H^{1/2}(\Gamma)$ and $\gamma_{\textup{N}}:H^{1}(\Delta,\Omega)\to H^{-1/2}(\Gamma)$, by

[TABLE]

We recall that if the Dirichlet and Neumann traces of a harmonic function are known, then the potentials 3 and 4 may be used to reconstruct the function in $\Omega$ using the following relation.

[TABLE]

It is also known [19, Lemma 6.6] that for all $\mu\in H^{-1/2}(\Gamma)$, the function

[TABLE]

satisfies $-\Delta u^{\mathcal{V}}_{\mu}=0$ and

[TABLE]

Similarly [19, Lemma 6.10], the function

[TABLE]

satisfies $-\Delta u^{\mathcal{K}}_{v}=0$ for all $v\in H^{1/2}(\Gamma)$ and

[TABLE]

We define $\left\{\gamma_{\textup{D}}f\right\}_{\Gamma}$ and $\left\{\gamma_{\textup{N}}f\right\}_{\Gamma}$ to be the averages of the interior and exterior Dirichlet and Neumann traces of $f$. We define the single layer, double layer, adjoint double layer, and hypersingular boundary integral operators, $\mathsf{V}:H^{-1/2}(\Gamma)\to H^{1/2}(\Gamma)$, $\mathsf{K}:H^{1/2}(\Gamma)\to H^{1/2}(\Gamma)$, ${\mathsf{K}^{\prime}}:H^{-1/2}(\Gamma)\to H^{-1/2}(\Gamma)$, and $\mathsf{W}:H^{1/2}(\Gamma)\to H^{-1/2}(\Gamma)$, by

[TABLE]

where ${\bm{x}}\in\Gamma$, $v\in H^{1/2}(\Gamma)$ and $\mu\in H^{-1/2}(\Gamma)$ [19, Chapter 6].

Next, we define the Calderón projector by

[TABLE]

where $\sigma$ is defined for ${\bm{x}}\in\Gamma$ by [19, Equation 6.11]

[TABLE]

Recall that if $u$ is a solution of 1 then it satisfies

[TABLE]

Taking the product of 15 with two test functions, and using the fact that $\sigma=\frac{1}{2}$ almost everywhere, we arrive at the following equations.

[TABLE]

For a more compact notation, we introduce $\lambda=\gamma_{\textup{N}}u$ and $u=\gamma_{\textup{D}}u$ and the Calderón form

[TABLE]

We may then rewrite 16 and 17 as

[TABLE]

We will also frequently use the multitrace form, defined by

[TABLE]

Using this, we may rewrite 19 as

[TABLE]

To quantify the two traces we introduce the product space

[TABLE]

and the associated norm

[TABLE]

The continuity and coercivity of $\mathcal{A}$ are immediate consequences of the properties of the operators $\mathsf{V}$, $\mathsf{K}$, ${\mathsf{K}^{\prime}}$ and $\mathsf{W}$:

Lemma 2.1** (Continuity & coercivity).**

There exists $C>0$ such that

[TABLE]

There exists $\alpha>0$ such that

[TABLE]

Proof 2.2**.**

See [2].

3 Discretisation and weak imposition of Dirichlet boundary conditions

In this section, we introduce the discrete spaces and review briefly how (non-homogeneous) Dirichlet boundary conditions can be imposed weakly within the variational formulations introduced above. For a detailed derivation, and for different boundary conditions, we refer to [2].

To reduce the number of constants that appear, we introduce the following notation.

• •

If $\exists C>0$ such that $a\leqslant Cb$, then we write $a\lesssim b$.

• •

If $a\lesssim b$ and $b\lesssim a$, then we write $a\eqsim b$.

We assume that $\Omega$ is a polygonal domain with faces denoted by $\{\Gamma_{i}\}_{i=1}^{M}$. We introduce a family of conforming, shape regular triangulations of $\Gamma$, $\{\mathcal{T}_{h}\}_{h>0}$, indexed by the largest element diameter of the mesh, $h$. We let $T_{1},..,T_{m}\in\mathcal{T}_{h}$ be the triangles of a triangulation.

We consider the following finite element spaces

[TABLE]

where $\mathbb{P}_{k}(T_{i})$ denotes the space of polynomials of order less than or equal to $k$ on the triangle $T_{i}$.

In addition, we consider the space $\textup{DUAL}_{h}^{0}(\Gamma)$ of piecewise constant functions on the barycentric dual grid, as shown in Figure 1. On non-smooth domains, these spaces have lower order approximation properties than the standard space $\textup{DP}_{h}^{0}(\Gamma)$, as given in the following lemma.

Lemma 3.1**.**

Let $\mu\in H^{s}(\Gamma)$. If $\Gamma$ consists of a finite number of smooth faces meeting at edges, then

[TABLE]

*where $\xi=\min(\tfrac{1}{2},s)$. If $\Gamma$ is smooth, then the same result holds with $\xi=\min(1,s)$. *

Proof 3.2**.**

*See [16, Appendix 2]. *

We observe that $\textup{P}^{k}_{h}(\Gamma)\subset H^{1/2}(\Gamma)$, $\textup{DP}^{l}_{h}(\Gamma)\subset L^{2}(\Gamma)$, $\widetilde{\textup{DP}}^{l}_{h}(\Gamma)\subset L^{2}(\Gamma)$, and $\textup{DUAL}_{h}^{0}(\Gamma)\subset L^{2}(\Gamma)$. We define the discrete product space

[TABLE]

where $\Lambda_{h}^{l}$ can be any of the spaces $\textup{DP}^{l}_{h}(\Gamma),\widetilde{\textup{DP}}^{l}_{h}(\Gamma)$ or $\textup{DUAL}_{h}^{0}(\Gamma)$.

3.1 Dirichlet boundary conditions

Let us for the moment assume that $\Gamma\equiv\Gamma_{\textup{D}}$. Then, the basic idea is to add the following suitably weighted boundary residual to the weak formulation.

[TABLE]

This is defined such that $R_{\Gamma_{\textup{D}}}(u_{h},\lambda_{h})=0$ is equivalent to the boundary condition 1b. We obtain an expression of the form

[TABLE]

or equivalently

[TABLE]

where $\beta_{1}$ and $\beta_{2}$ are problem dependent scaling operators that can be chosen as a function of the physical parameters in order to obtain robustness of the method.

For the Dirichlet problem, we choose $\beta_{1}=\beta_{\textup{D}}^{1/2}$, $\beta_{2}=\beta_{\textup{D}}^{-1/2}$, where different choices for $\beta_{\textup{D}}$ in the range $0\leqslant\beta_{\textup{D}}\lesssim h^{-1}$ are possible. Inserting this into 24, we obtain the formulation:

[TABLE]

By formally identifying $\lambda_{h}$ with $\partial_{\nu}u_{h}$ and $\mu_{h}$ with $\partial_{\nu}v_{h}$, we obtain the classical (non-symmetric) Nitsche’s method (up to the multiplicative factor $\tfrac{1}{2}$).

For a more compact notation, we introduce the boundary operator associated with the non-homogeneous Dirichlet condition

[TABLE]

the operator corresponding to the left-hand side

[TABLE]

and the operator associated with the right-hand side

[TABLE]

Using these and 25, we arrive at the following boundary element formulation: Find $(u_{h},\lambda_{h})\in\mathbb{V}_{h}$ such that

[TABLE]

We introduce the following $\mathcal{B}_{\textup{D}}$-norm

[TABLE]

and summarise the properties of the bilinear form $\mathcal{A}_{\textup{D}}$ in the following lemma.

Lemma 3.3** (Properties of the bilinear form).**

Let $\mathbb{W}$ be a product Hilbert space for the primal and flux variables, such that $\mathbb{V}\subset\mathbb{W}$. The bilinear form has the following properties:

Property 1 (Coercivity):* If $\beta_{\textup{D}}=0$ or if there exists $\beta_{\min}>0$ (independent of $h$) such that $\beta_{\textup{D}}>\beta_{\min}$, then there exists $\alpha>0$ such that $\forall(v,\mu)\in\mathbb{W}$*

[TABLE]

Property 2 (Continuity):* There exists $M>0$ such that*

$\forall(w,\eta),(v,\mu)\in\mathbb{W}$**

[TABLE]

Proof 3.4**.**

*See [2, Section 4.1]. *

4 Weak imposition of Signorini boundary conditions

Recently Chouly, Hild and Renard [7, 8] showed how contact problems can be treated in the context of Nitsche’s method. We will here show how we may use arguments similar to theirs in the present framework to integrate unilateral contact seamlessly. The result is a nonlinear system to which one may apply Newton’s method or a fixed-point iteration in a straightforward manner. We prove existence and uniqueness of solutions to the nonlinear system and optimal order error estimates.

For the derivation of the formulation on the contact boundary we will first omit the Dirichlet part, letting $\Gamma=\Gamma_{\textup{C}}$. To impose the contact conditions, we recall the following relations, introduced by Alart and Curnier [1], with $\left[x\right]_{\pm}:=\pm\max(0,\pm x)$.

[TABLE]

for all $\tau>0$. It is straighforward [7] to show that each of these two conditions is equivalent to the contact boundary conditions 1c and 1d.

To simplify the notation, we introduce the operators

[TABLE]

Using 30, we arrive at the following boundary term for the contact conditions

[TABLE]

Alternatively, by using 31, we arrive at the following boundary term

[TABLE]

By using the fact that $x=\left[x\right]_{+}+\left[x\right]_{-}$, it can be shown that 32 and 33 are equal.

Substituting 32 into 24, and using the weights $\beta_{1}=\tau$ and $\beta_{2}=1$, we obtain

[TABLE]

Using 33, we have

[TABLE]

We see that 35 is similar to the non-symmetric version of the method proposed in [8] and 34 is similar to the non-symmetric Nitsche formulation for contact discussed in [5]. As pointed out in the latter reference, the two formulations are equivalent, with the same solutions. In what follows, we focus exclusively on the variant 35.

Defining

[TABLE]

we arrive at the boundary element method formulation: Find $(u_{h},\lambda_{h})\in\mathbb{V}_{h}$ such that

[TABLE]

4.1 Mixed Dirichlet and contact boundary conditions

Combining the formulations for the Dirichlet and contact conditions, we arrive at the following boundary element method for the problem 1: Find $(u_{h},\lambda_{h})\in\mathbb{V}_{h}$ such that

[TABLE]

where $\mathcal{A}_{\textup{D}}$, $\mathcal{L}_{\textup{D}}$, $\mathcal{B}_{\textup{C}}$ and $\mathcal{L}_{\textup{C}}$ are defined in 27, 36, 28, and 37. For discretisation, we use the assumptions and spaces introduced in Section 3. Note that the formulation (40) is consistent, i.e. the continuous solution $(u,\lambda)$ to (1) fulfills (40) for all $(v_{h},\mu_{h})\in\mathbb{V}_{h}$.

5 Analysis

In this section, we prove the existence of unique solutions to the nonlinear system of equations 40 as well as optimal error estimates.

We assume that the solution $(u,\lambda)$ of 1 lies in $\mathbb{W}:=H^{1+\epsilon}(\Gamma)\times H^{\epsilon}(\tilde{\Gamma})$ for some $\epsilon\in(0,1/2]$, where $\tilde{\Gamma}=\cup_{i=1}^{M}\Gamma_{i}\setminus\partial\Gamma_{i}$ is the set of boundary points that lie in the interior of the faces $\Gamma_{i}$. As the normal vectors $\bm{\nu}_{{\bm{x}}}$ are discontinuous between faces, we can not expect a higher global regularity for $\lambda$.

We define the distance function $d_{\textup{C}}$ and norm $\|\cdot\|_{*}$, for $(v,\mu),(w,\eta)\in\mathbb{W}$, by

[TABLE]

We note that due the appearance of $\left[\cdot\right]_{+}$ in its second term, $d_{\textup{C}}$ is not a norm. $d_{\textup{C}}$ does provide a bound on the error however, as for all $(v,\mu)\in\mathbb{W}$, $d_{\textup{C}}\left((v,\mu),(0,0)\right)\geqslant\|(v,\mu)\|_{\mathcal{B}_{\textup{D}}}\geqslant\|(v,\mu)\|_{\mathbb{V}}$.

When proving this section’s results, we will use properties of the $\left[\cdot\right]_{+}$ function that are given in the following lemma.

Lemma 5.1**.**

For all $a,b\in\mathbb{R}$,

[TABLE]

Proof 5.2**.**

*For a proof of these well-known properties see e.g. [7]. *

We now prove a result analogous to the coercivity assumption in [2].

Lemma 5.3**.**

If there is $\beta_{\min}>0$, independent of $h$, such that $\beta_{\textup{D}}>\beta_{\min}$, then there is $\alpha>0$ such that for all $(v,\mu),(w,\eta)\in\mathbb{W}$,

[TABLE]

Proof 5.4**.**

From the analysis of the Dirichlet problem (Lemma 3.3) we know that when $\beta_{\textup{D}}>\beta_{\min}>0$,

[TABLE]

Introducing the notation $\delta P:=\left[P^{\tau}(v,\mu)\right]_{+}-\left[P^{\tau}(w,\eta)\right]_{+}$, we have

[TABLE]

To estimate the expression on the right-hand side, we use

[TABLE]

Using 43, this implies the bound

[TABLE]

Observing that $P_{0}^{\tau}(v-w,\mu-\eta)+2(\mu-\eta)=\tau(v-w)+\mu-\eta$, we infer that

[TABLE]

We conclude the proof by noting that

[TABLE]

*and applying 45 and 47. *

Next, we prove a result analagous to the discrete coercivity assumption in [2].

Lemma 5.5**.**

If there is $\beta_{\min}>0$, independent of $h$, such that $\beta_{\textup{D}}>\beta_{\min}$, then there is $\alpha>0$ such that for all $(v_{h},\mu_{h})\in\mathbb{V}_{h}$,

[TABLE]

Proof 5.6**.**

The proof is similar to that of Lemma 5.3, but with $\mu_{h}$ and $v_{h}$ instead of $\mu-\eta$ and $v-w$. The appearance of the data term in the right-hand side is due to the relation

[TABLE]

Using Lemmas 5.3 and 5.5, we may now prove that 40 is well-posed.

Theorem 5.7**.**

*The finite dimensional nonlinear system 40 admits a unique solution. *

Proof 5.8**.**

To prove the existence of a solution, we show the continuity and the positivity of the nonlinear operator $\mathcal{A}+\mathcal{B}_{\textup{D}}+\mathcal{B}_{\textup{C}}$. This allows us to apply Brouwer’s fixed point theorem, see eg [21, Chapter 2, Lemma 1.4].

We define $\mathsf{F}:\mathbb{V}_{h}\to\mathbb{V}_{h}$, for $(v_{h},\mu_{h})\in\mathbb{V}_{h}$, by

[TABLE]

for all $(w_{h},\eta_{h})\in\mathbb{V}_{h}$. We may write the non-linear system 40 as

[TABLE]

For fixed $h$, by the equivalance of norms on discrete spaces, there exist $c_{1},c_{2}>0$ such that for all $(v_{h},\mu_{h})\in\mathbb{V}_{h}$,

[TABLE]

To show positivity, we let $(v_{h},\mu_{h})\in\mathbb{V}_{h}$. Using Lemma 5.5, we see that

[TABLE]

Using the Cauchy–Schwarz inequality and an arithmetic-geometric inequality, we see that there exists $C_{g_{\textup{C}}\psi_{\textup{C}}}>0$ such that

[TABLE]

Using norm equivalence, we obtain

[TABLE]

for some $C^{\prime}>0$. We conclude that for all $(v_{h},\mu_{h})\in\mathbb{V}_{h}$ with

[TABLE]

there holds $\left\langle\mathsf{F}(v_{h},\mu_{h}),(v_{h},\mu_{h})\right\rangle_{\Gamma}>0$.

To show continuity, let $(v^{1}_{h},\mu^{1}_{h}),(v^{2}_{h},\mu^{2}_{h})\in\mathbb{V}_{h}$. We have for all $(w_{h},\eta_{h})\in\mathbb{V}_{h}$,

[TABLE]

where we have used 44. By norm equivalence, this means that

[TABLE]

showing that $\mathsf{F}$ is continuous.

It then follows by Brouwer’s fixed point theorem [21, Chapter 2, Lemma 1.4] that there exists a solution to 48 and hence also to 40.

Uniqueness is an immediate consequence of Lemma 5.3. Assume that $(u_{h}^{1},\lambda_{h}^{1})$ and $(u_{h}^{2},\lambda_{h}^{2})$ are solutions to 40. We immediately see that

[TABLE]

*and we conclude that the solution is unique. *

We now proceed to prove the following best approximation result.

Lemma 5.9**.**

Let $(u,\lambda)\in\mathbb{W}$ be the solution of 1 and $(u_{h},\lambda_{h})\in\mathbb{V}_{h}$ the solution of 40. Then there holds

[TABLE]

Proof 5.10**.**

Using Lemma 5.3 and Galerkin orthogonality, we see that, for arbitrary $(v_{h},\mu_{h})\in\mathbb{V}_{h}$,

[TABLE]

Next, we use

[TABLE]

to show that

[TABLE]

We estimate the three parts of the right-hand separately. For the first term, we use the continuity of $\mathcal{A}+\mathcal{B}_{\textup{D}}$ (Lemma 3.3) to obtain

[TABLE]

For the second line, we use $H^{1/2}(\Gamma)$–$H^{-1/2}(\Gamma)$ duality and the Cauchy–Schwarz inequality to obtain

[TABLE]

For the last term, we use the Cauchy–Schwarz inequality to get

[TABLE]

Collecting these bounds, we see that

[TABLE]

*Dividing through by $d_{\textup{C}}\left((u,\lambda),(u_{h},\lambda_{h})\right)$, and taking the infimum yields the desired result. *

We now prove the main result of this section, an a priori bound on the error of the solution of 40.

Theorem 5.11**.**

Let $(u,\lambda)\in H^{s}(\Gamma)\times H^{r}(\tilde{\Gamma})$ for some $s\geqslant 1,r\geqslant 0$ and $(u_{h},\lambda_{h})\in\textup{P}_{h}^{k}(\Gamma)\times\Lambda_{h}^{l}$ be the solutions of 1 and the discrete problem 40, respectively. If there is $\beta_{\min}>0$ such that $\beta_{\min}<\beta_{\textup{D}}\lesssim h^{-1}$ and $\tau\eqsim h^{-1}$, then

[TABLE]

where $\zeta=\min(k+1,s)$ and $\xi=\min(l+1,r)$ for $\Lambda_{h}^{l}\in\{\textup{DP}^{l}_{h}(\Gamma),\widetilde{\textup{DP}}^{l}_{h}(\Gamma)\}$ and $\zeta=\min(2,s)$ and $\xi=\min(\frac{1}{2},r)$ for $\Lambda_{h}^{l}=\textup{DUAL}^{0}_{h}(\Gamma)$. Additionally,

[TABLE]

*where $\tilde{u}$ and $\tilde{u}_{h}$ are the solutions in $\Omega$ defined by 7. *

Proof 5.12**.**

First, we observe that for all $(v,\mu)$ and $(w,\eta)$ in $\mathbb{W}$

[TABLE]

Using standard approximation results for $\Lambda_{h}^{l}\in\{\textup{DP}^{l}_{h}(\Gamma),\widetilde{\textup{DP}}^{l}_{h}(\Gamma)\}$ (see eg [19, chapter 10]) and Lemma 3.1 for $\Lambda_{h}^{l}=\textup{DUAL}^{0}_{h}(\Gamma)$, we see that

[TABLE]

Applying these to the definition of $\|\cdot\|_{*}$ gives

[TABLE]

*By means of Lemma 5.9 and the given choice of the parameters $\tau$ and $\beta_{\textup{D}}$ this proves the first assertion. The estimate in the domain $\Omega$ follows by using the relations 9 and 11. *

If $\lambda$ is smooth enough and $k=l$, the bounds on $\tau$ can be replaced with $h\lesssim\tau\lesssim h^{-1}$ without reducing the order of convergence.

6 Numerical results

We now demonstrate the theory with a series of numerical examples. In this section, we consider the following test problem. Let $\Omega=[0,1]\times[0,1]\times[0,1]$ be the unit cube, $\Gamma_{\textup{C}}:=\{(x,y,z)\in\Gamma:z=1\}$, and $\Gamma_{\textup{D}}:=\Gamma\setminus\Gamma_{\textup{C}}$. Let

[TABLE]

It can be shown that

[TABLE]

is the solution to 1 with these boundary conditions.

To solve the non-linear system 39, we will treat the nonlinear term explicitly. Therefore, we define

[TABLE]

Note that $\mathcal{B}^{\prime}_{\textup{C}}$ differs from $\mathcal{B}_{\textup{C}}$ only by the missing nonlinear term.

We pick initial guesses $(u_{0},\lambda_{0})\in\mathbb{V}_{h}$ and define $(u_{n+1},\lambda_{n+1})\in\mathbb{V}_{h}$, for $n\in\mathbb{N}$, to be the solution of

[TABLE]

This leads us to Algorithm 1, an iterative method for solving the contact problem.

In all the computations in this section, we preconditioned the GMRES solver using a mass matrix preconditioner applied blockwise from the left, as described in [3].

Inspired by the parameter choices in [2], we fix $\beta_{\textup{D}}=0.01$ and look for suitable values of the parameter $\tau$. Figure 2 shows how the error, number of outer iterations, and the average number of GMRES iterations inside each outer iteration change as the parameter $\tau$ is varied, for both $\mathbb{V}_{h}=\textup{P}_{h}^{1}(\Gamma)\times\textup{DUAL}_{h}^{0}(\Gamma)$ (left, blue) and $\mathbb{V}_{h}=\textup{P}_{h}^{1}(\Gamma)\times\textup{DP}_{h}^{0}(\Gamma)$ (right, orange). Here, we see that the error and number of outer iterations are lowest when $\tau$ is between around 1 and 10.

Motivated by Figure 2 and the bounds in Theorem 5.11, we take $\tau=0.5/h$, and look at the convergence as $h$ is decreased. Figure 3 shows how the error and iteration counts vary as $h$ is decreased when $\mathbb{V}_{h}=\textup{P}_{h}^{1}(\Gamma)\times\textup{DUAL}_{h}^{0}(\Gamma)$ (left, blue circles) and $\mathbb{V}_{h}=\textup{P}_{h}^{1}(\Gamma)\times\textup{DP}_{h}^{0}(\Gamma)$ (right, orange squares).

For $\mathbb{V}_{h}=\textup{P}_{h}^{1}(\Gamma)\times\textup{DUAL}_{h}^{0}(\Gamma)$, we observe slightly higher than the order 1 convergence predicted by Theorem 5.11. In this case, the mass matrix preconditioner is effective, as the number of GMRES iterations required inside each outer iteration is reasonably low, and only grows slowly as $h$ is decreased. We believe that the effectiveness of the preconditioner for this choice of spaces is due to the spaces $\textup{P}_{h}^{1}(\Gamma)$ and $\textup{DUAL}_{h}^{0}(\Gamma)$ forming an inf-sup stable pair [18, Lemma 3.1].

When $\mathbb{V}_{h}=\textup{P}_{h}^{1}(\Gamma)\times\textup{DP}_{h}^{0}(\Gamma)$, Theorem 5.11 tells us to expect order 1.5 convergence. However, we observe a slightly lower order. This appears to be due to the ill-conditioning of this system, and the mass matrix preconditioner being ineffective, leading to an inaccurate solution when using GMRES. In this case, the spaces $\textup{P}_{h}^{1}(\Gamma)$ and $\textup{DP}_{h}^{0}(\Gamma)$ do not form an inf-sup stable pair, and so the mass-matrix between them is not guaranteed to be invertible leading to a less effective preconditioner.

In order to obtain order 1.5 convergence with a well-conditioned system, we could look for $(u_{h},\lambda_{h})\in\textup{P}_{h}^{1}(\Gamma)\times\textup{DP}_{h}^{0}(\Gamma)$ and test with $(v_{h},\mu_{h})\in\textup{DUAL}_{h}^{1}(\Gamma)\times\textup{DUAL}_{h}^{0}(\Gamma)$, where $\textup{DUAL}_{h}^{1}(\Gamma)$ is the space of piecewise linear functions on the dual grid that forms an inf-sup stable pair with the space $\textup{DP}_{h}^{0}(\Gamma)$, as defined in [4]. With this choice of spaces, we obtain the higher order convergence as in Theorem 5.11, while having stable dual pairings and hence more effective mass matrix preconditioning.

For the problems discussed in [2], we have run numerical experiments using this space pairing and observe the full order $\tfrac{3}{2}$ convergence in a low number of iterations. A deeper investigation of this method using these dual spaces, and the adaption of the theory to this case, warrants future work.

7 Conclusions

Based on our work in [2], we have analysed and demonstrated the effectiveness of Nitsche type coupling methods for boundary element formulations of contact problems.

An open problem is preconditioning. While the iteration counts in the presented examples were already practically useful, for large and complex structures preconditioning is still essential. The hope is to use the properties of the Calderón projector to build effective operator preconditioning techniques for the presented Nitsche type frameworks.

Avenues of future research include looking at how this approach could be applied to problems in linear elasticity, and an extension of this method to problems involving friction.