Augmented Lagrangian finite element methods for contact problems

Augmented Lagrangian finite element methods for contact problems

Erik Burman Department of Mathematics, University College London, Gower Street, London, UK–WC1E 6BT, United Kingdom; (e.burman@ucl.ac.uk)    Peter Hansbo Department of Mechanical Engineering, Jönköping University, SE-55111 Jönköping, Sweden; (peter.hansbo@ju.se)    Mats G. Larson Department of Mathematics and Mathematical Statistics, Umeå University, SE-901 87 Umeå, Sweden; (mats.larson@math.umu.se)
Abstract

We propose two different Lagrange multiplier methods for contact problems derived from the augmented Lagrangian variational formulation. Both the obstacle problem, where a constraint on the solution is imposed in the bulk domain and the Signorini problem, where a lateral contact condition is imposed are considered. We consider both continuous and discontinuous approximation spaces for the Lagrange multiplier. In the latter case the method is unstable and a penalty on the jump of the multiplier must be applied for stability. We prove the existence and uniqueness of discrete solutions, best approximation estimates and convergence estimates that are optimal compared to the regularity of the solution.

1 Introduction

We consider the Signorini problem, find and such that

 −Δu=f in Ωu=0 on ΓDu≤0,λ≤0,uλ=0 on ΓC, \hb@xt@.01(1.1)

or the obstacle problem

 −Δu−λ=f in Ωu=0 on ∂Ωu≤0,λ≤0,uλ=0 in Ω. \hb@xt@.01(1.2)

Here , is a bounded polyhedral (polygonal) domain and . It is well known that these problems admit unique solutions . This follows from the theory of Stampacchia applied to the corresponding variational inequality (see for instance [19]).

From a mechanical point of view, these equations model the deflection of a membrane in isotropic tension under the load , assuming small deformations. The membrane is either in contact with an obstacle on part of the boundary, (1), or in the interior of the membrane, (1), preventing positive displacements . In both cases the Lagrange multiplier has the interpretation of a distributed reaction force enforcing the contact condition .

2 Finite element discretization

Our aim in this paper is to design a consistent penalty method for contact problems that can easily be included in a standard Lagrange-multiplier method, without having to resort to the solution of variational inequalities. We consider two different choices for the multiplier spaces, either a stable choice or an unstable choice where a stabilization term is needed to ensure the stability of the formulation. In the latter case we add a penalty on the jump of the multiplier over element faces in the spirit of [11, 10].

There exists a large body of litterature treating finite element methods for contact problems [8, 22, 17, 5, 4, 6, 28, 27, 14]. Discretization of (1) is usually performed on the variational inequality or using a penalty method. The first case 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 [21, 16]. The latter case, on the other hand leads to the usual consistency and conditioning problems of penalty methods. Another approach proposed by Hild and Renard [20] is to use a stabilized Lagrange-multiplier in the spirit of Barbosa and Hughes [3]. As a further development one may use the reformulation of the contact condition

 λ=−γ−1[u−γλ]+ \hb@xt@.01(2.1)

where , introduced by Alart and Curnier [1] in an augmented Lagrangian framework. Using the close relationship between the Barbosa–Hughes method and Nitsche’s method [23] discussed by Stenberg [25], this method was then further developed in the elegant Nitsche-type formulation for the Signorini problem introduced by Chouly, Hild and Renard [13, 15]. In these works optimal error estimates for the above model problem were obtained for the first time.

Using the notation for the inner product over we have in the case of the Signorini problem (1) that corresponds to , the boundary part where the contact conditions hold and

 ⟨u,v⟩C:=∫ΓCuv\leavevmode\nobreak ds,

while for the obstacle problem (1) and

 ⟨u,v⟩C:=∫Ωuv\leavevmode\nobreak dx.

Finally, we define . With this notation, the augmented Lagrangian multiplier seeks stationary points to the functional

 F(u,λ):=12a(u,u)+12γ∥[u−γλ]+∥2C−γ2∥λ∥2C, \hb@xt@.01(2.2)

cf. Alart and Curnier [1]. Observe that formally the stationary points are given by such that

 a(u,v)+⟨γ−1[u−γλ]+,v⟩C=(f,v)Ω⟨λ+γ−1[u−γλ]+,μ⟩C=0 \hb@xt@.01(2.3)

for all , or by substituting the second equation in the first

 a(u,v)−⟨λ,v⟩C=(f,v)Ω⟨γλ+[u−γλ]+,μ⟩C=0. \hb@xt@.01(2.4)

Observing now that the contact condition equally well can be written on the primal variable as we get by adding and subtracting in the second equation of (2.4)

 a(u,v)−⟨λ,v⟩C=(f,v)Ω⟨u+[γλ−u]+,μ⟩C=0. \hb@xt@.01(2.5)

In this paper we consider two different method, resulting from this approach. The first formulation is the straightforward discretization of (2.3) resulting in a method that gives the stationary points of the functional (2.2) over the discrete spaces. The second formulation is a discretization of (2.5) that is chosen for its closeness to the standard Lagrange multiplier method for the imposition of Dirichlet boundary conditions.

We consider discretization either with a choice of approximation spaces that results in a stable approximation, or a choice that is stable only with an added stabilizing term. Here we consider stabilization based on the interior penalty stabilized Lagrange multiplier method introduced by Burman and Hansbo [11] for solving elliptic interface problems. The appeal of this latter approach is that we may use the lowest order approximation spaces where the displacement is piecewise linear and the multiplier constant per element (or element side). When considering the Signorini problem (1) these spaces match the regularity of the physical problem perfectly and therefore in some sense is the most economical choice. Contact problems also present non trivial quadrature problems so that in practice it can be very difficult to integrate the terms of the formulation to a sufficient accuracy to get optimal accuracy when higher order interpolations are used. Herein we will assume that integration can be performed exactly on the interface between the contact and non-contact subdomain.

For an alternative stabilization method of Barbosa–Hughes type in the augmented Lagrangian setting, see Hansbo, Rashid, and Salomonsson [18].

We assume that is a family of quaisuniform meshes of , such that the mesh is fitted to the zone . That is is a subset of boundary element faces of simplices such that , , with for the Signorini problem. For the obstacle problem is defined by and hence and . Below we will denote the elements of by in both cases. We define to be the space of -conforming functions on , satisfying the homogeneous boundary condition of .

 Vkh:={vh∈H1(Ω):v|ΓD=0;v|K∈Pk(K),∀K∈T},

where denotes the set of polynnomials of order less than or equal to on the simplex . Whenever the superscript is dropped we refer to the generic space of order . For the multipliers we introduce the space defined as the space piecewise polynomials of order less than or equal to defined on .

 Λlh:={μh∈L2(C):μh|K∈Pl(K),∀K∈TC}.

Whenever the superscript is dropped. We will detail the case of discontinuous multipliers, but all arguments below are valid also in case the Lagrange multiplier is approximated in the space of continuous functions, , , in this case no stabilization is necessary. The differences in the analysis will be outlined.

Both formulations that we consider herein take the form: Find such that

 a(uh,vh)+b[(uh,λh);(vh,μh)]=(f,vh)Ω∀(vh,μh)∈Vh×Λh \hb@xt@.01(2.6)

where denotes the standard -inner product, and the methods are distinguished by the definition of the form that acts only in the zone where contact may occur. The stabilization will be included in the form . As already pointed out this term is necessary if the choice , does not satisfy the inf-sup condition. In our framework, this is the case where the multiplier is discontinuous over element faces. In this paper we will focus on a stabilization using a penalty on the jumps over element faces of the multiplier variable in the spirit of [11, 10],

 s(λh,μh):=∑F∈FCδγ∫Fh\llbracketλh\rrbracket\llbracketμh\rrbracket\leavevmode\nobreak ds, \hb@xt@.01(2.7)

where is a parameter, denotes the jump of the quantity over the face and denotes the set of interior element faces of the elements in . The semi-norm associated with the stabilization operator will be defined as .

We will also below use the compact notation

 Ah[(uh,λh),(vh,μh)]:=a(uh,vh)+b[(uh,λh);(vh,μh)]

and the associated formulation, find such that

 Ah[(uh,λh),(vh,μh)]=(f,vh)Ω, for all (vh,μh)∈Vh×Λh. \hb@xt@.01(2.8)

We will now specify two different choices of leading to two different Lagrange-multiplier methods.

FORMULATION 1

: In the first formulation we use the original formula for the contact condition proposed by Alart and Curnier,

 b[(uh,λh);(vh,μh)]:= ⟨γ−1[uh−γλh]+,vh⟩C +⟨γ−1[uh−γλh]+,γμh⟩C +⟨γλh,μh⟩C+s(λh,μh) \hb@xt@.01(2.9)

or, writing the nonlinearity as the derivative of a quadratic form, and using the notation

 b[(uh,λh);(vh,μh)]:= ⟨γ−1[Pγ+(uh,λh)]+,Pγ+(vh,μh)⟩C −⟨γμh,λh⟩C−s(λh,μh), \hb@xt@.01(2.10)

with a parameter to determine. In this case the finite element formulation corresponds to the approximate solutions of (2.3) in the finite element space.

FORMULATION 2

: In the second formulation we use a reformulation of the contact condition on the displacement variable, to obtain the semi-linear form

 b[(uh,λh);(vh,μh)]:= −⟨λh,vh⟩C+⟨μh,uh⟩C +⟨μh,[Pγ−(uh,λh)]+⟩C+s(λh,μh), \hb@xt@.01(2.11)

with a parameter to determine. In this case the finite element formulation corresponds to the approximate solutions of (2.5) in the finite element space.

2.1 Alternative formulations

In both formulation 1 and 2 above it is possible to derive an alternative formulation of the same method using the relation

 [Pγ−(uh,λh)]+=[Pγ+(uh,λh)]+−Pγ+(uh,λh).

Considering the form (2.10) and adding and subtracting in the nonlinear term we have the alternative form (omitting the stabilization term)

 b[(uh,λh);(vh,μh)]= ⟨γ−1[Pγ+(uh,λh)]+,Pγ+(vh,μh)⟩C−⟨γμh,λh⟩C = ⟨γ−1([Pγ+(uh,λh)]+−Pγ+(uh,λh)),Pγ+(vh,μh)⟩C +⟨γ−1Pγ+(uh,λh),Pγ+(vh,μh)⟩C−⟨γμh,λh⟩C = −⟨λh,vh⟩C+⟨μh,uh⟩C+γ−1⟨uh,vh⟩C +⟨γ−1([Pγ−(uh,λh)]+,Pγ+(vh,μh)⟩C. \hb@xt@.01(2.12)

Similarly for formulation 2 we obtain in (2.11) omitting for simplicity the stabilization term

 b[(uh,λh);(vh,μh)]= −⟨λh,vh⟩C+⟨μh,uh⟩C +⟨μh,[Pγ−(uh,λh)]++Pγ+(uh,λh)−Pγ+(uh,λh)⟩C = −⟨λh,vh⟩C+γ⟨μh,λh⟩C +⟨μh,[Pγ+(uh,λh)]+⟩C. \hb@xt@.01(2.13)

We see that this semi-linear form corresponds to a discretization of (2.4).

The methods defined by (2.11) and (2.13) or (2.10) and (2.12) respectively are equivalent, but if during the solution process the linear and nonlinear parts are separated in the nonlinear solver, one can expect the different formulations to have different behavior and give rise to different sequences of approximations in the iterative procedure.

3 Technical results

Here we will collects some useful elementary results. First recall the following inverse inequalities and trace inequalities (for a proof see, e.g., [2])

 ∥∇uh∥K≤Cih−1∥uh∥K,∀uh∈Vh \hb@xt@.01(3.1)
 ∥u∥∂K≤CT(h−12∥u∥K+h12∥∇u∥K),∀u∈H1(K) \hb@xt@.01(3.2)
 ∥uh∥∂K≤CTh−12∥uh∥K,∀uh∈Vh \hb@xt@.01(3.3)

Similar inequalities hold for functions in and we will use them without making any distinction between the two cases. We let denote the standard projection onto and we observe that there holds, by standard approximation properties of the projection onto constants (and a trace inequality in the case of lateral contact),

 ∥(1−π0)vh∥C≤c0hs∥∇vh∥Ω

with for the Obstacle problem where and for the Signorini problem where . Similarly we define and note that the corresponding inequality holds for

 ∥(1−π1)vh∥C≤c1hs∥∇vh∥Ω.

We also observe for future reference that in both cases.

For the analysis below it is useful to introduce an indicator function for the contact domain defined on the space . Let denote a finite element function such that with for nodes in , that is nodes outside the contact zone. For all other nodes with , , . The following bound is well known, see for instance [12]

 ∃cξ∈R+ such that cξ∥μh∥C≤∥ξ12hμh∥C,∀μh∈Λlh,l≥0. \hb@xt@.01(3.4)

Stability of the method will rely on the satisfaction of the following assumption:

Assumption 3.1

There exists such that for all there holds

 ∥(1−ξh)μh∥C≤cD∥μh∥C.

The assumption holds whenever there exists a quadrature rule on the simplex with positive weights and only interior quadrature points. This is easily shown by observing that

 ∥(1−ξh)μh∥2C =∑K∈TC∑i∈QK(1−ξh(xi))2μh(xi)2ωi ≤maxK∈TC(maxi∈QK(1−ξh(xi))2∑K∈T∑i∈QK(μh(xi))2ωi =c2D∥μh∥2C

where is a set of integers indexing the quadrature points in and

 cD≡maxK∈TC(maxi∈QK(1−ξh(xi))2.

Since is zero only on the boundary of and no points are on the boundary we conclude that .

This is a very mild condition, on triangles it has been showed to hold at least up to integration degree , see [26, 29]. It follows that for the Signorini problem in three dimensions and the obstacle problem in two space dimensions the analysis holds at least up to . For the lowest order case where the multipliers are constant per element it is straightforward to show that if and if .

Lemma 3.1

Let ; then there holds

 ([a]+−[b]+)2≤([a]+−[b]+)(a−b),
 |[a]+−[b]+|≤|a−b|.

Proof. Expanding 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 3.2

(Continuity of ) The forms (2.11) and (2.10) satisfy

 |b[(u1,λ1);(v,μ)]−b[(u2,λ2);(v,μ)]| ≤(γ−12∥(u1−u2)∥H1(Ω)+γ12∥λ1−λ2∥C)(γ−12∥v∥C+γ12∥μ∥C) +|λ1−λ2|s|μ|s.

Proof. Immediate by the definitions of , the second inequality of Lemma 3.1, the Cauchy-Schwarz inequality and the assumptions on .

Next we define the local averaging interpolation operator such that for every Lagrangian node

 Icfλh(xi)=κ−1i∑K:xi∈Kλh(xi),

where denotes the cardinality of the set . Observe that since , for any there are functions in such that . We recall the following interpolation result between discrete spaces:

Proposition 3.3

For all there holds

 ∥ξhμh−Icf(ξhμh)∥C≤cs∥h12\llbracketμh\rrbracket∥FC,∥Icfμh∥C≤ccf∥μh∥C

and

 |μh|2s≤Cδ∥μh∥2C.

Proof. For a proof of the first inequality we refer to [9, Lemma 5.3]. The second inequality is immediate by applying the trace inequality (3.3) to each term in the definition (2.7) of .

Lemma 3.4

Let , then there exists such that and , with when is a subset of and when is a subset of .

Proof. Define so that for all nodes in and for all other nodes in the mesh. First consider the case when is a subset of the bulk domain . Then, using an inverse inequality,

 ∥∇Rh∥Ω≤Cih−1∥rh∥Ω=Cih−1∥rh∥Ω=Cih−1∥rh∥C.

In the case is a subset of the boundary of we observe that

Using that is defined by the nodes in , combined with the shape regularity of the mesh, we may use the following inverse trace inequality [9, Lemma 3.1] on every ,

 ∥Rh∥K≤Ch12∥Rh∥∂K∩C.

It follows, since , that

 ∥∇Rh∥Ω≤Ch−1/2∥Rh∥C≤Ch−1/2∥rh∥C.

4 Existence of unique discrete solution

In the previous works on Nitsche’s method for contact problems [13, 15] 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 Brouwer’s fixed point theorem to establish existence and the monotonicity of the nonlinearity for uniqueness. To this end we introduce the finite dimensional nonlinear system corresponding to the formulation (2.6).

Let , where and denote the number of degrees of freedom of and respectively. Then define , where , , where and denote the vectors of unknowns associated to the basis functions of and respectively.

Consider the mapping defined by

 (G(U),V)RM:=Ah[(uh,λh),(vh,μh)]−(f,vh)Ω.

Existence and uniqueness of a solution to (2.6) is equivalent to showing that there exists a unique such that .

We start by showing some positivity results and a priori bounds

Lemma 4.1

There exists and an associated constant so that with the form defined by (2.10), and with there holds, for all

 ∥∇uh∥2Ω+γ∥λh+γ−1[Pγ+(uh,λh)]+∥2C+cα∥γ12λh∥2C≲Ah[(uh,λh),(uh−αRh,λh)], \hb@xt@.01(4.1)

where is defined in Lemma 3.4, such that .

There exists and an associated constant so that with the form defined by (2.11), and with , sufficiently large, and there holds, for all

 ∥∇uh∥2Ω+γ−1∥uh+[Pγ−(uh,λh)]+∥2C+cα∥γ12λh∥2C≲Ah[(uh,λh),(uh+αRh,λh+γ−1π0uh)], \hb@xt@.01(4.2)

with as before. In case (4.2) holds under the additional that .

Under the same conditions on the parameters as above, for both formulations there also holds, for solution of (2.8),

 ∥∇uh∥Ω+∥γ12λh∥C≲∥f∥Ω. \hb@xt@.01(4.3)

The hidden constants are independent of .

Remark 4.1

For and continuous multiplier space the parameter and the term can be dropped above.

Proof. First consider the claims for formulation 1. By testing in (2.8) with and and observing that

 ⟨γ−1[Pγ+(uh,λh)]+,γλh⟩C+∥γ12λh∥2C=∥γ12λh∥2C+⟨γ−1[Pγ+(uh,λh)]+,−γλh⟩C+2⟨γ−12[Pγ+(uh,λh)]+,γ12λh⟩C

we obtain the relation

 ∥∇uh∥2Ω+γ∥γ−12[Pγ+(uh,λh)]++λh∥2C+|λh|2s=Ah[(uh,λh),(uh,λh)]

and hence by using the Cauchy-Schwarz inequality and a Poincaré inequality in the right hand side

 12∥∇uh∥2Ω+γ∥γ−12[Pγ+(uh,λh)]++λh∥2C+|λh|2s≲∥f∥2Ω \hb@xt@.01(4.4)

Using now the first equation we have testing with , with and ,

 Ah[(uh,λh),(−Rh,0)] =a(uh,−Rh)+⟨γ−1[Pγ+(uh,λh)]+,−γIcf(ξhλh)⟩C =a(uh,−Rh)+⟨γ−1[Pγ+(uh,λh)]++λh,−γIcf(ξhλh)⟩C −⟨λh,γ(ξhλh−Icf(ξhλh))⟩C+⟨λh,γξhλh⟩C. \hb@xt@.01(4.5)

For the last term in the right hand side we have by the inequality (3.4), . The second to last term of the right hand side, which is zero for continuous multiplier spaces, can be bounded using Proposition 3.3

 (γλh,ξhλh−Icf(ξhλh))C≤c2ξ14∥γ12λh∥2C+c2sc−2ξδ−1|λh|2s. \hb@xt@.01(4.6)

The second term is bounded using a Cauchy-Schwarz inequality and the stability of ,

 ⟨γ−1[Pγ+(uλ,λh)]++λh,γIcf(ξhλh)⟩C≤12(ccfcξ)−2γ∥γ−1[Pγ+(uλ,λh)]++λh∥2C+14c2ξ∥γ12λh∥2C \hb@xt@.01(4.7)

for the first term we use the Caucy-Schwarz inequality followed by the stability of and of to obtain

 a(uh,Rh)≤C2Rh−2sγc2cfc−2ξ∥∇uh∥2Ω+c2ξ14∥γ12λh∥2C. \hb@xt@.01(4.8)

Applying the inequalities (4.6)-(4.8) to (4.5) we have

 c2ξα∥γ12λh∥2C−C2Rh−2sγc2cfc−2ξα∥∇uh∥Ω −12(ccfcξ)−2γα∥γ−1[Pγ+(uλ,λh)]++λh∥2C −c2sc−2ξδ−1α|λh|2s ≤Ah[(uh,λh),(−αRh,0)] \hb@xt@.01(4.9)

We conclude by observing that and by combining the bounds (3.4), (4.4) and (4.9) with small enough. The a priori estimate follows noting that for solution of (2.6) there holds using the Poincaré inequality and the properties of ,

 Ah[(uh,λh),(uh−αRh,0)]=(f,uh−αRh)≤C∥f∥Ω(∥∇uh∥Ω+∥γ12λh∥C).

To prove (4.2) we start by testing in the left hand side of (2.8) with and , where if and for . Observing this time that, using the definition (2.11) and adding and subtracting at suitable places the following equality holds

 b[(uh,λh),(uh,λh+γ−1πiuh)] =γ−1⟨uh,uh⟩C−γ−1∥πiuh−uh∥2C +γ−1∥[Pγ−(uh,λh)]+∥2C +2⟨γ−1uh,[Pγ−(uh,λ