Stabilized Mixed Finite Element Methods for Linear Elasticity on Simplicial Grids in ℝⁿ

1 Introduction

Assume that Ω⊂ℝn is a bounded polytope. Denote by 𝕊 the space of all symmetric n×n tensors. The Hellinger–Reissner mixed formulation of the linear elasticity under the load 𝒇∈𝑳2⁢(Ω;ℝn) is given as follows: Find (𝝈,𝒖)∈𝚺×𝑽:=𝑯⁢(div,Ω;𝕊)×𝑳2⁢(Ω;ℝn) such that

where

with 𝒜 being the compliance tensor of fourth order defined by

Here 𝜹:=(δi⁢j)n×n is the Kronecker tensor, tr is the trace operator, and positive constants λ and μ are the Lamé constants. It is arduous to design 𝑯⁢(div,Ω;𝕊) conforming finite element with polynomial shape functions due to the symmetry requirement of the stress tensor. Hence composite elements were one main choice to approximate the stress in the last century (cf. [7, 30, 46, 56]).

In the early years of this century, Arnold and Winther constructed the first 𝑯⁢(div,Ω;𝕊) conforming mixed finite element with polynomial shape functions in two dimensions in [10], which was extended to tetrahedral grids in three dimensions in [1, 4] and simplicial grids in any dimension in [42]. In those elements, the displacement space is approximated by 𝑳2⁢(Ω;ℝn)-Pk-1 while the stress space is approximated by the space of functions in 𝑯⁢(div,Ω;𝕊)-Pk+n-1 whose divergence is in 𝑳2⁢(Ω;ℝn)-Pk-1 for k≥2. Recently, Hu and Zhang showed that the more compact pair of 𝑯⁢(div,Ω;𝕊)-Pk and 𝑳2⁢(Ω;ℝn)-Pk-1 spaces is stable on triangular and tetrahedral grids for k≥n+1 with n=2,3 in [40, 41]. And Hu generalized those stable finite elements to simplicial grids in any dimension for k≥n+1 in [37]. One key observation there is that the divergence space of the 𝑯⁢(div) bubble function space of polynomials on each element is just the orthogonal complement space of the piecewise rigid motion space with respect to the discrete displacement space. Then the discrete inf-sup condition was proved for k≥n+1 through controlling the piecewise rigid motion space by 𝑯1⁢(Ω;𝕊)-Pk space. It is, however, troublesome to prove that the pair of 𝑯⁢(div,Ω;𝕊)-Pk and 𝑳2⁢(Ω;ℝn)-Pk-1 is still stable for 1≤k≤n. For such a reason, Hu and Zhang enriched the 𝑯⁢(div,Ω;𝕊)-Pk space with 𝑯⁢(div,Ω;𝕊)-Pn+1 face-bubble functions of piecewise polynomials for each n-1 dimensional simplex in [42]. Gong, Wu and Xu constructed two types of interior penalty mixed finite element methods by using nonconforming symmetric stress approximation in [31]. The stability of those nonconforming mixed methods is ensured by 𝑯⁢(div) nonconforming face-bubble spaces. An interior penalty mixed finite element method using Crouzeix–Raviart nonconforming linear element to approximate the stress was studied in [21]. To get rid of the vertex degrees of freedom appeared in the 𝑯⁢(div,Ω;𝕊)-conforming elements and make the resulting mixed finite element methods hybridizable, nonconforming mixed elements on triangular and tetrahedral grids were developed in [11, 5, 32, 57]. On rectangular grids, we refer to [3, 36, 38, 12, 23] for symmetric conforming mixed finite elements and [39, 47, 58, 59] for symmetric nonconforming mixed finite elements. For keeping the symmetry of the discrete stress space and relaxing the continuity across the interior faces of the triangulation, many discontinuous Galerkin methods were proposed in [20, 27, 24, 43, 44, 45], hybridizable discontinuous Galerkin methods in [51, 35], weak Galerkin methods in [55, 22], hybrid high-order method in [29]. For the weakly symmetric mixed finite element methods for linear elasticity, we refer to [2, 6, 53, 48, 8, 15, 49, 26, 34, 50, 33, 28, 9].

In this paper, we intend to design stable mixed finite element methods for the linear elasticity using as few global degrees of freedom as possible. To this end, two classes of stabilized mixed finite element methods on simplicial grids in any dimension are proposed. In the first one, we use the pair of 𝑯⁢(div,Ω;𝕊)-Pk and 𝑳2⁢(Ω;ℝn)-Pk-1 constructed in [37] to approximate the stress and displacement for 1≤k≤n. To simplify the notation, we shall use superscript (⋅)div for 𝑯⁢(div,Ω;𝕊) conforming elements, (⋅)-1 for discontinuous elements, and (⋅)0 for 𝑯1⁢(Ω;ℝn) or 𝑯1⁢(Ω;𝕊) continuous elements. Instead of enriching Pkdiv elements with face bubble functions of piecewise polynomials as in [42, 31], we include a jump stabilization term into the Hellinger–Reissner mixed formulation to make the discrete method stable, inspired by the discontinuous Galerkin methods constructed in [24] for the linear elasticity problem. The discrete inf-sup condition in a compact form is established with the help of a partial inf-sup condition (2.1) estabilished in [37, 40, 41] and a well-tailored interpolation operator for the stress. Then we show the a priori error analysis for the resulting stabilized Pkdiv-Pk-1-1 elements.

In the second class of stabilized mixed finite element methods, we adopt the stabilization technique suggested in [19] and use 𝑯01⁢(Ω;ℝn)-Pk to approximate the displacement space. The merit of this stabilization technique is that the coercivity condition for the bilinear form related to the stress holds automatically, thus we only need to focus on the discrete inf-sup condition. To recover the inf-sup condition, we first employ 𝑯1⁢(Ω;𝕊)-Pk enriched with (k+1)-st order 𝑯⁢(div) bubble function space of polynomials on each element to approximate the stress space. The discrete inf-sup condition is established by using another special interpolation operator for the stress. The a priori error estimate for the (Pk0+Bk+1div)-Pk0 is then derived by the standard theory of mixed finite element methods. The rate of convergence for the displacement in 𝑳2⁢(Ω;ℝn) norm is, however, suboptimal due to the coupling of stress error measured in 𝑯⁢(div)-norm. To remedy this, we use 𝑯⁢(div,Ω;𝕊)-Pk+1 to approximate the stress space instead. The resulting stable finite element pair is the Hood–Taylor type Pk+1div-Pk0. It was mentioned in [19] that it is not known if the Hood–Taylor element in [54, 13, 14] is stable for the linear elasticity. We solve this problem by enriching the Hood–Taylor element space Pk+10 for the stress with the same degree of 𝑯⁢(div) bubble function space of polynomials on each element.

Note that the key component in constructing the previous two interpolation operators is the 𝑯⁢(div) bubble function space of polynomials on each element. We would like to mention that for the stress we obtain the optimal error estimates in 𝑯⁢(div,Ω;𝕊) norm, whereas these error estimates in 𝑳2⁢(Ω;𝕊) norm are suboptimal. To the best of our knowledge, the global degrees of freedom of our methods for the lowest order case k=1 are fewer than those of any existing mixed-type symmetric finite element method for the linear elasticity in the literature. To be specific, the global degrees of freedom for the stress and the displacement for the stabilized mixed finite element methods (3.1)–(3.2), (4.1)–(4.2) and (4.9)–(4.10) with k=1 are, respectively,

Here |𝒱|,|ℰ|,|𝒯| are the numbers of vertices, edges and elements of the triangulation.

When adopting the same degree of polynomial spaces for displacement, the Hood–Taylor type elements Pk+1div-Pk0 and the stabilized elements Pk+1div-Pk-1 share the same convergence rate, which is one order higher than the stabilized elements (Pk0+Bk+1div)-Pk0. It is worth mentioning that to keep the same convergence rate, the stabilized elements Pk+1div-Pk-1 need larger global degrees of freedom than the Hood–Taylor type elements Pk+1div-Pk0.

The rest of this article is organized as follows. We present some notations and definitions in Section 2 for later uses. In Section 3, a stabilized mixed finite element method with discontinuous displacement for the linear elasticity is designed and analyzed. Then we propose a second class of stabilized mixed finite element methods with continuous displacement for the linear elasticity in Section 4. In Section 5, some numerical experiments are given to demonstrate the theoretical results.

2 Preliminaries

Given a bounded domain G⊂ℝn and a non-negative integer m, let Hm⁢(G) be the usual Sobolev space of functions on G, and 𝑯m⁢(G;𝕏) the usual Sobolev space of functions taking values in the finite-dimensional vector space 𝕏 for 𝕏 being 𝕊 or ℝn. The corresponding norm and seminorm are denoted respectively by ∥⋅∥m,G and |⋅|m,G. If G is Ω, we abbreviate them by ∥⋅∥m and |⋅|m, respectively. Let 𝑯0m⁢(G;ℝn) be the closure of 𝑪0∞⁢(G;ℝn) with respect to the norm ∥⋅∥m,G. Denote by 𝑯⁢(div,G;𝕊) the Sobolev space of square-integrable symmetric tensor fields with square-integrable divergence. For any 𝝉∈𝑯⁢(div,Ω;𝕊), we equip the norm

When 𝝉∈𝑯⁢(div,Ω;𝕊) satisfying ∫Ωtr⁡𝝉⁢d⁢x=0, it follows from [16, Proposition 9.1.1] that there exists a constant C>0 such that

which means ∥𝝉∥𝑯⁢(div,𝒜) and ∥𝝉∥𝑯⁢(div) are equivalent uniformly with respect to the Lamé constant λ. Hence the norm ∥⋅∥𝑯⁢(div,𝒜) presented in all of the estimates in this paper can be replaced by the norm ∥⋅∥𝑯⁢(div).

Suppose the domain Ω is subdivided by a family of shape regular simplicial grids 𝒯h (cf. [17, 25]) with h:=maxK∈𝒯h⁡hK and hK:=diam⁢(K). Let ℱh be the union of all n-1 dimensional faces of 𝒯h. For any F∈ℱh, denote by hF its diameter. Let Pm⁢(G) stand for the set of all polynomials in G with the total degree no more than m, and 𝑷m⁢(G;𝕏) denote the tensor or vector version of Pm⁢(G) for 𝕏 being 𝕊 or ℝn, respectively. Throughout this paper, we also use “≲⋯” to mean that “≤C⁢⋯”, where C is a generic positive constant independent of h and the Lamé constant λ, which may take different values at different appearances.

Consider two adjacent simplices K+ and K- sharing an interior face F. Denote by 𝝂+ and 𝝂- the unit outward normals to the common face F of the simplices K+ and K-, respectively. For a vector-valued function 𝒘, write 𝒘+:=𝒘|K+ and 𝒘-:=𝒘|K-. Then define a matrix-valued jump as

On a face F lying on the boundary ∂⁡Ω, the above term is defined by

For each K∈𝒯h, define an 𝑯⁢(div,K;𝕊) bubble function space of polynomials of degree k as

It is easy to check that 𝑩K,1 is merely the zero space. Denote the vertices of simplex K by 𝒙0,…,𝒙n. For any edge 𝒙i⁢𝒙j (i≠j) of element K, let 𝒕i,j be the associated tangent vectors and

It follows from [37] that, for k≥2,

where λi are the associated barycentric coordinates corresponding to 𝒙i for i=0,…,n. Some global finite element spaces are given by

with integer k≥1. It follows from [37, 40, 41] that

where the local rigid motion space and its orthogonal complement space with respect to 𝑷k-1⁢(K;ℝn) on each simplex K∈𝒯h are defined as (cf. [37])

with 𝜺⁢(𝒗):=(∇⁡𝒗+(∇⁡𝒗)T)/2 being the linearized strain tensor.

To introduce an elementwise 𝑯⁢(div) bubble function interpolation operator, we first present the degrees of freedom for 𝚺k,h which are slightly different from those given in [37]. For the ease of notation, we understand Pk=∅ for negative integers k.

It is easy to see that we have the same first set of degrees of freedom as those in [37], whereas the second set of degrees of freedom is different.

Now we present an elementwise 𝑯⁢(div) bubble function interpolation operator. Given 𝝉∈𝑳2⁢(Ω;𝕊), define 𝑰k,hb⁢𝝉∈𝚺k,h as follows: on each simplex K∈𝒯h,

1. for any degree of freedom D in the first set of degrees of freedom in Lemma 2.1,

   D⁢(𝑰k,hb⁢𝝉)=0,

2. for any 𝝇∈𝑷k-2⁢(K;𝕊),

   (2.2)∫K𝑰k,hb⁢𝝉:𝝇⁢d⁢x=∫K𝝉:𝝇⁢d⁢x.

Since the first set of degrees of freedom in Lemma 2.1 completely determines 𝝉⁢𝝂 on ∂⁡K for any 𝝉∈𝑷k⁢(K;𝕊) (cf. [37, Theorem 2.1]), thus 𝑰k,hb⁢𝝉∈𝑩k,h. Applying a scaling argument, we have, for any 𝝉∈𝑳2⁢(Ω;𝕊),

3 A Stabilized Mixed Finite Element Method with Discontinuous Displacement

In this section, we devise a stabilized mixed finite element method for linear elasticity. In [37], the pair of 𝑯⁢(div,Ω;𝕊)-Pk and 𝑳2⁢(Ω;ℝn)-Pk-1 is shown to be stable for k≥n+1. Here we consider the range 1≤k≤n.

With previous preparation, a stabilized mixed finite element method for linear elasticity is defined as follows: Find (𝝈h,𝒖h)∈𝚺k,h×𝑽k-1,h such that

where the jump stabilization term for the displacement is

With this jump stabilization term, a jump seminorm for 𝑽k-1,h+𝑯1⁢(Ω;ℝn) is defined as

We also define the following two norms:

Let 𝑸h be the L2 orthogonal projection from 𝑳2⁢(Ω;ℝn) onto 𝑽k-1,h. The following error estimate holds (cf. [25, 17]):

with integer m≥1. Let 𝑰hSZ be a tensorial or vectorial Scott–Zhang interpolation operator designed in [52], which possesses the following error estimate:

for any 𝝉∈𝑯m⁢(Ω;𝕊) with integer m≥1. Then for each 𝝉∈𝑯1⁢(Ω;𝕊), define

Apparently we have 𝑰h⁢𝝉∈𝚺k,h. And it follows from (2.2) that

To derive a discrete inf-sup condition for the stabilized mixed finite element method (3.1)–(3.2), we rewrite it in a compact way: Find (𝝈h,𝒖h)∈𝚺k,h×𝑽k-1,h such that

where

Similarly, problem (1.1)–(1.2) can be rewritten as

Obviously the bilinear form ℬ is continuous with respect to the norm ∥⋅∥𝑯⁢(div,𝒜)+∥⋅∥0,c. Now we present the following inf-sup condition for (3.6).

The unique solvability of the stabilized mixed finite element method (3.1)–(3.2) is the immediate result of the inf-sup condition (3.8).

Next we show the a priori error analysis for the stabilized mixed finite element method (3.1)–(3.2). Subtracting (3.6) from (3.7), we have the following error equation from the definition of 𝑸h:

for any (𝝉h,𝒗h)∈𝚺k,h×𝑽k-1,h.