and Raviart–Thomas-type mixed finite elements in elastostatics

Numerische Mathematik

c Springer-Verlag 1996

Electronic Edition

Symmetric coupling of boundary elements

and Raviart–Thomas-type mixed finite elements in elastostatics

Ulrich Brink¹, Carsten Carstensen², Erwin Stein¹

1 Institut f¨ur Baumechanik und Numerische Mechanik, Universit¨at Hannover, D-30167 Hannover, Germany

2 Mathematisches Seminar II, Christian-Albrechts-Universit¨at zu Kiel, D-24098 Kiel, Germany

Received February 21, 1995 / Revised version received December 21, 1995

Dedicated to Professor G. C. Hsiao on occasion of his 60th birthday.

Summary. Both mixed finite element methods and boundary integral methods are important tools in computational mechanics according to a good stress ap- proximation. Recently, even low order mixed methods of Raviart–Thomas-type became available for problems in elasticity. Since either methods are robust for critical Poisson ratios, it appears natural to couple the two methods as proposed in this paper. The symmetric coupling changes the elliptic part of the bilinear form only. Hence the convergence analysis of mixed finite element methods is applicable to the coupled problem as well. Specifically, we couple boundary el- ements with a family of mixed elements analyzed by Stenberg. The locking-free implementation is performed via Lagrange multipliers, numerical examples are included.

Mathematics Subject Classification (1991): 65N30, 65R20

1. Introduction

In the classical finite element approach, the displacements are the unknowns while the stresses are computed afterwards in lower accuracy. In many applications, the stresses rather than the displacements are of primary interest. In the boundary element method (BEM), displacements and stresses in the interior of the domain are approximated with the same order.

Regarding finite elements, the approximation of the stresses can be improved by mixed methods. Here we consider Raviart–Thomas-type finite elements due to Stenberg [11], which are an improvement of Arnold–Brezzi–Douglas’s PEERS element (plane elasticity element with reduced symmetry) [1]. The unknowns are three independent fields, namely the stress tensor, which is not a priori assumed to be symmetric, the rotation, which acts as a Lagrange multiplier to enforce the symmetry of the stress tensor in a weak form, and the displacement vector.

The analysis is based on the theory of mixed methods but refined by using mesh-dependent norms (Pitk¨aranta and Stenberg [9], [10]), by stability proofs on patches of elements and by weakening of the symmetry condition.

Another advantage of these elements is the capability of modeling nearly incompressible elasticity; more precisely, the relative error is independent of Poisson’s ratioν: there is no locking.

In many applications, for example, if nonlinearities in a bounded domain ΩF are present as well as homogeneous, isotropic linear elastic material in an unbounded, e.g., exterior, domain ΩB, one might combine the finite element method (inΩF) with the boundary integral method (acting on∂ΩBbut treating the problem inΩ_B). The symmetric coupling with boundary elements was proposed and analyzed by Costabel in [6] where the displacements inside the FEM domain are sought in the Sobolev space H¹(ΩF); traces on the interface are inserted into the boundary integral equations, while the equilibrium of tractions across the interface is satisfied in a weak sense only.

In contrast, in the mixed FEM under consideration here, the stresses σ are required to satisfyσ∈L²(ΩF) and divσ∈L²(ΩF), while the displacements are sought in L²(ΩF) only. Hence, the discrete tractions across interelement sides are continuous and, consequently, our coupled scheme is designed to yield continuous tractions across the interface between FEM and BEM while continuity of the displacements across the interelement sides and the interface is satisfied in a weak sense.

The construction of finite element spaces with continuous interelement trac- tions may be cumbersome and hence is enforced by using Lagrange multipliers.

Then all unknowns, except the Lagrange multipliers, are discontinuous across the interelement sides and can be eliminated on each element before assembling the global linear system. We extended this scheme to the coupled system.

A similar coupling is possible with other mixed finite element methods than those considered in this paper. The weakening of the symmetry of the stress tensor is a particular way to construct stable finite element spaces, but is by no means necessary for the coupling.

The paper is organized as follows: A model problem (as depicted in Fig. 1) is described in its strong form in Sect. 2 and rewritten in a weak form. We prove existence and uniqueness of solutions and a norm estimate by using Brezzi’s theory of mixed problems. The discretization is described in Sect. 3 where we state convergence estimates under some hypotheses on the mixed finite element methods. Proofs are given in Sect. 4 while in Sect. 5 Stenberg’s locking-free fam- ily of methods is discussed. We describe the implementation by using Lagrange multipliers in Sect. 6 and give some numerical examples in Sect. 7.

2. The coupled problem

Let Ω be a bounded polygonal (resp. polyhedral) domain in R^d, d = 2 (resp.

d = 3) with boundary ∂Ω = Γu ∪Γt with disjointΓu and Γt, either having a positive surface measure.

Γu ΩF

FEM

Γt

-n

n ΩB

BEM

Γ

Fig. 1. Notation

Throughout this paper, we consider the linear elasticity problem divσ = −f inΩ

σ = Cε(u) inΩ

u = 0 onΓ_u

σn = 0 onΓt. (2.1)

The displacement field is denoted by u, the related (linear Green) strain tensor is ε(u) = ¹₂( grad u + ( grad u)^T). The elasticity tensor C describes the stress–

strain relationship; in the simplest case we have σ= Cε= λtrεI + 2µε where λwhereµ are the Lam´e coefficients and I denotes the identity matrix. (In the two-dimensional case, this is the constitutive equation of plane strain).

As depicted in Fig. 1, the domain Ω is partitioned into Ω= ΩF∪Γ ∪ΩB; the outward unit normal ofΩ_Fis denoted by n.

OnΩ_Fa mixed finite element method based on the Hellinger–Reissner prin- ciple is applied. This means the displacement u : ΩF → R^d and the stress σ : ΩF → R^d×d are independent unknowns. Furthermore, the stress tensor σ is not a priori assumed to be symmetric. Symmetry will be enforced in a weak form by a Lagrange multiplier technique. To obtain a variational formulation, we choose test functions τ : ΩF → R^d×d satisfying τn = 0 on Γt, and gain from

(2.1)_b Z

ΩF

τ:C⁻¹σdΩ− Z

ΩF

τ :ε(u) dΩ= 0 so that integration by parts gives

Z

ΩF

τ:C⁻¹σdΩ+ Z

ΩF

divτ·u dΩ+ Z

ΩF

τ:γdΩ= Z

Γ

τn·u ds

where the rotation γ := ¹₂( grad u−( grad u)^T) is a new variable for the skew- symmetric part of the displacement gradient,

γ∈W :={η∈L²(ΩF)^d×d :η+η^T= 0}. The stresses and displacements are sought in

H := {τ∈L²(Ω_F)^d×d : divτ∈L²(Ω_F)^d, τn = 0 onΓ_t}, L := L²(Ω_F)^d.

H is a Hilbert space when endowed with the norm kτkdiv := (kτk²0,ΩF+ kdivτk²0,ΩF)^1/2.

In summary, the variational form of the problem inΩFreads: Given a body load f and a displacementϕonΓ, find (u, σ, γ)∈L ×H ×W satisfying

R

ΩFτ :C⁻¹σdΩ+R

ΩF divτ·u dΩ +R

ΩFτ:γdΩ = R

Γτn·ϕds ∀τ ∈H R

ΩF divσ·vdΩ = −R

ΩFf ·vdΩ ∀v∈L R

ΩFσ:ηdΩ = 0 ∀η∈W .

(2.2)

The identity (2.2)_a is derived above, (2.2)_b is a weak form of (2.1)_a, and (2.2)_c is a weak form of the symmetry ofσ.

We assume that the body load f ∈L²(Ω)^d vanishes onΩB for simplicity and that the Lam´e coefficients are constant onΩ_B. Then, at any point x ∈Ω_B, the displacement field can be represented by the Betti formula

u(x ) =− Z

Γ

G(x,y) T_yu(y) ds_y+ Z

Γ

T_yG(x,y)T

u(y) ds_y.

Here T_yu(y) =Cε(u(y)) n(y) is the traction corresponding to u at a point y ∈ Γ, and TyG(x,y) are the columnwise tractions of G(x,y) at y. G(x,y) is the fundamental solution and equals

λ+3µ 4πµ(λ+2µ)

n

log_|x−¹_y|I +_λ^λ₊₃^+µ_µ^(x−_|x^y)(x_−y⁻_|2^y)T

o if d = 2,

λ+3µ 8πµ(λ+2µ)

n 1

|x−y|I +_λ^λ₊₃^+µ_µ^(x−_|x^y)(x_−y⁻_|3^y)T

o if d = 3.

Letting x →Γ we obtain with the classical jump relations the boundary integral equation

1

2u =−Vt + Ku (2.3)

with t (y) = Tyu(y) and the integral operators (Vt )(x ) =

Z

Γ

G(x,y) t (y) ds_y, x ∈Γ,

(Ku)(x ) = Z

Γ

T_yG(x,y)T

u(y) ds_y, x∈Γ .

Applying the traction operator T_x we get another boundary integral equation

1

2t =−K⁰t−Wu (2.4)

where

(K⁰t )(x ) = Z

Γ

TxG(x,y) t (y) dsy, x ∈Γ , (Wu)(x ) = −Tx

Z

Γ

TyG(x,y)T

u(y) dsy, x∈Γ .

The symmetric coupling of (2.2) with the integral equations (2.3) and (2.4) is performed as follows: Introduce a new variable ϕ:= u|Γ belonging to the trace space

H^1/2:= H^1/2(Γ)^d :={u|_Γ : u∈H¹(Ω)^d},

and use continuity of the tractions onΓ, i.e., t =σn. Then, (2.3) reads ϕ=−V (σn) + (¹₂I + K )ϕ

and this is inserted in the right-hand side of (2.2)_a while (2.4) reads Wϕ+ (¹₂I + K⁰)(σn) = 0

and this is added in a weak form to (2.2). (This coupling is in a sense ‘dual’ to the more classical approach, see e.g. [4, 6, 7], where a new variable is introduced for σn on the interface while for u|Γ continuity is used.) The resulting weak formulation is rewritten in a saddle point structure: Find (σ, ϕ,u, γ) ∈ H × H^1/2×L ×W such that

a(σ, ϕ;τ, ψ) + b(τ; u, γ) = 0 b(σ;v, η) = −R

ΩFf ·vdΩ (2.5)

for all (τ, ψ , v, η)∈H×H^1/2×L ×W . Here, a(σ, ϕ;τ, ψ) :=

Z

ΩF

τ:C⁻¹σdΩ +hτn, V (σn)i −

τn,(¹₂I + K )ϕ

− hψ ,Wϕi −

ψ ,(¹₂I + K⁰)(σn) b(σ;v, η) :=

Z

ΩF

divσ·vdΩ+ Z

ΩF

σ:ηdΩ.

Throughout this paper, hϕ , ψi denotes the extension of the L²-scalar product R

Γϕ·ψds to the duality in H^−1/2×H^1/2; H^s := H^s(Γ)^d.

Remark 2.1. It is known that the mappings V : H^−1/2 → H^1/2, K : H^1/2 → H^1/2, K⁰ : H^−1/2 →H^−1/2, W : H^1/2 →H^−1/2 are well defined, linear and continuous (cf., e.g., [5]). V and W are symmetric; K⁰ is dual to K and W is positive semi-definite and ker W = kerε|Γ, i.e., the kernel of W consists of the (linearized) rigid body motions. W is positive definite on (H^1/2/kerε)². For d = 3, V is positive definite. For d = 2, V is positive definite when restricted to H₀^−1/2 where

H₀^s:={w∈H^s : Z

Γ

wds = 0} ≡H^s/R^d.

We refer to [7] for proofs in case d = 3 and mention that the proofs work verbatim in case d = 2 (provided the radiation condition gives a sufficiently strong decay which is guaranteed owing to the restriction on H₀^−1/2).

Remark 2.2. We note thatσn is defined in H^−1/2 via Green’s formula even if σ ∈ H as follows. Given v ∈ H^1/2 extend it to some v ∈ H¹(ΩF)^d with v|Γu = 0. Then, let

hv , σni:=

Z

ΩF

σ: gradvdΩ+ Z

ΩF

v·divσdΩ.

The right-hand side is well defined and depends linearly and continuously on v andσ. Furthermore,

kσnk₋1/2,Γ ≤Ckσkdiv. (2.6)

In view of the remarks, a is a symmetric and continuous bilinear form on (H ×H^1/2)², and b is continuous onH ×(L ×W ).

Theorem 2.1. For every f ∈L²(ΩF)^dthe saddle point problem (2.5) has a unique solution satisfying

kσk0,ΩF+ kdivσk0,ΩF+ kϕk1/2,Γ+ kuk0,ΩF+ kγk0,ΩF ≤Ckfk0,ΩF

(2.7)

with a positive constant C which is independent of f .

Proof. We apply the theory of saddle point problems, cf., e.g., [3, Sect. II.1]. It is sufficient to verify surjectivity of b and the inf–sup condition on a. As it is well known, the bilinear form b has the following surjectivity property: For all (v, η)∈L ×W there is someτ∈H satisfying

divτ=v and asτ :=¹₂(τ−τ^T) =η.

(2.8)

Therefore, it remains to prove that, for a constantα >0, inf

(σ,ϕ) sup

(τ,ψ)

a(σ, ϕ;τ, ψ)

(kσk0,ΩF+ kϕk1/2,Γ)(kτk0,ΩF+ kψk1/2,Γ)≥α (2.9)

where in inf(σ,ϕ) and sup_(τ,ψ) the nonzero arguments run through ker B , ker B :={τ ∈H : divτ= 0 and asτ= 0} ×H^1/2.

We will prove (2.9) in two steps partly arguing as in [4]: Given (σ, ϕ)∈ker B letϕ=ϕ0+ r0 whereϕ0∈H^1/2/kerεand r0∈kerε, kerεthe (linearized) rigid body motions. Let t_c∈R^d be defined by

Z

Γ

(σn−tc) ds = 0.

Then, letσ_c:=Cε(u_c) where u_c∈H¹(Ω_F)^d is the solution of divσ_c= 0 onΩ_F whileσcn = tconΓ,σcn = 0 onΓt and uc= 0 onΓu. Furthermore, letϕc∈R^d

satisfy

V (σn−tc) + (¹₂I + K )ϕ ,tc

− hσn,uci=htc, ϕci. (2.10)

(If t_c= 0, every ϕ_c ∈R^d solves (2.10).) Note thatϕ_c can be chosen such that kϕck1/2,Γ ≤C (kσk0,ΩF+ kϕk1/2,Γ). Thus,

kσ−σck0,ΩF+ kϕ−ϕck1/2,Γ ≤C (kσk0,ΩF+ kϕk1/2,Γ) (2.11)

where C >0 is independent ofσandϕ. Integration by parts shows Z

ΩF

σc:C⁻¹σdΩ=hσn,uci. (2.12)

Then, using (2.10) and Kϕc= ¹₂ϕc, a(σ, ϕ;σ−σc,−ϕ+ϕc) =

Z

ΩF

σ:C⁻¹σdΩ

+hσn−tc,V (σn−tc)i+hϕ , Wϕi. Since W defines a positive definite bilinear form on (H^1/2/kerε)² and since V is positive definite on (H₀^−1/2)² (recallσn−tc∈H₀^−1/2 by definition of tc) we obtain

a(σ, ϕ;σ−σc,−ϕ+ϕc)≥C1(kσk²0,ΩF+ kϕ0k²1/2,Γ) (2.13)

with C₁>0 depending only on C and W .

Now we show that for every rigid body motion r0 there is a stress field τ0

with (τ₀,0)∈ker B and

a(0,r0;τ0,0) =hr0,r0i. (2.14)

For example, let τ0 ∈H satisfy divτ0 = 0 in ΩF and τ0n =−r0 on Γ. Since Kr0= ¹₂r0, we conclude (2.14).

In the second step we assume for contradiction that (2.9) is false. Hence we may find a sequence (σj, ϕj) in ker B with kσjk0,ΩF+ kϕjk1/2,Γ = 1 and

sup

(τ,ψ)∈ker B\{0}

a(σ_j, ϕ_j;τ, ψ)/(kτk0,ΩF+ kψk1/2,Γ)<1/j (2.15)

for all j . Letϕ_j =ϕ_0j + r_j whereϕ_0j ∈H^1/2/kerεand r_j ∈kerε.

According to (2.11), (2.13) and (2.15) we get

0 = lim

j→∞(kσ_jk0,ΩF+ kϕ_0jk1/2,Γ)

and hence, since rj is a bounded sequence in a finite dimensional space, (σ_j, ϕ_j)→(0,r₀) as j → ∞

(2.16)

for a subsequence (which is not relabeled) and a rigid body motion r0. In particu- lar, kr0k1/2,Γ = 1. Letτ0satisfy (2.14) with r0as in (2.16). Since a is continuous, (2.16) shows

0 = lim

j→∞

a(σj, ϕj;τ0,0) kτ₀k0,ΩF

= a(0,r0;τ0,0) kτ₀k0,ΩF

>0

owing to (2.14) and r₀/= 0. This contradiction verifies (2.9); the proof is finished.

u t

3. Approximation and convergence

For the finite element method we consider a regular family of triangulationsTh of Ω¯F. Two different triangles (resp. tetrahedrons) inTh are either disjoint or have one common side or edge or vertex. Let the sides of finite elements in Th on the interfaceΓ define a partitionEh ofΓ and takeEh as boundary elements for simplicity. These partitions give rise to finite-dimensional subspacesHh, H_h^1/2, Lh andWh which are assumed to satisfy (H1)–(H3); examples are considered in Sect. 5.

(H1) (Conformity and approximation property) There holds Hh×H_h^1/2×Lh ×Wh ⊂ H×H^1/2×L ×W

andR^d ⊂H_h^1/2. For all (τ, ψ , v, η)∈H ×H^1/2×L ×W , 0 = lim

h→0

τhinf∈Hh

kτ−τ_hkdiv + inf

ψh∈H_h^1/2

kψ −ψ_hk1/2,Γ

+ inf

vh∈Lh

kv−vhk0,ΩF+ inf

ηh∈Wh

kη−ηhk0,ΩF

.

(H2) (Equilibrium condition) For eachτh ∈Hh, the condition b(τ_h;v_h, η_h) = 0 ∀(v_h, η_h)∈Lh ×Wh

implies divτ_h = 0.

Remark 3.1. Note that a piecewise polynomial function τh satisfies divτh ∈ L²(Ω_F)^d if and only if the tractionsτ_hn are continuous across interelement sides (which is thus implied byHh ⊂H).

(H3) (inf–sup condition) The bilinear form b satisfies the inf–sup condition in some mesh-dependent norm. That is there exist a norm k·k_H,hon some space H(h)⊂H withHh ⊂H(h)and

kτhk0,ΩF ≤ kτhk_H,h≤Ckτhkdiv

(3.1)

for allτ_h ∈Hh and a norm k(·,·)kh inL(h)×Wh withLh ⊂L(h)⊂L such that a and b are uniformly continuous (i.e., with h-independent bounds), and there is a positive constantβ such that for all (vh, ηh)∈Lh ×Wh

sup

τh∈H0h\{0}

b(τh;vh, ηh) kτhkH,h

≥βk(v_h, η_h)kh. (3.2)

where H0h :={τh ∈ Hh : τhn|Γ = 0}. Both C andβ are assumed to be independent of h.

Remark 3.2. To ensure the continuity of a, we need the estimate kτnk₋1/2,Γ ≤Ckτk_H,h ∀τ∈H(h)

(3.3)

Note that (3.3) does not affect (3.2) owing toH0h.

Example 3.1. Let Th be a finite element triangulation of ΩF and letSh denote the set of sides in the interior ofΩ_FandSh¯the set of all finite element sides. The value of the jump ofvacross an interelement side S is denoted by [v], while the diameter of S is h_S. We follow Stenberg [10, 11]: The first mesh-dependent norm

kτk²_H,h := kτk²0,ΩF+ X

S∈Sh¯

h_S Z

S

|τn|²ds

is defined on the ‘intermediate’ space

H(h):={τ∈H :τn∈L²(S )^d ∀S ∈Sh¯}, and the second

k(v, η)k²h := X

T∈Th

(kε(v)k²0,T + kη−ω(v)k²0,T)

+ X

S∈Sh

h_S⁻¹ Z

S

|[v]|²ds + X

S⊂Γu

h_S⁻¹ Z

S

|v|²ds

onL(h)×Wh withω(v) := ¹₂( gradv−( gradv)^T) and

L(h):={v∈L²(Ω)^d :v|T ∈H¹(T )^d ∀T ∈Th}.

On Hh, the norms k·kH,h and k·k0,ΩF are uniformly equivalent, i.e., for all τh∈Hh

kτhk0,ΩF ≤ kτhk_H,h ≤Ckτhk0,ΩF. (3.4)

In these norms, the bilinear forms a and b are continuous, i.e., there exists C >0 such that for all (σ, ϕ), (τ, ψ)∈H(h)×H^1/2

a(σ, ϕ;τ, ψ)≤C (kσk_H,h+ kϕk1/2,Γ)(kτk_H,h+ kψk1/2,Γ) (3.5)

and for all (τ, v, η)∈H(h)×L(h)×W

b(τ;v, η)≤Ckτk_H,hk(v, η)kh. (3.6)

To show continuity of a, it remains to verify (3.3). Recall kτnk₋1/2,Γ = sup

v∈H1/2(Γ)^d

hτn, vi kvk1/2,Γ.

Eachv∈H^1/2(Γ)^dcan be extended tov∈H¹(Ω_F)^dwithv|_Γu = 0 andkvk1,ΩF ≤ Ckvk1/2,Γ. Integration by parts yields

hτn, vi= Z

ΩF

τ: gradvdΩ+ b(τ;v,0).

Further, k(v,0)kh ≤ kvk1,ΩFsince the jumps [v] on the interelement sides vanish.

Thus, using the continuity of b in the mesh-dependent norms, we obtain hτn, vi ≤CkτkH,hkvk1,ΩF.

This proves (3.3).

Remark 3.3. We assume throughout this paper that the solution to (2.5) satisfies σ ∈ H(h) and u ∈ L(h) for all f ∈ L²(ΩF)^d, which is a regularity condition.

In the situation of Example 3.1, the condition u ∈ H^s(ΩF)^d with s > 3/2 is sufficient.

The discretized saddle point problem consists in finding (σh, ϕh,uh, γh) ∈ Hh×H_h^1/2×Lh×Whsuch that for all (τ_h, ψ_h, v_h, η_h)∈Hh×H_h^1/2×Lh×Wh

a(σh, ϕh;τh, ψh) + b(τh; uh, γh) = 0 b(σ_h;v_h, η_h) = −R

ΩFf ·v_hdΩ (3.7)

Theorem 3.1. Assuming (H1)–(H3), the discrete problem (3.7) has exactly one solution. There exists some h-independent constant C >0 such that

kσ−σ_hk_H,h+ kϕ−ϕ_hk1/2,Γ

≤ C

τhinf∈Hh

kσ−τhk_H,h+ inf

ψh∈H_h^1/2

kϕ−ψhk1/2,Γ+ inf

ηh∈Wh

kγ−ηhk0,ΩF

.

To derive L²-estimates for the displacements, we need a regularity–approxi- mation estimate and a mapping Ph :L →Lh×Wh.

(H4) (Regularity–approximation property) In the mesh-dependent norms we have the following estimate where%(h) is an h-dependent constant such that

inf

(τh,ψh,vh,ηh)∈Hh×H_h^1/2×Lh×Wh

k(σ−τh, ϕ−ψh,u−vh, γ−ηh)kh

≤ %(h)kfk0,ΩF

(3.8)

for all f ∈L²(ΩF)^d and (σ, ϕ,u, γ) solving (2.5) (cf. Theorem 2.1). The mesh- dependent norm k·kh in (3.8) is defined by

k(τ, ψ , v, η)kh := kτk_H,h+ kψk1/2,Γ+ k(v, η)kh. (3.9)

Lemma 3.1. There is a linear mapping Ph :L →Lh×Wh which satisfies b(τ_h;v,0) = b(τ_h; P_h(v))

for allτh ∈Hh and allv∈L.

Modifications of this mapping Ph are frequently used in the literature for a proof of the inf–sup condition; we conversely may construct it from (H3). In particular cases, P_h is known explicitly (see Sect. 5).

Theorem 3.2. Assuming (H1)–(H4), there exists some h-independent constant C >0 such that, with%(h) as in (H4),

kPh(u)−(uh, γh−γ)k0,ΩF ≤C·%(h)·

τhinf∈Hh

kσ−τhk_H,h

+ inf

ψh∈H_h^1/2

kϕ−ψ_hk1/2,Γ + inf

ηh∈Wh

kγ−η_hk0,ΩF

.

Remark 3.4. To obtain an estimate for ku−u_hk0,ΩF, let ( ˜u,γ) := P˜ _h(u) and observe that

ku−uhk0,ΩF ≤ ku−u˜k0,ΩF+ kPh(u)−(uh, γh−γ)k0,ΩF. (3.10)

The last term is estimated by Theorem 3.2, and, according to the specific P_h, an a priori estimate of ku−u˜k0,ΩF is available.

Proofs are given in Sect. 4 while examples are studied in Sect. 5.

4. Proofs

The proofs of Theorem 3.1 and 3.2 use several lemmas where Z := {τ∈H : divτ= 0 andτ =τ^T}

Zh := {τh ∈Hh : b(τh;vh, ηh) = 0 ∀(vh, ηh)∈Lh×Wh} ker Bh := Zh×H_h^1/2.

Throughout this section, C >0 denotes a generic h-independent constant.

Lemma 4.1. [11] For allτ_h ∈Z_h Z

ΩF

τ_h :C⁻¹τ_hdΩ≥Ckτ_hk²0,ΩF. Lemma 4.2. For all ¯σ∈Z there exists ¯σ_h ∈Z_h satisfying

kσ¯ −σ¯hkdiv ≤C inf

τh∈Hh

kσ¯−τhk0,ΩF.

Proof. The result follows as in [11]; for related results we refer to [3, Proposition II.2.5] and [2, Remark III.4.6]. We give a proof for completeness and to stress that (H1)–(H3) are sufficient (even for different situations on the boundary). Let σ˜ be the best approximant to ¯σ∈Z inHh with respect to the L²(ΩF)-norm. By Lemma 4.1 and (H3), the mixed finite element problem

Z

ΩF

σ¯h :C⁻¹τhdΩ+ b(τh; ¯uh,γ¯h) = L1(τh) b( ¯σh;vh, ηh) = L2(vh, ηh)

(with linear forms L₁,L₂) satisfies ellipticity and inf–sup conditions in mesh- dependent norms so that we have a unique solution ( ¯σh,u¯h,γ¯h)∈Hh×Lh×Wh

satisfying

Z

ΩF

( ¯σh−σ) :¯ C⁻¹τhdΩ+ b(τh; ¯uh,γ¯h) = 0 (4.1)

b( ¯σh−σ;˜ vh, ηh) = 0 (4.2)

for all (τh, vh, ηh)∈Hh×Lh×Wh. Note that ¯σh−σ˜ ∈Zh according to (4.2).

Hence, Lemma 4.1 yields Ckσ¯_h−σ˜k²0,ΩF ≤

Z

ΩF

( ¯σ_h−σ) :˜ C⁻¹( ¯σ_h−σ) dΩ˜

= Z

ΩF

( ¯σ−σ) :˜ C⁻¹( ¯σh−σ) dΩ˜ owing to (4.1). Thus, by Cauchy’s inequality,

kσ¯_h −σ˜k0,ΩF ≤Ckσ¯−σ˜k0,ΩF.

Then, the triangle inequality and (H2) finish the proof of the lemma. ut Lemma 4.3. There is a constantα >0 such that

inf

(σh,ϕh) sup

(τh,ψh)

a(σh, ϕh;τh, ψh)

(kσ_hk_H,h+ kϕ_hk1/2,Γ)(kτ_hk_H,h+ kψ_hk1/2,Γ) ≥α (4.3)

where the nonzero arguments in inf(σh,ϕh)and sup_(τ

h,ψh) belong to ker Bh.

Proof. According to the equilibrium condition (H2) it is sufficient to consider the norm k·k0,ΩF instead of k·k_H,h (because, by (3.1), kτhk0,ΩF ≤ kτhk_H,h ≤ Ckτ_hkdiv = Ckτ_hk0,ΩF for allτ_h ∈Z_h). We proceed as in the proof of Theorem 2.1. Assuming that (4.3) is false we find a sequence of meshes, discrete spaces and discrete functions (σj, ϕj) in ker Bh_j with kσjk0,ΩF+ kϕjk1/2,Γ = 1 and

sup

(τ,ψ)∈ker B_hj\{0}

a(σj, ϕj;τ, ψ)/(kτk0,ΩF+ kψk1/2,Γ)<1/j (4.4)

for all j . Letϕj =ϕ0j + rj whereϕ0j ∈H^1/2/kerεand rj ∈kerε.

Since (σj, ϕj) are bounded inH×H^1/2we may extract a weakly convergent subsequence (not relabeled); so assume,

(σ_j, ϕ_0j)*(σ, ϕ₀) (weakly) inH ×H^1/2/kerε

for some (σ, ϕ0) ∈ H ×H^1/2/kerε. By (H2) we have divσj = 0 and so divσ= 0. Define t_j ∈R^d by

Z

Γ

(σjn−tj) ds = 0.

The weak convergence of (σ_j) causes strong convergence of t_j inR^d so that lim

j→∞t_j =: ¯t∈R^d. (4.5)

Let ¯σ ∈Z satisfy ¯σn = ¯t on Γ. Assume ¯t = 0 first. Then, for any t/ j /= 0, we choose ¯σj ∈Zh_j as in Lemma 4.2. If ¯t = 0 we choose ¯σj = 0 = ¯σ. Note ¯σj →σ¯ strongly inH as j→ ∞ because of Lemma 4.2, (H1) and (4.5).

By R^d ⊂H_h^1/2 (cf., (H1)) we have −ϕj + cj ∈H_h^1/2 for each cj ∈R^d. We define cj ∈R^d with minimal Euclid norm in R^d such that

hcj, tji = htj,V (2σjn−tj)i (4.6)

−

σ¯jn,V (σjn)−(¹₂I + K )ϕj

− Z

ΩF

σ¯j :C⁻¹σjdΩ.

Note that ( ¯σj,cj) are bounded inH ×H^1/2. Using (4.6) we compute a(σj, ϕj;σj−σ¯j,−ϕj+ cj) =

Z

ΩF

σj :C⁻¹σjdΩ +hσjn−tj,V (σjn−tj)i+hϕj,Wϕji

≥ C

kσjk²0,ΩF+ kϕ0jk²1/2,Γ

where C >0 depends on W and the constant in Lemma 4.1 only. Since (4.4) and kσj−σ¯jk0,ΩF+ k−ϕj + cjk1/2,Γ ≤C the above estimate proves

(σj, ϕj)→(0,r0) (strongly) inH ×H^1/2, (σ, ϕ0) = 0; r0∈kerε. Thus, kr0k1/2,Γ = 1.

Let ˆσ∈Z satisfy ˆσn =−r₀ onΓ. By Lemma 4.2 we find a discrete stress field ˆσj ∈ Zh_j which, by (H1), converges towards ˆσ in (H,k·kdiv). Letting (τ, ψ) := ( ˆσ_j,0)∈ker B_hj \ {0}in (4.4),

a(σ_j, ϕ_j; ˆσ_j,0)<1

j kσˆ_jk0,ΩF.

By strong convergence, for j → ∞, a(0,r0; ˆσ,0) ≤0. But, by construction of σˆ and because kr0k1/2,Γ = 1, a(0,r0; ˆσ,0) = hr0,r0i > 0. This contradiction proves the lemma. ut

Proof of Lemma 3.1. Let Z_h⁰denote the set of functionalsΦonHh withΦ(τh) = 0 for allτ_h ∈Z_h. By (H3), the mapping

Lh ×Wh → Z_h⁰ (v_h, η_h) 7→ b(·;v_h, η_h)

is an isomorphism [2, 3] (assuming Lh ×Wh be endowed with the mesh- dependent norm). For all v ∈ L, (H2) implies b(·;v,0)|_Hh ∈ Z_h⁰. Hence, for eachv∈L, there exists P_h(v) := (v_h, η_h)∈Lh×Wh satisfying

b(·;vh, ηh)|Hh = b(·;v,0)|Hh. ut

Proof of Theorem 3.1. We follow [11, Proof of Theorem 3.1] and let ( ˜σ,ϕ,˜ γ) be˜ the best approximant to (σ, ϕ, γ) inHh×H_h^1/2×Wh. LetHh×H_h^1/2×Lh×Wh

be endowed with the norm (3.9). Because of (H3) and Lemma 4.3 we get the inf–

sup condition for the discrete spaces. With the theory of saddle point problems there exists (τh, ψh, vh, ηh) inHh×H_h^1/2×Lh×Wh with

k(τh, ψh, vh, ηh)kh ≤C and

kσ˜−σhkH,h + kϕ˜−ϕhk1/2,Γ+ kPh(u)−(uh, γh −γ)˜ kh

≤ a( ˜σ−σh,ϕ˜−ϕh;τh, ψh)

+ b(τh; P_h(u)−(u_h, γh−γ)) + b( ˜˜ σ−σh;vh, ηh)

= a( ˜σ−σ,ϕ˜−ϕ;τ_h, ψ_h)

+ b( ˜σ−σ;v_h, η_h)−b(τ_h; 0, γ−γ)˜

where we used the discrete equations and b(τh; Ph(u)−(u,0)) = 0 according to Lemma 3.1. Since a and b are bounded (in the discrete norms)

kσ˜ −σ_hk_H,h+ kϕ˜−ϕ_hk1/2,Γ + kP_h(u)−(u_h, γ_h−γ)˜ kh

≤ C

kσ˜−σk_H,h+ kϕ˜−ϕk1/2,Γ + kγ˜−γk0,ΩF

.

This and the triangle inequality prove

kσ−σ_hk_H,h+ kϕ−ϕ_hk1/2,Γ+ kP_h(u)−(u_h, γ_h−γ)kh

(4.7)

≤ C

τhinf∈Hh

kσ−τhk_H,h+ inf

ψh∈H_h^1/2

kϕ−ψhk1/2,Γ+ inf

ηh∈Wh

kγ−ηhk0,ΩF

.

The proof of Theorem 3.1 is finished. ut

Proof of Theorem 3.2. We follow [11], assume (H4) and continue in the notations of the proof of Theorem 3.1. With Theorem 2.1 we find (Π, Ψ,z, µ) inH × H^1/2×L ×W satisfying, for all (τ, ψ , v, η) inH×H^1/2×L ×W ,

a(Π, Ψ;τ, ψ) + b(τ; z, µ) = 0 b(Π;v, η) =

Z

ΩF

Ph(u)−(uh, γh −γ)

·(v, η) dΩ.

Taking (τ, ψ , v, η) = (σ−σh, ϕ−ϕh,Ph(u)−(uh, γh−γ)) we obtain kPh(u)−(uh, γh−γ)k²0,ΩF = b(Π; Ph(u)−(uh, γh−γ))

= a(Π, Ψ;σ−σh, ϕ−ϕh)

+ b(σ−σh; z, µ) + b(Π; Ph(u)−(uh, γh−γ))

= a(σ−σ_h, ϕ−ϕ_h;Π−Π, Ψ˜ −Ψ˜)

+ b(σ−σ_h; z −˜z, µ−µ) + b(Π˜ −Π; P˜ _h(u)−(u_h, γ_h−γ)) using Galerkin equations and the definition of Ph in Lemma 3.1 for the best approximants ( ˜Π,Ψ,˜ ˜z,µ) in˜ Hh×H_h^1/2×Lh×Wh to (Π, Ψ,z, µ). Considering norms and using (4.7) forσ−σ_h,ϕ−ϕ_h and P_h(u)−(u_h, γ_h−γ) we gain

kP_h(u)−(u_h, γ_h−γ)k²0,ΩF ≤ C · k(Π−Π, Ψ˜ −Ψ,˜ z−˜z, µ−µ)˜ kh

·(kσ˜−σk_H,h+ kϕ˜−ϕk1/2,Γ+ kγ˜−γk0,ΩF).

By (H4), we have an a priori estimate of the form

k(Π−Π, Ψ˜ −Ψ,˜ z−˜z, µ−µ)˜ kh ≤%(h)kPh(u)−(uh, γh−γ)k0,ΩF. Using this in the former estimate and dividing by kPh(u)−(uh, γh−γ)k0,ΩF we conclude the claimed estimate. ut

5. A family of elements

In the sequel we describe a family of finite element spaces due to Stenberg [11].

For each tetrahedron (resp. triangle) T ∈Th we define a bubble function bT by b_T(x ) =

Yd

i =0

λ_i(x ),

whereλ0, . . . , λd are the barycentric coordinates in T . By Pk(T ) we denote the space of polynomials of degree≤k on T . For d = 3 we define

Bl(T ) := {(τij) : (τi 1, . . . , τi 3) = curl (bTwi 1, . . . ,bTwi 3), wij ∈Pl(T ), i,j = 1,2,3}