A Reduced Local Discontinuous Galerkin Method for Nearly Incompressible Linear Elasticity

Abstract

A reduced local discontinuous Galerkin (RLDG) method for nearly incompressible linear elasticity is proposed in this paper, which is locking-free. RLDG method can be formally regarded as a special case of LDG method with 𝐶₁₁=0. However, RLDG method is actually not covered by LDG method, where 𝐶₁₁ must be chosen to be positive to ensure the stability of LDG method. RLDG method can also be considered as the localization of some symmetric nonconforming mixed finite element method. The implementation of RLDG method is discussed. By introducing a lifting operator as LDG method, RLDG method can be rewritten as primal formulation with unknown displacement only. Next, we obtain that the convergence rates of the approximation to stress tensor in energy norm and displacement in 𝐿²-norm are 𝑂(ℎ^𝑘) and 𝑂(ℎ^𝑘+1), respectively, which are both uniform with respect to 𝜆. Moreover, we obtain a 𝐻(div)-conforming displacement by projecting the displacement and corresponding numerical trace of RLDG method into the Raviart-Thomas element space. And then we analyze the error estimates of this postprocessed displacement in 𝐻(div)-seminorm and 𝐿²-norm, which are also uniform with respect to 𝜆. Finally, some numerical results are shown to demonstrate the theoretical results.

1. Introduction

Volume locking phenomenon will happen when solving nearly incompressible linear elasticity problem by standard displacement method using continuous lower order finite elements (cf. [1–3]); that is, the convergence rates for these lower order conforming finite element methods will deteriorate as the Lamé constant 𝜆 tends to +∞. For this reason, robust numerical methods with respect to 𝜆 are necessary to deal with volume locking phenomenon, and these methods are also called locking-free methods. As for standard displacement method, Babuška and Suri showed detailedly in [3] that it is locking-free whenever the polynomial order 𝑝 of the local shape space is no less than 4 and locked for 𝑝<4 on triangular meshes, and locking phenomenon always takes place for any 𝑝≥1 on quadrilateral meshes. Besides higher order ℎ-version conforming finite element methods, 𝑝-version finite element method in [4], nonconforming finite element methods in [5–8], and mixed finite element methods in [9–11] are also suggested to overcome the volume locking phenomenon.

On the other hand, discontinuous Galerkin (DG) methods which possess the flexibility of constructing feasible local shape function spaces and the advantage of capturing non-smooth or oscillatory solutions effectively have been applied to solve nearly incompressible linear elasticity. Hansbo and Larson analyzed symmetric interior penalty discontinuous Galerkin (SIPDG) methods in [12, 13], whereas a dimensionless penalty parameter appeared in the stabilized term cannot be quantified in theory and can only be determined by numerical experience. A locking-free nonsymmetric interior penalty discontinuous Galerkin (NIPDG) method was developed by Wihler in [14]. And Wihler also designed an adaptive algorithm for NIPDG method in [15], whose robust a posteriori error analysis was given by Bridgeman and Wihler in [15, 16]. A local discontinuous Galerkin (LDG) method based on a first-order system with three mathematical equations was presented and discussed by Cockburn et al. in [17], and a post-processed 𝐻(div)-conforming displacement was achieved by projecting the displacement and corresponding trace of LDG method into the Raviart-Thomas element space. Using symmetric normal-normal continuous element for stress and tangential continuous Nédélec element for displacement, Pechstein and Schöberl proposed a locking-free mixed finite element method in [18], which is in fact some kind of discontinuous Galerkin method. Through introducing a DG derivative, a family of discontinuous Galerkin mixed methods were designed by Shen and Lew in [19] which were proved to be locking-free as well. They also observed that locking might appear if Dirichlet boundary conditions were imposed strongly rather than weakly.

In this paper, we intend to propose a reduced local discontinuous Galerkin (RLDG) method for nearly incompressible linear elasticity. RLDG method can be formally regarded as a special case of LDG method proposed in [20] with 𝐶₁₁=0. However, RLDG method is actually not covered by LDG method, where 𝐶₁₁ must be chosen to be positive to ensure the stability of LDG method, whereas RLDG method drops this requirement by choosing large enough finite element space for stress to match finite element space for displacement. On the other hand, RLDG method can also be considered as the localization of the symmetric nonconforming mixed finite element method given in [21]. Based on this observation, the unique solvability of RLDG method is proved using the similar argument as the method in [21]. Since there is not any continuity across the interior faces for the finite element space for stress, Schur complement is easy to obtain, and thus implementation of RLDG method is not difficult. On the other hand, RLDG method can be rewritten as primal formulation with unknown displacement only through introducing a lifting operator as LDG method in [20]. According to the error equations, an interpolation operator designed in [21] and BDM interpolation operator, we obtain that the convergence rate of the approximation to stress tensor and displacement in energy norm is 𝑂(ℎ^𝑘). And by duality argument, the convergence rate of the approximation to displacement in 𝐿²-norm is shown to be 𝑂(ℎ^𝑘+1). These two convergence rates are both uniform with respect to 𝜆; hence our RLDG method is locking-free. Moreover, we obtain a 𝐻(div)-conforming displacement by projecting the displacement and corresponding numerical trace of RLDG method into the Raviart-Thomas element space, just as the postprocessing technique in [17]. And then we analyze the error estimates of this post-processed displacement in 𝐻(div)-seminorm and 𝐿²-norm, which are also uniform with respect to 𝜆. Finally, some numerical results are shown to demonstrate the theoretical results.

The rest of this paper is organized as follows. The reduced local discontinuous Galerkin method for nearly incompressible linear elasticity is presented in Section 2. In Section 3, we give a robust error analysis for RLDG method. A postprocessing of the displacement of RLDG method is discussed in Section 4. And in Section 5, some numerical results are included to confirm our theoretical convergence orders. In the final section, we give some concluding statements.

2. The RLDG Method for Linear Elasticity

Assume that Ω⊂ℝ^𝑑(𝑑=2,3) is a bounded polygon or polyhedron. Let 𝝈=(𝜎_𝑖𝑗)_𝑑×𝑑 be the stress, let 𝐮=(𝑢₁,…,𝑢_𝑑)^𝑡 be the displacement, and let 𝐟=(𝑓₁,…,𝑓_𝑑)^𝑡 be the body force. Denote by 𝜺(𝐮)∶=(𝜀_𝑖𝑗(𝐮))_𝑑×𝑑 the linearized strain tensor with 𝜀_𝑖j(𝐮)∶=(𝜕𝑢_𝑖/𝜕𝑥_𝑗+𝜕𝑢_𝑗/𝜕𝑥_𝑖)/2, tr the trace operator, and 𝛁⋅ the divergence operator. Consider linear elasticity in the stress-displacement formulation: 𝒜𝝈−𝜺(𝐮)=𝟎inΩ,−𝛁⋅𝝈=𝐟inΩ,𝐮=𝟎on𝜕Ω, (1) where 𝒜 is the compliance tensor of the fourth order defined by 1𝒜𝝈=𝜆2𝜇𝝈−()𝜹.𝑑𝜆+2𝜇tr𝝈 (2) Here, tr is the trace operation of a second-order tensor (a matrix), 𝜹∶=(𝛿_𝑖𝑗)_𝑑×𝑑 the usual Kronecker delta, and the positive constants 𝜆∈(0,+∞) and 𝜇∈[𝜇₀,𝜇₁] are the Lamé coefficients of the elastic material under consideration with 0<𝜇₀<𝜇₁<+∞.

To define a reduced LDG method for solving problem (1), we first introduce some notations frequently used later on. Denote the space of all symmetric 𝑑×𝑑 tensor by 𝕊. Given a bounded domain 𝐺⊂ℝ^𝑑 and a nonnegative integer 𝑚, let 𝐻^𝑚(𝐺) be the usual Sobolev space of functions on 𝐺, and let 𝐇^𝑚(𝐺,𝕏) be the usual Sobolev space of functions taking values in the finite-dimensional vector space 𝕏 for 𝕏 being 𝕊 or ℝ^𝑑. The corresponding norm and seminorm are denoted, respectively, by ‖⋅‖_𝑚,𝐺 and |⋅|_𝑚,𝐺. If 𝐺 is Ω, we abbreviate them by ‖⋅‖_𝑚 and |⋅|_𝑚, respectively. Let 𝐻^𝑚₀(𝐺) be the closure of 𝐶^∞₀(𝐺) with respect to the norm ‖⋅‖_𝑚,𝐺. We also denote by 𝐇(div,𝐺,𝕊) the Sobolev space consisting of all 𝐋²(𝐺,𝕊) functions whose divergence are square integrable. Let 𝐇₀(div,Ω)∶=𝐯∈𝐋²Ω,ℝ^𝑑;𝛁⋅𝐯∈𝐿²(Ω),𝐯⋅𝐧|_𝜕Ω=0. (3)𝑃_𝑚(𝐺) stands for the set of all polynomials in 𝐺 with the total degree no more than 𝑚, and 𝐏_𝑚(𝐺,𝕏) denotes the tensor or vector version of 𝑃_𝑚(𝐺) for 𝕏 being 𝕊 or ℝ^𝑑, respectively.

Let {𝒯_ℎ}_ℎ>0 be a regular family of triangulations of Ω (cf. [22, 23]). For each 𝐾∈𝒯_ℎ, we denote by ℎ_𝐾 the diameter of 𝐾 and by 𝜌_𝐾 the diameter of the biggest ball included in 𝐾. Let ℰ_ℎ be the union of all faces of the triangulation 𝒯_ℎ and let ℰ^𝑖_ℎ be the union of all interior faces of the triangulation 𝒯_ℎ. For a face 𝑒∈ℰ_ℎ, ℎ_𝑒 is set to be the diameter of 𝑒. The triangulations we consider can have hanging nodes, but they have to be regular. Moreover, we restrict the ratio of the sizes of neighbor element domains. To formally state this property, we need to introduce the set ⟨𝐾,𝐾^′⟩ defined as follows: 𝐾,𝐾^′∶=∅,ifmeas𝜕𝐾∩𝜕𝐾^′=0,interiorof𝜕𝐾∩𝜕𝐾^′,otherwise. (4) Thus we assume that there exists a positive constant 𝐶₂<1 such that, for each element 𝐾∈𝒯_ℎ, 𝐶₂≤ℎ_𝐾′ℎ_𝐾≤1𝐶₂∀𝐾^′satisfying𝐾,𝐾^′≠∅. (5) This assumption forbids the situation where the mesh is indefinitely refined in only one of two adjacent subdomains. Denote by 𝛁_ℎ and 𝜺_ℎ the elementwise counterparts of 𝛁 and 𝜺 related to 𝒯_ℎ, respectively.

Let 𝐾⁺ and 𝐾⁻ be two adjacent elements of 𝒯_ℎ. Let 𝐱 be an arbitrary point of the set 𝑒^′=⟨𝐾⁺,𝐾⁻⟩, and let 𝐧⁺ and 𝐧⁻ be the corresponding outward unit normals at that point. For a vector-valued function 𝐯 smooth inside each element 𝐾^±, let us denote by 𝐯^± the trace of 𝐯 on 𝑒^′ from the interior of 𝐾^±. Then we define averages and jumps at 𝐱∈𝑒^′ as follows: {𝐯}=12𝐯⁺+𝐯⁻,[𝐯]=𝐯⁺⋅𝐧⁺+𝐯⁻⋅𝐧⁻,{𝝉}=12𝝉⁺+𝝉⁻,[𝝉]=𝝉⁺⋅𝐧⁺+𝝉⁻⋅𝐧⁻. (6) If 𝐱 is on a face 𝑒 lying on the boundary 𝜕Ω, the above terms are defined by {𝐯}[𝐯]{𝝉}[𝝉]=𝐯,=𝐯⋅𝐧,=𝝉,=𝝉⋅𝐧, (7) where 𝐧 is the unit outward normal vector on 𝜕Ω. We define a matrix valued jump ⟦⋅⟧ of a vector 𝐯 as follows: ⟦𝐯⟧=12𝐯⁺⊗𝐧⁺+𝐧⁺⊗𝐯⁺+𝐯⁻⊗𝐧⁻+𝐧⁻⊗𝐯⁻,if𝐱∈𝑒∈ℰ^𝑖_ℎ,⟦𝐯⟧=12(𝐯⊗𝐧+𝐧⊗𝐯),if𝐱∈𝑒∈ℰ_ℎ∩𝜕Ω, (8) where 𝐯⊗𝐧 denote the matrix whose (𝑖,𝑗)th component is 𝑣_𝑖𝑛_𝑗 for two vectors 𝐯 and 𝐧.

Based on the triangulation 𝒯_ℎ, let 𝚺∶=𝝉∈𝐋²(Ω,𝕊);𝝉|_𝐾∈𝐇(div,𝐾,𝕊)∀𝐾∈𝒯_ℎ,𝐕∶=𝐯∈𝐋²Ω,ℝ^𝑑;𝐯|_𝐾∈𝐇¹𝐾,ℝ^𝑑∀𝐾∈𝒯_ℎ. (9) The corresponding finite element spaces are given by 𝚺_ℎ∶=𝝉∈𝐋²(Ω,𝕊);𝝉|_𝐾∈𝓢_𝐾∀𝐾∈𝒯_ℎ,𝐕_ℎ∶=𝐯∈𝐋²Ω,ℝ^𝑑;𝐯|_𝐾∈𝐏_𝑘𝐾,ℝ^𝑑∀𝐾∈𝒯_ℎ, (10) where, for each 𝐾∈𝒯_ℎ, 𝓢_𝐾 is a finite-dimensional symmetric tensor space of polynomials in 𝐾 containing 𝐏_𝑘+1(𝐾,𝕊), with 𝑘≥1. Let 𝐖_ℎ∶=𝝉∈𝐋²(Ω,𝕊);𝝉|_𝐾∈𝐏_𝑘+1(𝐾,𝕊)∀𝐾∈𝒯_ℎ,and_𝑒[𝝉]⋅𝐯𝑑𝑠=0∀𝐯∈𝐏_𝑘𝑒,ℝ^𝑑foreach𝑒∈ℰ^𝑖_ℎ, (11) which was used to construct a family of symmetric nonconforming mixed finite elements for linear elasticity in [21]. It is obvious that 𝐖_ℎ⊂𝚺_ℎ.

With former notations, we propose a reduced local discontinuous Galerkin (RLDG) method for nearly incompressible linear elasticity problem (1). Find (𝝈_ℎ,𝐮_ℎ)∈𝚺_ℎ×𝐕_ℎ such that 𝑎𝝈_ℎ,𝝉+𝑏𝝉,𝐮_ℎ=0∀𝝉∈𝚺_ℎ,𝑏𝝈_ℎ,𝐯=−_Ω𝐟⋅𝐯𝑑𝑥∀𝐯∈𝐕_ℎ, (12)(13) where 𝑎(𝝈,𝝉)∶=_Ω𝑏()𝒜𝝈∶𝝉𝑑𝑥,𝝉,𝐯∶=_{𝐾∈𝒯ℎ𝐾}(𝛁⋅𝝉)⋅𝐯𝑑𝑥−_ℰ𝑖ℎ[𝝉]⋅{𝐯}𝑑𝑠. (14) For comparison, let us recall the LDG method proposed in [20] in the following. Find (𝝈_ℎ,𝐮_ℎ)∈𝚺_ℎ×𝐕_ℎ such that 𝑎𝝈_ℎ,𝝉+𝑏𝝉,𝐮_ℎ=0∀𝝉∈𝚺_ℎ,𝑏𝝈_ℎ,𝐯−𝑐𝐮_ℎ,𝐯=−_Ω𝐟⋅𝐯𝑑𝑥∀𝐯∈𝐕_ℎ, (15) where 𝑐(𝐮,𝐯)=_ℰℎ𝐶₁₁⟦𝐮⟧∶⟦𝐯⟧𝑑𝑠 (16) with 𝐶₁₁>0. Formally, RLDG method (12)-(13) can be viewed as a special case of LDG method (15) with 𝐶₁₁=0. These two methods are actually different. LDG method requires 𝐶₁₁ to be positive to ensure the stability, whereas RLDG method drops this requirement by choosing large enough finite element space 𝚺_ℎ to match 𝐕_ℎ in the sense that 𝚺_ℎ⊃𝐖_ℎ.

On the other hand, according to the definition of 𝐖_ℎ, the following holds: 𝑏(𝝉,𝐯)=_{𝐾∈𝒯ℎ𝐾}(𝛁⋅𝝉)⋅𝐯𝑑𝑥∀𝝉∈𝐖_ℎ,𝐯∈𝐕_ℎ. (17) With (17), RLDG method (12)-(13) will be reduced to the symmetric nonconforming mixed finite element method designed in [21] by taking 𝚺_ℎ=𝐖_ℎ. Thus RLDG method can be also regarded as the discretization of the symmetric nonconforming mixed finite element method when 𝓢_𝐾=𝐏_𝑘+1(𝐾,𝕊), for there is not any continuity of finite element space 𝚺_ℎ across interior faces. The merit of RLDG method over the symmetric nonconforming mixed finite element method is that we do not need to solve a painful saddle problem in RLDG method. By the way, the symmetric nonconforming mixed finite element method is solved by the hybridized technique in [21].

Theorem 1. Reduced local discontinuous Galerkin method (12)-(13) is uniquely solvable.

Proof. Since 𝚺_ℎ and 𝐕_ℎ are both finite-dimensional, it is enough to show that if 𝐟=𝟎, then 𝝈_ℎ=𝟎, 𝐮_ℎ=𝟎. Taking 𝝉=𝝈_ℎ in (12) and 𝐯=𝐮_ℎ in (13) and subtracting (13) from (12), it follows that 𝑎(𝝈_ℎ,𝝈_ℎ)=0. Thus 𝝈_ℎ=0 from the definition of 𝑎(⋅,⋅). Then (12) is reduced to 𝑏𝝉,𝐮_ℎ=0∀𝝉∈𝚺_ℎ. (18) Since 𝐖_ℎ⊂𝚺_ℎ and (17), thus _{𝐾∈𝒯ℎ𝐾}(𝛁⋅𝝉)⋅𝐮_ℎ𝑑𝑥=0∀𝝉∈𝐖_ℎ. (19) Then we can get 𝐮_ℎ=𝟎 by adopting the same argument used in the proof of Corollary 2.4 in [21], which ends the proof.

Next, let us consider the implementation of RLDG method (12)-(13). Let {𝝋_𝑖}^𝑁1_𝑖=1 and {𝝓_𝑗}^𝑁2_𝑗=1 be the bases of 𝚺_ℎ and 𝐕_ℎ, respectively. Write 𝝈_ℎ=^𝑁1_𝑖=1𝜎_𝑖𝝋_𝑖,𝐮_ℎ=^𝑁2_𝑗=1𝑢_𝑗𝝓_𝑗. (20) And denote 𝐴_𝑖𝑗∶=𝑎𝝋_𝑖,𝝋_𝑗,𝐵_𝑖𝑗∶=𝑏𝝋_𝑗,𝝓_𝑖,𝐹_𝑖∶=_Ω𝑓𝝓_𝑖𝑑𝑥. (21) Then RLDG method (12)-(13) has the matrix representation 𝐴𝐵^𝑇𝑌𝑋=0𝐵0−𝐹 (22) with 𝐴=(𝐴_𝑖𝑗)_𝑁1×𝑁1, 𝐵=(𝐵_𝑖𝑗)_𝑁2×𝑁1, 𝐹=(𝐹_𝑖)_𝑁2×1, 𝑋=(𝑢_𝑖)_𝑁2×1, and 𝑌=(𝜎_𝑖)_𝑁1×1. Since 𝚺_ℎ is a piecewise polynomial space without any continuity across faces in ℰ^𝑖_ℎ, we can rearrange the basis {𝝋_𝑖}^𝑁1_𝑖=1 such that 𝐴 is a block diagonal matrix associated with triangulation 𝒯_ℎ. Thus it is trivial to compute the inverse of matrix 𝐴. Then the last symmetric indefinite linear system can be reduced to 𝐵𝐴⁻¹𝐵^𝑇𝑋=𝐹. (23) The new linear system is much easier to solve compared to original system, for new linear system has only unknown 𝑋 and 𝐵𝐴⁻¹𝐵^𝑇 is symmetric and positive definite. After obtaining 𝑋 by solving the new linear system, we can get 𝑌 from 𝑌=−𝐴⁻¹𝐵^𝑇𝑋. (24) Therefore, RLDG method (12)-(13) can be implemented efficiently.

In the end of this section, we will intend to drive the primal formulation of RLDG method (12)-(13) using the argument taken in [20]. For this, assume that 𝒜𝚺_ℎ⊂𝚺_ℎ (25) in this paragraph. Define a global lifting operator 𝐫∶𝐋²(ℰ_ℎ,𝕊)→𝚺_ℎ by (cf. [20]) _Ω𝐫(𝝓)∶𝝉𝑑𝑥=−_ℰℎ𝝓∶{𝝉}𝑑𝑠∀𝝉∈𝚺_ℎ,𝝓∈𝐋²ℰ_ℎ,𝕊. (26) It holds the following inequality for lifting operator (cf. [20]): ‖𝐫(𝝓)‖²₀≲_{𝑒∈ℰℎ}ℎ_𝑒⁻¹‖𝝓‖²_0,𝑒∀𝝓∈𝐋²ℰ_ℎ,𝕊. (27) From (12), by integration by parts and definition of lifting operator 𝐫 we have _Ω𝒜𝝈_ℎ∶𝝉𝑑𝑥=_{𝐾∈𝒯ℎ𝐾}𝜺𝐮_ℎ∶𝝉𝑑𝑥−_ℰℎ{𝝉}∶𝐮_ℎ𝑑𝑠=_{𝐾∈𝒯ℎ𝐾}𝜺𝐮_ℎ∶𝝉𝑑𝑥+_Ω𝐫𝐮_ℎ∶𝝉𝑑𝑥∀𝝉∈𝚺_ℎ. (28) Thus we have from (25) 𝝈_ℎ=𝒜⁻¹𝜺_ℎ𝐮_ℎ+𝐫𝐮_ℎ. (29) Substituting (29) into (13) and using integration by parts and definition of lifting operator 𝐫, we gain the primal formulation of RLDG method (12)-(13) as follows. Find 𝐮_ℎ∈𝐕_ℎ such that _Ω𝒜⁻¹𝜺_ℎ𝐮_ℎ+𝐫𝐮_ℎ∶𝜺_ℎ(𝐯)𝐯=+𝐫(⟦⟧)𝑑𝑥_Ω𝐟⋅𝐯𝑑𝑥∀𝐯∈𝐕_ℎ. (30)

3. A Priori Error Analysis

In this section, we provide an a priori error analysis for RLDG method (12)-(13). For a function 𝑣∈𝐿²(Ω) with 𝑣|_𝐾∈𝐻^𝑚(𝐾) for all 𝐾∈𝒯_ℎ, let ‖𝑣‖_𝑚,ℎ and |𝑣|_𝑚,ℎ be the usual broken 𝐻^𝑚-type norm and seminorm of 𝑣: ‖𝑣‖_𝑚,ℎ∶=_{𝐾∈𝒯ℎ}‖𝑣‖²_𝑚,𝐾^1/2,|𝑣|_𝑚,ℎ∶=_{𝐾∈𝒯ℎ}|𝑣|²_𝑚,𝐾^1/2. (31) If 𝑣 is a vector-valued or tensor-valued function, the corresponding ‖⋅‖_𝑚,ℎ and |⋅|_𝑚,ℎ are defined in the similar manners. For a vector or tensor 𝐯, its length |𝐯| is (𝐯⋅𝐯)^1/2 or (𝐯∶𝐯)^1/2. Here the symbol ∶ denotes the double dot product operation of tensors. Throughout this paper, we use the notation “≲⋯” to mean that “≤𝐶⋯,” where 𝐶 is a generic positive constant independent of local element sizes and 𝜆, which may take different values at different appearances.

According to the definition of 𝒜, we have, for any second-order tensor 𝝉, 𝒜⁻¹𝝉=2𝜇𝝉+𝜆(tr𝝉)𝜹. (32) Given a bounded domain 𝐺⊂ℝ^𝑑, we define two energy-type norms for 𝐋²(𝐺,𝕊) in the following. For any 𝝉∈𝐋²(𝐺,𝕊), let ‖𝝉‖²_𝒜,𝐺∶=_𝐺𝒜𝝉∶𝝉𝑑𝑥,‖𝝉‖²_{𝒜−1,𝐺}∶=_𝐺𝒜⁻¹𝝉∶𝝉𝑑𝑥. (33) By the definition of 𝒜, it is obvious that ‖𝝉‖_𝒜,𝐺≤1‖𝝉‖2𝜇^0,𝐺∀𝝉∈𝐋²(𝐺,𝕊). (34) Then we denote by (⋅,⋅)_𝒜,𝐺 the Hilbert inner product corresponding to energy norm ‖⋅‖_𝒜,𝐺. If 𝐺 is Ω, we abbreviate ‖𝝉‖_𝒜,Ω, ‖𝝉‖_𝒜⁻¹_,Ω, and (⋅,⋅)_𝒜,Ω by ‖𝝉‖_𝒜, ‖𝝉‖_𝒜⁻¹, and (⋅,⋅)_𝒜, respectively. For any 𝝉₁,𝝉₂∈𝐋²(𝐺,𝕊), we obtain from Cauchy-Schwarz inequality _Ω𝝉₁∶𝝉₂𝑑𝑥=(𝝉₁,𝒜⁻¹𝝉₂)_𝒜≤𝝉_1𝒜𝒜⁻¹𝝉_2𝒜=𝝉_1𝒜𝝉_2𝒜−1. (35) For later uses, we introduce two interpolation operators. Let 𝚷_ℎ be the linear operator from 𝐇(div,Ω,𝕊)∩𝐋^𝑝(Ω,𝕊) for any 𝑝>2 to 𝐖_ℎ defined by (3.12) in [21], and let 𝐐_ℎ be the 𝐿² projection operator from 𝐕 onto the finite element space 𝐕_ℎ. From (3.5) in [21], it holds for all 𝝉∈𝐇(div,Ω,𝕊)∩𝐋^𝑝(Ω,𝕊)𝛁_ℎ⋅𝚷_ℎ𝝉=𝐐_ℎ(𝛁⋅𝝉). (36)

Lemma 2. For all 𝝉∈𝐇(div,Ω,𝕊)∩𝐋^𝑝(Ω,𝕊) and 𝐯∈𝐕_ℎ, it follows 𝑏𝝉−𝚷_ℎ𝝉,𝐯=0. (37)

Proof. Since [𝝉]|_𝑒=𝟎 for any 𝑒∈ℰ^𝑖_ℎ and 𝚷_ℎ𝝉∈𝐖_ℎ, we have by (17) 𝑏𝝉−𝚷_ℎ𝝉,𝐯=_Ω(𝛁⋅𝝉)⋅𝐯𝑑𝑥−_Ω𝛁_ℎ⋅𝚷_ℎ𝝉⋅𝐯𝑑𝑥. (38) According to (36) and the definition of 𝐐_ℎ, 𝑏𝝉−𝚷_ℎ𝝉,𝐯=_Ω(𝛁⋅𝝉)⋅𝐯𝑑𝑥−_Ω𝐐_ℎ(𝛁⋅𝝉)⋅𝐯𝑑𝑥=0. (39) This ends the proof.

Next, let us study the Galerkin orthogonality of RLDG method (12)-(13). Using integration by parts, we can get from problem (1) that 𝑎()()𝑏()𝝈,𝝉+𝑏𝝉,𝐮=0∀𝝉∈𝚺,𝝈,𝐯=−_Ω𝐟⋅𝐯𝑑𝑥∀𝐯∈𝐕. (40) Subtracting (12)-(13) from (40), we obtain the following error equations: 𝑎𝝈−𝝈_ℎ,𝝉+𝑏𝝉,𝐮−𝐮_ℎ=0∀𝝉∈𝚺_ℎ,𝑏𝝈−𝝈_ℎ,𝐯=0∀𝐯∈𝐕_ℎ. (41)(42)

To derive the a priori error estimates of RLDG method (12)-(13), we require as in [17] some regularity result for the solution (𝝈,𝐮) of problem (1). We say the solution (𝝈,𝐮) of problem (1) is (𝑠+2)-regular, where 𝑠 is a nonnegative integer whenever 𝝈∈𝐇^𝑠+1(Ω,𝕊),𝐮∈𝐇^𝑠+2(Ω,ℝ^𝑑) and there holds (cf. [4, 5, 17, 24]) ‖𝝈‖_𝑠+1+‖𝐮‖_𝑠+2+𝜆‖𝛁⋅𝐮‖_𝑠+1≲‖𝐟‖_𝑠 (43) for 𝐟∈𝐇^𝑠(Ω,ℝ^𝑑). We also assume in this paper that tr𝚺_ℎ𝜹⊂𝚺_ℎ. (44) Moreover, introduce the usual BDM element space (cf. [25, 26]) 𝐕^BDM_ℎ∶=𝐯∈𝐕_ℎ;[𝐯]|_𝑒=0∀𝑒∈ℰ_ℎ. (45) Let 𝐈_ℎ∶(𝐻¹₀(Ω))^𝑑→𝐕^BDM_ℎ be the BDM nodal interpolation operator. For simplicity, we still write 𝚷_ℎ and 𝐈_ℎ for 𝚷_ℎ|_𝐾 and 𝐈_ℎ|_𝐾. Using Lemma 2, error equations (41)-(42) can be rewritten as 𝑎𝝈−𝝈_ℎ,𝝉+𝑏𝝉,𝐈_ℎ𝐮−𝐮_ℎ=𝑏𝝉,𝐈_ℎ𝐮−𝐮∀𝝉∈𝚺_ℎ,𝑏𝚷_ℎ𝝈−𝝈_ℎ,𝐯=0∀𝐯∈𝐕_ℎ. (46)(47) The error estimates for interpolation operators 𝚷_ℎ and 𝐈_ℎ are summarized in the following lemma (cf. [21, 26]).

Lemma 3. For all 𝐯∈𝐇^𝑚+1(𝐾,ℝ^𝑑), 𝝉∈𝐇^𝑚(Ω,𝕊) with 𝑚 a positive integer and all 𝐾∈𝒯_ℎ, one has the estimates 𝐯−𝐈_ℎ𝐯_0,𝐾+ℎ_𝐾𝐯−𝐈_ℎ𝐯_1,𝐾+ℎ_𝐾^1/2𝐯−𝐈_ℎ𝐯_0,𝜕𝐾≲ℎ_𝐾^{min{𝑚,𝑘}+1}‖𝐯‖_𝑚+1,𝐾,𝛁⋅𝐯−𝐈_ℎ𝐯_0,𝐾≲ℎ_𝐾^{min{𝑚,𝑘}}‖𝛁⋅𝐯‖_𝑚,𝐾,𝝉−𝚷_ℎ𝝉_0,𝐾+ℎ_𝐾^1/2𝝉−𝚷_ℎ𝝉_0,𝜕𝐾+ℎ_𝐾𝝉−𝚷_ℎ𝝉_1,𝐾≲ℎ_𝐾^{min{𝑚,𝑘+2}}‖𝝉‖_{𝑚,𝜔𝐾}, (48) where 𝜔_𝐾∶=⋃{𝐾^′∈𝒯_ℎ;𝐾^′∩𝐾≠∅}.

Lemma 4. For all 𝝉∈𝚺_ℎ and 𝐯∈𝐇₀(div,Ω)∩𝐕, it follows 𝑏(𝝉,𝐯)≤‖𝝉‖_𝒜‖𝜺(𝐯)‖_𝒜⁻¹+2𝜇‖𝐫(⟦𝐯⟧)‖₀. (49)

Proof . Applying integration by parts, 𝑏(𝝉,𝐯)=−_Ω𝝉∶𝜺(𝐯)𝑑𝑥+_ℰℎ{𝝉}∶⟦𝐯⟧𝐯𝑑𝑠=−_Ω𝝉∶(𝜺(𝐯)+𝐫(⟦𝐯⟧))𝑑𝑥. (50) Then using (35) and triangle inequality, 𝑏(𝝉,𝐯)≤‖𝝉‖_𝒜‖𝜺(𝐯)‖_𝒜⁻¹+‖𝐫(⟦𝐯⟧)‖_𝒜⁻¹. (51) On the other hand, according to (44) and the definitions of trace operator tr and lifting operator 𝐫, ‖tr(𝐫(⟦𝐯⟧))‖²₀=_Ωtr(𝐫(⟦𝐯⟧))𝜹∶𝐫(⟦𝐯⟧)𝑑𝑥=−_ℰℎ{tr(𝐫(⟦𝐯⟧))𝜹}∶⟦𝐯⟧𝑑𝑠. (52) Since 𝐯∈𝐇₀(div,Ω) implies [𝐯]=0 on ℰ_ℎ, ‖tr(𝐫(⟦𝐯⟧))‖²₀=−_ℰℎtr({𝐫(⟦𝐯⟧)})𝜹∶⟦𝐯⟧𝑑𝑠=−_ℰℎtr({𝐫(⟦𝐯⟧)})[𝐯]𝑑𝑠=0, (53) which means tr(𝐫(⟦𝐯⟧))=0 in Ω. Thus by the definition of 𝒜⁻¹, ‖𝐫(⟦𝐯⟧)‖²_𝒜−1=2𝜇‖𝐫(⟦𝐯⟧)‖²₀+𝜆‖tr(𝐫(⟦𝐯⟧))𝜹‖²₀=2𝜇‖𝐫(⟦𝐯⟧)‖²₀. (54) Therefore, we can finish the proof by combining last inequality and (51).

Lemma 5. Assume that the solution (𝝈,𝐮) of problem (1) is (𝑠+2)-regular, where 𝑠 is a nonnegative integer. For all 𝝉∈𝚺_ℎ, it follows 𝑏𝝉,𝐈_ℎ𝐮−𝐮≲ℎ^{min{𝑘,𝑠+1}}‖𝝉‖_𝒜‖𝐟‖_𝑠. (55)

Proof. It follows from the definition of 𝒜⁻¹, (27), Lemma 3, and regularity (43) that 𝜺𝐮−𝐈_ℎ𝐮²_𝒜−1+2𝜇𝐫𝐮−𝐈_ℎ𝐮²₀=2𝜇𝜺𝐮−𝐈_ℎ𝐮²₀+𝐫𝐮−𝐈_ℎ𝐮²₀+𝜆𝛁⋅𝐮−𝐈_ℎ𝐮²₀≲2𝜇𝜺𝐮−𝐈_ℎ𝐮²₀+_{𝑒∈ℰℎ}ℎ_𝑒⁻¹𝐮−𝐈_ℎ𝐮²_0,𝑒+𝜆𝛁⋅𝐮−𝐈_ℎ𝐮²₀≲ℎ^{2min{𝑘,𝑠+1}}‖𝐮‖²_𝑠+2+𝜆‖𝛁⋅𝐮‖²_𝑠+1≲ℎ^{2min{𝑘,𝑠+1}}‖𝐟‖²_𝑠. (56) Noting the fact that 𝐈_ℎ𝐮−𝐮∈𝐇₀(div,Ω)∩𝐕, we have from Lemma 4𝑏𝝉,𝐈_ℎ𝐮−𝐮≤‖𝝉‖_𝒜𝜺𝐮−𝐈_ℎ𝐮_𝒜−1+2𝜇𝐫𝐮−𝐈_ℎ𝐮₀≲ℎ^{min{𝑘,𝑠+1}}‖𝝉‖_𝒜‖𝐟‖_𝑠, (57) which completes the proof.

With previous preparation, we are ready to give the a priori error estimates.

Theorem 6. Assume that the solution (𝝈,𝐮) of problem (1) is (𝑠+2)-regular, where 𝑠 is a nonnegative integer. Let (𝝈_ℎ,𝐮_ℎ) be the solution of RLDG method (12)-(13). Then 𝝈−𝝈_ℎ𝒜≲ℎ^{min{𝑘,𝑠+1}}‖𝐟‖_𝑠. (58)

Proof . Choosing 𝝉=𝚷_ℎ𝝈−𝝈_ℎ in (46) and 𝐯=𝐈_ℎ𝐮−𝐮_ℎ in (47) and subtracting (47) from (46), we get 𝑎𝝈−𝝈_ℎ,𝚷_ℎ𝝈−𝝈_ℎ=𝑏𝚷_ℎ𝝈−𝝈_ℎ,𝐈_ℎ𝐮−𝐮. (59) Then we have 𝚷_ℎ𝝈−𝝈_ℎ²_𝒜=𝑎𝚷_ℎ𝝈−𝝈_ℎ,𝚷_ℎ𝝈−𝝈_ℎ=𝑎𝚷_ℎ𝝈−𝝈,𝚷_ℎ𝝈−𝝈_ℎ+𝑎𝝈−𝝈_ℎ,𝚷_ℎ𝝈−𝝈_ℎ=𝑎𝚷_ℎ𝝈−𝝈,𝚷_ℎ𝝈−𝝈_ℎ+𝑏𝚷_ℎ𝝈−𝝈_ℎ,𝐈_ℎ𝐮−𝐮. (60) Applying Cauchy-Schwarz inequality and Lemma 5, it follows 𝚷_ℎ𝝈−𝝈_ℎ²_𝒜≲𝚷_ℎ𝝈−𝝈_𝒜𝚷_ℎ𝝈−𝝈_ℎ𝒜+ℎ^{min{𝑘,𝑠+1}}𝚷_ℎ𝝈−𝝈_ℎ𝒜‖𝐟‖_𝑠, (61) that is, 𝚷_ℎ𝝈−𝝈_ℎ𝒜≲𝝈−𝚷_ℎ𝝈_𝒜+ℎ^{min{𝑘,𝑠+1}}‖𝐟‖_𝑠. (62) By triangle inequality and (34), we get 𝝈−𝝈_ℎ𝒜≤𝝈−𝚷_ℎ𝝈_𝒜+𝚷_ℎ𝝈−𝝈_ℎ𝒜≲𝝈−𝚷_ℎ𝝈₀+ℎ^{min{𝑘,𝑠+1}}‖𝐟‖_𝑠. (63) The proof is finished from Lemma 3.

Under the assumption (25), we can get the following error estimate for 𝐮_ℎ in energy norm from Theorem 6 and (29).

Corollary 7. Assume that (25) is true and the solution (𝝈,𝐮) of problem (1) is (𝑠+2)-regular, where 𝑠 is a nonnegative integer. Let (𝝈_ℎ,𝐮_ℎ) be the solution of RLDG method (12)-(13). Then 𝜺_ℎ𝐮−𝐮_ℎ+𝐫𝐮−𝐮_ℎ𝒜−1=𝝈−𝝈_ℎ𝒜≲ℎ^{min{𝑘,𝑠+1}}‖𝐟‖_𝑠. (64)

Using the usual duality argument, we can additionally derive an optimal order error estimate between 𝐮_ℎ and 𝐮 in 𝐋²(Ω,ℝ^𝑑)-norm.

Theorem 8. Let Ω be a convex bounded polygonal domain. Assume that the solution (𝝈,𝐮) of problem (1) is (𝑠+2)-regular, where 𝑠 is a nonnegative integer. Let (𝝈_ℎ,𝐮_ℎ) be the solution of RLDG method (12)-(13). Then 𝐮−𝐮_ℎ0≲ℎ^{min{𝑘+1,𝑠+2}}‖𝐟‖_𝑠. (65)

Proof . Let (𝝈^∼,𝐮^∼) be the solution of the auxiliary problem it follows: 𝒜𝝈^∼−𝜺𝐮^∼=𝟎inΩ,−𝛁⋅𝝈^∼=𝐮−𝐮_ℎinΩ,_∼𝐮=𝟎on𝜕Ω. (66) For Ω is a convex bounded polygon or polyhedron, the solution (𝝈^∼,𝐮^∼) of problem (66) has 2-regularity estimate (cf. [24]): 𝝈^∼₁+𝐮^∼₂+𝜆𝛁⋅𝐮^∼₁≲𝐮−𝐮_ℎ0. (67) Multiplying 𝐮−𝐮_ℎ on both sides of the second equation of (66) and integrating on Ω, we have 𝐮−𝐮_ℎ²₀=−𝑏𝝈^∼,𝐮−𝐮_ℎ=−𝑏𝝈^∼−𝚷𝝈_ℎ^∼,𝐮−𝐮_ℎ−𝑏𝚷𝝈_ℎ^∼,𝐮−𝐮_ℎ. (68) By Lemma 2, we get 𝑏𝝈^∼−𝚷𝝈_ℎ^∼,𝐮−𝐮_ℎ=𝑏𝝈^∼−𝚷𝝈_ℎ^∼,𝐮−𝐈_ℎ𝐮. (69) Choosing 𝝉=𝚷𝝈_ℎ^∼ in error equation (41), it holds −𝑏𝚷𝝈_ℎ^∼,𝐮−𝐮_ℎ=𝑎𝝈−𝝈_ℎ,𝚷𝝈_ℎ^∼=𝑎𝝈−𝝈_ℎ,𝚷𝝈_ℎ^∼−𝝈^∼+𝑎𝝈−𝝈_ℎ,𝝈^∼. (70) Then it follows from the last three equalities that 𝐮−𝐮_ℎ²₀=−𝑏𝝈^∼−𝚷𝝈_ℎ^∼,𝐮−𝐈_ℎ𝐮+𝑎𝝈−𝝈_ℎ,𝚷𝝈_ℎ^∼−𝝈^∼+𝑎𝝈−𝝈_ℎ,𝝈^∼. (71) For the first term of the right hand side in (71), applying Cauchy-Schwarz inequality, Lemma 3, and regularity (67), it holds −𝑏𝝈^∼−𝚷𝝈_ℎ^∼,𝐮−𝐈_ℎ𝐮=−_{𝐾∈𝒯ℎ𝐾}𝛁⋅𝝈^∼−𝚷𝝈_ℎ^∼⋅𝐮−𝐈_ℎ𝐮+𝑑𝑥_{ℰ𝑖ℎ∼}𝝈−𝚷𝝈_ℎ^∼⋅𝐮−𝐈_ℎ𝐮𝑑𝑠≲ℎ𝝈^{min{𝑘+1,𝑠+2}∼}₁‖𝐮‖_𝑠+2≲ℎ^{min{𝑘+1,𝑠+2}}‖𝐟‖_𝑠𝐮−𝐮_ℎ0. (72) For the second term of the right hand side in (71), we get from Cauchy-Schwarz inequality, (34), Theorem 6, Lemma 3, and regularity (67) 𝑎𝝈−𝝈_ℎ,𝚷𝝈_ℎ^∼−𝝈^∼≤𝝈−𝝈_ℎ𝒜𝚷𝝈_ℎ^∼−𝝈^∼_𝒜≤𝝈−𝝈_ℎ𝒜𝚷𝝈_ℎ^∼−𝝈^∼₀≲ℎ^{min{𝑘+1,𝑠+2}}‖𝐟‖𝝈_𝑠^∼₁≲ℎ^{min{𝑘+1,𝑠+2}}‖𝐟‖_𝑠𝐮−𝐮_ℎ0. (73) For the third term of the right hand side in (71), by the first equation of (66), integration by parts, and using error equation (42), we get 𝑎𝝈−𝝈_ℎ,𝝈^∼=_Ω𝝈−𝝈_ℎ∶𝜺𝐮^∼𝑑𝑥=−𝑏𝝈−𝝈_ℎ,𝐮^∼=−𝑏𝝈−𝝈_ℎ,𝐮^∼−𝐈𝐮_ℎ^∼=−𝑏𝝈−𝚷_ℎ𝝈,𝐮^∼−𝐈𝐮_ℎ^∼−𝑏𝚷_ℎ𝝈−𝝈_ℎ,𝐮^∼−𝐈𝐮_ℎ^∼. (74) According to Lemma 5 with regularity (67), triangle inequality, (34), Theorem 6, and Lemma 3, we have −𝑏𝚷_ℎ𝝈−𝝈_ℎ,𝐮^∼−𝐈𝐮_ℎ^∼≲ℎ𝚷_ℎ𝝈−𝝈_ℎ𝒜𝐮−𝐮_ℎ0≲ℎ𝝈_ℎ−𝝈_𝒜+𝝈−𝚷_ℎ𝝈₀𝐮−𝐮_ℎ0≲ℎ^{min{𝑘+1,𝑠+2}}‖𝐟‖_𝑠+‖𝝈‖_𝑠+1𝐮−𝐮_ℎ0≲ℎ^{min{𝑘+1,𝑠+2}}‖𝐟‖_𝑠𝐮−𝐮_ℎ0. (75) On the other hand, using Cauchy-Schwarz inequality, Lemma 3, and regularity (67), −𝑏𝝈−𝚷_ℎ𝝈,𝐮^∼−𝐈𝐮_ℎ^∼=−_{𝐾∈𝒯ℎ𝐾}𝛁⋅𝝈−𝚷_ℎ𝝈⋅𝐮^∼−𝐈𝐮_ℎ^∼+𝑑𝑥_ℰ𝑖ℎ𝝈−𝚷_ℎ𝝈⋅𝐮^∼−𝐈𝐮_ℎ^∼𝑑𝑠≲ℎ^{min{𝑘+1,𝑠+2}}‖𝝈‖𝐮_𝑠+1^∼₂≲ℎ^{min{𝑘+1,𝑠+2}}‖𝐟‖_𝑠𝐮−𝐮_ℎ0. (76) Thus we know from the last two inequalities that 𝑎𝝈−𝝈_ℎ,𝝈^∼≲ℎ^{min{𝑘+1,𝑠+2}}‖𝐟‖_𝑠𝐮−𝐮_ℎ0. (77) Combining (71)–(77) gives 𝐮−𝐮_ℎ²₀≲ℎ^{min{𝑘+1,𝑠+2}}‖𝐟‖_𝑠𝐮−𝐮_ℎ0. (78) Finally, we can conclude the result by cancelling ‖𝐮−𝐮_ℎ‖₀ on both sides.

4. Postprocessing of RLDG Solution

In this section, we will focus on reconstructing an approximate displacement in the Raviart-Thomas element space via the ideas in [17]. Assume that triangulation 𝒯_ℎ is conforming in this section; that is, there is no hanging nodes in 𝒯_ℎ. The Raviart-Thomas element space is defined by (cf. [26, 27]) 𝐑𝐓_ℎ∶=𝐯∈𝐇₀();𝐯|div,Ω_𝐾∈𝐏_𝑘𝐾,ℝ^𝑑+𝐱𝑃_𝑘(𝐾)∀𝐾∈𝒯_ℎ. (79) With this space, the postprocessing can be described as follows. Find 𝐮_ℎ∈𝐑𝐓_ℎ such that for any 𝐾∈𝒯_ℎ, _𝐾𝐮_ℎ−𝐮_ℎ⋅𝐯𝑑𝑥=0∀𝐯∈𝐏_𝑘−1𝐾,ℝ^𝑑,_𝑒𝐮_ℎ−𝐮_ℎ⋅𝐧_𝐾𝑣𝑑𝑠=0∀𝐯∈𝐏_𝑘𝑒,ℝ^𝑑,∀face𝑒⊂𝜕𝐾⧵𝜕Ω, (80) where 𝐧_𝐾 is the unit outward normal on 𝜕𝐾, pointing at outside of 𝐾. Since (80) are defined locally, the function 𝐮_ℎ can be obtained very conveniently and in parallel. To show its approximation property, we first establish a relation between the divergence of the post-processed displacement 𝐮_ℎ and the RLDG displacement 𝐮_ℎ. For this, let 𝑉_ℎ∶=𝑣∈𝐿²(Ω);𝑣|_𝐾∈𝑃_𝑘(𝐾)∀𝐾∈𝒯_ℎ. (81) And 𝑄_ℎ denote the 𝐿²-projection operator from 𝐿²(Ω) onto the finite element space 𝑉_ℎ. We first recall two results which have been proved in [28] 𝛁⋅𝐮_ℎ=𝛁_ℎ⋅𝐮_ℎ+𝑄_ℎtr𝐫𝐮_ℎ,𝐮_ℎ−𝐮_ℎ0,𝐾≲ℎ_𝐾^1/2tr𝐮_{ℎ0,𝜕𝐾}=ℎ_𝐾^1/2[𝐮_ℎ]_0,𝜕𝐾∀𝐾∈𝒯_ℎ. (82)(83)

Now we can give the error estimate for 𝐮_ℎ as follows.

Theorem 9. Let (𝝈_ℎ,𝐮_ℎ) be the solution of RLDG method (12)-(13) and let 𝐮_ℎ be the solution of (80). Under the same assumption as in Corollary 7, one has √𝜆𝛁⋅(𝐮−𝐮_ℎ)₀≲ℎ^{min{𝑘,𝑠+1}}‖𝐟‖_𝑠. (84) Moreover, if Ω is a convex bounded polygon or polyhedron, 𝐮−𝐮_ℎ0≲ℎ^{min{𝑘+1,𝑠+2}}‖𝐟‖_𝑠. (85)

Proof . From (82), and the definition of 𝑄_ℎ, it follows that √𝜆𝛁⋅𝐮−𝐮_ℎ0=√𝜆𝛁_ℎ⋅𝐮−𝐮_ℎ−𝑄_ℎtr𝐫𝐮_ℎ0=√𝜆𝛁_ℎ⋅𝐮−𝐮_ℎ+𝑄_ℎtr𝐫𝐮−𝐮_ℎ0=√𝜆𝛁⋅𝐮−𝑄_ℎ(𝛁⋅𝐮)+𝑄_ℎ𝛁_ℎ⋅𝐮−𝐮_ℎ+tr𝐫𝐮−𝐮_ℎ0. (86) Then by triangle inequality, and the boundedness of 𝑄_ℎ (cf. [22, 23]), it follows that √𝜆𝛁⋅𝐮−𝐮_ℎ0=√𝜆𝛁⋅𝐮−𝑄_ℎ(𝛁⋅𝐮)₀+√𝜆𝑄_ℎ𝛁_ℎ⋅𝐮−𝐮_ℎ+tr𝐫𝐮−𝐮_ℎ0=√𝜆𝛁⋅𝐮−𝑄_ℎ(𝛁⋅𝐮)₀+√𝜆𝛁_ℎ⋅𝐮−𝐮_ℎ+tr𝐫𝐮−𝐮_ℎ0≤√𝜆𝛁⋅𝐮−𝑄_ℎ(𝛁⋅𝐮)₀+𝜺_ℎ(𝐮−𝐮_ℎ)+𝐫𝐮−𝐮_ℎ𝒜−1. (87) On the other hand, it is well known that 𝐿²-projection operator 𝑄_ℎ has the following estimate (cf. [22, 23]): 𝛁⋅𝐮−𝑄_ℎ(𝛁⋅𝐮)₀≲ℎ^{min{𝑘+1,𝑠+1}}‖𝛁⋅𝐮‖_𝑠+1. (88) Hence, we obtain from the last two inequalities and Corollary 7 that √𝜆𝛁⋅𝐮−𝐮_ℎ0≲√𝜆ℎ^{min{𝑘+1,𝑠+1}}‖𝛁⋅𝐮‖_𝑠+1+ℎ^{min{𝑘,𝑠+1}}‖𝐟‖_𝑠≲ℎ^{min{𝑘+1,𝑠+1}}‖𝐟‖_𝑠+ℎ^{min{𝑘,𝑠+1}}‖𝐟‖_𝑠≲ℎ^{min{𝑘,𝑠+1}}‖𝐟‖_𝑠. (89) This completes the proof of (84).
As for (85), noting the fact that [𝐈_ℎ𝐮]=0 for any 𝑒∈ℰ_ℎ, we get from (83) and inverse inequality 𝐮_ℎ−𝐮_ℎ0,𝐾≲ℎ_𝐾^1/2𝐈_ℎ𝐮−𝐮_{ℎ0,𝜕𝐾}≲𝐈_ℎ𝐮−𝐮_{ℎ0,𝜔𝐾}. (90) Summing over all elements in 𝒯_ℎ, it gives 𝐮_ℎ−𝐮_ℎ0≲𝐈_ℎ𝐮−𝐮_ℎ0. (91) Then using triangle inequality, 𝐮−𝐮_ℎ0≤𝐮−𝐮_ℎ0+𝐮_ℎ−𝐮_ℎ0≲𝐮−𝐮_ℎ0+𝐈_ℎ𝐮−𝐮_ℎ0≤𝐮−𝐮_ℎ0+𝐮−𝐈_ℎ𝐮₀. (92) Therefore, we can obtain (85) from Theorem 8 and Lemma 3.

5. Numerical Results

In this section, some numerical results are shown to validate the theoretical convergence rates of RLDG method. Let Ω=(−1,1)×(−1,1), 𝜇=1, and 𝐟𝑥₁,𝑥₂=−8𝑥₁+𝑥₂3𝑥₁𝑥₂−2𝑥²₁+𝑥²₂+5(𝑥₁𝑥₂−1)²−2𝑥²₁𝑥²₂−8𝑥₁−𝑥₂3𝑥₁𝑥₂+2𝑥²₁+𝑥²₂−5(𝑥₁𝑥₂+1)²+2𝑥²₁𝑥²₂. (93) It is not difficult to check that the exact solution of (1) for this load is 𝐮𝑥₁,𝑥₂=−4𝑥₂1−𝑥²₂(1−𝑥²₁)²4𝑥₁1−𝑥²₁(1−𝑥²₂)²+12+𝜆−4𝑥₁1−𝑥²₁(1−𝑥²₂)²−4𝑥₂1−𝑥²₂(1−𝑥²₁)²,𝝈𝑥₁,𝑥₂=2𝜺(𝐮)+𝜆(𝛁⋅𝐮)𝜹. (94)

We make a uniform triangulation 𝒯_ℎ on Ω. Set 𝓢_𝐾=𝐏_𝑘+1(𝐾,𝕊) for all 𝐾∈𝒯_ℎ and 𝑘=1.

First, let us consider the convergence of RLDG method for fixed Lamé coefficient 𝜆. Errors between (𝝈_ℎ,𝐮_ℎ) and (𝝈,𝐮) for 𝜆=1 are listed in Table 1, from which we can see that the convergence rates of ‖𝐮−𝐮_ℎ‖₀ and ‖𝝈−𝝈_ℎ‖_𝒜 are 𝑂(ℎ²) and 𝑂(ℎ), respectively. These numerical results agree with theoretical results in Theorems 6 and 8. We also present the numerical convergence rates of postprocessing solution 𝐮_ℎ in Table 2 for 𝜆=1. As Theorem 9 indicates that the convergence rates of ‖𝐮−𝐮_ℎ‖₀ and √𝜆‖𝛁⋅(𝐮−𝐮_ℎ)‖₀ are 𝑂(ℎ²) and 𝑂(ℎ), respectively. To illustrate the robustness of RLDG method with respect to 𝜆, we give numerical results for 𝜆=10⁹ in Table 3, which is very close to the incompressible case. It can be observed from Table 3 that the convergence rates of ‖𝐮−𝐮_ℎ‖₀ and ‖𝝈−𝝈_ℎ‖_𝒜 being again 𝑂(ℎ²) and 𝑂(ℎ) do not deteriorate even for such large 𝜆; hence our RLDG method is locking-free.

Table 1 

Errors between (𝝈_ℎ, 𝐮_ℎ) and (𝝈, 𝐮) for different ℎ when 𝜆=1.

Table 2 

Errors between 𝐮_ℎ and 𝐮 for different ℎ when 𝜆=1.

Table 3 

Errors between (𝝈_ℎ, 𝐮_ℎ) and (𝝈, 𝐮) for different ℎ when 𝜆=10⁹.

Then, we fix ℎ=2⁻⁴ and let 𝜆 vary from 1 to 10⁹. The corresponding errors ‖𝐮−𝐮_ℎ‖₀, ‖𝝈−𝝈_ℎ‖_𝒜, ‖𝐮−𝐮_ℎ‖₀, and √𝜆‖𝛁⋅(𝐮−𝐮_ℎ)‖₀ are shown in Table 4. From Table 4, we note that the errors ‖𝐮−𝐮_ℎ‖₀, ‖𝝈−𝝈_ℎ‖_𝒜, and ‖𝐮−𝐮_ℎ‖₀ are insensitive to the choice of 𝜆. And √𝜆‖𝛁⋅(𝐮−𝐮_ℎ)‖₀ becomes smaller as 𝜆 grows larger; thus it is of course bounded. In fact, the quantity 𝜆‖𝛁⋅(𝐮−𝐮_ℎ)‖₀ is insensitive to the choice of 𝜆 by an in-depth observation. In a word, all the errors in Table 4 are uniform with respect to 𝜆, which coincide with the theoretical results.

Table 4 

Numerical results for different 𝜆 when ℎ=2⁻⁴.

In the end, let us compare RLDG method and LDG method. Set 𝜆=1 again. Errors ‖𝐮−𝐮_ℎ‖₀ and ‖𝝈−𝝈_ℎ‖_𝒜 of RLDG and LDG methods with various choices of 𝐶₁₁ are shown in Tables 5 and 6, from which we can see that errors ‖𝐮−𝐮_ℎ‖₀ and ‖𝝈−𝝈_ℎ‖_𝒜 increase as 𝐶₁₁ becomes larger. To be more precise, errors ‖𝐮−𝐮_ℎ‖₀ and ‖𝝈−𝝈_ℎ‖_𝒜 grow slowly with respect to 𝐶₁₁ when 𝐶₁₁|_𝑒≲ℎ_𝑒⁻¹ for all 𝑒∈ℰ_ℎ; however, they become much larger when 𝐶₁₁ reaches 𝑂(ℎ⁻²) especially for finer mesh. We use PCG method with diagonal preconditioner to solve the resulting linear systems, and the stop criterion requires that the relative residual is less than 10⁻¹⁰. The numbers of diagonal-PCG iterations for RLDG and LDG methods are given in Table 7, which show similar variation as those in Tables 5 and 6. Therefore, we can conclude that RLDG method performs better than LDG method, especially for large 𝐶₁₁. The performance of LDG method is almost robust with respect to 𝐶₁₁ when 𝐶₁₁≲ℎ⁻¹ and deteriorates when 𝐶₁₁=𝑂(ℎ⁻²). In addition, if stress tensor and displacement are both approximated by piecewise linear polynomials, RLDG method will be ill-posed, whereas LDG method is still well-posed whose numerical results can be found in [20].

Table 5 

Error ‖𝐮−𝐮_ℎ‖₀ of RLDG and LDG methods with various choices of 𝐶₁₁.

Table 6 

Error ‖𝝈−𝝈_ℎ‖_𝒜 of RLDG and LDG methods with various choices of 𝐶₁₁.

Table 7 

The numbers of diagonal-PCG iterations for RLDG and LDG methods with various choices of 𝐶₁₁.

6. Conclusions

In this paper, we develop a locking-free RLDG method for nearly incompressible linear elasticity. RLDG method can be formally regarded as a special case of LDG method with 𝐶₁₁=0 proposed in [20] and can also be viewed as the localization of symmetric nonconforming mixed finite element method presented in [21]. Then we discuss the implementation of RLDG method and obtain robust error estimates in energy norm and 𝐿²-norm which are uniform with respect to 𝜆. Error estimates of the postprocessing 𝐻(div)-nonconforming displacement in 𝐻(div)-norm and 𝐿²-norm are achieved, which are also uniform with respect to 𝜆. According to numerical results, RLDG method performs a little better than LDG method for small 𝐶₁₁ and much better than LDG method for large 𝐶₁₁. The performance of LDG method is almost robust with respect to 𝐶₁₁ when 𝐶₁₁≲ℎ⁻¹ and deteriorates when 𝐶₁₁=𝑂(ℎ⁻²).

Acknowledgments

The author thanks the referee for valuable comments and suggestions. This work was partly supported by the NSFC (Grant nos. 11126226, 11171257) and Zhejiang Provincial Natural Science Foundation of China (Y6110240, LY12A01015).