A lowest-order weak Galerkin method for linear elasticity

Accepted Manuscript A lowest-order weak Galerkin method for linear elasticity Son-Young Yi

PII: DOI: Reference:

S0377-0427(18)30623-X https://doi.org/10.1016/j.cam.2018.10.016 CAM 11961

To appear in:

Journal of Computational and Applied Mathematics

Received date : 18 December 2017 Revised date : 17 September 2018 Please cite this article as: S.-Y. Yi, A lowest-order weak Galerkin method for linear elasticity, Journal of Computational and Applied Mathematics (2018), https://doi.org/10.1016/j.cam.2018.10.016 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

A lowest-order weak Galerkin method for linear elasticity Son-Young Yia,∗ a

Department of Mathematical Sciences, The University of Texas at El Paso, El Paso, TX 79968-0506, USA

Abstract The lowest-order weak Galerkin (WG) method is considered for linear elasticity based on the displacement formulation. The new method approximates the displacement using piecewise constant vector functions both inside and on the boundary of each mesh element and its bilinear form does not require a stabilization term for the existence and uniqueness of the solution. A-priori error estimates of optimal order in the discrete H 1 - and L2 -norms for the displacement are proved when the solution is smooth. The error estimates are independent of the Lam´e constant λ, thus the performance of the new method does not deteriorate as the elastic material becomes incompressible. Further, a simple post-processing technique to obtain a numerical approximation of the stress is presented. A careful error analysis reveals that the L2 -norm error in the stress is also optimal and independent of λ. Several numerical experiments confirm the locking-free property and optimal convergence rates of the new method. Keywords: linear elasticity, weak Galerkin method, finite element method, locking-free 1. Introduction In this paper, we present and analyze a robust and efficient weak Galerkin method for the equations of linear elasticity. It is well known that standard continuous Galerkin methods, such as continuous piecewise linear or bilinear elements, yield poor approximations to the displacement as the material ∗

Corresponding author. Email address: [email protected] (Son-Young Yi)

Preprint submitted to Journal of Computational and Applied MathematicsOctober 23, 2018

becomes incompressible or equivalently, the Lam´e constant λ → ∞. In this case, either a finite element method produces a significantly smaller displacement than it should or numerical solutions exhibit oscillations in the stress variable. This phenomenon is known as “Poisson locking.” In order to overcome locking, various approaches have been proposed, most of which are characterized as mixed finite element methods [1, 2, 3, 4, 5], nonconforming finite element methods [6, 7], or discontinuous Galerkin methods [8, 9]. Mixed finite element methods for linear elasticity, which approximate both the displacement and stress simultaneously, can be free of locking and offer an accurate approximation of the stress. However, it is well known that a mixed finite element method requires a careful choice of mixed finite element spaces that satisfy the stability conditions for saddle point problems [10]. Discontinuous Galerkin methods based on the primal formulation do not require the stability conditions, but they usually require lots of degrees of freedom, hence a solution of a large linear system [8, 9]. Very recently, weak Galerkin methods have been also proposed for linear elasticity [11, 12, 13]. The weak Galerkin method refers to a class of finite element methods that employ weakly defined differential operators to discretize partial differential equations. It was first introduced by Wang and Ye [14] for the second-order elliptic problem and later applied to various other partial differential equations, including biharmonic equations [15, 16, 17, 18], parabolic equations [19], Stokes equations [20], Brinkman equations [21], Maxwell’s equations [22], and Helmholtz equation [23]. In the papers by Zhang [11] and Chen and Xie [13], weak Galerkin mixed finite element methods were proposed for linear elasiticy, in which numerical approximations of the stress and the displacement are found simultaneously based on the Hellinger-Reissner principle and Hu-Washizu principle, respectively. Both methods also employ a bilinear form with a stabilization term, without which the existence and uniqueness of the solution cannot be proved. On the other hand, Wang et al. [12] presented a family of locking-free weak Galerkin methods based on the primal formulation, though their convergence analysis was performed for an equivalent weak Galerkin mixed finite element method. The method by Wang et al. [12] requires a stabilization term for the existence and uniqueness of the solution like the methods by Zhang [11] and Chen and Xie [13] do. The lowest-order method in the family approximates the displacement using a piecewise linear vector function in the interior of each element and using a trace of a rigid motion on the boundary of each element. 2

The main goal of this paper is to propose and analyze a locking-free and lowest-order weak Galerkin finite element method for the equations of linear elasticity based on the displacement formulation. The numerical solution of the displacement is approximated by piecewise constant vector functions in the interior and on the boundary of each element, whereas the discrete weak gradient and discrete weak divergence operators are approximated by a local Raviart-Thomas element and a piecewise constant function, respectively. Unlike the previously proposed methods by Zhang [11], Wang et al. [12], and Chen and Xie [13], the new lowest-order method does not require a stabilization term for the existence and uniqueness of the solution. We prove a-priori error estimates of optimal order in the discrete H 1 - and L2 -norms for the displacement when the solution is smooth. They are independent of the Lam´e constant λ, thus the performance of the new method does not deteriorate as the material becomes incompressible. Moreover, our L2 -norm error estimate does not use the usual duality argument, but instead it uses a direct relationship between the errors in the L2 - and the discrete H 1 -norms. Another aspect in which the current work distinguishes itself from the work by Wang et al. [12] is that we present a straightforward and computationally inexpensive post-processing technique to recover a numerical approximation of the stress from the approximation of the displacement. The technique is rather trivial, but our careful error analysis shows that the L2 -norm error estimate for the stress is of optimal order and independent of λ. It is worth mentioning that the method and error analysis are presented assuming the finite element mesh is either a simplicial or a rectangular mesh in two and three dimensions. However, owing to the recent work by Chen and Wang [24] on H(div)-conforming finite element spaces, we can extend our method to a general polygonal mesh in two dimensions or to a certain type of polyhedral mesh in three dimensions. The rest of this paper is organized as follows. In Section 2, we introduce the governing equations of linear elasticity, some notation, and preliminary results for linear elasticity. In Section 3, we review the definitions of weak gradient and weak divergence operators. In Sections 4 and 5, we describe our weak finite element spaces and the new weak Galerkin method and prove the existence and uniqueness of the solution. Section 6 is devoted to a-priori error estimates and Section 7 is concerned with how to recover an approximate stress solution from the approximate displacement solution. Then, we discuss extensions of the method to more general meshes in Section 8. Finally, in Section 9, we present some numerical results. 3

2. Governing equations and some notation Consider an isotropic elastic body in the configuration space Ω ⊂ Rd , d = 2, 3, and assume that Ω is an open bounded and connected domain with Lipschitz continuous boundary ∂Ω. The kinematic model for linear elasticity then reads as follows. Given an exterior body force f , find the displacement u ∈ Rd that satisfies −∇ · (2µε(u) + λ(∇ · u)I) = f u = gD

in Ω, on ∂Ω.

(1a) (1b)

Here, ε(u) = 21 (∇u + (∇u)T ) is the linear strain tensor, µ and λ are Lam´e constants, satisfying µ ∈ [µ1 , µ2 ] with 0 < µ1 ≤ µ2 < ∞ and 0 < λ < ∞, and I is a d × d identity tensor. For linear plane strain, the Lam´e constants are given by E Eν ,µ= , λ= (1 + ν)(1 − 2ν) 2(1 + ν) where E is Young’s modulus and ν is Poisson’s ratio. Also, the stress tensor, denoted by σ, is defined by σ = 2µε(u) + λ(∇ · u)I.

(2)

We now introduce some notation for Hilbert spaces and their associated inner products and norms. For any open bounded domain T ⊂ Rd with Lipschitz continuous boundary, let H m (T ) be the usual Sobolev space defined by H m (T ) = {v ∈ L2 (T ) : Dα v ∈ L2 (T ) ∀|α| ≤ m}, where Dα v =

∂ |α| v α α1 ∂x1 ···∂xd d

with |α| = α1 + · · · αd , and let H0m (T ) be a closed

subspace of H m (T ) such that

H0m (T ) = {v ∈ H m (T ) : v = 0 on ∂Ω}. Then, denote by | · |m,T and k · km,T , respectively, the semi-norm and norm in the Sobolev spaces H m (T ), [H m (T )]d , and [H m (T )]d×d . When m = 0, H m (T ) coincides with L2 (T ). In this case, denote by (·, ·)T and k·kT , respectively, the inner product and the norm. The subscript T will be dropped if T = Ω. Also, let h·, ·iT denote the L2 (∂T ) inner product or duality pairing. In addition, define H(div; T) = {τ ∈ [L2 (T)]d×d : ∇ · τ ∈ (L2 (T))d } 4

1

with its natural norm kvkH(div;T) = (kvk2T + k∇ · vk2T ) 2 . We are now in a position to derive a weak formulation of the problem. By using ∇ · (∇u)T = ∇ · ((∇ · u)I), the equation of linear elasticity (1a) can be rewritten as −µ∆u − (µ + λ)∇(∇ · u) = f . (3) Multiplying (3) by a test function v ∈ (H10 (Ω))d and using integration by parts, we arrive at the weak formulation µ(∇u, ∇v) + (µ + λ)(∇ · u, ∇ · v) = (f , v) ∀v ∈ (H10 (Ω))d .

(4)

Throughout this paper, we assume the H 2 -regularity estimate for the solution [7, 25]. kuk2 + λk∇ · uk1 ≤ Ckf k, (5) where C is independent of λ. 3. Weak gradient and divergence operators Weak Galerkin methods employ weakly defined differential operators to discretize partial differential equations. The goal of this section is to review the definitions of weak gradient and weak divergence operators defined on a space of generalized functions. Let T be any polygonal domain with boundary ∂T and e ⊂ ∂T be an edge (in two dimensions) or a face (in three dimensions) of T . In order to define the weak divergence operator, we introduce a space of weak vector-valued functions on T as follows: 1

V(T ) = {v = {v0 , vb } : v0 ∈ [L2 (T )]d , vb · n ∈ [H − 2 (∂T )]d }, where n is the unit outward normal vector on ∂T . Here, v0 and vb represent the value of v in T and on the boundary ∂T , respectively. However, it should be noted that vb is not necessarily related to the trace of v0 on ∂T , should a trace be well defined. Definition 3.1. (weak divergence) For any weak vector-valued function v ∈ V(T ), the weak divergence, denoted by ∇w,T · v, is defined as a linear functional in the Sobolev space H 1 (T ): (∇w,T · v, φ)T = −(v0 , ∇φ)T + hvb · n, φi∂T where n is the unit outward normal vector on ∂T . 5

∀φ ∈ H 1 (T ),

We now introduce another space of weak vector-valued functions on T to define the weak gradient operator: 1

W(T ) = {v = {v0 , vb } : v0 ∈ [L2 (T )]d , vb ∈ [H 2 (∂T )]d }, Definition 3.2. (weak gradient) For any weak vector-valued function v ∈ W(T ), the weak gradient is defined as a linear functional ∇w v in the Sobolev space of [H(div; T )]d : (∇w v, ψ)T = −(v0 , ∇ · ψ)T + hvb , ψni∂T

∀ψ ∈ [H(div; T)]d ,

where n is the unit outward normal vector on ∂T . 4. A weak Galerkin finite element method Let Th denote a shape regular triangulation of Ω into simplicial or rectangular meshes and let Eh be the union of all edges/faces of the mesh elements. For each T ∈ Th , hT denotes the diameter of T and the mesh size of Th is defined as h = maxT ∈Th hT . In order to define our discrete weak gradient and discrete weak divergence operators, denote by Pk (T ) and Pk (e) the set of polynomials of degree no greater than k on T and e ⊂ ∂T , respectively. Here, e represents an edge in two dimensions and a face in three dimensions. It will also be useful to define Pk1 ,k2 (T ) (in two dimensions), the space of polynomials of degree ≤ k1 in x1 and ≤ k2 in x2 . When k1 = k2 = k, we define Qk (T ) = Pk1 ,k2 (T ). Analogously, we can define Pk1 ,k2 ,k3 (T ) and Qk (T ) in three dimensions. Our local weak finite element space on T ∈ Th is defined by V (T ) = {v = {v0 , vb } : v0 ∈ [P0 (T )]d , vb ∈ [P0 (e)]d

∀e ⊂ ∂T }.

Then, the global weak finite element space Vh is assembled by the local space V (T ) through a continuity condition of the component vb on each edge/face e ∈ Eh . That is, Vh = {v = {v0 , vb } : v0 |T ∈ [P0 (T )]d ∀T ∈ Th , vb |e ∈ [P0 (e)]d ∀e ∈ Eh }. Note that Vh0 = {v = {v0 , vb } ∈ Vh : vb = 0 on ∂Ω} is a subspace of Vh . 6

In order to define our discrete weak gradient operator, denote the standard Raviart-Thomas space of index k by RTk and RT[k] on a simplicial mesh and a rectangular mesh, respectively. Then, the local RTk and RT[k] are defined as follows [10]: RTk (T ) = [Pk (T )]d + xPk (T ),

d = 2, 3

and ( Pk+1,k (T ) × Pk,k+1 (T ) in 2d, RT[k] (T ) = Pk+1,k,k (T ) × Pk,k+1,k (T ) × Pk,k,k+1 (T ) in 3d. Let Σ(T ) denote a tensor space each row of which belongs to RT0 (T ) on a simplicial mesh and RT[0] (T ) on a rectangular mesh. That is, Σ(T ) = {τ ∈ (L2 (T ))d×d : τ j ∈ RT0 (T ) or RT[0] (T ), j = 1, · · · , d}, where τ j denotes the jth row of τ . Then, let Σh = {τ ∈ (L2 (Ω))d×d : τ |T ∈ Σ(T )}. Note that Σh is not a subspace of H(div; Ω). Now, we define the discrete weak divergence and discrete weak gradient operators as follows: Definition 4.1. (discrete weak divergence) The discrete weak divergence of v ∈ V (T ), denoted by ∇w,T · v, is the unique polynomial in P0 (T ) that satisfies (∇w,T · v, φ)T = −(v0 , ∇φ)T + hvb · n, φi∂T

∀φ ∈ P0 (T ),

(6)

where n is the unit outward normal vector on ∂T . Definition 4.2. (discrete weak gradient) The discrete weak gradient of v ∈ V (T ), denoted by ∇w,T v, is the unique polynomial in Σ(T ) that satisfies (∇w,T v, ψ)T = −(v0 , ∇ · ψ)T + hvb , ψni∂T where n is the unit outward normal vector on ∂T .

7

∀ψ ∈ Σ(T ),

(7)

For each v ∈ Vh , the discrete weak divergence ∇w · v and the discrete weak gradient ∇w v are computed locally using (6) and (7), respectively: (∇w · v)|T = ∇w,T · (v|T ) ∀T ∈ Th , (∇w v)|T = ∇w,T (v|T ) ∀T ∈ Th . Now, we turn to present the lowest-order weak Galerkin method. The weak Galerkin method to solve (1) is to find uh = {uh0 , uhb } ∈ Vh with uhb = Qb gD on ∂Ω such that ∀v = {v0 , vb } ∈ Vh0 ,

aw (uh , v) = (f , v0 )

(8)

where the bilinear form aw (·, ·) is defined by X X (∇w · w, ∇w · v)T (∇w w, ∇w v)T + (µ + λ) aw (w, v) = µ T ∈Th

T ∈Th

and Qb is the local L2 -projection onto a set of all piecewise constant vector functions on Eh as defined in the next section. 5. Existence and uniqueness of the solution In the weak Galerkin finite element space Vh , we introduce the following semi-norm ! 21 X |||v||| := k∇w vk2T . (9) T ∈Th

Lemma 5.1. The semi-norm |||·||| defined in (9) is a norm in the linear space Vh0 . Proof. Let |||v||| = 0. On each T ∈ Th , we can use integration by parts to rewrite the definition (7) of the discrete weak gradient of v as follows: (∇w v, τ )T = −(v0 , ∇ · τ )T + hvb , τ ni∂T = (∇v0 , τ )T − hv0 − vb , τ ni∂T = −hv0 − vb , τ ni∂T ∀τ ∈ Σ(T ). Since ∇w v = 0 on each T ∈ Th , we have −hv0 − vb , τ ni∂T = 0 ∀τ ∈ Σ(T ). 8

(10)

Then, as v0 − vb ∈ P0 (∂T ), we can find a unique element τ ∗ ∈ Σ(T ) such that τ ∗ n = −(v0 − vb ) on ∂T [10], hence kv0 − vb k∂T = 0. This implies that v0 = vb on ∂T for any T ∈ Th . Since vb is a single-valued function on the boundaries of the elements, v0 = vb ≡ c on Ω for some constant c. On the other hand, vb = 0 on ∂Ω, hence v0 = vb ≡ 0 on Ω. That is, v ≡ 0 on Ω. It is now trivial to prove the coercivity of the bilinear form aw with respect to the norm defined in (9). Lemma 5.2. There exists a positive constant α such that α|||v|||2 ≤ aw (v, v)

∀v ∈ Vh0 .

(11)

Theorem 5.1. There exists a unique solution to the weak Galerkin finite element method defined in (8). Proof. Due to the finite dimensionality of the solution space, uniqueness and existence are equivalent. Therefore, it suffices to prove the uniqueness of the solution. Suppose there are two solutions, uh,1 and uh,2 , to the discrete problem (8). Let eu = uh,1 − uh,2 . Then, it is trivial to see that aw (eu , v) = 0 ∀v ∈ Vh0 . As eu ∈ Vh0 , we have aw (eu , eu ) = 0. Then, in light of Lemmas 5.2 and 5.1, one can obtain that eu = 0, which implies the uniqueness of the solution. 6. Error analysis 6.1. Projection operators We begin this section by defining several projection operators that will be useful in the error analysis. On each element T ∈ Th , let Q0 be the L2 projection from [L2 (T )]d onto [P0 (T )]d . Also, on each edge or face e ∈ Eh , let Qb be the L2 -projection from [L2 (e)]d onto [P0 (e)]d . Then, define the L2 -projection Qh onto the weak Galerkin finite element space Vh by Qh := {Q0 , Qb }. In addition, Ph denotes the local L2 -projection from L2 (T ) onto P0 (T ). Finally, let Ph : [L2 (Ω)]d×d → Σh be the local L2 -projections onto Σh . Then, using a standard scaling argument, we can prove the following approximation estimates [26]: X kQ0 v − vk2T ≤ Ch2 kvk21 ∀v ∈ [H 1 (Ω)d , (12a) T ∈Th

9

X

T ∈Th

X

T ∈Th

kPh τ − τ k2T ≤ Ch2m kτ k2m kPh w − wk2T ≤ Ch2 kwk21

∀τ ∈ [H 1 (Ω)]d×d , m = 0, 1, ∀w ∈ H 1 (Ω).

(12b) (12c)

Lemma 6.1. For any w ∈ (H 1 (Ω))d , the following commutativity properties hold true: ∇w · (Qh w) = Ph (∇ · w), ∇w (Qh w) = Ph (∇w).

(13a) (13b)

Proof. We shall prove only (13a) as (13b) can be proved in a similar fashion. For any φ ∈ P0 (T ), use (6), integration by parts, and the definition of the projections to see (∇w · (Qh w), φ)T = −(Q0 w, ∇φ)T + hQb w · n, φi∂T = hw · n, φi∂T = (∇ · w, φ)T = (Ph (∇ · w), φ)T . Since both ∇w ·(Qh w) and Ph (∇·w) are constants on T , (13a) holds true. 6.2. Discrete H 1 - and L2 -norm error estimates of the displacement We will first derive an error equation for the weak Galerkin approximation uh and Qh u. Lemma 6.2. For any v = {v0 , vb } ∈ Vh0 , aw (Qh u − uh , v) = l1 (v) + l2 (v), where l1 (v) = µ

X

T ∈Th

h(∇u − Ph ∇u) · n, v0 − vb i∂T ,

l2 (v) = (µ + λ)

X

T ∈Th

h∇ · u − Ph (∇ · u), (v0 − vb ) · ni∂T .

10

(14)

Proof. Let vh = {v0 , vb } be any test function in Vh0 . Then, multiply (3) by v0 , integrate both sides, and sum over T ∈ Th to get X − (∇ · (µ∇u + (µ + λ)(∇ · u)I), v0 )T = (f , v0 ) = aw (uh , v). T ∈Th

On the other hand, by applying integration by parts multiple times and invoking (7), (6), (13a), and (13b), the left-hand side can be rewritten as X − (∇ · (µ∇u + (µ + λ)(∇ · u)I), v0 )T T ∈Th

=µ

X

T ∈Th

−µ =µ

X

T ∈Th

X

T ∈Th

−µ =µ

T ∈Th

T ∈Th

= −µ +µ

T ∈Th

X

T ∈Th

T ∈Th

−µ =µ

h∇u · n, v0 i∂T − (µ + λ)

(∇ · u, ∇ · v0 )T X

T ∈Th

X

T ∈Th

h∇u · n, v0 i∂T − (µ + λ)

h∇u · n, v0 i∂T − (µ + λ)

T ∈Th

X

T ∈Th

X

T ∈Th

T ∈Th

T ∈Th

−µ

h∇u · n, v0 i∂T − (µ + λ)

(∇w · Qh u, ∇ · v0 )T h∇ · u, v0 · ni∂T X

T ∈Th

X

T ∈Th

X

T ∈Th

(∇w Qh u, ∇w v)T + (µ + λ)

X

h∇ · u, v0 · ni∂T

(∇ · (∇w Qh u), v0 )T − (µ + λ)

X

T ∈Th

h∇ · u, v0 · ni∂T

(Ph ∇ · u, ∇ · v0 )T

X

h∇w Qh u · n, v0 i∂T + (µ + λ)

X

+µ

T ∈Th

(∇w Qh u, ∇v0 )T + (µ + λ)

X

X

X

(Ph ∇u, ∇v0 )T + (µ + λ)

X

X

−µ

(∇u, ∇v0 )T + (µ + λ)

h∇w · Qh u, v0 · ni∂T

h∇ · u, v0 · ni∂T

X

T ∈Th

(∇w · Qh u, ∇w · v)T

h∇w Qh u · n, v0 − vb i∂T + (µ + λ)

X

T ∈Th

h∇u · n, v0 − vb i∂T − (µ + λ) 11

(∇∇w · Qh u, v0 )T

X

T ∈Th

X

T ∈Th

h∇w · Qh u, (v0 − vb ) · ni∂T

h∇ · u, (v0 − vb ) · ni∂T

=µ

X

T ∈Th

−µ

(∇w Qh u, ∇w v)T + (µ + λ)

X

T ∈Th

X

T ∈Th

(∇w · Qh u, ∇w · v)T

h(∇u − Ph ∇u) · n, v0 − vb i∂T

− (µ + λ)

X

T ∈Th

= aw (Qh u, v) − µ − (µ + λ)

X

T ∈Th

h∇ · u − Ph ∇ · u, (v0 − vb ) · ni∂T X

T ∈Th

h(∇u − Ph ∇u) · n, v0 − vb i∂T

h∇ · u − Ph ∇ · u, (v0 − vb ) · ni∂T

= aw (Qh u, v) − l1 (v) − l2 (v). Therefore, (14) follows. The following lemma will be useful for the subsequent error analysis. Lemma 6.3. X

T ∈Th

2 2 h−1 T kv0 − vb k∂T ≤ C|||v|||

∀v ∈ Vh0 .

(15)

Proof. As discussed in the proof of Lemma 5.1, one can choose a unique τ ∗ ∈ Σ(T ) such that τ ∗ n = −(v0 − vb ) on ∂T on each T ∈ Th . Then, taking τ = τ ∗ in (10) and using the Cauchy-Schwarz inequality, one can see kv0 − vb k2∂T ≤ k∇w vkT kτ ∗ kT . Therefore, to prove (15), it is sufficient to prove that kτ ∗ k2T ≤ ChT kv0 − vb k2∂T ,

(16)

which requires more specific information on τ ∗ . Let {qi , i = 1, · · · , ne } be a set of local basis functions of Σ(T ) such that qi ·nj = δij on ej , j = 1, · · · , ne , where nj is an outward normal vector on the edge or face ej ⊂ ∂T and ne is the number of edges or faces of T . See [24] for an explicit expression of qi . Then, it is easy to see that τ ∗ can be defined as follows: τ

∗,j

=

ne X i=1

−(v0j − vbj )|ei qi , j = 1, · · · , d, 12

where τ ∗,j , j = 1, 2, denotes the jth row of τ ∗ and v0j and vbj denotes the jth component of v0 and vb , respectively. In [24], it was proved when T is a triangle or a rectangle that kqi k2T ≤ C(ne )|ei |2 , 1 ≤ i ≤ ne ,

(17)

where |ei | is the measure of ei . A similar argument can be used to obtain an analogous result in three dimensions: kqi k2T ≤ C(ne )hT |ei |, 1 ≤ i ≤ ne .

(18)

Note that (17) and (18) imply that kτ ∗ k2T = ≤

d X j=1

kτ ∗,j k2T

d X ne X j=1 i=1

((v0j − vbj )|ei )2 kqi k2T

ne X 2 X ≤C ((v0j − vbj )|ei )2 hT |ei | i=1 j=1

≤ ChT kv0 − vb k2∂T .

Lemma 6.4. For any v ∈ Vh0 , aw (Qh u − uh , v) ≤ Chkf k|||v|||.

(19)

Proof. In light of (14), we have aw (Qh u − uh , v) ≤ |l1 (v)| + |l2 (v)|. Therefore, we need to bound |l1 (v)| and |l2 (v)|. Using the inverse inequality, (12b), and (15), one can show that X |l1 (v)| = |µ h(∇u − Ph ∇u) · n, v0 − vb i∂T | T ∈Th

≤µ

X

T ∈Th

k(∇u − Ph ∇u) · nk∂T kv0 − vb k∂T 13

≤C

X

T ∈Th

≤C

−1

(k∇u − Ph ∇ukT )(hT 2 kv0 − vb k∂T )

X

T ∈Th

hT k∇uk21,T

≤ Chkuk2 |||v||| ∀v ∈

! 12 Vh0 .

X

T ∈Th

2 h−1 T kv0 − vb k∂T

! 21

In a similar way, |l2 (v)| = |(µ + λ)

X

T ∈Th

h∇ · u − Ph ∇ · u, (v0 − vb ) · ni∂T |

≤ C(µ + λ)hk∇ · uk1 |||v||| ∀v ∈ Vh0 . Finally, (19) can be obtained by using the H 2 -regularity (5) of the solution. Lemma 6.5. |||Qh u − uh ||| ≤ Chkf k.

(20)

Proof. Taking v = Qh u − uh ∈ Vh0 in (19) and using Lemma 5.2, we can easily obtain the result. Now, in preparation for the L2 -norm error estimate, we prove the following lemma. Lemma 6.6. kv0 k ≤ C|||v||| v ∈ Vh0 . Proof. Since v0 is a piecewise constant vector function, one can find τ ∈ Σh ∩ (H(div; Ω))d [10] such that ∇ · τ = −v0

and kτ kH(div) ≤ Ckv0 k.

Take ψ = τ in (7), sum over T ∈ Th , and use (21) and the fact that X hvb , τ ni∂T = 0 T ∈Th

to see that kv0 k2 =

X

T ∈Th

kv0 k2T =

X

T ∈Th

(∇w v, τ )T ≤

X

T ∈Th

k∇w vkT kτ kT .

Then, the result follows from the Cauchy-Schwarz inequality and (21). 14

(21)

We now have the following auxiliary error estimates. Corollary 6.1. aw (Qh u − uh , Qh u − uh ) ≤ Ch2 kf k2

(22)

kQ0 u − uh0 k ≤ Chkf k.

(23)

and Proof. These are immediate results of (19) and Lemma 6.6 with v = Qh u − uh , respectively, combined with Lemma 6.5. Finally, we shall present our main theorem that shows the optimal error estimates of the displacement in the broken H 1 - and L2 -norms. Theorem 6.1. Let u and uh be the solution of (4) and (8), respectively. Then, we have the following error estimates. k∇u − ∇w uh k ≤ Chkf k, ku − uh0 k ≤ Chkf k,

(24a) (24b)

where C is independent of h and λ. Proof. First, to prove the broken H 1 -norm error estimate (24a), note that ∇u−∇w uh = (∇u−∇w Qh u)+∇w (Qh u−uh ) = (∇u−Ph ∇u)+∇w (Qh u−uh ) thanks to (13b). Then, (24a) follows from Young’s inequality, (12b), and (5). The L2 error estimate (24b) follows from a similar argument using the L2 -projection Q0 and (23). 7. Stress recovery and an L2 -norm error estimate When a standard Galerkin method based on the primal formulation is used for linear elasticity, the stress tensor can be approximated by post processing the numerical solution of the displacement. A straightforward post processing involves numerical differentiations of the displacement and one loses accuracy while doing so. In this section, we present a simple, and rather trivial, post processing technique applied to the numerical displacement uh to obtain an approximation to the stress.

15

Let uh be the solution to the weak Galerkin method (8) and define the approximate stress σh by σh = 2µεw (uh ) + λ(∇w · uh )I,

(25)

where εw (uh ) = 12 (∇w uh +(∇w uh )T ). From (22), we can deduce the following error estimate for ∇w · (Qh u − uh ): √ λk∇w · (Qh u − uh )k ≤ Chkf k. However, a careful error analysis can prove an improved error estimate for ∇w · (Qh u − uh ), which will play a critical role in the λ-independent error estimate of the stress. Lemma 7.1. λk∇w · (Qh u − uh )k ≤ Chkf k,

(26)

where C is independent of λ and h.

Proof. Using the definition of the weak divergence (6), one can R easily check that ∇w · (Qh u − uh ) ∈ L20 (Ω), where L20 (Ω) = {v ∈ L2 (Ω) : Ω v dx = 0.}. Therefore, there exists z ∈ [H01 (Ω)]d [27] such that ∇ · z = ∇w · (Qh u − uh ) and kzk1 ≤ Ck∇w · (Qh u − uh )k. Then, invoking (27), the definition of Ph , (13b), (12b), and (20), λk∇w · (Qh u − uh )k2 X =λ (∇w · (Qh u − uh ), ∇w · (Qh u − uh ))T T ∈Th

≤ (µ + λ)

= (µ + λ)

X

T ∈Th

X

T ∈Th

= (µ + λ)

X

T ∈Th

= (µ + λ)

X

T ∈Th

(∇w · (Qh u − uh ), ∇w · (Qh u − uh ))T

(∇w · (Qh u − uh ), ∇ · z)T (∇w · (Qh u − uh ), Ph (∇ · z))T (∇w · (Qh u − uh ), ∇w · Qh z)T

= aw (Qh u − uh , Qh z) − µ

X

T ∈Th

16

(∇w (Qh u − uh ), ∇w Qh z)T

(27)

≤ C(hkf k + |||Qh u − uh |||) |||Qh z||| ≤ Chkf k kPh (∇z)k ≤ Chkf k kzk1 ≤ Chkf k k∇w · (Qh u − uh )k. Thus, (26) follows. We shall now present an L2 error estimate for the stress. Theorem 7.1. Let σ be the stress tensor defined in (2) and σh be the approximate stress tensor defined in (25). Then, we have kσ − σh k ≤ Chkf k,

(28)

where C is independent of h and λ. Proof. We use a similar argument in the proof of Theorem 6.1. Note that kσ − σh k ≤ 2µkε(u) − εw (uh ))k + λk∇ · u − ∇w · uh k ≤ C(k∇u − ∇w uh k + λk∇ · u − ∇w · uh k). Since we already proved (24a), it suffices to prove that the second term λk∇ · u − ∇w · uh k is also bounded by Chkf k. To end this, note that ∇·u−∇w ·uh = (∇·u−∇w ·Qh u)+∇w ·(Qh u−uh ) = (∇·u−Ph (∇·u))+∇w ·(Qh u−uh ). Then, we can obtain (28) by using the Young’s inequality, (12c), (5), and (26). 8. Extensions of the method to more general meshes The new method that was presented and analyzed in the previous sections assumed a simplicial or rectangular mesh in two and three dimensions. In this section, we briefly discuss how we can extend the method to more general polygonal or polyhedral meshes. The weak finite element spaces for the displacement and the discrete weak divergence operator can trivially be extended to any polygonal or polyhedral meshes. So, the main issue here would be to define a function space for the discrete weak gradient operator on more general meshes that can reduce to the standard, lowest-order Raviart-Thomas element on a simplicial or rectangular mesh in two and 17

three dimensions. Very recently, Chen and Wang [24] constructed H(div)conforming finite elements with minimal possible degrees of freedom on arbitrary convex polygons and certain polyhedra. They showed that their new finite elements coincide with the lowest-order Raviart-Thomas elements when the mesh is a simplicial or a rectangular mesh in two and three dimensions. Therefore, we can define the discrete weak gradient operator using their local basis functions on more general meshes. Then, the error analysis with some modifications can go through without any problem. 9. Numerical Experiments In this section, we present some numerical results to validate the accuracy and efficiency of the new methods we proposed in previous sections. Here we solve the problems on both triangular and rectangular meshes in two dimensions. 9.1. Smooth solution In order to confirm the optimal convergence rates that we proved in Section 6, we first consider a pure displacement problem with a homogeneous boundary condition. Let Ω be the unit square and λ = 1.0, µ = 0.5. The external body force f in (1) is chosen so that the exact displacement is u = [ e(x−y) x(1 − x)y(1 − y), sin (πx) sin (πy) ].

(29)

Tables 1 and 2 summarize the errors and convergence rates on triangular and rectangular meshes, respectively. Optimal convergence rates for the displacement in the discrete H 1 - and L2 -norms and for the stress in the L2 -norm are observed on both type of meshes. 9.2. Locking-free property The second test problem is intended to show not only the optimal convergence rates but also the locking-free property of the proposed method. Here we choose the body force, f , so that the exact solution on the computational domain Ω = (0, 1)2 is the following. 1 1 u = [ sin x · sin y + x, cos x · cos y + y ]. λ λ

(30)

We impose a Dirichlet boundary condition that is calculated from the known solution. Note that the solution is designed to satisfy ∇ · u = λ2 → 0 as 18

1/h 8 16 32 64 128 256

ku − uh0 k Rate k∇u − ∇w uh k 7.0657e-02 8.6489e-01 3.3626e-02 1.07 4.5312e-01 1.6520e-02 1.03 2.3010e-01 8.2191e-03 1.01 1.1559e-01 4.1043e-03 1.00 5.7877e-02 2.0515e-03 1.00 2.8950e-02

Rate 0.93 0.98 0.99 1.00 1.00

kσ − σh k Rate 7.8708e-01 4.1064e-01 0.94 2.0807e-01 0.98 1.0445e-01 0.99 5.2284e-02 1.00 2.6151e-02 1.00

Table 1: Convergence study for a homogeneous Dirichlet boundary problem (29) with λ = 1.0, µ = 0.5 on uniform triangular meshes.

1/h 8 16 32 64 128 256

ku − uh0 k Rate k∇u − ∇w uh k 8.2069e-02 7.7773e-01 4.0433e-02 1.02 3.9848e-01 2.0129e-02 1.01 2.0065e-01 1.0053e-02 1.00 1.0052e-01 5.0251e-03 1.00 5.0290e-02 2.5124e-03 1.00 2.5149e-02

Rate 0.96 0.99 1.00 1.00 1.00

kσ − σh k Rate 6.4266e-01 3.2778e-01 0.97 1.6486e-01 0.99 8.2572e-02 1.00 4.1306e-02 1.00 2.0656e-02 1.00

Table 2: Convergence study for a homogeneous Dirichlet boundary problem (29) with λ = 1.0, µ = 0.5 on uniform rectangular meshes.

λ → ∞ and f is independent of λ. In this case, standard linear or bilinear finite elements become over constrained by the divergence-free condition as λ → ∞, hence suffer from Poisson locking. Two different values of λ, one small and one large, are chosen to investigate the locking-free property of the new method. Tables 3 and 4 show the errors and convergence rates for the displacement and stress when λ = 1.0 and λ = 106 , respectively, on triangular meshes and Tables 5 and 6 summarize the results on rectangular meshes. It is evident from these results that the performance of the new method does not deteriorate even for a very large value of λ, as predicted by the mathematical analysis in Section 6. 9.3. Solution with a corner singularity For the last test problem, we consider a linear elasticity problem with a known analytic solution [28] in an L-shaped domain in two dimensions. The L-shaped domain Ω can be obtained by removing a section [0, 1] × [−1, 0] from a square domain [−1, 1] × [−1, 1]. See Figure 1 for an illustration of Ω and a coarse rectangular mesh Th . Let κ = 3 − 4ν. Then, the exact solution 19

1/h 8 16 32 64 128 256

ku − uh0 k 4.6961e-02 2.3488e-02 1.1745e-02 5.8725e-03 2.9362e-03 1.4681e-03

Rate 1.00 1.00 1.00 1.00 1.00

k∇u − ∇w uh k 3.7157e-02 1.8620e-02 9.3166e-03 4.6593e-03 2.3298e-03 1.1649e-03

Rate 1.00 1.00 1.00 1.00 1.00

kσ − σh k 3.1935e-02 1.6007e-02 8.0096e-03 4.0057e-03 2.0030e-03 1.0015e-03

Rate 1.00 1.00 1.00 1.00 1.00

Table 3: Convergence study for the locking susceptible solution (30) with λ = 1.0, µ = 0.5 on uniform triangular meshes.

1/h 8 16 32 64 128 256

ku − uh0 k 2.1662e-02 1.0846e-02 5.4251e-03 2.7128e-03 1.3564e-03 6.7821e-04

Rate 1.00 1.00 1.00 1.00 1.00

k∇u − ∇w uh k 3.7072e-02 1.8587e-02 9.3027e-03 4.6528e-03 2.3266e-03 1.1634e-03

Rate 1.00 1.00 1.00 1.00 1.00

kσ − σh k 3.2204e-02 1.6121e-02 8.0610e-03 4.0303e-03 2.0150e-03 1.0075e-03

Rate 1.00 1.00 1.00 1.00 1.00

Table 4: Convergence study for the locking susceptible solution (30) with λ = 106 , µ = 0.5 on uniform triangular meshes.

in polar coordinates (r, θ) is u = (u1 , u2 ), where 1 α r [(κ − Q(α + 1)) cos (αθ) − α cos ((α − 2)θ)], 2µ 1 α r [(κ + Q(α + 1)) sin (αθ) + α sin ((α − 2)θ)]. u2 = 2µ u1 =

Here, α is the solution of the equation     3π 3π sin α + α sin =0 2 2 and Q is given by


cos (α − 1) 3π 4
Q=− . cos (α + 1) 3π 4

(31a) (31b)

(32)

Then, the displacement u defined in (31) satisfies the linear elasticity equations (1) with the body force f = 0. The Dirichlet boundary condition was found by using the exact displacement. Numerical experiments are carried 20

1/h 8 16 32 64 128 256

ku − uh0 k 6.0291e-02 3.0152e-02 1.5077e-02 7.5385e-03 3.7693e-03 1.8846e-03

Rate 1.00 1.00 1.00 1.00 1.00

k∇u − ∇w uh k 3.9639e-02 1.9819e-02 9.9097e-03 4.9548e-03 2.4774e-03 1.2387e-03

Rate 1.00 1.00 1.00 1.00 1.00

kσ − σh k 3.4326e-02 1.7164e-02. 8.5820e-03 4.2910e-03 2.1455e-03 1.0728e-03

Rate 1.00 1.00 1.00 1.00 1.00 1.00

Table 5: Convergence study for the locking susceptible solution (30)with λ = 1.0, µ = 0.5 on uniform rectangular meshes.

1/h 8 16 32 64 128 256

ku − uh0 k 3.2106e-02 1.6066e-02 8.0343e-03 4.0174e-03 2.0087e-03 1.0044e-03

Rate 1.00 1.00 1.00 1.00 1.00

k∇u − ∇w uh k 3.9639e-02 1.9819e-02 9.9097e-03 4.9548e-03 2.4774e-03 1.2387e-03

Rate 1.00 1.00 1.00 1.00 1.00

kσ − σh k 3.4335e-02 1.7165e-02 8.5822e-03 4.2910e-03 2.1455e-03 1.0728e-03

Rate 1.00 1.00 1.00 1.00 1.00

Table 6: Convergence study for the locking susceptible solution (30) with λ = 106 , µ = 0.5 on uniform rectangular meshes.

out with Young’s modulus E = 105 and Poisson’s ratio ν = 0.4. Also, the solution α to the equation (32) was numerically approximated and the resulting values for α and Q are α = 0.5444837367825 and Q = 0.5430755788367. Since α < 1, all the stress components have a singularity in a neighborhood of the origin, while the displacement components are continuous. Indeed, u ∈ (H 1+α− (Ω))2 and σ ∈ (H α− (Ω))2×2 for  > 0 [29]. Therefore, one can expect a linear convergence for the displacement in the L2 -norm, but only a suboptimal convergence rate of α for the stress when the proposed weak Galerkin method is used. Table 7 confirms these theoretical expectations. References [1] D. Arnold, R. Winther, Mixed finite elements for elasticity, Numer. Math. 92 (2002) 401–419. [2] D. Arnold, G. Awanou, Rectangular mixed finite elements for elasticity, Math. Models Methods Appl. Sci. 15 (2005) 1417–1429.

21

(-1, 1)

(1, 1)

(0,0)

(-1, -1)

Figure 1: L-shaped domain and a coarse mesh Th .

1/h 8 16 32 64 128 256

ku − uh0 k Rate kσ − σh k Rate 3.6402e-06 6.7953e-01 1.8350e-06 0.99 4.7258e-01 0.52 9.1981e-07 1.00 3.2633e-01 0.53 4.5994e-07 1.00 2.2456e-01 0.54 2.0087e-03 1.00 1.5427e-01 0.54 1.1476e-07 1.00 1.0588e-01 0.54

Table 7: Convergence study for the singular solution (31) in an L-shaped domain.

[3] S.-Y. Yi, Nonconforming mixed finite element methods for linear elasticity using rectangular elements in two and three dimensions, Calcolo 42 (2005) 115–133. [4] S.-Y. Yi, A new nonconforming mixed finite element method for linear elasticity, Math. Models Methods Appl. Sci. 16 (2006) 979–999. [5] D. N. Arnold, R. S. Falk, R. Winther, Mixed finite element methods for linear elasticity with weakly imposed symmetry, Math. Comp. 76 (2007) 1699–1723. [6] R. S. Falk, Nonconforming finite element methods for the equations of linear elasticity, Math. Comp. 57 (1991) 529–550. [7] S. C. Brenner, L.-Y. Sung, Linear finite element methods for planar linear elasticity, Math. Comp. 59 (1992) 321–338.

22

[8] B. Rivi`ere, S. Shaw, M. F. Wheeler, J. R. Whiteman, Discontinuous Galerkin finite element methods for linear elasticity and quasistatic linear viscoelasticity, Numer. Math. 95 (2003) 347–376. [9] D. A. Di Pietro, S. Nicaise, A locking-free discontinuous Galerkin method for linear elasticity in locally nearly incompressible heterogeneous media, Appl. Numer. Math. 63 (2013) 105–116. [10] F. Brezzi, M. Fortin, Mixed and hybrid finite element methods, volume 15 of Springer Series in Computational Mathematics, SpringerVerlag, New York, 1991. [11] Y. Zhang, A weak Galerkin mixed finite element method for linear elasticity equations, ProQuest LLC, Ann Arbor, MI, 2015. Thesis (Ph.D.)– Oklahoma State University. [12] C. Wang, J. Wang, R. Wang, R. Zhang, A locking-free weak Galerkin finite element method for elasticity problems in the primal formulation, J. Comput. Appl. Math. 307 (2016) 346–366. [13] G. Chen, X. Xie, A robust weak Galerkin finite element method for linear elasticity with strong symmetric stresses, Comput. Methods Appl. Math. 16 (2016) 389–408. [14] J. Wang, X. Ye, A weak Galerkin finite element method for second-order elliptic problems, J. Comput. Appl. Math. 241 (2013) 103–115. [15] L. Mu, J. Wang, Y. Wang, X. Ye, A weak Galerkin mixed finite element method for biharmonic equations, in: Numerical solution of partial differential equations: theory, algorithms, and their applications, volume 45 of Springer Proc. Math. Stat., Springer, New York, 2013, pp. 247–277. [16] L. Mu, J. Wang, X. Ye, Weak Galerkin finite element methods for the biharmonic equation on polytopal meshes, Numer. Methods Partial Differential Equations 30 (2014) 1003–1029. [17] C. Wang, J. Wang, An efficient numerical scheme for the biharmonic equation by weak Galerkin finite element methods on polygonal or polyhedral meshes, Comput. Math. Appl. 68 (2014) 2314–2330. 23

[18] C. Wang, J. Wang, A hybridized weak Galerkin finite element method for the biharmonic equation, Int. J. Numer. Anal. Model. 12 (2015) 302–317. [19] Q. H. Li, J. Wang, Weak Galerkin finite element methods for parabolic equations, Numer. Methods Partial Differential Equations 29 (2013) 2004–2024. [20] J. Wang, X. Ye, A weak Galerkin finite element method for the stokes equations, Adv. Comput. Math. 42 (2016) 155–174. [21] L. Mu, J. Wang, X. Ye, A stable numerical algorithm for the Brinkman equations by weak Galerkin finite element methods, J. Comput. Phys. 273 (2014) 327–342. [22] L. Mu, J. Wang, X. Ye, S. Zhang, A weak Galerkin finite element method for the Maxwell equations, J. Sci. Comput. 65 (2015) 363–386. [23] L. Mu, J. Wang, X. Ye, A new weak Galerkin finite element method for the Helmholtz equation, IMA J. Numer. Anal. 35 (2015) 1228–1255. [24] W. Chen, Y. Wang, Minimal degree H(curl) and H(div) conforming finite elements on polytopal meshes, Math. Comp. 86 (2017) 2053–2087. [25] S. C. Brenner, L. R. Scott, The mathematical theory of finite element methods, volume 15 of Texts in Applied Mathematics, Springer-Verlag, New York, 1994. [26] P. Ciarlet, The Finite Element Method for Elliptic Problems, Studies in Mathematics and its Applications, Elsevier Science, 1978. [27] V. Girault, P.-A. Raviart, Finite element methods for Navier-Stokes equations, volume 5 of Springer Series in Computational Mathematics, Springer-Verlag, Berlin, 1986. Theory and algorithms. [28] M. Ainsworth, B. Senior, Aspects of an adaptive hp-finite element method: adaptive strategy, conforming approximation and efficient solvers, Comput. Methods Appl. Mech. Engrg. 150 (1997) 65–87. Symposium on Advances in Computational Mechanics, Vol. 2 (Austin, TX, 1997).

24

[29] I. Babuska, M. Suri, The h-p version of the finite element method with quasi-uniform meshes, RAIRO Mod´el. Math. Anal. Num´er. 21 (1987) 199–238.

25