A Mixed Discontinuous Galerkin Method for Linear Elasticity with Strongly Imposed Symmetry

TL;DR

This paper introduces a mixed discontinuous Galerkin method for linear elasticity that achieves optimal error estimates and demonstrates sharp convergence orders through numerical validation.

Contribution

The paper develops a novel mixed DG scheme with arbitrary order discontinuous finite element spaces for linear elasticity, providing rigorous stability and error analysis.

Findings

Optimal error estimates for stress and displacement are established.

Numerical results confirm the sharpness of convergence orders.

The method is well-posed for any polynomial degree k ≥ 0.

Abstract

In this paper, we study a mixed discontinuous Galerkin (MDG) method to solve linear elasticity problem with arbitrary order discontinuous finite element spaces in $d$-dimension ($d=2,3$). This method uses polynomials of degree $k+1$ for the stress and of degree $k$ for the displacement ($k\geq 0$). The mixed DG scheme is proved to be well-posed under proper norms. Specifically, we prove that, for any $k \geq 0$, the $H({\rm div})$-like error estimate for the stress and $L^2$ error estimate for the displacement are optimal. We further establish the optimal $L^2$ error estimate for the stress provided that the $\mathcal{P}_{k+2}-\mathcal{P}_{k+1}^{-1}$ Stokes pair is stable and $k \geq d$. We also provide numerical results of MDG showing that the orders of convergence are actually sharp.

Tables

Table 1. (a) Linear elasticity: 𝒫 1 − 1 − 𝒫 0 − 1 superscript subscript 𝒫 1 1 superscript subscript 𝒫 0 1 \mathcal{P}_{1}^{-1}-\mathcal{P}_{0}^{-1} , 2D uniform grids

$1/h$	${\Vert u-u_h\Vert }_{0,\Omega }$	$h^n$	${\Vert 𝝈-{𝝈}_h\Vert }_{0,\Omega }$	$h^n$	${\Vert {\mathrm{div}}_h\left(𝝈-{𝝈}_h\right)\Vert }_{0,\Omega }$	$h^n$
4	0.135877	—	0.445892	—	3.839803	—
8	0.067302	1.01	0.177473	1.33	1.936584	0.99
16	0.033543	1.00	0.080752	1.14	0.970346	1.00
32	0.016757	1.00	0.039257	1.04	0.485431	1.00

Table 2. (a) Linear elasticity: 𝒫 1 − 1 − 𝒫 0 − 1 superscript subscript 𝒫 1 1 superscript subscript 𝒫 0 1 \mathcal{P}_{1}^{-1}-\mathcal{P}_{0}^{-1} , 2D uniform grids

$1/h$	${\Vert u-u_h\Vert }_{0,\Omega }$	$h^n$	${\Vert 𝝈-{𝝈}_h\Vert }_{0,\Omega }$	$h^n$	${\Vert {\mathrm{div}}_h\left(𝝈-{𝝈}_h\right)\Vert }_{0,\Omega }$	$h^n$
4	0.135877	—	0.445892	—	3.839803	—
8	0.067302	1.01	0.177473	1.33	1.936584	0.99
16	0.033543	1.00	0.080752	1.14	0.970346	1.00
32	0.016757	1.00	0.039257	1.04	0.485431	1.00

Table 3. (b) Linear elasticity: 𝒫 2 − 1 − 𝒫 1 − 1 superscript subscript 𝒫 2 1 superscript subscript 𝒫 1 1 \mathcal{P}_{2}^{-1}-\mathcal{P}_{1}^{-1} , 2D uniform grids

$1/h$	${\Vert u-u_h\Vert }_{0,\Omega }$	$h^n$	${\Vert 𝝈-{𝝈}_h\Vert }_{0,\Omega }$	$h^n$	${\Vert {\mathrm{div}}_h\left(𝝈-{𝝈}_h\right)\Vert }_{0,\Omega }$	$h^n$
4	0.0198206	—	0.0425699	—	0.5850957	—
8	0.0050264	1.98	0.0079777	2.42	0.1483264	1.98
16	0.0012616	1.99	0.0017692	2.17	0.0372321	1.99
32	0.0003158	2.00	0.0004284	2.05	0.0093191	2.00

Table 4. (c) Linear elasticity: 𝒫 3 − 1 − 𝒫 2 − 1 superscript subscript 𝒫 3 1 superscript subscript 𝒫 2 1 \mathcal{P}_{3}^{-1}-\mathcal{P}_{2}^{-1} , 2D uniform grids

$1/h$	${\Vert u-u_h\Vert }_{0,\Omega }$	$h^n$	${\Vert 𝝈-{𝝈}_h\Vert }_{0,\Omega }$	$h^n$	${\Vert {\mathrm{div}}_h\left(𝝈-{𝝈}_h\right)\Vert }_{0,\Omega }$	$h^n$
4	0.00217252	—	0.00341919	—	0.06370927	—
8	0.00027548	2.98	0.00024533	3.80	0.00805005	2.98
16	0.00003456	2.99	0.00001627	3.91	0.00100892	3.00
32	0.00000432	3.00	0.00000104	3.96	0.00012620	3.00

Table 5. (a) Linear elasticity: 𝒫 1 − 1 − 𝒫 0 − 1 superscript subscript 𝒫 1 1 superscript subscript 𝒫 0 1 \mathcal{P}_{1}^{-1}-\mathcal{P}_{0}^{-1} , 3D uniform grids

$1/h$	${\Vert u-u_h\Vert }_{0,\Omega }$	$h^n$	${\Vert 𝝈-{𝝈}_h\Vert }_{0,\Omega }$	$h^n$	${\Vert {\mathrm{div}}_h\left(𝝈-{𝝈}_h\right)\Vert }_{0,\Omega }$	$h^n$
2	0.235741	—	1.221265	—	7.534218	—
4	0.127481	0.89	0.536012	1.19	4.420875	0.77
8	0.063704	1.00	0.210303	1.35	2.294909	0.95

Table 6. (a) Linear elasticity: 𝒫 1 − 1 − 𝒫 0 − 1 superscript subscript 𝒫 1 1 superscript subscript 𝒫 0 1 \mathcal{P}_{1}^{-1}-\mathcal{P}_{0}^{-1} , 3D uniform grids

$1/h$	${\Vert u-u_h\Vert }_{0,\Omega }$	$h^n$	${\Vert 𝝈-{𝝈}_h\Vert }_{0,\Omega }$	$h^n$	${\Vert {\mathrm{div}}_h\left(𝝈-{𝝈}_h\right)\Vert }_{0,\Omega }$	$h^n$
2	0.235741	—	1.221265	—	7.534218	—
4	0.127481	0.89	0.536012	1.19	4.420875	0.77
8	0.063704	1.00	0.210303	1.35	2.294909	0.95

Table 7. (b) Linear elasticity: 𝒫 2 − 1 − 𝒫 1 − 1 superscript subscript 𝒫 2 1 superscript subscript 𝒫 1 1 \mathcal{P}_{2}^{-1}-\mathcal{P}_{1}^{-1} , 3D uniform grids

$1/h$	${\Vert u-u_h\Vert }_{0,\Omega }$	$h^n$	${\Vert 𝝈-{𝝈}_h\Vert }_{0,\Omega }$	$h^n$	${\Vert {\mathrm{div}}_h\left(𝝈-{𝝈}_h\right)\Vert }_{0,\Omega }$	$h^n$
2	0.0831048	—	0.3641751	—	2.8564400	—
4	0.0227446	1.87	0.0664638	2.45	0.7833919	1.87
8	0.0058207	1.97	0.0123827	2.42	0.2007023	1.96

Equations

$$\{\begin{array}{rlr}\mathcal{A}\mathbf{\sigma }-\mathbf{\epsilon }(u)& =0& \text{in }\Omega ,\\ \mathrm{div}\mathbf{\sigma }& =f& \text{in }\Omega ,\\ u& =0& \text{on }\partial \Omega ,\end{array}$$

$$\hat{\mathbf{p}}$$

$$\hat{u}$$

$$\{\begin{array}{rlr}(\mathcal{A}\mathbf{\sigma },\mathbf{\tau }{)}_{\Omega }+(u,\mathrm{div}\mathbf{\tau }{)}_{\Omega }& =0& \forall \mathbf{\tau }\in \mathbf{\Sigma },\\ (\mathrm{div}\mathbf{\sigma },v{)}_{\Omega }& =(f,v{)}_{\Omega }& \forall v\in V.\end{array}$$

$$\Vert \mathbf{\tau }{\Vert }_{\mathrm{div},\Omega }^2:=\Vert \mathbf{\tau }{\Vert }_{0,\Omega }^2+\Vert \mathrm{div}\mathbf{\tau }{\Vert }_{0,\Omega }^2\forall \mathbf{\tau }\in H(\mathrm{div},\Omega ;\mathbb{S}).$$

$$u\odot v:=\frac{1}{2}(u\otimes v+v\otimes u)\in \mathbb{S}\forall u,v\in {\mathbb{R}}^d,$$

$$\{v\}=\frac{1}{2}(v^{+}+v^{-}),$$

$$[\mathbf{\tau }]={\mathbf{\tau }}^{+}{n}^{+}+{\mathbf{\tau }}^{-}{n}^{-},$$

$$\{\mathbf{\tau }\}=\mathbf{\tau },$$

$${\int }_K\mathrm{div}\mathbf{\tau }\cdot v\mathrm{d}x=-{\int }_K\mathbf{\tau }:\mathbf{\epsilon }(v)\mathrm{d}x+{\int }_{\partial K}(\mathbf{\tau }{n}_K)\cdot v\mathrm{d}s,$$

$$\sum _{K\in {\mathcal{T}}_h}{\int }_{\partial K}(\mathbf{\tau }{n}_K)\cdot v\mathrm{d}s={\int }_{{\mathcal{E}}_h}\{\mathbf{\tau }\}:[[v]]\mathrm{d}s+{\int }_{{\mathcal{E}}_h^i}[\mathbf{\tau }]\cdot \{v\}\mathrm{d}s.$$

$${\nabla }_{h}v{|}_K=\nabla v{|}_K,{\mathrm{div}}_h\mathbf{\tau }{|}_K=\mathrm{div}\mathbf{\tau }{|}_K,$$

$$\{\begin{array}{rlr}\sum _{K\in {\mathcal{T}}_h}(\mathcal{A}\mathbf{\sigma },{\mathbf{\tau }}_h{)}_K+\sum _{K\in {\mathcal{T}}_h}(u,\mathrm{div}{\mathbf{\tau }}_h{)}_K-\sum _{K\in {\mathcal{T}}_h}⟨u,{\mathbf{\tau }}_h{n}_K{⟩}_{\partial K}& =0& \forall {\mathbf{\tau }}_h\in {\mathbf{\Sigma }}_h,\\ -\sum _{K\in {\mathcal{T}}_h}(\mathbf{\sigma },{\mathbf{\epsilon }}_h(v_h){)}_K+\sum _{K\in {\mathcal{T}}_h}⟨\mathbf{\sigma }{n}_K,v_h{⟩}_{\partial K}& =\sum _{K\in {\mathcal{T}}_h}(f,v_h{)}_K& \forall v_h\in V_h.\end{array}$$

$$\{\begin{array}{rlr}\sum _{K\in {\mathcal{T}}_h}(\mathcal{A}{\mathbf{\sigma }}_h,{\mathbf{\tau }}_h{)}_K+\sum _{K\in {\mathcal{T}}_h}(u_h,\mathrm{div}{\mathbf{\tau }}_h{)}_K-\sum _{K\in {\mathcal{T}}_h}⟨{\hat{u}}_h,{\mathbf{\tau }}_h{n}_K{⟩}_{\partial K}& =0& \forall {\mathbf{\tau }}_h\in {\mathbf{\Sigma }}_h,\\ -\sum _{K\in {\mathcal{T}}_h}({\mathbf{\sigma }}_h,{\mathbf{\epsilon }}_h(v_h){)}_K+\sum _{K\in {\mathcal{T}}_h}⟨{\hat{\mathbf{\sigma }}}_h{n}_K,v_h{⟩}_{\partial K}& =\sum _{K\in {\mathcal{T}}_h}(f,v_h{)}_K& \forall v_h\in V_h,\end{array}$$

$$\{\begin{array}{rlr}{\int }_{\Omega }\mathcal{A}{\mathbf{\sigma }}_h:{\mathbf{\tau }}_h\mathrm{d}x+{\int }_{\Omega }{u}_h\cdot {\mathrm{div}}_h{\mathbf{\tau }}_h\mathrm{d}x-{\int }_{{\mathcal{E}}_h}[[{\hat{u}}_h]]:\{{\mathbf{\tau }}_h\}\mathrm{d}s-{\int }_{{\mathcal{E}}_h^i}\{{\hat{u}}_h\}\cdot [{\mathbf{\tau }}_h]\mathrm{d}s& =0& \forall {\mathbf{\tau }}_h\in {\mathbf{\Sigma }}_h,\\ {\int }_{\Omega }{\mathrm{div}}_h{\mathbf{\sigma }}_h\cdot v_h\mathrm{d}x+{\int }_{{\mathcal{E}}_h}\{{\hat{\mathbf{\sigma }}}_h-{\mathbf{\sigma }}_h\}:[[v_h]]\mathrm{d}s+{\int }_{{\mathcal{E}}_h^i}[{\hat{\mathbf{\sigma }}}_h-{\mathbf{\sigma }}_h]\cdot \{v_h\}\mathrm{d}s& ={\int }_{\Omega }f\cdot v_h\mathrm{d}x& \forall v_h\in V_h.\end{array}$$

$$\{\begin{array}{rl}{\hat{u}}_h& =\{u_h\}-\eta [{\mathbf{\sigma }}_h]\mathrm{on}{\mathcal{E}}_h^i,{\hat{u}}_h=0\mathrm{on}{\mathcal{E}}_h^{\partial },\\ {\hat{\mathbf{\sigma }}}_h& =\{{\mathbf{\sigma }}_h\}\mathrm{on}{\mathcal{E}}_h.\end{array}$$

$$[[{\hat{u}}_h]]=0,[{\hat{\mathbf{\sigma }}}_h]=0,\mathrm{and}\{{\hat{\mathbf{\sigma }}}_h-{\mathbf{\sigma }}_h\}=0.$$

$$\{\begin{array}{rlr}{a}_h({\mathbf{\sigma }}_h,{\mathbf{\tau }}_h)+b_h({\mathbf{\tau }}_h,u_h)& =0& \forall {\mathbf{\tau }}_h\in {\mathbf{\Sigma }}_h,\\ b_h({\mathbf{\sigma }}_h,v_h)& =(f,v_h{)}_{\Omega }& \forall v_h\in V_h.\end{array}$$

$$a_h(\mathbf{\sigma },\mathbf{\tau })$$

$$b_h(\mathbf{\tau },v)$$

$$\Vert \mathbf{\tau }{\Vert }_{\ast ,\Omega }^2:={\int }_{\Omega }(|\mathbf{\tau }{|}^2+|{\mathrm{div}}_h\mathbf{\tau }{|}^2)\mathrm{d}x+{\int }_{{\mathcal{E}}_h^i}{\eta }_e{h}_e^{-1}|[\mathbf{\tau }]{|}^2\mathrm{d}s\forall \mathbf{\tau }\in {\mathbf{\Sigma }}_h\cup \mathbf{\Sigma }.$$

$$V_h=V_h^k$$

$${\mathbf{\Sigma }}_h={\mathbf{\Sigma }}_h^{k+1}$$

$$a_h(\mathbf{\sigma },\mathbf{\tau })\lesssim \Vert \mathbf{\sigma }{\Vert }_{\ast ,\Omega }\Vert \mathbf{\tau }{\Vert }_{\ast ,\Omega }\forall \mathbf{\sigma },\mathbf{\tau }\in {\mathbf{\Sigma }}_h\cup \mathbf{\Sigma }.$$

$${\int }_{\Omega }{r}_e(w)\cdot v_h\mathrm{d}x=-{\int }_{e}w\cdot \{v_h\}\mathrm{d}s\forall v_h\in V_h.$$

$$\Vert r_e(w){\Vert }_{0,\Omega }\lesssim h_e^{-1/2}\Vert w{\Vert }_{0,e}.$$

$$\Vert r_e(w){\Vert }_{0,\Omega }^2\le \frac{1}{2}\Vert w{\Vert }_{0,e}(\Vert r_e(w{)}^{+}{\Vert }_{0,e}+\Vert r_e(w{)}^{-}{\Vert }_{0,e})\lesssim h_e^{-1/2}\Vert w{\Vert }_{0,e}\Vert r_e(w){\Vert }_{0,\Omega },$$

$$b_h(\mathbf{\tau },v_h)$$

$$b_h(\mathbf{\tau },v)$$

$$b_h(\mathbf{\tau },v_h)$$

$$\lesssim \Vert v_h{\Vert }_{0,\Omega }{\left(\Vert {\mathrm{div}}_h\mathbf{\tau }{\Vert }_{0,\Omega }^2+{\int }_{{\mathcal{E}}_h^i}{h}_e^{-1}|[\mathbf{\tau }]{|}^2\mathrm{d}s\right)}^{1/2}$$

$$\le \Vert v_h{\Vert }_{0,\Omega }\Vert \mathbf{\tau }{\Vert }_{\ast ,\Omega }.$$

$$b_h(\mathbf{\tau },v)$$

$$a_h({\mathbf{\tau }}_h,{\mathbf{\tau }}_h)\ge C\Vert {\mathbf{\tau }}_h{\Vert }_{\ast ,\Omega }^2\forall {\mathbf{\tau }}_h\in Z_h,$$

Full text

Keywords. Mixed DG method, linear elasticity, well-posedness, a priori error analysis

Mathematics Subject Classification. 65N30, 65M60

1 Introduction

In this paper, we present a mixed discontinuous Galerkin (MDG) method for the following linear elasticity problem:

[TABLE]

where $u:\Omega\mapsto\mathbb{R}^{d}$ and $\boldsymbol{\sigma}:\Omega\mapsto\mathbb{S}$, denote displacement and stress, respectively. Here, $\mathbb{S}$ represents the space of real symmetric matrices of order $d\times d$. The tensor $\mathcal{A}:\mathbb{S}\mapsto\mathbb{S}$ is assumed to be bounded and symmetric positive definite, and the linearized strain tensor is denoted by $\boldsymbol{\varepsilon}(u)=(\nabla u+(\nabla u)^{t})/2$.

For the mixed methods for linear elasticity problem (1.1), it is very challenging to develop the stable mixed finite element methods because the stress tensor needs to be symmetric. One approach to circumvent this difficulty is to introduce the antisymmetric part of $\nabla u$ as a new variable, and hence, to enforce stress symmetry weakly [2, 6, 11, 22, 26, 37, 29]. Another approach is to use the composite element for the stress [36, 5]. The first stable non-composite conforming mixed finite element method for plane elasticity was proposed by Arnold and Winther in 2002 [7], and analogs of the results in the 3D case were reported in [1, 3]. In this class of elements, the displacement is discretized by discontinuous piecewise $\mathcal{P}_{k}^{-1}$ ($k\geq 1$) polynomial, while the stress is discretized by the conforming $\mathcal{P}_{k+2}$ tensors whose divergence is $\mathcal{P}_{k}$ vector on each triangle. In recent years, Hu and Zhang [33, 34] and Hu [32] proposed a family of conforming mixed elements for $\mathbb{R}^{d}$ that apply the $\mathcal{P}_{k+1}-\mathcal{P}_{k}$ pair for the stress and displacement when $k\geq d$. These elements also admit a unified theory and a relatively easy implementation. The lower order conforming approximations of stress were also considered in [35], and a simpler stress element with jump stabilization term for the displacement [19].

Because of the lack of suitable conforming mixed elasticity elements, several authors have resorted to the nonconforming elements [8, 4, 28], where the optimal convergence order for the displacement can be proved under the full elliptic regularity assumption but the convergence order of $L^{2}$ error for stress is still suboptimal. To improve the convergence order for stress, an interior penalty mixed finite element method using Crouzeix-Raviart nonconforming linear element to approximate each component of the symmetric stress was studied in [17]. In [40], Wu, Gong, and Xu proposed two classes of interior penalty mixed finite elements for linear elasticity of arbitrary order in arbitrary dimension, where the stability is guaranteed by introducing the nonconforming face-bubble spaces based on the local decomposition of discrete symmetric tensors.

Discontinuous Galerkin (DG) methods have been applied to solve various differential equations due to their flexibility in constructing feasible local shape function spaces and the advantage to capture non-smooth or oscillatory solutions effectively. The DG methods are attracting the interest of many applied mathematicians and engineers because they discretize the equations in an element-by-element fashion, and glue each element through numerical traces, which can give rise to locally conservative methods. In [9], Arnold, Brezzi, Cockburn, and Marini proposed a unified framework for the devising and analysis of most DG methods for second-order elliptic equations. The LDG method, which is introduced in [24], is one of several discontinuous Galerkin methods which are being vigorously studied [18, 9, 21, 23]. As proposed in [18, Equ. (2.4)], the numerical traces for second-order elliptic equations have the general expressions as

[TABLE]

where $u$ and $\boldsymbol{p}$ are the approximations of primal variable and flux, respectively. In most literature, the parameter $C_{22}$ is taken as [math] or $\mathcal{O}(h)$ so that the resulting scheme is of the category of primal DG method. When taking $C_{22}$ as $\mathcal{O}(h^{-1})$, the penalty term on the jump of $\boldsymbol{p}$ leads to a mixed DG scheme [31].

For linear elasticity problem, a primal LDG method was studied in [20], where the discontinuous $\mathcal{P}_{k}^{-1}-\mathcal{P}_{k+1}^{-1}$ pairs were used to approximate the stress and displacement for $k\geq 0$. In the weak formulation, two penalty terms for stress and displacement are adopted, however, the error analysis was only given for the case when the penalty term of the stress vanishes, i.e. $C_{22}=0$.

In this paper, we study the mixed LDG method for solving linear elasticity by discontinuous $\mathcal{P}_{k+1}^{-1}-\mathcal{P}_{k}^{-1}$ finite element pairs for the stress and displacement with $k\geq 0$ for any spatial dimension in a unified fashion. Our contributions are twofold. First, by introducing a mesh-dependent norm for the stress, we give a prior error analysis, which shows that optimal $L^{2}$-error estimate for displacement and optimal $H_{h}({\rm div})$ error estimate for stress. Second, when the $\mathcal{P}_{k+2}-\mathcal{P}_{k+1}^{-1}$ Stokes pair is stable and $k\geq d$, we prove the optimal $L^{2}$ error estimate for the stress by the BDM projection [14] and a symmetrization technique.

The rest of the paper is organized as follows. In Section 2, we derive the mixed DG scheme to solve the linear elasticity problem. Then based on Brezzi theory, we prove the well-posedness of the scheme in Section 3, and the optimal convergence rates are obtained for both stress and displacement variables in Section 4. In addition, the optimal $L^{2}$ error estimate for the stress is shown in Section 5. In Section 6, numerical tests are given for solving the linear elasticity problems by the mixed LDG methods, and the numerical results verify the theoretical error analysis. Finally, we give several concluding remarks in the last section.

2 Mixed DG method for linear elasticity problem

In this section, we study a mixed discontinuous Galerkin method for the linear elasticity problem (1.1), whose weak formulation reads: Find $(\boldsymbol{\sigma},u)\in\boldsymbol{\Sigma}\times V$ such that

[TABLE]

Here, $V=L^{2}(\Omega;\mathbb{R}^{d})$ denotes the space of vector-valued functions which are square-integrable with the $L^{2}$ norm, and $\boldsymbol{\Sigma}=H({\rm div},\Omega;\mathbb{S})$ consists of square-integrable symmetric matrix fields with square-integrable divergence, and the corresponding norm is defined by

[TABLE]

For the symmetric tensor space $\mathbb{S}$, we define the inner products by $\boldsymbol{\sigma}:\boldsymbol{\tau}=\sum_{i,j=1}^{d}\sigma_{ij}\tau_{ij}$ for any $\boldsymbol{\sigma},\boldsymbol{\tau}\in\mathbb{S}$. Further, we define the symmetric tensor product $\odot$ as

[TABLE]

where $u\otimes v$ is a tensor with $u_{i}v_{j}$ as its $(i,j)$-th entry.

2.1 DG notation

We introduce some notation before presenting the mixed DG scheme. Given a bounded domain $D\subset\mathbb{R}^{d}$ and a positive integer $m$, $H^{m}(D)$ is the Sobolev space with the corresponding usual norm and semi-norm, which are denoted respectively by $\|\cdot\|_{m,D}$ and $|\cdot|_{m,D}$. We abbreviate them by $\|\cdot\|_{m}$ and $|\cdot|_{m}$, respectively, when $D$ is chosen as $\Omega$. The $L^{2}$-inner product on $D$ and $\partial D$ are denoted by $(\cdot,\cdot)_{D}$ and $\langle\cdot,\cdot\rangle_{\partial D}$, respectively. $\|\cdot\|_{D}$ and $\|\cdot\|_{\partial D}$ are the norms of Lebesgue spaces $L^{2}(D)$ and $L^{2}(\partial D)$, respectively. We assume $\Omega$ is a polygonal domain and denote by $\{\mathcal{T}_{h}\}_{h}$ a family of triangulations of $\overline{\Omega}$, with the minimal angle condition satisfied. Let $h_{K}={\rm diam}(K)$ and $h=\max\{h_{K}:K\in\mathcal{T}_{h}\}$. Denote by ${\cal E}_{h}$ the union of the boundaries of the elements $K$ of $\mathcal{T}_{h}$, ${\cal E}_{h}^{i}$ is the set of interior edges and ${\cal E}_{h}^{\partial}={\cal E}_{h}\backslash{\cal E}_{h}^{i}$ is the set of boundary edges. Let $e$ be the common edge of two elements $K^{+}$ and $K^{-}$, and $\boldsymbol{n}^{i}$ = $\boldsymbol{n}|_{\partial K^{i}}$ be the unit outward normal vector on $\partial K^{i}$ with $i=+,-$. For any vector-valued function $v$ and tensor-valued function $\boldsymbol{\tau}$, let $v^{\pm}$ = $v|_{\partial K^{\pm}}$, $\boldsymbol{\tau}^{\pm}$ = $\boldsymbol{\tau}|_{\partial K^{\pm}}$. Then, we define the average $\{\cdot\}$, jump $[\cdot]$ and tensor jump $\llbracket\cdot\rrbracket$ as follows:

[TABLE]

where $n$ is the outward unit normal vector on $\partial\Omega$. Let us give the following identities which are used often in this section. For any vector-valued function $v$ and tensor-valued function $\boldsymbol{\tau}$, all being continuously differentiable over $K$, we have the following integration by parts formula:

[TABLE]

and the following identity:

[TABLE]

Throughout this paper, we shall use letter $C$ to denote a generic positive constant independent of $h$ which may stand for different values at its different occurrences. The notation $x\lesssim y$ means $x\leq Cy$. For piecewise smooth vector-valued function $v$ and tensor-valued function $\boldsymbol{\tau}$, let $\nabla_{h}$ and ${\rm div}_{h}$ be defined by the relation

[TABLE]

on any element $K\in{\cal T}_{h}$, respectively.

2.2 Mixed LDG scheme

Now, let us introduce the mixed LDG formulation for (1.1). We denote the piecewise vector and symmetric matrix valued discrete spaces by $V_{h}$ and $\boldsymbol{\Sigma}_{h}$, respectively. We multiply (1.1) by arbitrary test functions $\boldsymbol{\tau}_{h}\in\boldsymbol{\Sigma}_{h}$ and $v_{h}\in V_{h}$, respectively, and integration by parts over the element $K\in\mathcal{T}_{h}$ to obtain

[TABLE]

Let $\widehat{V}_{h}$ and $\widehat{\boldsymbol{\Sigma}}_{h}$ be the piecewise vector and symmetric matrix valued discrete spaces on $\mathcal{E}_{h}$, respectively. The approximate solution $(\boldsymbol{\sigma}_{h},u_{h})$ is then defined by using the weak formulation (2.5), namely

[TABLE]

where the numerical traces $\widehat{u}_{h}\in\widehat{V}_{h}$ and $\widehat{\boldsymbol{\sigma}}_{h}\in\widehat{\boldsymbol{\Sigma}}_{h}$ need to be suitably defined to ensure the stability of the method and to enhance its accuracy. By the identity (2.4) and integration by parts (2.3), we get from (2.6) that

[TABLE]

Similar to the discussion for Poisson problem in [31], we choose mixed LDG numerical traces as follow:

[TABLE]

In such choice, it is easy to see that the numerical traces are single valued. Further, we can see that if $u_{h}$ and $\boldsymbol{\sigma}_{h}$ are replaced by the exact solution $u$ and $\boldsymbol{\sigma}$, then $\widehat{u}_{h}=u|_{{\cal E}_{h}}$ and $\widehat{\boldsymbol{\sigma}}_{h}=\boldsymbol{\sigma}|_{{\cal E}_{h}}$ on ${\cal E}_{h}$. That is, the numerical traces are consistent. Moreover, we have

[TABLE]

Then, we obtain the mixed LDG formulation for (1.1): Find $(\boldsymbol{\sigma}_{h},u_{h})\in\boldsymbol{\Sigma}_{h}\times V_{h}$ such that

[TABLE]

Here, we choose $\eta=\eta_{e}h_{e}^{-1},\eta_{e}=\mathcal{O}(1)$, and define

[TABLE]

Moreover, we define the following star norm

[TABLE]

In the following subsections, we prove the boundedness, stability and consistency of the mixed LDG formulation (2.9) when choosing

[TABLE]

for $k\geq 0$, which lead to the optimal order of convergence.

3 Well-posedness of the mixed LDG method

The well-posedness of the mixed LDG methods (2.9) comes from the boundedness and the stability.

Boundedness.

It is easy to check by Cauchy-Schwarz inequality that $a_{h}(\cdot,\cdot)$ satisfies

[TABLE]

The remaining task is the boundedness of $b_{h}(\cdot,\cdot)$. To this end, let us recall the lifting operator $r_{e}:(L^{2}({\cal E}_{h}))^{d}\rightarrow V_{h}$ defined by

[TABLE]

Then, we have the following lemma (see also [9, 16]).

Lemma 3.1

For any edge $e\in\partial K$, it holds

[TABLE]

Proof. By taking $v_{h}=r_{e}(w)$ in (3.2) and applying the inverse inequality, we obtain

[TABLE]

which gives rise to (3.3).

Lemma 3.2

It holds that

[TABLE]

Proof. In light of Lemma 3.1, we have for any $v_{h}\in V_{h}$

[TABLE]

Furthermore, for any $v\in V\cap H^{1}(\Omega;\mathbb{R}^{d})$,

[TABLE]

Here, we use the trace inequality in the last step.

Stability.

According to the theory of mixed method, the stability of the saddle point problem (2.9) is the corollary of the following two conditions [13, 15]:

K-ellipticity: There exists a constant $C>0$, independent of the grid size such that

[TABLE]

where $Z_{h}=\{\boldsymbol{\tau}_{h}\in\boldsymbol{\Sigma}_{h}~{}|~{}b_{h}(\boldsymbol{\tau}_{h},v_{h})=0\;\forall v_{h}\in V_{h}\}$. 2. 2.

The discrete inf-sup condition: There exists a constant $C>0$, independent of the grid size such that

[TABLE]

First, we prove the inf-sup condition (3.7) in the following lemma.

Lemma 3.3** (Inf-sup condition)**

When choosing $\boldsymbol{\Sigma}_{h}\times V_{h}=\boldsymbol{\Sigma}_{h}^{k+1}\times V_{h}^{k}$ for $k\geq 0$, the discrete inf-sup condition (3.7) holds true for mixed LDG method (2.9) of linear elasticity problem.

Proof. In [40], Wu, Gong, and Xu introduced a class of nonconforming finite element spaces for $k\geq 0$ that

[TABLE]

Thanks to the Lemma 3.3 and Lemma 4.1 in [40], we know that for any $v_{h}\in V_{h}$, there exists a $\bar{\boldsymbol{\tau}}_{h}\in\Sigma_{k+1,h}^{(1)}$ such that

[TABLE]

Note that $\Sigma_{k+1,h}^{(1)}\subset\boldsymbol{\Sigma}_{h}^{k+1}$ and the property of $\Sigma_{k+1,h}^{(1)}$ implies that

[TABLE]

Here, we use the fact that $\{v_{h}\}$ is of degree $k$ on the edge. Therefore, for any $v_{h}\in V_{h}^{k}$

[TABLE]

Then, we finish the proof.

Theorem 3.4

The mixed LDG scheme (2.9) is well-posed for $(\boldsymbol{\Sigma}_{h}^{k+1},\|\cdot\|_{*,\Omega})$ and $(V_{h}^{k},\|\cdot\|_{0,\Omega})$.

Proof. In light of the boundedness and Lemma 3.3, we only need to prove the K-ellipticity (3.6). By the definition of lifting operator (3.2), we have

[TABLE]

which implies that

[TABLE]

With the help of the Lemma 3.1, we see that

[TABLE]

Let $\eta_{0}=\inf_{e\in{\cal E}_{h}^{i}}\eta_{e}$ be a positive constant that independent of the grid size. Then,

[TABLE]

Then, we finish the proof.

Remark 3.5

From Lemma 3.1, we, Gong, can see that the penalty term $\int_{{\cal E}_{h}^{i}}\eta_{e}h_{e}^{-1}[\boldsymbol{\sigma}_{h}]\cdot[\boldsymbol{\tau}_{h}]\,\mathrm{d}s$ can be replaced by $\sum_{e\in{\cal E}_{h}^{i}}\int_{\Omega}\eta_{e}r_{e}([\boldsymbol{\sigma}_{h}])\cdot r_{e}([\boldsymbol{\tau}_{h}])\,\mathrm{d}x$, and the well-posedness of the corresponding scheme can be proved similarly with a modified norm $\|\boldsymbol{\tau}\|_{*,\Omega}^{2}:=\int_{\Omega}(|\boldsymbol{\tau}|^{2}+|{\rm div}_{h}\boldsymbol{\tau}|^{2}+\sum_{e\in{\cal E}_{h}^{i}}|r_{e}([\boldsymbol{\tau}])|^{2})\,\mathrm{d}x$.

4 A priori error estimates in energy norms

Lemma 4.1

Assume the solution $(\boldsymbol{\sigma},u)\in\Sigma\times H^{1}(\Omega;\mathbb{R}^{d})$, we have

[TABLE]

Proof. It can be seen that $[\boldsymbol{\sigma}]=0$ and $\llbracket u\rrbracket=0$ on ${\cal E}_{h}^{i}$ as $(\boldsymbol{\sigma},u)\in\Sigma\times H^{1}(\Omega;\mathbb{R}^{d})$. Therefore,

[TABLE]

Hence, we prove the first equality in (4.1). On the other hand,

[TABLE]

which implies the second equality in the lemma.

By combining Lemma 4.1 and the well-posedness of mixed LDG formulation (2.9), we have the following a priori error estimates.

Theorem 4.2

Let $(\boldsymbol{\sigma}_{h},u_{h})$ be the solution of the mixed LDG problem (2.9), and $(\boldsymbol{\sigma},u)\in\Sigma\times H^{1}(\Omega;\mathbb{R}^{d})$ be the solution of (1.1). Then,

[TABLE]

Proof. Define

[TABLE]

which satisfies discrete inf-sup condition based on the well-posedness of (2.9). In the light of Lemma 4.1 and the boundedness (3.1), (3.4) and (3.5), we have for any $(\boldsymbol{\tau}_{h},v_{h})\in\boldsymbol{\Sigma}_{h}^{k+1}\times V_{h}^{k}$,

[TABLE]

By triangle inequality, we finish the proof.

For $(\boldsymbol{\sigma},u)\in H^{k+2}(\Omega;\mathbb{S})\times H^{k+1}(\Omega;\mathbb{R}^{d})$, it is well-known that the Scott-Zhang interpolation [39] $I_{h}^{r}$ satisfies:

[TABLE]

Hence, we have the following theorem.

Theorem 4.3

Assume that the solution of (1.1) satisfies $(\boldsymbol{\sigma},u)\in H^{k+2}(\Omega;\mathbb{S})\times H^{k+1}(\Omega;\mathbb{R}^{d})$. Then, the solution of the mixed LDG problem (2.9) satisfies

[TABLE]

5 $L^{2}$ error estimate of stress

In this section, we prove the optimal $L^{2}$ error estimate of $\boldsymbol{\sigma}$ provided that the Stokes pair $\mathcal{P}_{k+2}-\mathcal{P}_{k+1}^{-1}$ is stable and $k\geq d$.

First, we recall the definition of classical BDM projection $\Pi_{h}^{\rm BDM}$ [14]. Given a function $q\in H({\rm div},\Omega;\mathbb{R}^{d})$, the restriction of $\Pi_{h}^{\rm BDM}$ to $K$ is defined as the element of $\mathcal{P}_{k+1}(K;\mathbb{R}^{d})$ such that

[TABLE]

where

[TABLE]

Let $\mathbb{M}$ be the space of real matrices of order $d\times d$. In light of the BDM projection (5.1), on each $K\in\mathcal{T}_{h}$, we first define a matrix-valued function $\widetilde{\boldsymbol{\sigma}}_{h}$ as the only element of $\mathcal{P}_{k+1}(K;\mathbb{M})$ through the numerical solution $\boldsymbol{\sigma}_{h}$ and $\widehat{\boldsymbol{\sigma}}_{h}$ in (2.8):

[TABLE]

where

[TABLE]

Here, the $\boldsymbol{\nabla}$ is regarded as the row-wise operator, i.e.,

[TABLE]

Define the following space

[TABLE]

Then, we have the following lemma.

Lemma 5.1

The $\widetilde{\boldsymbol{\sigma}}_{h}$ in (5.2) is well-defined, and

[TABLE]

Proof. Since (5.2) can be viewed as the row-wise BDM projection, then the well-posedness and (5.3a) follows directly by the definition of $\Pi_{h}^{\rm BDM}$, and by the fact that the normal component of the numerical trace for the flux is single-valued. Let $\boldsymbol{\delta}=\widetilde{\boldsymbol{\sigma}}_{h}-\boldsymbol{\sigma}_{h}$, then

[TABLE]

Then, (5.3b) follows easily by the standard scaling argument; see [12].

Next, we symmetrize $\widetilde{\boldsymbol{\sigma}}_{h}$ by the Stokes pair $\mathcal{P}_{k+2}-\mathcal{P}_{k+1}^{-1}$. A similar technique can be found in [25, 29, 27].

Lemma 5.2

Suppose that the Stokes pair $\mathcal{P}_{k+2}-\mathcal{P}_{k+1}^{-1}$ is stable on the grid $\mathcal{T}_{h}$. Having $\widetilde{\boldsymbol{\sigma}}_{h}$ defined in (5.2), there exists a matrix-valued function $\widetilde{\boldsymbol{\tau}}_{h}\in{\rm BDM}_{k+1}^{d\times d}$ such that $\boldsymbol{\sigma}_{h}^{\star}:=\widetilde{\boldsymbol{\sigma}}_{h}+\widetilde{\boldsymbol{\tau}}_{h}\in H({\rm div},\Omega;\mathbb{S})$, and

[TABLE]

Proof. We construct a divergence-free term $\widetilde{\boldsymbol{\tau}}_{h}=\boldsymbol{\rm curl}\rho_{h}$ where $\rho_{h}$ satisfies

For $d=2$: $\rho_{h}\in H^{1}(\Omega;\mathbb{R}^{2})$ is a vector-valued function and $\rho_{h}|_{K}\in\mathcal{P}_{k+2}(K;\mathbb{R}^{2})$; 2. 2.

For $d=3$: $\rho_{h}\in H^{1}(\Omega;\mathbb{M})$ is a matrix-valued function and $\rho_{h}|_{K}\in\mathcal{P}_{k+2}(K;\mathbb{M})$.

For the 2D case, the $\mathrm{curl}$ operator is a rotation of the operator $\nabla$ (i.e., $\mathrm{curl}=(-\partial_{y},\partial_{x})$) and applies on each entry of the vector $\rho_{h}$. For the 3D case, the $\mathrm{curl}$ operator applies on each row of the matrix $\rho_{h}$. By direct calculation, the symmetry of $\widetilde{\boldsymbol{\sigma}}_{h}+\widetilde{\boldsymbol{\tau}}_{h}$ is equivalent to the following equation,

[TABLE]

where $\boldsymbol{\rm skw}\boldsymbol{\tau}:=(\boldsymbol{\tau}-\boldsymbol{\tau}^{T})/2$. For a scalar function $v$ or a vector-valued function $v=(v_{1},v_{2},v_{3})^{T}$, we further define

[TABLE]

Then, the proof can be divided into the following two cases:

For $n=2$: from [10], we have $\boldsymbol{\rm skw}(\boldsymbol{\rm curl}\rho_{h})=\frac{1}{2}\boldsymbol{\rm Skw}_{2}(\mathrm{div}\rho_{h})$. Thus, (5.5) can be written as:

[TABLE]

The stability of Stokes pair $\mathcal{P}_{k+2}-\mathcal{P}_{k+1}^{-1}$ then implies that there exists a $\rho_{h}\in\{v\in H^{1}(\Omega;\mathbb{R}^{2}):~{}v|_{K}\in\mathcal{P}_{k+2}(K;\mathbb{R}^{2})\}$ satisfying (5.6) and

[TABLE] 2. 2.

For $n=3$: from [10], we have $\boldsymbol{\rm skw}(\boldsymbol{\rm curl}\rho_{h})=-\frac{1}{2}\boldsymbol{\rm Skw}_{3}(\mathrm{div}~{}\Xi\rho_{h})$, where $\Xi$ is an algebraic operator defined as ${\Xi}\rho_{h}=\rho_{h}^{T}-\mathrm{tr}(\rho_{h})\boldsymbol{I}$. Denoting $\eta_{h}=\Xi\rho_{h}$, it is obvious that $\rho_{h}=\Xi^{-1}\eta_{h}=\eta_{h}^{T}-\frac{1}{2}\mathrm{tr}(\eta_{h})\boldsymbol{I}$. Thus, (5.5) can be written as:

[TABLE]

Again, there exists a $\eta_{h}\in\{\tau\in H^{1}(\Omega;\mathbb{M}):~{}\tau|_{K}\in\mathcal{P}_{k+2}(K;\mathbb{M})\}$ satisfying (5.7) and

[TABLE]

To summarize, we obtain the desired $\widetilde{\boldsymbol{\tau}}_{h}=\boldsymbol{\rm curl}\rho_{h}$ that satisfies (5.4). This completes the proof.

We are now in the position to prove the optimal $L^{2}$ error estimate.

Theorem 5.3

Assume that the Stokes pair $\mathcal{P}_{k+2}-\mathcal{P}_{k+1}^{-1}$ is stable on $\mathcal{T}_{h}$ and $k\geq d$. Assume further that the solution of (1.1) satisfies $(\boldsymbol{\sigma},u)\in H^{k+2}(\Omega;\mathbb{S})\times H^{k+1}(\Omega;\mathbb{R}^{d})$. Then, the solution of the mixed LDG problem (2.9) satisfies

[TABLE]

where $\|\boldsymbol{\sigma}\|_{\mathcal{A},\Omega}^{2}:=(\mathcal{A}\boldsymbol{\sigma},\boldsymbol{\sigma})_{\Omega}^{1/2}$.

Proof. By (2.6), (5.2) and Lemma 5.1, we have that for any $v_{h}\in V_{h}$,

[TABLE]

By Lemma 5.2, the symmetrized variable $\boldsymbol{\sigma}_{h}^{\star}=\widetilde{\boldsymbol{\sigma}}_{h}+\widetilde{\boldsymbol{\tau}}_{h}$ is piecewise $\mathcal{P}_{k+1}(K;\mathbb{S})$ and belongs to $H({\rm div},\Omega;\mathbb{S})$. Further, the divergence-free of $\widetilde{\boldsymbol{\tau}}_{h}$ implies that

[TABLE]

In [33, 34], Hu and Zhang constructed the conforming $\mathcal{P}_{k+1}-\mathcal{P}_{k}^{-1}$ mixed methods for linear elasticity on simplicial grids when $k\geq d$. Hu also show that (cf. [32, Remark 3.1]), when $k\geq d$, there exists a projection $\Pi_{h}^{c}$ such that,

[TABLE]

By (5.9) and (5.10a), we have

[TABLE]

Taking $\boldsymbol{\tau}_{h}=\boldsymbol{\sigma}_{h}^{\star}-\Pi_{h}^{c}\boldsymbol{\sigma}$ in the error equation (4.1), we immediately have the $\mathcal{A}$-orthogonality condition:

[TABLE]

Hence, by the energy estimate (4.3), (5.3b) and (5.10b),

[TABLE]

This completes the proof.

Remark 5.4

In the 2D case, the Scott-Vogelius elements $\mathcal{P}_{k+2}-\mathcal{P}_{k+1}^{-1}$ are stable when $k\geq 2$ and the grid does not contain singular vertices (cf. [38, 30]). Hence, in the 2D case, we have the optimal $L^{2}$ estimate when $k\geq 2$ with some mild constrain pertaining to the grids.

6 Numerical examples

In this section, we present some numerical results of the mixed LDG method for linear elasticity problem. The compliance tensor is given by

[TABLE]

where $\boldsymbol{I}_{d}$ is the $d\times d$ identity matrix. In the computation, the Lamé constants are set to be $\mu=1/2$ and $\lambda=1$. The parameter in (2.10a) is chosen as $\eta_{e}=1$ on all $e\in\mathcal{E}_{h}^{i}$.

2D example.

The 2D problem is computed on the unit square $\Omega=(0,1)^{2}$ with a homogeneous boundary condition that $u=0$ on $\partial\Omega$. Let the exact solution be

[TABLE]

The exact stress function $\boldsymbol{\sigma}$ and the load function $f$ can be analytically derived from (1.1) and for a given $u$. Uniform grids with different grid sizes are adopted in the computation.

We list the errors and the rates of convergence of the computed solution in Table 1. The $(k+1)$-th order convergence is observed for both the $L_{2}$ error of $u$ and the $H_{h}({\rm div})$ error of $\boldsymbol{\sigma}$, which is in agreement with Theorem 4.3. Further, we see from Table 1(c) that $\|\boldsymbol{\sigma}-\boldsymbol{\sigma}_{h}\|_{0,\Omega}=\mathcal{O}(h^{4})$ when $k=2$. This convergence rate coincides with the statements in Theorem 5.3, which is also shown sharp from the $L^{2}$ errors of stress in Table 1(a)-1(b).

3D example.

Let the exact solution on the unit cube be

[TABLE]

Again, the true stress function $\boldsymbol{\sigma}$ and the load function $f$ are defined by the relations in (1.1), for the given solution $u$. In Table 2, the errors and the convergence order in various norms are listed when $k=0,1$. The optimal orders of convergence are achieved respectively under the $H_{h}({\rm div})$ norm for the stress and $L^{2}$ norm for the displacement, which confirms Theorem 4.3.

7 Concluding remarks

In this paper, we present the first a priori error analysis of mixed DG method for solving the linear elasticity problem. We provide numerical evidence indicating the sharpness of our estimates, namely, the convergence order of $k+1$ both stress in $H_{h}({\rm div})$-norm and displacement in $L^{2}$-norm with the elements pair $(\boldsymbol{\sigma}_{h},u_{h})\in\boldsymbol{\Sigma}_{h}^{k+1}\times V_{h}^{k}$. The estimate holds for any $k\geq 0$ in arbitrary dimension, making the MDG more meaningful for the linear elasticity as the lower order conforming $\mathcal{P}_{k+1}$-$\mathcal{P}_{k}^{-1}$ elasticity element does not exist on general simplicial grids [40]. Since there is a close connection between elasticity elements and Stokes elements, we also prove the optimal $L^{2}$ error estimate for the stress provided that the $\mathcal{P}_{k+2}-\mathcal{P}_{k+1}^{-1}$ Stokes pair is stable and $k\geq d$.