Residual-Based A Posteriori Error Estimates for Symmetric Conforming Mixed Finite Elements for Linear Elasticity Problems

TL;DR

This paper develops and verifies residual-based a posteriori error estimators for symmetric mixed finite element methods applied to linear elasticity problems, ensuring their stability and efficiency through theoretical proofs and numerical validation.

Contribution

It introduces new a posteriori error estimators for symmetric mixed finite elements in linear elasticity, with proven stability and efficiency, supported by numerical experiments.

Findings

Estimators are stable and efficient.

Numerical examples confirm theoretical results.

Error bounds are effectively estimated.

Abstract

A posteriori error estimators for the symmetric mixed finite element methods for linear elasticity problems of Dirichlet and mixed boundary conditions are proposed. Stability and efficiency of the estimators are proved. Finally, we provide numerical examples to verify the theoretical results.

Figures

Tables

Table 1. Table 6.1. Numerical errors for the first example when λ = 10 𝜆 10 \lambda=10

$h$	${\Vert \sigma -{\sigma }_h\Vert }_A$	order	${\Vert {\nabla }_h\left(u-u_h^{\ast }\right)\Vert }_0$	order	$\eta ({\sigma }_h,{𝒯}_h)$	order	${\Vert A{\sigma }_h-{\epsilon }_h\left(u_h^{\ast }\right)\Vert }_0$	order
$2^{-1}$	6.6998E-01	$-$	7.9544E-01	$-$	1.6615E+01	$-$	4.0073E-02	$-$
$2^{-2}$	5.2451E-02	3.68	6.0585E-02	3.71	1.3585E+00	3.61	9.3899E-03	2.09
$2^{-3}$	3.6139E-03	3.86	4.5839E-03	3.72	1.0918E-01	3.64	7.1387E-04	3.72
$2^{-4}$	2.2714E-04	3.99	3.0676E-04	3.90	7.4510E-03	3.87	4.5925E-05	3.96
$2^{-5}$	1.4193E-05	4.00	1.9600E-05	3.97	4.7919E-04	3.96	2.8824E-06	3.99
$2^{-6}$	8.8742E-07	4.00	1.2347E-06	3.99	3.0263E-05	3.99	1.8040E-07	4.00
$2^{-7}$	5.5567E-08	4.00	7.7435E-08	3.99	1.8992E-06	3.99	1.1306E-08	4.00

Table 2. Table 6.2. Numerical errors for the first example when λ = 10000 𝜆 10000 \lambda=10000

$h$	${\Vert \sigma -{\sigma }_h\Vert }_A$	order	${\Vert {\nabla }_h\left(u-u_h^{\ast }\right)\Vert }_0$	order	$\eta ({\sigma }_h,{𝒯}_h)$	order	${\Vert A{\sigma }_h-{\epsilon }_h\left(u_h^{\ast }\right)\Vert }_0$	order
$2^{-1}$	6.6096E-01	$-$	7.7905E-01	$-$	1.6050E+01	$-$	4.3292E-02	$-$
$2^{-2}$	5.1630E-02	3.68	5.8762E-02	3.73	1.3066E+00	3.62	9.0182E-03	2.26
$2^{-3}$	3.5430E-03	3.87	4.3977E-03	3.74	1.0508E-01	3.64	6.8780E-04	3.71
$2^{-4}$	2.2220E-04	4.00	2.9277E-04	3.91	7.1542E-03	3.88	4.4330E-05	3.96
$2^{-5}$	1.3873E-05	4.00	1.8668E-05	3.97	4.5947E-04	3.96	2.7853E-06	3.99
$2^{-6}$	8.6708E-07	4.00	1.1751E-06	3.99	2.8998E-05	3.99	1.7442E-07	4.00
$2^{-7}$	5.4210E-08	4.00	7.3695E-08	4.00	1.8195E-06	3.99	1.0922E-08	4.00

Equations

$$\{\begin{array}{rlrl}(A\sigma ,\tau )+(\mathrm{div}\tau ,u)& =0& & for\; all\tau \in \Sigma ,\\ (\mathrm{div}\sigma ,v)& =(f,v)& & for\; allv\in V,\end{array}$$

$$H(\mathrm{div},\Omega ;\mathbb{S})$$

$$L^2(\Omega ;{\mathbb{R}}^2)$$

$${\eta }^2({\sigma }_h,{\mathcal{T}}_h):=\sum _{K\in {\mathcal{T}}_h}{\eta }_K^2({\sigma }_h)+\sum _{e\in {\mathcal{E}}_h}{\eta }_e^2({\sigma }_h)$$

$${\eta }_K^2({\sigma }_h):=h_K^4\Vert \mathrm{curl}\mathrm{curl}(A{\sigma }_h){\Vert }_{0,K}^2,{\eta }_e^2({\sigma }_h):=h_e\Vert {\mathcal{J}}_{e,1}{\Vert }_{0,e}^2+h_e^3\Vert {\mathcal{J}}_{e,2}{\Vert }_{0,e}^2,$$

$$\begin{array}[]{ll}\displaystyle\mathcal{J}_{e,1}:=&\left\{\begin{array}[]{lr}\displaystyle\big{[}(A\sigma_{h})t_{e}\cdot t_{e}\big{]}_{e}&\quad\quad\quad\quad\quad\quad\quad\ \ {\rm if~{}}e\in\mathcal{E}_{h}(\Omega),\\ \displaystyle\big{(}(A\sigma_{h})t_{e}\cdot t_{e}\big{)}|_{e}&~{}~{}~{}~{}{\rm if~{}}e\in\mathcal{E}_{h}(\Gamma),\end{array}\right.\\ \\ \displaystyle\mathcal{J}_{e,2}:=&\left\{\begin{array}[]{lr}\displaystyle\big{[}{\rm curl}(A\sigma_{h})\cdot t_{e}\big{]}_{e}&~{}~{}~{}~{}{\rm if~{}}e\in\mathcal{E}_{h}(\Omega),\\ \displaystyle\big{(}{\rm curl}(A\sigma_{h})\cdot t_{e}-\partial_{t_{e}}\big{(}(A\sigma_{h})t_{e}\cdot\nu_{e}\big{)}\big{)}|_{e}&~{}~{}~{}~{}{\rm if~{}}e\in\mathcal{E}_{h}(\Gamma),\end{array}\right.\end{array}$$

$${\mathrm{osc}}^2(f,{\mathcal{T}}_h):=\sum _{K\in {\mathcal{T}}_h}{h}_K^2\Vert f-Q_{h}f{\Vert }_{0,K}^2,$$

$$\Vert \sigma -{\sigma }_h{\Vert }_A\le C_1(\eta ({\sigma }_h,{\mathcal{T}}_h)+\mathrm{osc}(f,{\mathcal{T}}_h)).$$

$$C_2\eta ({\sigma }_h,{\mathcal{T}}_h)\le \Vert \sigma -{\sigma }_h{\Vert }_A$$

$${\tilde{\eta }}^2({\sigma }_h,{\mathcal{T}}_h)$$

$$+\sum _{K\in \mathcal{T}}{h}_K^2\Vert \mathrm{curl}(A{\sigma }_h+{\gamma }_h){\Vert }_{0,K}^2$$

$$+\sum _{e\in {\mathcal{E}}_h(\Omega )}{h}_e\Vert [A{\sigma }_h+{\gamma }_h{]}_e{\tau }_e{\Vert }_{0,e}^2$$

$$+\sum _{e\in {\mathcal{E}}_h({\Gamma }_D)}{h}_e\Vert (A{\sigma }_h+{\gamma }_h-\nabla u_D){\tau }_e{\Vert }_{0,e}^2.$$

$${\rm\bf{Curl}}\phi:=\left(-\partial\phi/\partial x_{2},\ \partial\phi/\partial x_{1}\right),\quad\quad\quad{\rm\bf{Curl}}v:=\left(\begin{array}[]{cc}-\partial v_{1}/\partial x_{2}&\partial v_{1}/\partial x_{1}\\ -\partial v_{2}/\partial x_{2}&\partial v_{2}/\partial x_{1}\end{array}\right).$$

$${\rm curl}\tau:=\left(\begin{array}[]{c}\partial\tau_{12}/\partial x_{1}-\partial\tau_{11}/\partial x_{2}\\ \partial\tau_{22}/\partial x_{1}-\partial\tau_{21}/\partial x_{2}\end{array}\right),\qquad{\rm div}\tau:=\left(\begin{array}[]{c}\partial\tau_{11}/\partial x_{1}+\partial\tau_{12}/\partial x_{2}\\ \partial\tau_{21}/\partial x_{1}+\partial\tau_{22}/\partial x_{2}\end{array}\right).$$

$$[w{]}_e:=(w{|}_{K_{+}}{)}_e-(w{|}_{K_{-}}{)}_e.$$

$$\left\{\begin{array}[]{rcl}(A\sigma_{h},\tau_{h})+({\rm div}\tau_{h},u_{h})&&=0\quad\quad\quad\quad\;{\rm for~{}all~{}}\tau_{h}\in\Sigma_{h},\\ ({\rm div}\sigma_{h},v_{h})&&=(f,v_{h})\quad\quad{\rm for~{}all~{}}v_{h}\in V_{h}.\end{array}\right.$$

$$P_k(K;\mathbb{S}):=\{\tau \in L^2(K;{\mathbb{R}}^{2\times 2})|{\tau }_{i,j}\in P_k(K),{\tau }_{ij}={\tau }_{ji},1\le i\le 2,1\le j\le 2\},$$

$$P_k(K;{\mathbb{R}}^2):=\{v\in L^2(K;{\mathbb{R}}^2)|v_i\in P_k(K),1\le i\le 2\},$$

$$B_{K,k}:=\left\{\tau \in P_k(K;\mathbb{S}):\tau \nu {|}_{\partial K}=0\right\}.$$

$${\Sigma }_h$$

$$V_h$$

$$B_{k,h}$$

$${\tilde{\Sigma }}_{k,h}$$

$$\Vert \sigma -{\sigma }_h{\Vert }_0$$

$$\Vert \mathrm{div}(\sigma -{\sigma }_h){\Vert }_0$$

$$\Vert u-u_h{\Vert }_0$$

$$P_1(\Omega )⟶H^2(\Omega )\stackrel{\mathbf{Curl}\mathbf{Curl}}{⟶}H(\mathrm{div},\Omega ;\mathbb{S})\stackrel{\mathrm{div}}{⟶}{L}^2(\Omega ,{\mathbb{R}}^2)$$

$$P_1(\Omega )⟶{\Phi }_h\stackrel{\mathbf{Curl}\mathbf{Curl}}{⟶}{\Sigma }_h\stackrel{\mathrm{div}}{⟶}{V}_h.$$

$$\tau =\mathcal{C}\epsilon (v)+\mathbf{Curl}\mathbf{Curl}\varphi ,$$

$$\Vert \tau {\Vert }_A^2=\Vert \epsilon (v){\Vert }_{A^{-1}}^2+\Vert \mathbf{Curl}\mathbf{Curl}\varphi {\Vert }_A^2,$$

$$(A^{-1}A\tau ,\tau )=(\tau ,\tau )=(A(A^{-1}\tau ),\tau ),$$

$$C_2\Vert \tau {\Vert }_A^2=C_2(A\tau ,\tau )\le \Vert \tau {\Vert }_0^2\le C_1(A^{-1}\tau ,\tau )=C_1\Vert \tau {\Vert }_{A^{-1}}^2.$$

$${\eta }^2({\tau }_h,{\mathcal{T}}_h):=\sum _{K\in {\mathcal{T}}_h}{\eta }_K^2({\tau }_h)+\sum _{e\in {\mathcal{E}}_h}{\eta }_e^2({\tau }_h),$$

Full text

Keywords. symmetric mixed finite element, linear elasticity problems, a posteriori error estimator, adaptive method.

AMS subject classifications. 65N30, 73C02.

The first author was supported by NSF Grant DMS-1418934. This work was finished when L. Chen visited Peking University in the fall of 2015. He would like to thank Peking University for the support and hospitality, as well as for their exciting research atmosphere.

The second author was supported by the NSFC Projects 11625101, 91430213 and 11421101.

The third author was supported by the NSFC Projects 11301396 and 11671304, Zhejiang Provincial Natural Science Foundation of China Projects LY17A010010, LY15A010015 and LY15A010016, and Wenzhou Science and Technology Plan Project G20160019.

The last author was supported by the NSFC Project 11401026. She would like to thank the support of the China Scholarship Council and the university of California, Irvine during her visit to UC Irvine from 2014 to 2015.

1. Introduction

In this paper, we are concerned with the development of residual-based a posteriori error estimators for the symmetric mixed finite element methods for planar linear elasticity problems. Let $\Omega\subset\mathbb{R}^{2}$ be a bounded polygonal domain with boundary $\Gamma:=\partial\Omega$, based on the Hellinger-Reissner principle, the linear elasticity problem with homogeneous Dirichlet boundary condition within a stress-displacement form reads: Find $(\sigma,u)\in\Sigma\times V:=H({\rm div},\Omega;\mathbb{S})\times L^{2}(\Omega;\mathbb{R}^{2})$, such that

[TABLE]

where $\mathbb{S}\subset\mathbb{R}^{2\times 2}$ is the space of symmetric matrices, and the symmetric tensor space for stress and the space for vector displacement are, respectively,

[TABLE]

Throughout the paper, the compliance tensor $A:\mathbb{S}~{}\rightarrow~{}\mathbb{S}$, characterizing the properties of the material, is bounded and symmetric positive definite. In the homogeneous isotropic case, the compliance tensor is given by $A\tau=(\tau-\lambda/(2\mu+2\lambda){\rm tr}\tau~{}{\rm I})/(2\mu)$, where $\mu>0,\lambda\geq 0$ are the Lamé constants, ${\rm I}$ is the identity matrix, ${\rm tr}\tau=\tau_{11}+\tau_{22}$ is the trace of the matrix $\tau$. For simplicity, we assume $A$ is a constant matrix in this paper and comment on the generalization to the piecewise constant matrix case.

Because of the symmetry constraint on the stress tensor, it is extremely difficult to construct stable conforming finite elements of (1.1) even for 2D problems, as stated in the plenary presentation to the 2002 International Congress of Mathematicians by Arnold [3]. An important progress in this direction is the work of Arnold and Winther [6] and Arnold, Awanou, and Winther [5]. In particular, a sufficient condition of the discrete stable method is proposed in these two papers, which states that a discrete exact sequence guarantees the stability of the mixed method. Based on such a condition, conforming mixed finite elements on the simplical and rectangular meshes are developed for both 2D and 3D [1, 4, 7, 18, 25]. Recently, based on a crucial structure of symmetric matrix valued piecewise polynomial $H(\operatorname{div})$ space and two basic algebraic results, Hu and Zhang developed a new framework to design and analyze the mixed finite element of elasticity problems. As a result, on both simplicial and tensor product grids, several families of both symmetric and optimal mixed elements with polynomial shape functions in any space dimension are constructed, see more details in [23, 24, 26, 27, 28]. Theoretical and numerical analysis show that symmetric mixed finite element method is a popular choice for a robust stress approximation [12, 14].

Computation with adaptive grid refinement has proved to be a useful and efficient tool in scientific computing over the last several decades. When the domain contains a re-entering corner, the stress has a singularity at that corner, non-uniform mesh is necessary to catch the singularity. Adaptive finite element methods based on local mesh refinement can recovery the optimal rate of convergence. The key behind this technique is to design a good a posteriori error estimator that provides a guidance on how and where grids should be refined. The residual-based a posteriori error estimators provide indicators for refining and coarsening the mesh and allow to control whether the error is below a given threshold. Various error estimators for mixed finite element discretizations of the Poisson equation have been obtained in [2, 10, 16, 19, 22, 30, 32]. Extension to the mixed finite element for linear elasticity is, however, very limited. In [11, 29, 31], the authors gave the a posteriori error estimators for the nonsymmetric mixed finite elements only.

The symmetry of the stress tensor brings essential difficulty to the a posteriori error analysis. Since only the symmetric part is approximated and not the full gradient, the approach of a posteriori error analysis developed in [11, 15, 29, 31] cannot be applied directly. In order to overcome this difficulty, Carstensen and Gedicke propose to generalize the framework of the a posteriori analysis for nonsymmetric mixed finite elements to the case of symmetric elements by decomposing the stress into the gradient and the asymmetric part of the gradient. A robust residual-based a posteriori error estimator for Arnold-Winther’s symmetric element was proposed in [13], but an arbitrary asymmetric approximation $\gamma_{h}$ of the asymmetric part of the gradient skew$(Du)=(Du-D^{T}u)/2$ was involved in this estimator. Furthermore $\gamma_{h}$ was chosen as the asymmetric gradient of a post-processed displacement to ensure the efficiency of the estimator. More details can be found below.

The goal of this paper is to present an a posteriori error estimator together with a theoretical upper and lower bounds, for the conforming and symmetric mixed finite element solutions developed in [6, 26]. We shall follow the guide principle in [6]: use the continuous and discrete linear elasticity complex, c.f. (2.2) and (2.3).

Given an approximation $\sigma_{h}$ on the triangulation $\mathcal{T}_{h}$ consisting of triangles, we construct the following a posteriori error estimator, denoted by $\eta$,

[TABLE]

where

[TABLE]

[TABLE]

with $\mathcal{E}_{h}$ being the collection of all edges of $\mathcal{T}_{h}$. We write $\mathcal{E}_{h}=\mathcal{E}_{h}(\Omega)\bigcup\mathcal{E}_{h}(\Gamma)$, where $\mathcal{E}_{h}(\Omega)$ is the collection of interior edges and $\mathcal{E}_{h}(\Gamma)$ is the collection of all element edges on the boundary. For any edge $e\in\mathcal{E}_{h}$, let $t_{e}=(-n_{2},n_{1})^{T}$ be the unit tangential vector along edge $e$ for the unit outward normal $\nu_{e}=(n_{1},n_{2})^{T}$. Let $h_{K}$ be the diameter of the element $K$ and $h_{e}$ be the length of edge $e$. The data oscillation is defined as

[TABLE]

where $Q_{h}$ is the $L^{2}$ orthogonal projection operator onto the discrete displacement space.

Using the Helmholtz decomposition induced from the linear elasticity complex [6, 11], we establish the following reliability

[TABLE]

In addition, we will prove the following efficiency estimate

[TABLE]

by following the approach from  [2].

We also generalize the above results to the mixed boundary problems, for which the error estimator is modified on the Dirichlet boundary edges. Reliability and efficiency of the modified error estimator can be proved similarly.

In [17], a superconvergent approximate displacement $u_{h}^{*}$ was constructed by a postprocessing of $(\sigma_{h},u_{h})$ . Using this result and the a posteriori error estimation of the stress, we also give the a posteriori error estimation for the displacement $\|u-u_{h}^{*}\|_{1,h}$ in a mesh dependent norm.

In order to compare with the a posteriori error estimator in [13], we present their estimator as follows:

[TABLE]

(The estimator is rewritten in our notation and the details of the standard notation can be found below.) To ensure the efficiency of the estimator, a sufficiently accurate polynomial asymmetric approximation $\gamma_{h}$ of the asymmetric gradient ${\rm skew}(Du):=(Du-D^{T}u)/2$ is involved in the above estimator. Since the global approximation or even minimization may be too costly, Carstensen and Gedicke compute the sufficiently accurate approximation $\gamma_{h}={\rm skew}(Du_{h}^{*})$ by the post-processed displacement $u_{h}^{*}$ in the spirit of Stenberg [35]. As we can see, this estimator is totally different to ours. The estimators we propose use the symmetric stress directly and do not need any estimation of the asymmetric part. Therefore it is more computationally efficient.

The remaining parts of the paper is organized as follows. Section 2 presents the notations and the discrete finite element problems. Section 3 proposes an a posteriori error estimator for the stress and proves the reliability and efficiency of the estimator. Section 4 generalizes the results of section 3 to mixed boundary problems. Section 5 gives a posteriori error estimation for the displacement. Section 6 presents numerical experiments to show the effectiveness of the estimator. Throughout this paper, we use “$\lesssim\cdots$” to mean that “$\leq C\cdots$”, where $C$ is a generic positive constant independent of $h$ and the Lamé constant $\lambda$, which may take different values at different appearances.

2. Notations and preliminaries

Standard notations on Sobolev spaces and norms are adopted throughout this paper and, for brevity, $\|\cdot\|:=\|\cdot\|_{L^{2}(\Omega)}$ denotes the $L^{2}$ norm. $(\cdot,\cdot)_{K}$ represents, as usual, the $L^{2}$ inner product on the domain $K$, the subscript $K$ is omitted when $K=\Omega$. $\langle\cdot,\cdot\rangle_{\Gamma}$ represents the $L^{2}$ inner product on the boundary $\Gamma$. For brevity, let $\displaystyle\partial_{x_{i}}:=\partial/\partial x_{i}$ and $\partial^{2}_{x_{i}x_{j}}:=\partial^{2}/\partial x_{i}\partial x_{j},j=1,2,$ $\partial_{\nu}:=\partial/\partial\nu$, $\partial_{t}:=\partial/\partial t$. For $\phi\in H^{1}(\Omega;\mathbb{R})$, $v=(v_{1},v_{2})^{T}\in H^{1}(\Omega;\mathbb{R}^{2})$, set

[TABLE]

For $\tau=(\tau_{i,j})_{2\times 2}\in H^{1}(\Omega;\mathbb{R}^{2\times 2})$, set

[TABLE]

Namely the differential operators $\operatorname{curl}$ and $\operatorname{div}$ are applied rowwise for tensors.

Let $\mathcal{T}_{h}$ be a shape-regular triangulation of $\bar{\Omega}$ into triangles with the set of edges $\mathcal{E}_{h}$. Denote by $\mathcal{E}_{h}(\Omega)$ the collection of all interior element edges in $\mathcal{T}_{h}$ and $\mathcal{E}_{h}(\Gamma)$ the collection of all element edges on the boundary. For any triangle $K\in\mathcal{T}_{h}$, let $\mathcal{E}(K)$ be the set of its edges. For any edge $e\in\mathcal{E}(K)$, let $t_{e}=(-n_{2},n_{1})^{T}$ be the unit tangential vector along edge $e$ for the unit outward normal vector $\nu_{e}=(n_{1},n_{2})^{T}$, $h_{K}$ be the diameter of the element $K$ and $h_{e}$ be the length of the edge $e$, $h=\max\limits_{K\in\mathcal{T}_{h}}\{h_{K}\}$ be the diameter of the partition $\mathcal{T}_{h}$. The jump $[w]_{e}$ of $w$ across edge $e=\bar{K}_{+}\cap\bar{K}_{-}$ reads

[TABLE]

Particularly, if $e\in\mathcal{E}_{h}(\Gamma),\ [w]_{e}:=w|_{e}$.

Let $\Sigma_{h}\times V_{h}\subseteq\Sigma\times V$ be a symmetric conforming mixed element defined on the mesh $\mathcal{T}_{h}$, then the discrete mixed formulation for the problem (1.1) is: find $(\sigma_{h},u_{h})\in\Sigma_{h}\times V_{h}$, such that

[TABLE]

In the sequel, we briefly introduce Hu-Zhang element  [23, 26, 28]. For each $K\in\mathcal{T}_{h}$, let $P_{k}(K)$ be the space of polynomials of total degree at most $k$ on $K$ and

[TABLE]

[TABLE]

define an $H(\operatorname{div},K;\mathbb{S})$ bubble function as

[TABLE]

The Hu-Zhang element space is given by

[TABLE]

with integer $k\geq 3$, where

[TABLE]

For the above elements, the following *a priori * error estimate holds.

Theorem 2.1** (A priori error estimate [23, 26, 28]).**

The exact solution $(\sigma,u)$ of problem (1.1) and the approximate solution $(\sigma_{h},u_{h})$ of problem (2.1) satisfy

[TABLE]

In the continuous case, the following exact sequence

[TABLE]

holds for linear elasticity [6]. In the discrete case, the exact sequence holds similarly

[TABLE]

As stated in [6], the space $\Phi_{h}$ for the Arnold-Winther element is precisely the space of $C^{1}$ piecewise polynomials which are $C^{2}$ at the vertices, that is, the well-known high-order Hermite or Argyris finite element. The Hu-Zhang element is an enrichment of the Arnold-Winther element, adding all the piecewise polynomial matrices of degree $k$ which are not divergence-free on each element and belong to $H(\operatorname{div},\Omega;\mathbb{S})$ globally. So the space $\Phi_{h}$ for the Hu-Zhang element is the same as the one for the Arnold-Winther element.

Lemma 2.1** (Helmholtz-type decomposition [6, 11]).**

For any $\tau\in L^{2}(\Omega;\mathbb{S})$, there exists $v\in H_{0}^{1}(\Omega;\mathbb{R}^{2})$ and $\phi\in H^{2}(\Omega)/P_{1}(\Omega)$, such that

[TABLE]

and the decomposition is orthogonal in the weighted $L^{2}$-inner product $(\mathcal{C}^{-1}\cdot,\cdot):=(A~{}\cdot,\cdot)$, i.e.,

[TABLE]

where $P_{1}(\Omega)$ is the linear polynomial space on $\Omega$, the norm $\|\cdot\|_{A}=(A~{}\cdot,\cdot)$.

Since

[TABLE]

by the boundedness and coerciveness of the operator $A$, we obtain the following relationship of the norms: for any $\tau\in\Sigma$, there exist positive constants $C_{1}$ and $C_{2}$, which are independent of the Lamé constant $\lambda$, such that

[TABLE]

It is the goal of this paper to present a posterior error estimate of $\sigma-\sigma_{h}$ for the Hu-Zhang element method. It is worth mentioning that the a posterior error estimator designed in this paper can be easily extended to the Arnold-Winther element [6].

3. A posteriori Error Estimation for Stress

In this section, we shall prove the reliability and efficiency of the error estimator. The main observation is that: although it is a saddle point problem, the error of stress $\sigma-\sigma_{h}$ is orthogonal to the divergence-free subspace, while the part of the error that is not divergence- free can be bounded by the data oscillation using the stability of the discretization.

For any $\tau_{h}\in\Sigma_{h}$, the error estimator is defined as

[TABLE]

where

[TABLE]

[TABLE]

The data oscillation is defined as

[TABLE]

where $Q_{h}$ is the $L^{2}$ orthogonal projection operator onto the discrete displacement space $V_{h}$.

3.1. Stability result

For the easy of exposition, we write the mixed formulation for linear elasticity as $\mathcal{L}(\sigma,u)=f$. The natural stability of the operator is $\|\sigma\|_{H(\operatorname{div})}+\|u\|\lesssim\|f\|.$ However, a stronger stability can be proved for a special perturbation of the data.

Lemma 3.1**.**

Let $f_{h}$ be the $L^{2}$ projection of $f$ onto $V_{h}$ and let $(\sigma,u)=\mathcal{L}^{-1}f$ and $(\tilde{\sigma},\tilde{u})=\mathcal{L}^{-1}f_{h}$. Then we have

[TABLE]

Proof.

Use the first equation of (1.1) and let $v=u-\tilde{u}$

[TABLE]

where the Korn’s inequality is used and the symmetric gradient $\varepsilon(v)=\frac{1}{2}(\nabla v+(\nabla v)^{T})$. Since $\varepsilon(v)=A(\sigma-\tilde{\sigma})$, by (2.6), $\|\varepsilon(v)\|_{0}\lesssim\|\sigma-\tilde{\sigma}\|_{A}$. We acquire the desirable stability result. ∎

The oscillation ${\rm osc}(f,\mathcal{T}_{h})$ is an upper bound of $\|f-f_{h}\|_{-1}$ and is of high order comparing with the error estimator.

3.2. Orthogonality

For any $\phi\in H^{2}(\Omega)$, ${\rm\bf{Curl}\,\bf{Curl}\,}\phi\in H({\rm div},\Omega;\mathbb{S})$, we can use the exact sequence property $\operatorname{div}{\rm\bf{Curl}\,\bf{Curl}\,}=0$ to get

[TABLE]

Similarly

[TABLE]

for any $\phi_{h}\in\Phi_{h}$. Therefore we have a partial orthogonality

[TABLE]

3.3. Upper bound

Let $S_{h}^{5}$ denote the Argyris finite element space, which consists of $C^{1}$ piecewise polynomials of degree less than or equal to $5$

[TABLE]

Following [34, 20], we can define a quasi-interpolation operator $I_{h}:~{}H^{2}(\Omega)~{}\rightarrow S_{h}^{5}$, which preserves the values of the function on all vertices of $\mathcal{T}_{h}$. On each element $K\in\mathcal{T}_{h}$, for any $v\in H^{2}(\Omega)$, $I_{h}v|_{K}\in P_{5}(K)$ and it satisfies

• •

$I_{h}v|_{K}(a_{i,K})=v(a_{i,K}),\quad\quad 1\leq i\leq 3;$

• •

$\partial_{x_{j}}(I_{h}v|_{K})(a_{i,K})=\frac{1}{\mathcal{N}_{h}(a_{i,K})}\sum\limits_{K^{\prime}\in S(a_{i,K})}\partial_{x_{j}}(P_{h}v|_{K^{\prime}})(a_{i,K}),\ 1\leq i\leq 3,\\ \ j=1,2;$

• •

$\partial^{2}_{x_{j}x_{l}}(I_{h}v|_{K})(a_{i,K})=\frac{1}{\mathcal{N}_{h}(a_{i,K})}\sum\limits_{K^{\prime}\in S(a_{i,K})}\partial^{2}_{x_{j}x_{l}}(P_{h}v|_{K^{\prime}})(a_{i,K}),\ 1\leq i\leq 3,\ 1\leq j\leq l\leq 2;$

• •

$\partial_{\nu}(I_{h}v|_{K})(a_{3+i,K})=\frac{1}{\mathcal{N}_{h}(a_{3+i,K})}\sum\limits_{K^{\prime}\in S(a_{3+i,K})}\partial_{\nu}(P_{h}v|_{K^{\prime}})(a_{3+i,K}),\ 1\leq i\leq 3;$

where $a_{i,K},1\leq i\leq 3,$ are the vertices of $K$, $a_{3+i,K},1\leq i\leq 3,$ are the edge midpoints of $K$, $\nu$ is the edge outer normal of the element $K$ on the edge midpoint, $S(a_{i,K})=\bigcup\{K\in\mathcal{T}_{h}:a_{i,K}\in K\}$ and $\mathcal{N}_{h}(a_{i,K})=\text{card}\{K:K\in S(a_{i,K})\}$, $P_{h}$ is the projection operator from $L^{2}(\Omega)$ onto the piecewise linear polynomial finite element space on $\mathcal{T}_{h}$. It is obvious that the interpolation operator $I_{h}$ is uniquely determined by the above degrees of freedom. Furthermore, $I_{h}$ is a projection, i.e.

[TABLE]

and it preserves the value of the function on vertices for any $v\in H^{2}(\Omega)$, i.e.

[TABLE]

A similar scaling argument as in [34, 20] gives the following interpolation estimates

[TABLE]

[TABLE]

where $S_{K}=\bigcup\{K_{i}\in\mathcal{T}_{h}:K_{i}\bigcap\bar{K}\neq\varnothing\}$, $S_{e}=\bigcup\{K_{i}\in\mathcal{T}_{h}:K_{i}\bigcap e\neq\varnothing\}.$

Applying the Helmholtz decomposition to the error $\tilde{\sigma}-\sigma_{h}$, we have

[TABLE]

and

[TABLE]

where $v\in H_{0}^{1}(\Omega;\mathbb{R}^{2})$ and $\phi\in H^{2}(\Omega)/P_{1}(\Omega)$. By this orthogonal decomposition and the fact $\operatorname{div}(\tilde{\sigma}-\sigma_{h})=0$,

[TABLE]

Since ${\rm\bf{Curl}\,\bf{Curl}\,}(I_{h}\phi)\in\Sigma_{h}$, by the orthogonality (3.4) and the equation (3.3),

[TABLE]

An integration by parts gives

[TABLE]

The second term of the right hand side can be rewritten as

[TABLE]

Since the compliance tensor $A$ is symmetric and continuous, $(A\sigma_{h}t)\cdot\nu=(A\sigma_{h}\nu)\cdot t$ and $(A\sigma_{h}t)\cdot\nu$ is continuous across the interior element edge. This implies

[TABLE]

where the fact $(\phi-I_{h}\phi)$ vanishing at the boundary vertices (3.6) is used. So

[TABLE]

Substituting it into (3.11), we get

[TABLE]

Then applying the Cauchy-Schwarz inequality, the error estimate of the quasi-interpolation (3.7), (3.8), we have

[TABLE]

By [11], the $\phi$ defined in (3.9) satisfies that ${\rm div}(\bf{Curl}~{}\bf{Curl}~{}\phi)=0$ and

[TABLE]

Using Proposition 9.1.1 in [8], we get

[TABLE]

where the constant $C$ is independent of the Lamé constant $\lambda$. Combining this with (3.10), (3.3), we obtain

[TABLE]

Together with the triangle inequality and perturbation result (3.2), we get the desired error bound

[TABLE]

In summary, we obtain the following upper bound estimation.

Theorem 3.1** (Reliability of the error estimator).**

Let $(\sigma,u)$ be the solution of the mixed formulation (1.1) and $(\sigma_{h},u_{h})$ be the solution of the mixed finite element method (2.1). If the compliance tensor $A$ is continuous, there exists positive constant $C_{1}$ depending only on the shape-regularity of the triangulation and the polynomial degree $k$ such that

[TABLE]

Remark 3.1**.**

When $A$ is discontinuous, we can modify $\eta(\sigma_{h},\mathcal{T}_{h})$ as follows:

[TABLE]

Compared to the case of continuous coefficient $A$, this estimator includes an additional term, the jump of $\partial_{t_{e}}((A\sigma_{h})t_{e}\cdot\nu_{e})$ on all interior edges, owing to the discontinuity of the matrix $A$. Similarly, we can prove the reliability of the estimator

[TABLE]

Remark 3.2**.**

By Proposition 9.1.1 in [8], it holds

[TABLE]

where $\hat{\Sigma}:=\{\tau\in\Sigma:(\mathrm{tr}\tau,1)=0\}$ with $\mathrm{tr}$ being the trace operator of matrix. Then we also have from (3.13) and the fact that $\|f-f_{h}\|_{-1}\lesssim{\rm osc}(f,\mathcal{T}_{h})$

[TABLE]

That is we can control the $L^{2}$ norm of the stress with constant independent of the Lamé constant $\lambda$. $\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\;\Box$

3.4. Lower bound

We shall follow Alonso [2] to prove the efficiency of the error estimator defined in (3.1). Similar to [2], we need the following lemma.

Lemma 3.2**.**

For any $K\in\mathcal{T}_{h}$, given $p_{K}\in L^{2}(K),~{}q_{e}\in L^{2}(e),~{}r_{e}\in L^{2}(e),~{}e\in\partial K,$ there exists a unique $\psi_{K}\in P_{k+4}(K)$ satisfying that

[TABLE]

where $P_{k}(e)$ denotes the spaces of polynomial of degree less than or equal to $k$ on edge $e$. Moreover it holds that

[TABLE]

Proof.

Similar as in [33], such a function $\psi_{K}$ is determined uniquely by the above degrees of freedoms. A standard homogeneity argument gives (3.15). ∎

Theorem 3.2** (Efficiency of the error estimator).**

Let $(\sigma,u)$ be the solution of the mixed formulation (1.1) and $(\sigma_{h},u_{h})$ be the solution of the mixed finite element method (2.1). If the compliance tensor $A$ is continuous, there exists positive constant $C_{2}$ depending only on the shape-regularity of the triangulations and the polynomial degree $k$ such that

[TABLE]

Proof.

The estimator $\eta^{2}(\sigma_{h},{\mathcal{T}_{h}})$ can be rewritten as

[TABLE]

On each element $K\in\mathcal{T}_{h}$, we apply Lemma 3.2 for $p_{K}=h_{K}^{4}{\rm curl\,curl\,}(A\sigma_{h})|_{K}$, $q_{e}=-h_{e}^{3}\mathcal{J}_{e,2}$, $r_{e}=h_{e}\mathcal{J}_{e,1}$ for each edge $e\in\partial K$. Let $\psi|_{K}=\psi_{K}$, such a defined $\psi$ is in the high-order Argyris finite element space of degree $k+4$, hence $\psi\in H^{2}(\Omega)$. Using (3.15), it follows that

[TABLE]

This, in conjunction with (3.14), yields

[TABLE]

Since $(A\sigma_{h})t_{e}\cdot\nu_{e}$ is continuous across the element edge $e$, $[A\sigma_{h}t_{e}\cdot\nu_{e}]_{e}=0$ on interior edges. Noting that $\psi\in H^{2}(\Omega)$ and vanishes at the mesh vertices,

[TABLE]

Hence the last two terms of (3.18) become

[TABLE]

Substituting (3.20) into (3.18) leads to

[TABLE]

Integrating the first term by parts twice,

[TABLE]

where ${\rm\bf{Curl}\,\bf{Curl}}~{}\psi\in\Sigma$ and the inverse inequality are used. By (3.17),

[TABLE]

Combining the above two inequalities, we have that

[TABLE]

∎

Remark 3.3**.**

For discontinuous $A$ and the modified error estimator in Remark 3.1, efficiency can be also proved using a similar argument.

4. A posteriori error estimation for mixed boundary problems

The a posteriori error estimation for the linear elasticity problems with the homogeneous Dirichlet boundary condition can be generalized to problems with mixed boundary conditions. In this section, we will discuss the following linear elasticity problems with mixed boundary conditions. Let $\Omega\subset\mathbb{R}^{2}$ be a bounded polygonal domain with boundary $\Gamma:=\partial\Omega=\Gamma_{D}\cup\Gamma_{N}$, $\Gamma_{D}\cap\Gamma_{N}=\emptyset$, $\Gamma_{N}\neq\emptyset$. Given data $f\in L^{2}(\Omega;\mathbb{R}^{2})$, $u_{D}\in H^{1}(\Omega;\mathbb{R}^{2})$, and $g\in L^{2}(\Gamma_{N};\mathbb{R}^{2})$, seek the solution $(\sigma,u)\in\Sigma_{g}\times V$, such that

[TABLE]

where

[TABLE]

where $\mathcal{D}$ denotes the space of test functions. Let $\Sigma_{0,h}:=\Sigma_{0}\bigcap\Sigma_{h},~{}~{}~{}~{}\Sigma_{g,h}:=\Sigma_{g}\bigcap\Sigma_{h},$ the mixed finite element method seeks $(\sigma_{h},u_{h})\in\Sigma_{g,h}\times V_{h}$, such that

[TABLE]

We modify the a posterior error estimator defined in Section 3 as the following:

[TABLE]

where

[TABLE]

[TABLE]

where $\mathcal{E}_{h}(\Gamma_{D})$, $\mathcal{E}_{h}(\Gamma_{N})$ are the collection of element edges for Dirichlet boundary and Neumann boundary respectively.

Similar to Section 3, we can prove the reliability and efficiency of this a posteriori error estimator.

Theorem 4.1** (Reliability and efficiency of the error estimator).**

Let $(\sigma,u)$ be the solution of the mixed formulation (4.1) and $(\sigma_{h},u_{h})$ be the solution of the mixed finite element method (4.2). If the compliance tensor $A$ is continuous, there exist positive constant $C_{3}$ and $C_{4}$ depending only on the shape-regularity of the triangulation and the polynomial degree $k$ such that

[TABLE]

and

[TABLE]

where the data oscillations for the Dirichlet boundary $u_{D}$ and the Neumann boundary condition $g$ are defined as

[TABLE]

$g_{h}$* is the piecewise $L^{2}$ projection of $g$ onto $P_{k}(\mathcal{E}_{h}(\Gamma_{N}),\mathbb{R}^{2})$ and $u_{D,h}$ is the piecewise $L^{2}$ projection of $u_{D}$ onto $P_{k}(\mathcal{E}_{h}(\Gamma_{D}),\mathbb{R}^{2})$.*

5. A Posteriori Error Estimation for Displacement

In this section, we shall discuss the a posteriori error estimate for a superconvergent postprocessed displacement recently constructed in [17]. The key points of the theoretical analysis involve the discrete inf-sup condition and the norm equivalence on $H^{1}(\mathcal{T}_{h};\mathbb{R}^{2})$ developed in [17], and the a posteriori error estimates (3.13) and (3.16). Here the broken space

[TABLE]

For any $v\in H^{1}(\mathcal{T}_{h};\mathbb{R}^{2})$, define mesh dependent norm

[TABLE]

We first recall the superconvergent postprocessed displacement from $(\sigma_{h},u_{h})$ developed in [17]. To this end, let

[TABLE]

Then a postprocessed displacement is defined as follows [17, 32, 9]: Find $u_{h}^{\ast}\in V_{h}^{\ast}$ such that

[TABLE]

[TABLE]

for any $K\in\mathcal{T}_{h}$.

We recall the following two useful results [17]: the discrete inf-sup condition

[TABLE]

and norm equivalence

[TABLE]

Theorem 5.1**.**

Let $(\sigma,u)$ be the solution of the mixed formulation (1.1), $(\sigma_{h},u_{h})$ be the solution of the mixed finite element method (2.1), and $u_{h}^{\ast}$ be the postprocessed displacement defined by (5.1)-(5.2). Then we have

[TABLE]

[TABLE]

Proof.

Using the discrete inf-sup condition (5.3) with $v_{h}=Q_{h}(u-u_{h}^{\ast})$, (5.1), the first equations of (1.1) and (2.1), we get

[TABLE]

Choosing $v=u-u_{h}^{\ast}$ in (5.4),

[TABLE]

Then it follows from the last two inequalities that

[TABLE]

which combined with (3.13) implies (5.5).

Next we prove the efficiency (5.6). By the triangle inequality,

[TABLE]

Therefore we can end the proof by using (3.16). ∎

6. Numerical experiments

We will testify the a posteriori error estimator by some numerical examples in this section.

In the first example, let $\Omega=(0,1)^{2}$, $k=3$, $\mu=1$, the right-hand side

[TABLE]

and the exact solution [13, Section 5.2]

[TABLE]

We subdivide $\Omega$ by a uniform triangular mesh. The a priori and a posteriori error estimates for $\lambda=10$ and $\lambda=10000$ are listed in Tables 6.1-6.2, from which we can see that the convergence rates of $\|\sigma-\sigma_{h}\|_{A}$, $\|\nabla_{h}(u-u_{h}^{\ast})\|_{0}$, $\eta(\sigma_{h},\mathcal{T}_{h})$ and $\|A\sigma_{h}-\varepsilon_{h}(u_{h}^{\ast})\|_{0}$ are all $O(h^{4})$. Hence the a posteriori error estimators $\eta(\sigma_{h},\mathcal{T}_{h})$ and $\eta(\sigma_{h},\mathcal{T}_{h})+\|A\sigma_{h}-\varepsilon_{h}(u_{h}^{\ast})\|_{0}$ are both uniformly reliable and efficient with respect to the mesh size $h$ and $\lambda$ for smooth solutions.

Next we use the a posteriori error estimator $\eta(\sigma_{h},\mathcal{T}_{h})$ to design an adaptive mixed finite element method, i.e. Algorithm 1. The approximate block factorization preconditioner with GMRES [17] is adopted in the SOLVE part of Algorithm 1, which is verified to be highly efficient and robust even on adaptive meshes by our numerical examples.

Now we construct a problem with singularity in the solution to test Algorithm 1. Set L-shaped domain $\Omega=(-1,1)\times(-1,1)\backslash[0,1)\times(-1,0]$. Let

[TABLE]

[TABLE]

[TABLE]

where $z\in(0,1)$ is a real root of $(\lambda+3\mu)^{2}\sin^{2}(z\omega)=(\lambda+\mu)^{2}z^{2}\sin^{2}\omega$ with $\omega=3\pi/2$. The exact singular solution in polar coordinates is taken as [21, Section 4.6]

[TABLE]

It can be computed that $z=0.561586549334359$ for $\lambda=10$, and $z=0.544505718203590$ for $\lambda=10000$. We also take $k=3$ and $\mu=1$.

Some meshes generated by Algorithm 1 for different bulk parameter $\vartheta$ and Lamé constant $\lambda$ are shown in Figure 6.1, where $\#\textrm{dofs}$ is the number of degrees of freedom. The adaptive Algorithm 1 captures the singularity of the exact solution on the corner $(0,0)$ very well. The histories of the adaptive Algorithm 1 for $\vartheta=0.1,0.2$ and $\lambda=10,10000$ are presented in Figures 6.2-6.3. We can see from Figures 6.2-6.3 that the convergence rates of errors $\|\sigma-\sigma_{h}\|_{A}$ and $\eta(\sigma_{h},\mathcal{T}_{h})$ are both $O((\#\textrm{dofs})^{-2})$ no matter $\lambda=10$ or $\lambda=10000$, which demonstrates the theoretical results. For uniform grid, $\#(\textrm{dofs})^{-2}\cong h^{4}$, this means that the errors $\|\sigma-\sigma_{h}\|_{A}$ and $\eta(\sigma_{h},\mathcal{T}_{h})$ converge with an optimal rate.

The third example considers the L-shape benchmark problem with general boundary conditions testified in [13, section 5.3] on the rotated L-shaped domain with the initial mesh as depicted in Figure 6.4. We impose the Neumann boundary condition on the boundary $x^{2}=y^{2}$ and the Dirichlet boundary condition on the rest boundary of $\Omega$. The exact solution in the polar coordinates is given as follows

[TABLE]

The constants are $C_{1}:=-\cos((\alpha+1)\omega)/\cos((\alpha-1)\omega)$ and $C_{2}:=-2(\lambda+2\mu)/(\lambda+\mu)$, where $\alpha=0.544483736782$ is the positive solution of $\alpha\sin(2\omega)+\sin(2\omega\alpha)=0$ for $\omega=3\pi/4$. The Lamé parameters

[TABLE]

with the elasticity modulus $E=10^{5}$ and the Poisson ratio $\nu=0.4999$. The volume force $f(x,y)$ and the Neumann boundary data vanish, and the Dirichlet boundary condition is taken from the exact solution. The histories of Algorithm 1 for $k=3,4,5$ and $\vartheta=0.1$ are presented in Figures 6.5-6.6, which indicate that the convergence rates of errors $\|\sigma-\sigma_{h}\|_{A}$ and $\eta(\sigma_{h},\mathcal{T}_{h})$ are both $O((\#\textrm{dofs})^{-(k+1)/2})$.