Galerkin least squares finite element method for the obstacle problem

Galerkin least squares finite element method for the obstacle problem

Erik Burman, Peter Hansbo, Mats G. Larson, Rolf Stenberg Department of Mathematics, University College London, London, UK–WC1E 6BT, UK Department of Mechanical Engineering, Jönköping University, SE-55111 Jönköping, Sweden Department of Mathematics and Mathematical Statistics, Umeå University, SE-901 87 Umeå, Sweden Institute of Mathematics, Helsinki University of Technology, P. O. Box 1100, 02015 TKK, Finland
Abstract

We construct a consistent multiplier free method for the finite element solution of the obstacle problem. The method is based on an augmented Lagrangian formulation in which we eliminate the multiplier by use of its definition in a discrete setting. We prove existence and uniqueness of discrete solutions and optimal order a priori error estimates for smooth exact solutions. Using a saturation assumption we also prove an a posteriori error estimate. Numerical examples show the performance of the method and of an adaptive algorithm for the control of the discretization error.

keywords:
Obstacle problem, augmented Lagrangian method, a priori error estimate, a posteriori error estimate, adaptive method

1 Introduction

Our aim in this paper is to design a simple consistent penalty method for contact problems that avoids the solution of variational inequalities. We eliminate the need for Lagrange multipliers to enforce the contact conditions by using its definition in a discrete setting, following an idea of Chouly and Hild [4] used for elastic contact.

1.1 The model problem

We consider the obstacle problem of finding the displacement of a membrane constrained to stay above an obstacle given by (with at ):

 −Δu−f≥0 in Ω⊂R2u≥ψ in Ω(u−ψ)(f+Δu)=0 in Ωu=0 on ∂Ω (1)

where is a convex polygon. It is well known that this problem admits a unique solutions . This follows from the theory of Stampacchia applied to the corresponding variational inequality (see for instance [8]).

1.2 The finite element method

There exists a large body of literature treating finite element methods for unilateral problems in general and obstacle problems in particular, e.g., [11, 13, 6, 10, 8, 3, 16, 18, 2, 17]. Discretization of (1) is usually performed directly starting from the variational inequality or using a penalty method. The first approach however leads to some nontrivial choices in the construction of the discretization spaces in order to satisfy the nonpenetration condition and associated inf-sup conditions and until recently it has proved difficult to obtain optimal error estimates [9, 5]. The latter approach, o n the other hand leads to the usual consistency and conditioning problems of penalty methods.

An alternative is to use the augmented Lagrangian method. We introduce the Lagrange multiplier such that

 −Δu+λ=f in Ωu=0 on ∂Ω (2)

under the Kuhn-Tucker side conditions

 ψ−u≤0 in Ωλ≤0 in Ω(ψ−u)λ=0 on Ω. (3)

Using the standard trick of rewriting the Kuhn-Tucker conditions as

 λ=−1γ[ψ−u−γλ]+ (4)

where and , cf., e.g., Chouly and Hild [4], we can formulate the augmented Lagrangian problem of finding that are stationary points to the functional

 F(u,λ):= 12∫Ω|∇u|2dΩ+∫Ω12γ[ψ−u−γλ]2+dΩ −12∫Ωγλ2dΩ−∫ΩfudΩ, (5)

cf. Alart and Curnier [1], leading to seeking such that

 (6)

and

 ∫Ω1γ[ψ−u−γλ]+μdΩ+∫ΩλμdΩ=0∀μ∈L2(Ω). (7)

For our discrete method, we assume that is a family of conforming shape regular meshes on , consisting of triangles and define as the space of –conforming piecewise polynomial functions on , satisfying the homogeneous boundary condition of .

 Vh:={vh∈H10(Ω):vh|T∈Pk(T),∀T∈T}, for k≥2.

We then formally replace element–wise by to obtain a discrete minimization problem: seek such that

 uh=argminv∈VhFh(v) (8)

where

 Fh(v):= 12∫Ω|∇v|2dΩ+∑T∈T∫T12γ[ψ−v−γ(Δv+f)]2+dΩ −12∑T∈T∫Tγ(Δv+f)2dΩ−∫ΩfvdΩ. (9)

The Euler–Lagrange equations corresponding to (9) take the form: Find such that

 a(uh,vh)+b(uh,ψ,f;vh)=(f,vh)Ω∀vh∈Vh (10)

where denotes the standard -inner product, and

 b(uh,ψ,f;vh):= ⟨γ−1[ψ−uh−γ(Δuh+f)]+,−vh−γΔvh⟩h −⟨γ(Δuh+f),Δvh⟩h (11)

where

 ⟨xh,yh⟩h:=∑T∈T∫Txhyh\leavevmode\nobreak dx

and, for use below,

 ∥xh∥h:=⟨xh,xh⟩1/2h.

To simplify the notation below we introduce and

 b(uh,ψ,f;vh):=⟨γ−1[ψ−γf−Pγ(uh)]+,−Pγ(vh)⟩h−⟨γ(Δuh+f),Δvh⟩h.

We will also omit and from the argument of below, and use the notation so that

 b(uh;vh):=⟨γ−1[Ψ−Pγ(uh)]+,−Pγ(vh)⟩h−⟨γ(Δuh+f),Δvh⟩h.

Note that the form can be interpreted as a nonlinear consistent least squares penalty term for the imposition of the contact condition. A similar method was proposed in Stenberg et al. [7] in the framework of variational inequalities.

We will below alternatively use the compact notation

 Ah(uh,vh):=a(uh,vh)+b(uh;vh)

and the associated formulation, find such that

 Ah(uh;vh)=(f,vh)Ω, for all vh∈Vh. (12)

1.3 Summary of main results and outline

In Section 2 we recall some technical results, in Section 3 we derive an existence result for the discrete solution using Brouwer’s fixed point theorem and we prove uniqueness of the solution using monotinicity of the the nonlinearity, in Section 4 we prove an a priori error estimate and using a saturation assumption we also prove an a posteriori error estimate, finally in Section 5 we present numerical results confirming our theoretical results and illustrating the performance of an adaptive algorithm based on our a posteriori error estimate.

2 Technical results

Below we will use the notation for where is a constant independent of , but not of the local mesh geometry.

First we recall the following inverse inequality,

 ∥∇vh∥T≤Cih−1T∥vh∥T,T∈T (13)

see Thomée [15].

We will use the Scott-Zhang interpolant preserving boundary conditions, denoted . This operator is -stable, and the following interpolation error estimate is known to hold [14],

 ∥u−ihu∥Ω+h∥u−ihu∥H1(Ω)+h2∥Δ(u−ihu)∥h≲hk+1|u|Hk+1(Ω). (14)

The essential properties of the nonlinearity are collected in the following lemmas.

Lemma 1

Let then there holds

 ([a]+−[b]+)2≤([a]+−[b]+)(a−b),
 |[a]+−[b]+|≤|a−b|.
{@proof}

[Proof] Developing the left hand side of the expression we have

 [a]2++[b]2+−2[a]+[b]+≤[a]+a+[b]+b−a[b]+−[a]+b=([a]+−[b]+)(a−b).

For the proof of the second claim, this is trivially true in case both and are positive or negative. If is negative and positive then

 |[a]+−[b]+|=|b|≤|b−a|

and similarly if is negative and positive

 |[a]+−[b]+|=|a|≤|b−a|.
Lemma 2

(Continuity of ) For all , the form (11) satisfies

 |b(u1;v)−b(u2;v)|≲γ−1(∥(u1−u2)∥Ω+γh−1∥∇(u1−u2)∥Ω)(∥v∥Ω+γh−1∥∇vh∥Ω). (15)
{@proof}

[Proof]

 b(u1;vh)−b(u2;vh)= γ−1⟨[Ψ−Pγ(u1)]+−[Ψ−Pγ(u2)]+,−Pγ(vh)⟩h −⟨γΔ(u1−u2),Δvh⟩h.

Using the second inequality of Lemma 1 we see that the nonlinearity satisfies

 γ−1|⟨[Ψ−Pγ(u1)]+−[Ψ−Pγ(u2)]+,−Pγ(vh)⟩h|≤γ−1∥Pγ(u1)+Ψ−Pγ(u2)−Ψ∥C,f∥Pγ(vh)∥h. (16)

By the inverse inequality (13) we have

 ⟨γΔ(u1−u2),Δvh⟩h≤C2iγh−2∥∇(u1−u2)∥Ω∥∇vh∥Ω (17)

and

 ∥Pγ(vh)∥h≤∥vh∥Ω+Ciγh−1∥∇vh∥Ω. (18)

Collecting (16),(17) and (18) we have

 |b(u1;vh)−b(u2;vh)|≤ (∥u1−u2∥Ω+Ciγh−1∥∇(u1−u2)∥Ω)(∥vh∥Ω+Ciγh−1∥∇vh∥Ω) +C2iγh−2∥∇(u1−u2)∥Ω∥∇vh∥Ω

and the claim follows.

3 Existence of unique discrete solution

In the previous works on Nitsche’s method existence and uniqueness has been proven by using the monotonicity and hemi-continuity of the operator. Here we propose a different approach where we use the Brouwer’s fixed point theorem to establish existence and the monotonicity of the nonlinearity for uniqueness. We start by showing some positivity results and a priori bounds. Since we are interested in existence and uniqueness for a fixed mesh parameter , we do not require that the bounds in this section are uniform in .

Lemma 3

Let and assume that

 γ

then there holds

 α2∥u1−u2∥2H1(Ω)+γ−1∥[Ψ−Pγ(u1)]+−[Ψ−Pγ(u2)]+∥2h≤Ah(u1;u1−u2)−Ah(u2;u1−u2)

and

 α4∥u1∥2H1(Ω)≤Ah(u1;u1)+Cα−1γ−2∥[Ψ]+∥2Ω.
{@proof}

[Proof] First observe consider the form ,

 b(u1;vh)−b(u2;vh)=γ−1⟨[Ψ−Pγ(u1)]+−[Ψ−Pγ(u2)]+,−vh−γΔvh+Ψ−Ψ⟩h−⟨γΔ(u1−u2),Δvh⟩h.

Using the monotonicity of Lemma 1 we may write

 b(u1;u1−u2)−b(u2;u1−u2)≥ γ−1∥[Ψ−Pγ(u1)]+−[Ψ−Pγ(u2)]+∥2h −γ∥Δ(u1−u2)∥2h.

Observe that using an inverse inequality (13) we have

 γ∥Δ(u1−u2)∥2h≤Ciγh−2∥∇(u1−u2)∥2Ω

. We may then write

 (1−C2ih−2γ)∥∇(u1−u2)∥2Ω+γ−1∥[Ψ−Pγ(u1)]+−[Ψ−Pγ(u2)]+∥2h≤Ah(u1;u1−u2)−Ah(u2;u1−u2)

It follows that choosing and applying the Poincaré inequality

 α12∥u∥H1(Ω)≤∥∇u∥Ω,∀u∈H10(Ω) (20)

there holds

 α2∥u1−u2∥2H1(Ω)+γ−1∥[Ψ−Pγ(u1)]+−[Ψ−Pγ(u2)]+∥2h≤Ah(u1;u1−u2)−Ah(u2;u1−u2).

The second inequality follows by taking above and noting that then

 α2∥u1∥2H1(Ω)+γ−1∥[Ψ−Pγ(u1)]+−[Ψ]+∥2h ≤Ah(u1;u1)−γ−1⟨[Ψ]+,−Pγ(u1)⟩h ≤Ah(u1;u1)+(1+Ciγh−1)α−12γ−1∥[Ψ]+∥Ωα12∥u1∥H1(Ω)

where we used (18) in the last step. Considering the condition on and using an arithmetic-geometric inequality we may conclude.

Proposition 4

Assume that saisfy (19). Then formulation (12) using the contact operator (11), admits a unique solution.

{@proof}

[Proof] The uniqueness is an immediate consequence of Lemma 3. If and both are solution to (12), then

 Ah(u1;u1−u2)−Ah(u2;u1−u2)=(f,u1−u2)Ω−(f,u1−u2)Ω=0

and we conclude that and hence . Let denote the number of degrees of freedom of .

Consider the mapping defined by

 (G(U),V)RNV:=Ah(uh,vh)−(f,vh)Ω,

where , , where and denotes the vectors of unknown associated to the basis functions of .

By the second claim of Lemma 3, there holds

 α4∥uh∥2H1(Ω)−Cα−1γ−2∥[Ψ]+∥2Ω−(f,uh)Ω≤Ah(uh,uh)−(f,uh)Ω=(G(U),U)RNV.

Since

 α4∥uh∥2H1(Ω)−(f,uh)Ω≥α8∥uh∥2H1(Ω)−C1α∥f∥2Ω

we have that for any fixed the following positivity holds for sufficiently large

 0<α8∥uh∥2H1(Ω)−Cα(γ−2∥[Ψ]+∥2Ω+∥f∥2Ω)≤(G(U),U)RNV.

Assume that this positivity holds whenever . Denote by the (closed) ball in with radius and assume that there is no such that . Define the function

 ϕ(U)=−qG(U)/|G(U)|.

Then , is continuous by Lemma 2 and the assumption that for all . Hence there exists a fixed point such that

 X=ϕ(X).

It follows that

 |X|2=−q(G(X),X)/|G(X)|,

but since , by assumption , which leads to a contradiction, since .

4 Error estimates

Theorem 5

(A priori error estimate) Assume that with is the solution of (1) and that is the solution to (10) with (11) and , where , . then there holds for all

 α∥u−uh∥2H1(Ω)+γ−1∥[Ψ−Pγ(uh)]+−[Ψ−Pγ(u)]+∥2h≲1α∥u−vh∥2H1(Ω)+∥γ−12(u−vh)∥Ω+∥γ12Δ(u−vh)∥2h. (21)

If in addition then there holds

 α∥u−uh∥H1(Ω)+γ−1/2∥[Ψ−Pγ(uh)]+−[Ψ−Pγ(u)]+∥h≲hk|u|Hk+1(Ω). (22)
{@proof}

[Proof] Using the definition of we may write

 ∥∇(u−uh)∥2Ω ≤a(u−uh,u−uh) =a(u−uh,u−vh)+a(u−uh,vh−uh) ≤α4∥u−uh∥2H1(Ω)+1α∥u−vh∥2H1(Ω)+a(u−uh,vh−uh).

Observe that

 a(u,vh−uh)= ⟨−Δu−f+f,vh−uh⟩V′,V = ⟨γ−1[Ψ−Pγ(u)]+,(vh−uh)⟩V′,V+(f,vh−uh)Ω. (23)

If we may also write

 ⟨Δu+f,γΔ(vh−uh)⟩h+⟨γ−1[Ψ−Pγ(u)]+,γΔ(vh−uh)⟩h=0.

It follows that

 a(u,vh−uh)= (f,vh−uh)Ω−⟨γ−1[Ψ−Pγ(u)]+,−Pγ(vh−uh)⟩h +⟨Δu+f,γΔ(vh−uh)⟩h = (f,vh−uh)Ω−b(u;vh−uh). (24)

As a consequence we have the following property reminiscent of Galerkin orthogonality,

 a(u−uh,vh−uh) =b(uh;vh−uh)−b(u;vh−uh) =⟨γ−1[Ψ−Pγ(uh)]+−γ−1[Ψ−Pγ(u)]+,−Pγ(vh−uh)⟩h (25) −γ⟨Δ(uh−u),Δ(vh−uh)⟩h

First observe that

 γ⟨Δ(uh−u),Δ(vh−uh)⟩h ≤∥γ12(Δuh−Δvh)∥2h+∥γ12(Δvh−Δu)∥h∥γ12(Δvh−Δuh)∥h ≤12C2ih−2γ∥∇(uh−vh)∥2Ω+∥γ12(Δvh−Δu)∥2h ≤C2ih−2γ∥∇(uh−u)∥2Ω+C2ih−2γ∥∇(vh−u)∥2Ω+∥γ12(Δvh−Δu)∥2h.

Considering the first term in the right hand side of equation (25) we may write

 ⟨γ−1[Ψ−Pγ(uh)]+−γ−1[Ψ−Pγ(u)]+,−Pγ(vh−uh)⟩h =⟨γ−1[Ψ−Pγ(uh)]+−γ−1[Ψ−Pγ(u)]+,−Pγ(vh−u)⟩hI +⟨γ−1[Ψ−Pγ(uh)]+−γ−1[Ψ−Pγ(u)]+,−Pγ(u−uh)⟩hII =I+II.

The term may be bounded using the Cauchy-Schwarz inequality followed by the arithmetic geometric inequality

 I≤ϵγ−1∥[Ψ−Pγ(uh)]+−[Ψ−Pγ(u)]+∥2h+14ϵ∥γ−12Pγ(vh−u)∥2h.

For the term we use the monotonicity property , with and so that

 ([a]+−[b]+)(b−a)=([Ψ−Pγ(uh)]+−γ−1[Ψ−Pγ(u)]+)(Ψ−Pγ(u)−Ψ+Pγ(uh))

to deduce that

 II≤−γ−1∥[Ψ−Pγ(uh)]+−[Ψ−Pγ(u)]+∥2h.

Collecting the above bounds and using the Poincaré inequality (20) we find,

 α(34−C2ih−2γ)∥u−uh∥2H1(Ω)+(1−ϵ)γ−1∥[Ψ−Pγ(uh)]+−[Ψ−Pγ(u)]+∥2h+(1−ϵ)∥γ12Δ(u−uh)∥2h≤1α∥u−vh∥2H1(Ω)+14ϵ∥γ−12Pγ(u−vh)∥2h+∥γ12Δ(u−vh)∥2h (26)

Fixing , and fixing sufficiently small so that then there holds

 α∥u−uh∥2H1(Ω)γ−1∥[Ψ−Pγ(uh)]+−[Ψ−Pγ(u)]+∥2h≲1α∥u−vh∥2H1(Ω)+∥γ−12(u−vh)∥2h+∥γ12Δ(u−vh)∥2h. (27)

This concludes the proof of (21). The error estimate (22) then follows by choosing to be the interpolant, , applying the approximation error estimate (14) on the form

 ∥u−ihu∥H1(Ω)+∥γ−12(u−ihu)∥h+∥γ12Δ(u−ihu)∥h ≲(hk+γ−1/2hk+1+γ1/2hk−1)|u|Hk+1(Ω)

and using the bound on . Assumption: (Saturation) We assume that there exists a constant such that

 ∥Δ(u−uh)∥h≤Csh−1∥∇(u−uh)∥Ω. (28)
Theorem 6

(A posteriori error estimate) Assume that with is the solution of (1) and the solution of (10) satisfying (28) and with the parameter satisfying , then

 α∥u−uh∥H1(Ω)+γ−1∥[Ψ−Pγ(uh)]+−[Ψ−Pγ(u)]+∥h≲E(h,γ,uh,f), (29)

where

 E(h,γ,uh,f):=h∥f+Δuh+γ−1[Ψ−Pγ(uh)]+∥h+∥h12\llbracket∂nuh\rrbracket∥F.
{@proof}

[Proof] Let then, under the assumption (28) and using (3) we may write

 (1−C2sh−2γ)∥∇e∥2Ω+γ−1∥[Ψ−Pγ(u)]+−[Ψ−Pγ(uh)]+∥2h ≤(∇(u−uh),∇e)Ω (30) +γ−1⟨[Ψ−Pγ(u)]+−[Ψ−Pγ(uh)]+,−Pγ(e)⟩h −γ⟨Δ(u−uh),Δe⟩h

Now, using similar arguments as in Theorem 5 we deduce

 (∇u,∇e)Ω−γ−1⟨[Ψ−Pγ(u)]+,Pγ(e)⟩h−γ⟨Δu,Δe⟩h=⟨f,Pγ(e)⟩h.

Therefore the bound (30) may be written, under the assumption ,

 12∥∇e∥2Ω+γ−1∥[Ψ−Pγ(u)]+−[Ψ−Pγ(uh)]+∥2h ≤⟨f,Pγ(e)⟩h−(∇uh,∇e)Ω (31) +γ−1⟨[Ψ−Pγ(uh)]+,Pγ