A posteriori error estimation for planar linear elasticity by stress reconstruction

TL;DR

This paper extends a posteriori error estimation techniques for linear elasticity, ensuring guaranteed reliability and efficiency, including in the incompressible limit, through stress reconstruction and conforming approximation modifications.

Contribution

It introduces modifications to stress reconstruction and conforming approximation methods for reliable a posteriori error estimation in linear elasticity, applicable even in incompressible cases.

Findings

Guaranteed error bounds with known constants achieved.

Local efficiency of the estimator demonstrated.

Effective adaptive computations across Lamé parameters.

Abstract

The nonconforming triangular piecewise quadratic finite element space by Fortin and Soulie can be used for the displacement approximation and its combination with discontinuous piecewise linear pressure elements is known to constitute a stable combination for incompressible linear elasticity computations. In this contribution, we extend the stress reconstruction procedure and resulting guaranteed a posteriori error estimator developed by Ainsworth, Allendes, Barrenechea and Rankin \cite{AinAllBarRan:12} and by Kim \cite{Kim:12a} to linear elasticity. In order to get a guaranteed reliability bound with respect to the energy norm involving only known constants, two modifications are carried out: (i) the stress reconstruction in next-to-lowest order Raviart-Thomas spaces is modified in such way that its anti-symmetric part vanishes in average on each element; (ii) the auxiliary conforming approximation is constructed under the constraint that its divergence coincides with the one for the nonconforming approximation. An important aspect of our construction is that all results hold uniformly in the incompressible limit. Local efficiency is also shown and the effectiveness is illustrated by adaptive computations involving different Lam\'e parameters including the incompressible limit case.

Equations

$$\begin{array}{rl}{V}_h& =\{v_h\in L^2(\Omega ):{v_h|}_T\in P_2(T)\text{mbox}forallT\in {\mathcal{T}}_h,\\ & ⟨[[v_h{]]}_E,z{⟩}_{L^2(E)}=0\text{mbox}forallz\in P_1(E),E\in {\mathcal{E}}_h\cap \Omega ,\\ & ⟨v_h,z{⟩}_{L^2(E)}=0\text{mbox}forallz\in P_1(E),E\in {\mathcal{E}}_h\cap {\Gamma }_D\},\end{array}$$

$$\begin{split}2\mu\>(\mbox{\boldmath$\varepsilon$}({\bf u}),\mbox{\boldmath$\varepsilon$}({\bf v}))+(p,{\rm div}\>{\bf v})&=({\bf f},{\bf v})+\langle{\bf g},{\bf v}\rangle_{L^{2}(\Gamma_{N})}\>,\\ ({\rm div}\>{\bf u},q)-\frac{1}{\lambda}(p,q)&=0\end{split}$$

$$\begin{split}2\mu(\mbox{\boldmath$\varepsilon$}({\bf u}_{h}),\mbox{\boldmath$\varepsilon$}({\bf v}_{h}))_{h}+(p_{h},{\rm div}\>{\bf v}_{h})_{h}&=({\cal P}_{h}{\bf f},{\bf v}_{h})+\langle{\cal P}_{h}^{\Gamma_{N}}{\bf g},{\bf v}_{h}\rangle_{\Gamma_{N}}\\ ({\rm div}\>{\bf u}_{h},q_{h})_{h}-\frac{1}{\lambda}(p_{h},q_{h})&=0\end{split}$$

$$\begin{array}{rl}{\mathbf{B}}_h^{NC}& =\{{\mathbf{b}}_h^{NC}\in L^2(\Omega {)}^2:{{\mathbf{b}}_h^{NC}|}_T\in P_2(T{)}^2\text{mbox}forallT\in {\mathcal{T}}_h,\\ & ⟨{\mathbf{b}}_h^{NC},\mathbf{z}{⟩}_{L^2(E)}=0\text{mbox}forall\mathbf{z}\in P_1(E{)}^2,E\in {\mathcal{E}}_h\},\end{array}$$

$$\|\mbox{\boldmath$\nabla$}{\bf v}\|_{h}\leq C_{K}\|\mbox{\boldmath$\varepsilon$}({\bf v})\|_{h}\mbox{ for all }{\bf v}\in H_{\Gamma_{D}}^{1}(\Omega)^{2}+{\bf V}_{h}\>,$$

$$\mbox{\boldmath$\sigma$}_{h}^{R}\cdot{\bf n}=\{\!\{\mbox{\boldmath$\sigma$}_{h}\cdot{\bf n}\}\!\}_{E}{\color[rgb]{0,0,1}=\left\{\begin{array}[]{lcl}{\cal P}_{h}^{\Gamma_{N}}{\bf g}&,&E\subset\Gamma_{N}\>,\\ \left.\mbox{\boldmath$\sigma$}_{h}\cdot{\bf n}\right|_{T-}&,&E\subset\Gamma_{D}\>,\\ \left(\left.\mbox{\boldmath$\sigma$}_{h}\cdot{\bf n}\right|_{T-}+\left.\mbox{\boldmath$\sigma$}_{h}\cdot{\bf n}\right|_{T+}\right)/2&,&\mbox{ otherwise}\>,\end{array}\right.}$$

$$({\rm div}\>\mbox{\boldmath$\sigma$}_{h}^{R}+{\cal P}_{h}{\bf f},\widehat{{\bf q}}_{h})_{L^{2}(T)}=0$$

$$({\rm div}\>\mbox{\boldmath$\sigma$}_{h}^{R}+{\cal P}_{h}{\bf f},1)_{L^{2}(T)}=0$$

$$\begin{split}\mbox{\boldmath$\Sigma$}_{h}^{R,1}&=\{\mbox{\boldmath$\sigma$}_{h}\in H_{\Gamma_{N}}({\rm div},\Omega)^{2}:\left.\mbox{\boldmath$\sigma$}_{h}\right|_{T}\in P_{1}(T)^{2\times 2}\mbox{ for all }T\in{\cal T}_{h}\>,\\ &\;\;\;\;\;\;\;\langle{\bf n}_{E}\cdot(\mbox{\boldmath$\sigma$}_{h}\cdot{\bf n}_{E}),1\rangle_{E}=\langle{\bf t}_{E}\cdot(\mbox{\boldmath$\sigma$}_{h}\cdot{\bf n}_{E}),1\rangle_{E}=0\mbox{ for all }E\in{\cal E}_{h}\}\>,\\ \mbox{\boldmath$\Sigma$}_{h}^{R,\Delta}&=\{\mbox{\boldmath$\sigma$}_{h}\in{\color[rgb]{0,0,1}\mbox{\boldmath$\Sigma$}_{h}^{R}}:\left.{\rm div}\>\mbox{\boldmath$\sigma$}_{h}\right|_{T}\in P_{1}(T)^{2}\mbox{ for all }T\in{\cal T}_{h}\>,\\ &\;\;\;\;\;\;\;\mbox{\boldmath$\sigma$}_{h}\cdot{\bf n}_{E}={\bf 0}\mbox{ on }E\mbox{ for all }E\in{\cal E}_{h}\}\>.\end{split}$$

$$\|\mbox{\boldmath$\tau$}\|\leq C_{A}\left(\|{\bf dev}\>\mbox{\boldmath$\tau$}\|+\|{\rm div}\>\mbox{\boldmath$\tau$}\|\right)\mbox{ for all }\mbox{\boldmath$\tau$}\in H_{\Gamma_{N}}({\rm div},\Omega)^{2}\>,$$

$$\begin{split}|||({\bf u}-{\bf u}_{h},p-p_{h})|||&=(2\mu\mbox{\boldmath$\varepsilon$}({\bf u}-{\bf u}_{h})+(p-p_{h}){\bf I},\mbox{\boldmath$\varepsilon$}({\bf u}-{\bf u}_{h}))_{h}^{1/2}\\ &=\left(2\mu\|\mbox{\boldmath$\varepsilon$}({\bf u}-{\bf u}_{h})\|_{h}^{2}+\frac{1}{\lambda}\|p-p_{h}\|^{2}\right)^{1/2}\>.\end{split}$$

$${\rm tr}\>\mbox{\boldmath$\sigma$}=2\mu{\rm div}\>{\bf u}+2p=2(\mu+\lambda){\rm div}\>{\bf u}\mbox{ and }{\rm tr}\>\mbox{\boldmath$\sigma$}_{h}=2\mu{\rm div}\>{\bf u}_{h}+2p_{h}=2(\mu+\lambda){\rm div}\>{\bf u}_{h}\>,$$

$$\mbox{\boldmath$\varepsilon$}({\bf u})=\frac{1}{2\mu}\left(\mbox{\boldmath$\sigma$}-p{\bf I}\right)=\frac{1}{2\mu}\left(\mbox{\boldmath$\sigma$}-\frac{\lambda}{2(\mu+\lambda)}({\rm tr}\>\mbox{\boldmath$\sigma$}){\bf I}\right)=:{\cal A}\mbox{\boldmath$\sigma$}$$

$$\mbox{\boldmath$\varepsilon$}({\bf u}_{h})=\frac{1}{2\mu}\left(\mbox{\boldmath$\sigma$}_{h}-p_{h}{\bf I}\right)=\frac{1}{2\mu}\left(\mbox{\boldmath$\sigma$}_{h}-\frac{\lambda}{2(\mu+\lambda)}({\rm tr}\>\mbox{\boldmath$\sigma$}_{h}){\bf I}\right)={\cal A}\mbox{\boldmath$\sigma$}_{h}\>.$$

$$\begin{split}\|\mbox{\boldmath$\tau$}\|_{\cal A}^{2}=({\cal A}\mbox{\boldmath$\tau$},\mbox{\boldmath$\tau$})&=\frac{1}{2\mu}(\mbox{\boldmath$\tau$}-\frac{\lambda}{2(\mu+\lambda)}({\rm tr}\>\mbox{\boldmath$\tau$}){\bf I},\mbox{\boldmath$\tau$})\\ &=\frac{1}{2\mu}({\bf dev}\>\mbox{\boldmath$\tau$},\mbox{\boldmath$\tau$})+\frac{1}{4(\mu+\lambda)}(({\rm tr}\>\mbox{\boldmath$\tau$}){\bf I},\mbox{\boldmath$\tau$})\\ &=\frac{1}{2\mu}\|{\bf dev}\>\mbox{\boldmath$\tau$}\|^{2}+\frac{1}{4(\mu+\lambda)}\|{\rm tr}\>\mbox{\boldmath$\tau$}\|^{2}\geq\frac{1}{2\mu}\|{\bf dev}\>\mbox{\boldmath$\tau$}\|^{2}\>,\end{split}$$

$$\sqrt{2\mu}C_{A}\|\mbox{\boldmath$\tau$}\|_{\cal A}\geq\|\mbox{\boldmath$\tau$}\|\mbox{ for all }\mbox{\boldmath$\tau$}\in H_{\Gamma_{N}}({\rm div},\Omega)^{2}\mbox{ with }{\rm div}\>\mbox{\boldmath$\tau$}={\bf 0}\>,$$

$$\begin{split}|||({\bf u}-{\bf u}_{h},p-p_{h})|||_{T}^{2}&=2\mu\|\mbox{\boldmath$\varepsilon$}({\bf u}-{\bf u}_{h})\|_{L^{2}(T)}^{2}+\frac{1}{\lambda}\|p-p_{h}\|_{L^{2}(T)}^{2}\\ &=\|2\mu\mbox{\boldmath$\varepsilon$}({\bf u}-{\bf u}_{h})-(p-p_{h}){\bf I}\|_{{\cal A},T}^{2}\end{split}$$

$$\|2\mu\mbox{\boldmath$\varepsilon$}({\bf u}-{\bf u}_{h})-(p-p_{h}){\bf I}\|_{{\cal A},h}^{2}=|||({\bf u}-{\bf u}_{h},p-p_{h})|||^{2}\>.$$

$$\begin{split}\|\mbox{\boldmath$\sigma$}_{h}^{R}&-2\mu\mbox{\boldmath$\varepsilon$}({\bf u}_{h})-p_{h}{\bf I}\|_{{\cal A},h}^{2}=\|\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h}^{R}-2\mu\mbox{\boldmath$\varepsilon$}({\bf u}-{\bf u}_{h})-(p-p_{h}){\bf I}\|_{{\cal A},h}^{2}\\ &=\|\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h}^{R}\|_{\cal A}^{2}+\|2\mu\mbox{\boldmath$\varepsilon$}({\bf u}-{\bf u}_{h})-(p-p_{h}){\bf I}\|_{{\cal A},h}^{2}\\ &\;\;\;-2(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h}^{R},{\cal A}(2\mu\mbox{\boldmath$\varepsilon$}({\bf u}-{\bf u}_{h})-(p-p_{h}){\bf I}))_{h}\\ &=\|\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h}^{R}\|_{\cal A}^{2}+\|2\mu\mbox{\boldmath$\varepsilon$}({\bf u}-{\bf u}_{h})-(p-p_{h}){\bf I}\|_{{\cal A},h}^{2}-2(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h}^{R},\mbox{\boldmath$\varepsilon$}({\bf u}-{\bf u}_{h}))_{h}\>,\end{split}$$

$$\begin{split}(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h}^{R},\>&\mbox{\boldmath$\varepsilon$}({\bf u}-{\bf u}_{h}))_{h}=(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h}^{R},\mbox{\boldmath$\nabla$}({\bf u}-{\bf u}_{h}))_{h}-(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h}^{R},{\rm as}\>\mbox{\boldmath$\nabla$}({\bf u}-{\bf u}_{h}))_{h}\\ &=\sum_{E\in{\cal E}_{h}}\langle(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h}^{R})\cdot{\bf n},\llbracket{\bf u}-{\bf u}_{h}\rrbracket_{E}\rangle_{E}-(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h}^{R},{\rm as}\>\mbox{\boldmath$\nabla$}({\bf u}-{\bf u}_{h}))_{h}\\ &=\sum_{E\in{\cal E}_{h}}\langle(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h}^{R})\cdot{\bf n},\llbracket{\bf u}_{h}^{C}-{\bf u}_{h}\rrbracket_{E}\rangle_{E}-(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h}^{R},{\rm as}\>\mbox{\boldmath$\nabla$}({\bf u}-{\bf u}_{h}))_{h}\\ &=(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h}^{R},\mbox{\boldmath$\nabla$}({\bf u}_{h}^{C}-{\bf u}_{h}))_{h}-(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h}^{R},{\rm as}\>\mbox{\boldmath$\nabla$}({\bf u}-{\bf u}_{h}))_{h}\\ &=(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h}^{R},\mbox{\boldmath$\varepsilon$}({\bf u}_{h}^{C}-{\bf u}_{h}))_{h}-(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h}^{R},{\rm as}\>\mbox{\boldmath$\nabla$}({\bf u}-{\bf u}_{h}^{C}))_{h}\\ &=(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h}^{R},\mbox{\boldmath$\varepsilon$}({\bf u}_{h}^{C}-{\bf u}_{h}))_{h}+({\rm as}\>\mbox{\boldmath$\sigma$}_{h}^{R},{\rm as}\>\mbox{\boldmath$\nabla$}({\bf u}-{\bf u}_{h}^{C}))_{h}\\ &=(\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h}^{R},\mbox{\boldmath$\varepsilon$}({\bf u}_{h}^{C}-{\bf u}_{h}))_{h}+({\rm as}\>\mbox{\boldmath$\sigma$}_{h}^{R},\mbox{\boldmath$\nabla$}({\bf u}-{\bf u}_{h}^{C}))_{h}\end{split}$$

$$\displaystyle 2(\mbox{\boldmath$\sigma$}-\>$$

$$\displaystyle\leq 2\sqrt{2\mu}C_{A}\|\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h}^{R}\|_{\cal A}\|\mbox{\boldmath$\varepsilon$}({\bf u}_{h}^{C}-{\bf u}_{h})\|_{h}+2C_{K}\|{\rm as}\>\mbox{\boldmath$\sigma$}_{h}^{R}\|\>\|\mbox{\boldmath$\varepsilon$}({\bf u}-{\bf u}_{h}^{C})\|_{h}$$

$$\displaystyle\leq\|\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h}^{R}\|_{\cal A}^{2}+2\mu C_{A}^{2}\|\mbox{\boldmath$\varepsilon$}({\bf u}_{h}^{C}-{\bf u}_{h})\|_{h}^{2}+\frac{2}{\mu}C_{K}^{2}\|{\rm as}\>\mbox{\boldmath$\sigma$}_{h}^{R}\|^{2}+\frac{\mu}{2}\|\mbox{\boldmath$\varepsilon$}({\bf u}-{\bf u}_{h}^{C})\|_{h}^{2}$$

$$\displaystyle\leq\|\mbox{\boldmath$\sigma$}-\mbox{\boldmath$\sigma$}_{h}^{R}\|_{\cal A}^{2}+(2C_{A}^{2}+1)\mu\|\mbox{\boldmath$\varepsilon$}({\bf u}_{h}^{C}-{\bf u}_{h})\|_{h}^{2}+\frac{2}{\mu}C_{K}^{2}\|{\rm as}\>\mbox{\boldmath$\sigma$}_{h}^{R}\|^{2}+\mu\|\mbox{\boldmath$\varepsilon$}({\bf u}-{\bf u}_{h})\|_{h}^{2}\>.$$

$$\displaystyle\|2\mu\mbox{\boldmath$\varepsilon$}({\bf u}-{\bf u}_{h})$$

$$\displaystyle+(2C_{A}^{2}+1)\mu\|\mbox{\boldmath$\varepsilon$}({\bf u}_{h}^{C}-{\bf u}_{h})\|_{h}^{2}+\frac{2}{\mu}C_{K}^{2}\|{\rm as}\>\mbox{\boldmath$\sigma$}_{h}^{R}\|^{2}+\mu\|\mbox{\boldmath$\varepsilon$}({\bf u}-{\bf u}_{h})\|_{h}^{2}\>,$$

$$\begin{split}|||({\bf u}-{\bf u}_{h},p-p_{h})|||^{2}&\leq 2\big{(}\|\mbox{\boldmath$\sigma$}_{h}^{R}-2\mu\mbox{\boldmath$\varepsilon$}({\bf u}_{h})-p_{h}{\bf I}\|_{{\cal A},h}^{2}\\ &\;\;\;+(2C_{A}^{2}+1)\mu\|\mbox{\boldmath$\varepsilon$}({\bf u}_{h}^{C}-{\bf u}_{h})\|_{h}^{2}+\frac{2}{\mu}C_{K}^{2}\|{\rm as}\>\mbox{\boldmath$\sigma$}_{h}^{R}\|^{2}\big{)}\>.\end{split}$$

$$\begin{split}2\mu&\|\mbox{\boldmath$\varepsilon$}(\widetilde{{\bf u}}-{\bf u})\|^{2}+\frac{1}{\lambda}\|\tilde{p}-p\|^{2}=({\cal P}_{h}{\bf f}-{\bf f},\widetilde{{\bf u}}-{\bf u})+\langle{\cal P}_{h}^{\Gamma_{N}}{\bf g}-{\bf g},\widetilde{{\bf u}}-{\bf u}\rangle_{L^{2}(\Gamma_{N})}\\ &=({\cal P}_{h}{\bf f}-{\bf f},\widetilde{{\bf u}}-{\bf u}-{\cal P}_{h}(\widetilde{{\bf u}}-{\bf u}))+\langle{\cal P}_{h}^{\Gamma_{N}}{\bf g}-{\bf g},\widetilde{{\bf u}}-{\bf u}-{\cal P}_{h}^{\Gamma_{N}}(\widetilde{{\bf u}}-{\bf u})\rangle_{L^{2}(\Gamma_{N})}\>.\end{split}$$

$$\begin{split}2\mu\|\mbox{\boldmath$\varepsilon$}(\widetilde{{\bf u}}-{\bf u})\|^{2}+\frac{1}{\lambda}\|\tilde{p}-p&\|^{2}\leq\sum_{T\in{\cal T}_{h}}\|{\bf f}-{\cal P}_{h}{\bf f}\|_{L^{2}(T)}h_{T}\|\nabla(\widetilde{{\bf u}}-{\bf u})\|_{L^{2}(T)}\\ +&\sum_{E\subset\Gamma_{N}}\|{\bf g}-{\cal P}_{h}^{\Gamma_{N}}{\bf g}\|_{L^{2}(E)}h_{E}^{1/2}\|\widetilde{{\bf u}}-{\bf u}\|_{H^{1/2}(E)}\\ \leq&\left(\sum_{T\in{\cal T}_{h}}h_{T}^{2}\|{\bf f}-{\cal P}_{h}{\bf f}\|_{L^{2}(T)}^{2}\right)^{1/2}\|\nabla(\widetilde{{\bf u}}-{\bf u})\|\\ +&\left(\sum_{E\subset\Gamma_{N}}h_{E}\|{\bf g}-{\cal P}_{h}^{\Gamma_{N}}{\bf g}\|_{L^{2}(E)}^{2}\right)^{1/2}\|\widetilde{{\bf u}}-{\bf u}\|_{H^{1/2}(\Gamma_{N})}\\ \leq C\left(\sum_{T\in{\cal T}_{h}}h_{T}^{2}\|{\bf f}-{\cal P}_{h}{\bf f}\right.&\|_{L^{2}(T)}^{2}\left.+\sum_{E\subset\Gamma_{N}}h_{E}\|{\bf g}-{\cal P}_{h}^{\Gamma_{N}}{\bf g}\|_{L^{2}(E)}^{2}\right)^{1/2}\|\mbox{\boldmath$\varepsilon$}(\widetilde{{\bf u}}-{\bf u})\|\end{split}$$

$$\begin{array}{rl}|||(\tilde{\mathbf{u}}-\mathbf{u},& \tilde{p}-p)|||\\ & \le C{\left(\sum _{T\in {\mathcal{T}}_h}{h}_T^2\Vert \mathbf{f}-{\mathcal{P}}_h\mathbf{f}{\Vert }_{L^2(T)}^2+\sum _{E\subset {\Gamma }_N}{h}_E\Vert \mathbf{g}-{\mathcal{P}}_h^{{\Gamma }_N}\mathbf{g}{\Vert }_{L^2(E)}^2\right)}^{1/2}\end{array}$$

$$({\rm as}\>\mbox{\boldmath$\sigma$}_{h}^{S},{\bf J})_{L^{2}(T)}=0\mbox{ with }{\bf J}=\begin{pmatrix}0&1\\ -1&0\end{pmatrix}\mbox{ for all }T\in{\cal T}_{h}\>.$$

$$\begin{split}({\rm as}\>\mbox{\boldmath$\sigma$}_{h}^{S},\nabla({\bf u}-{\bf u}_{h}^{C}))_{h}&=\sum_{T\in{\cal T}_{h}}({\rm as}\>\mbox{\boldmath$\sigma$}_{h}^{S},\nabla({\bf u}-{\bf u}_{h}^{C}))_{L^{2}(T)}\\ &=\sum_{T\in{\cal T}_{h}}({\rm as}\>\mbox{\boldmath$\sigma$}_{h}^{S},\nabla({\bf u}-{\bf u}_{h}^{C})-\alpha_{T}{\bf J})_{L^{2}(T)}\end{split}$$

Full text

{AMS}

65N30, 65N50

\headers

A posteriori error estimation by stress reconstruction F. Bertrand, M. Moldenhauer, and G. Starke

1 Introduction

This paper is concerned with the nonconforming triangular finite element space of piecewise quadratic functions applied to linear elasticity in two space dimensions. This finite element space possesses some peculiarities due to the existence of a non-zero quadratic polynomial which vanishes at both Gauss points on all three boundary edges. A suitable basis for this space was constructed by Fortin and Soulie in [13] and the corresponding generalization to three dimensions by Fortin in [12]. The special structure of this nonconforming space leads to some advantageous properties of the piecewise gradients which have been exploited for the purpose of flux or stress reconstruction by Ainsworth and co-workers in [3] and [2] and by Kim in [18]. Roughly speaking, their reconstruction algorithm based on the quadratic nonconforming finite element space remains more local and requires less computational work than similar approaches for more general finite element spaces.

Our contribution with this work consists in the modification of this approach to the stress reconstruction associated with incompressible linear elasticity and the corresponding guaranteed a posteriori error estimation. The applicability of the procedure in [2] and [18] is limited to bilinear forms involving the full gradient including the Stokes system which is equivalent to incompressible linear elasticity if Dirichlet conditions are prescribed on the entire boundary. For incompressible linear elasticity with traction forces prescribed on some part of the boundary, the symmetric gradient needs to be used instead and this leads to complications associated with the anti-symmetric part of the stress reconstruction. In order to keep the constants associated with the anti-symmetric stress part under control, a modification like the one presented in this work needs to be done. Of course, one could perform the stress reconstruction in a symmetric $H({\rm div})$-conforming stress space as it is done in [21] and [4] based on the Arnold-Winther elements [6]. But this complicates the stress reconstruction procedure significantly compared to the Raviart-Thomas elements of next-to-lowest order used here. This is particularly true in three dimensions where the symmetric $H({\rm div})$-conforming finite element space from [5] involves polynomials of degree 4 with 162 degrees of freedom per tetrahedron.

Although we restrict our investigation to two space dimensions in this paper, the treatment of three-dimensional elasticity problems is our ultimate goal. The ingredients of our approach can be generalized to the three-dimensional case on tetrahedral elements, in some aspects in a straightforward way, in other aspects with complications. We are convinced that this is one of the most promising routes towards an effective guaranteed a posteriori error estimator for three-dimensional incompressible linear elasticity. We will therefore remark on generalizations to three dimensions at the end of our paper and also point out where the generalization is not so straightforward.

With respect to a triangulation ${\cal T}_{h}$ with the corresponding set of edges denoted by ${\cal E}_{h}$, the nonconforming finite element space of degree 2 is defined by

[TABLE]

where $\llbracket\>\cdot\>\rrbracket_{E}$ denotes the jump across the side $E$. In contrast to the lowest-order case (the nonconforming Crouzeix-Raviart elements), the implementation is not straightforward due to the existence of non-trivial finite element functions whose support is restricted to one element. On the other hand, the quadratic nonconforming finite element space has the remarkable property that the continuity equation is satisfied in average on each element and that the jump of the associated directional derivative in normal direction is zero in average across sides. These properties and the construction of a suitable basis were already contained in the landmark paper by Fortin and Soulie [13] (and generalized to the three-dimensional case in Fortin [12]).

One of the strengths of the quadratic nonconforming finite element space is that it provides a stable mixed method for the approximation of Stokes or incompressible linear elasticity if it is combined with the space of discontinuous piecewise linear functions for the pressure. Due to the fact that it satisfies a discrete Korn’s inequality also in the presence of traction boundary conditions, it is among the popular mixed finite element approaches in common use (cf. [7, Sects. 8.6.2 and 8.7.2]). It received increased attention in recent years in the context of a posteriori error estimation based on reconstructed fluxes or stresses (cf. [3], [2] and [18]). The special properties associated with average element-wise and side-wise conservation of mass or momentum mentioned above also lead to simplified flux or stress reconstruction algorithms. As will be explained in detail in Section 3, the direct application of the approach from [2] and [18] to the displacement-pressure formulation of incompressible linear elasticity leaves two global constants in the reliability bound which remain dependent on the geometry of the domain. In order to overcome this dependence two modifications of the error estimator are performed. Firstly, we construct a correction to the stress reconstruction such that the element-wise average of its anti-symmetric part vanishes. This enables us to multiply the anti-symmetry term in the estimator by a constant associated with an element-wise Korn’s inequality which can be computed from the element shape. Such a shape-dependent Korn’s inequality was used by Kim [17] in association with a posteriori error estimation for a stress-based mixed finite element approach. Secondly, the auxiliary conforming approximation is computed under the additional constraint that its divergence coincides with the nonconforming one. Both steps can be carried out in a local fashion using a vertex-patch based partition of unity. This procedure ensures that the properties of the resulting guaranteed a posteriori error estimator remain uniformly valid in the incompressible limit.

A posteriori error estimation based on stress reconstruction has a long history with ideas dating back at least as far as [19] and [22]. Recently, a unified framework for a posteriori error estimation based on stress reconstruction for the Stokes system was carried out in [15] (see also [11] for polynomial-degree robust estimates). These two references include the treatment of nonconforming methods and both of them contain a historical perspective with a long list of relevant references.

The outline of this paper is as follows. In the next section we review the displacement-pressure formulation for planar linear elasticity, its approximation using quadratic nonconforming finite elements and the associated stress reconstruction procedure in Raviart-Thomas spaces. Based on this, a preliminary version of our a posteriori error estimator is derived in Section 3. Section 4 provides an improved and guaranteed a posteriori estimator where all the constants in the reliability bound are known and depend only on the shape-regularity of the triangulation. To this end, reconstructed stresses with element-wise average symmetry are required. A procedure for this construction is described in Section 5. Section 6 presents an approach to the construction of a conforming approximation with divergence constraint which is also needed for the error estimator of Section 4. Global efficiency is established in Section 7. Finally, Section 8 presents the computational results for the well-known Cook’s membrane problem for different Lamé parameters including the incompressible limit.

2 Displacement-pressure formulation for incompressible linear elasticity and weakly symmetric stress reconstruction

On a bounded domain $\Omega\subset\mathds{R}^{2}$, assumed to be polygonally bounded such that the union of elements in the triangulation ${\cal T}_{h}$ coincides with $\Omega$, the boundary is split into $\Gamma_{D}$ (of positive length) and $\Gamma_{N}=\partial\Omega\backslash\Gamma_{D}$. The boundary value problem of (possibly) incompressible linear elasticity consists in the saddle-point problem of finding ${\bf u}\in H_{\Gamma_{D}}^{1}(\Omega)^{2}$ and $p\in L^{2}(\Omega)$ such that

[TABLE]

holds for all ${\bf v}\in H_{\Gamma_{D}}^{1}(\Omega)^{2}$ and $q\in L^{2}(\Omega)$. ${\bf f}\in L^{2}(\Omega)^{2}$ and ${\bf g}\in L^{2}(\Gamma_{N})^{2}$ are prescribed volume and surface traction forces, respectively. $\mu$ and $\lambda$ are the Lamé parameters characterizing the material properties, where we may assume $\mu$ to be on the order of one and our particular interest lies in large values of $\lambda$ associated with near-incompressibility. In the limit $\lambda\rightarrow\infty$, (2) turns to the Stokes system modelling incompressible fluid flow. For $\Gamma_{D}=\partial\Omega$, the Stokes system may then be restated with the full gradients \mbox{\boldmath\nabla}{\bf u} and \mbox{\boldmath\nabla}{\bf v} and the approaches from [2], [18] may be applied directly component-wise. In general, however, the symmetric part of the gradient appearing in (2) cannot be avoided which complicates the derivation of an a posteriori error estimator as we shall see below.

The resulting finite-dimensional saddle-point problem is then to find ${\bf u}_{h}\in{\bf V}_{h}$ and $p_{h}\in Q_{h}$ such that

[TABLE]

holds for all ${\bf v}_{h}\in{\bf V}_{h}$ and $q_{h}\in Q_{h}$. The $L^{2}(\Omega)$-orthogonal projection ${\cal P}_{h}$ is meant component-wise and ${\cal P}_{h}^{\Gamma_{N}}$ stands for the (component-wise) $L^{2}(\Gamma_{N})$-orthogonal projection onto piecewise linear functions without continuity restrictions. The notation $(\>\cdot\>,\>\cdot\>)_{h}$ stands for the piecewise $L^{2}(T)$ inner product, summed over all elements $T\in{\cal T}_{h}$. The finite element combination of using nonconforming piecewise quadratic elements for (each component of) the displacements in ${\bf V}_{h}$ with piecewise linear discontinuous elements for the pressure in $Q_{h}$ is particularly attractive. Similarly to the Taylor-Hood pair of using conforming quadratic finite elements for ${\bf V}_{h}$ combined with continuous linear finite elements for $Q_{h}$, it achieves the same optimal approximation order for both variables. This leads to the same quality of the stress approximation \mbox{\boldmath\sigma}_{h}({\bf u}_{h},p_{h})=2\mu\mbox{\boldmath\varepsilon}({\bf u}_{h})+p_{h}\>{\bf I} with respect to the $L^{2}(\Omega)$ norm. At the expense of an increased number of degrees of freedom compared to the Taylor-Hood elements (about 1.6 times in 2D, almost 4 times in 3D) on the same triangulation, the nonconforming approach offers advantages for the stress reconstruction which justifies its use.

For the actual computation of the nonconforming finite element approximation $u_{h}$, a basis of the space ${\bf V}_{h}$ is required. Such a basis, consisting of functions with local support, was derived in [13] and [12] for the two- and three-dimensional situation, respectively. The construction uses the fact that ${\bf V}_{h}={\bf V}_{h}^{C}+{\bf B}_{h}^{NC}$ with the conforming subspace ${\bf V}_{h}^{C}$ of continuous piecewise quadratic functions and a nonconforming bubble space ${\bf B}_{h}^{NC}$. In the two-dimensional case, the nonconforming bubble space is given by

[TABLE]

i.e., there is exactly one nonconforming bubble function in ${\bf B}_{h}^{NC}$ for each displacement component per triangle. We denote the corresponding space restricted to an element $T\in{\cal T}_{h}$ by ${\bf B}_{h}^{NC}(T)$. It should be kept in mind that the representation ${\bf V}_{h}={\bf V}_{h}^{C}+{\bf B}_{h}^{NC}$ is not a direct sum, in general. For example, globally constant functions can be expressed in two different ways in these subspaces, in general. In any case, if $\Gamma_{D}$ is a connected subset of positive length of $\partial\Omega$, a basis of ${\bf V}_{h}$ consisting of nonconforming bubble functions in ${\bf B}_{h}^{NC}$ and conforming nodal basis functions in ${\bf V}_{h}^{C}$ can be selected.

The well-posedness of the discrete linear elasticity system (3) relies on the discrete Korn’s inequality

[TABLE]

which is satisfied with some constant $C_{K}$ by the nonconforming quadratic elements (under our assumption that $\Gamma_{D}$ is a subset of the boundary with positive length) due to [8, Thm. 3.1]. It is well-known that the linear nonconforming elements by Crouzeix-Raviart do not satisfy the discrete Korn’s inequality, in general, if $\Gamma_{N}\neq\emptyset$. A second ingredient is the discrete inf-sup stability which can already be found in the original paper [13] (cf. [12] for the three-dimensional case). The a posteriori error estimator will be based on an approximation to the stress tensor \mbox{\boldmath\sigma}=2\mu\mbox{\boldmath\varepsilon}({\bf u})+p{\bf I} in conforming and nonconforming subspaces of $H({\rm div},\Omega)^{2}$.

From the discontinuous piecewise linear approximation \mbox{\boldmath\sigma}_{h}=2\mu\mbox{\boldmath\varepsilon}({\bf u}_{h})+p_{h}{\bf I} obtained from the solution $({\bf u}_{h},p_{h})\in{\bf V}_{h}\times Q_{h}$ of the discrete saddle point problem (3), we first reconstruct an $H({\rm div})$-conforming nonsymmetric stress tensor \mbox{\boldmath\sigma}_{h}^{R} in the Raviart-Thomas space (componentwise) \mbox{\boldmath\Sigma}_{h}^{R} of next-to-lowest order usually denoted by $RT_{1}$ (see, e.g., [7, Sect. 2.3.1]). We will work with the corresponding subspace \mbox{\boldmath\Sigma}_{h}^{R}\subset H_{\Gamma_{N}}({\rm div},\Omega)^{2}, where the normal flux is set to zero on $\Gamma_{N}=\partial\Omega\backslash\Gamma_{D}$. From now on, we also assume that the families of triangulations ${\cal T}_{h}$ are shape-regular to make sure that the constants appearing in our estimates remain independent of $h$. The diameter of an element $T\in{\cal T}_{h}$ is then denoted by $h_{T}$.

We use a stress reconstruction \mbox{\boldmath\sigma}_{h}^{R} in the next-to-lowest order Raviart-Thomas space \mbox{\boldmath\Sigma}_{h}^{R} which is constructed by the following procedure which is equivalent to the algorithms in [2] and [18]. On each edge $E\in{\cal E}_{h}$ we set

[TABLE]

where $T-$ and $T+$ denote the left and right adjacent triangles to $E$, respectively. The continuity of \langle\mbox{\boldmath\sigma}_{h}\cdot{\bf n}_{E},1\rangle_{E} and the identity ({\rm div}\>\mbox{\boldmath\sigma}_{h}+{\bf f},1)_{L^{2}(T)}=0 hold for Fortin-Soulie elements, see [13, Thm. 1]. The remaining four degrees of freedom per element are chosen such that, on each element $T\in{\cal T}_{h}$,

[TABLE]

holds for all $\widehat{{\bf q}}_{h}$ polynomial of degree 1 satisfying $(\widehat{{\bf q}}_{h},1)_{L^{2}(T)}=0$. The same line of proof as in [13, Thm. 1] for the conservation properties fulfilled by the quadratic nonconforming finite elements leads to

[TABLE]

for all $T\in{\cal T}_{h}$ (note that ${\cal P}_{h}{\bf f}$ may be replaced by its average over $T$, cf. [18, Lemma 1]). Combined with (7), {\rm div}\>\mbox{\boldmath\sigma}_{h}^{R}+{\cal P}_{h}{\bf f}={\bf 0} holds if \mbox{\boldmath\sigma}_{h}^{R} is computed from the nonconforming Galerkin approximation.

The above construction may be divided into elementary substeps using the following decomposition of the space \mbox{\boldmath\Sigma}_{h}^{R}: \mbox{\boldmath\Sigma}_{h}^{R}=\mbox{\boldmath\Sigma}_{h}^{R,0}\oplus\mbox{\boldmath\Sigma}_{h}^{R,1}\oplus\mbox{\boldmath\Sigma}_{h}^{R,\Delta} with \mbox{\boldmath\Sigma}_{h}^{R,0} and \mbox{\boldmath\Sigma}_{h}^{R,0}\oplus\mbox{\boldmath\Sigma}_{h}^{R,1} being the corresponding subspaces of the lowest-order Raviart-Thomas elements $RT_{0}$ and Brezzi-Douglas-Marini elements $BDM_{1}$, respectively. In particular,

[TABLE]

The above stress reconstruction procedure can be rewritten in the following steps:

Step 1. Compute \mbox{\boldmath\sigma}_{h}^{R,0}\in\mbox{\boldmath\Sigma}_{h}^{R,0} such that \langle{\bf n}_{E}\cdot(\mbox{\boldmath\sigma}_{h}^{R,0}\cdot{\bf n}_{E}),1\rangle_{E}=\langle{\bf n}_{E}\cdot(\mbox{\boldmath\sigma}_{h}\cdot{\bf n}_{E}),1\rangle_{E}

Step 1. and \langle{\bf t}_{E}\cdot(\mbox{\boldmath\sigma}_{h}^{R,0}\cdot{\bf n}_{E}),1\rangle_{E}=\langle{\bf t}_{E}\cdot(\mbox{\boldmath\sigma}_{h}\cdot{\bf n}_{E}),1\rangle_{E} holds for all $E\in{\cal E}_{h}$.

Step 2. Compute \mbox{\boldmath\sigma}_{h}^{R,1}\in\mbox{\boldmath\Sigma}_{h}^{R,1} s.t. \langle{\bf n}_{E}\cdot(\mbox{\boldmath\sigma}_{h}^{R,1}\cdot{\bf n}_{E}),q_{E}\rangle_{E}=\langle\{\!\{{\bf n}_{E}\cdot(\mbox{\boldmath\sigma}_{h}\cdot{\bf n}_{E})\}\!\},q_{E}\rangle_{E}

Step 2. and \langle{\bf t}_{E}\cdot(\mbox{\boldmath\sigma}_{h}^{R,1}\cdot{\bf n}_{E}),q_{E}\rangle_{E}=\langle\{\!\{{\bf t}_{E}\cdot(\mbox{\boldmath\sigma}_{h}\cdot{\bf n}_{E})\}\!\},q_{E}\rangle_{E} holds for all $q_{E}\in P_{1}(E)$

Step 2. with $\langle q_{E},1\rangle_{E}=0$ for all $E\in{\cal E}_{h}$.

Step 3. Compute \mbox{\boldmath\sigma}_{h}^{R,\Delta}\in\mbox{\boldmath\Sigma}_{h}^{R,\Delta} such that ({\rm div}\>\mbox{\boldmath\sigma}_{h}^{R,\Delta},{\bf q}_{T})_{L^{2}(T)}=({\cal P}_{h}{\bf f},{\bf q}_{T})_{L^{2}(T)} holds Step 3. for all ${\bf q}_{T}\in P_{1}(T)^{2}$ with $({\bf q}_{T},1)_{L^{2}(\Omega)}=0$ for all $T\in{\cal T}_{h}$.

Set \mbox{\boldmath\sigma}_{h}^{R}=\mbox{\boldmath\sigma}_{h}^{R,0}+\mbox{\boldmath\sigma}_{h}^{R,1}+\mbox{\boldmath\sigma}_{h}^{R,\Delta}.

3 A preliminary version of our a posteriori error estimator

In this section we present an a posteriori error estimator based on the stress reconstruction \mbox{\boldmath\sigma}_{h}^{R}. It will be shown to be robust in the incompressible limit $\lambda\rightarrow\infty$ but its reliability bound contains two constants which are generally not known and may become rather large depending on the shapes of $\Omega$, $\Gamma_{D}$ and $\Gamma_{N}$. Therefore, we will modify this a posteriori error estimator in subsequent sections in order to get a guaranteed upper bound for the error which involves only constants that are at our disposal. The derivation of the more straightforward error estimator in this section does, however, contribute to the understanding of the situation and already contains some of the crucial steps of the analysis. One of the unknown constants appearing in the reliability bound is $C_{K}$, the constant from the Korn inequality (5) which, roughly speaking, becomes big when $\Gamma_{D}$ is small. The other constant comes from the following dev-div inequality which was proved under different assumptions in [7, Prop. 9.1.1], [10] and [9]:

[TABLE]

where ${\bf dev}$ denotes the trace-free part given by {\bf dev}\>\mbox{\boldmath\tau}=\mbox{\boldmath\tau}-({\rm tr}\>\mbox{\boldmath\tau})/2. Note that the constant $C_{A}$ may become big if $\Gamma_{N}$ is small (cf., for example, [20, Sect. 2]); the case $\Gamma_{N}=\emptyset$ requires an additional constraint ({\rm tr}\>\mbox{\boldmath\tau},1)=0.

Our aim is to estimate the error in the energy norm, expressed in terms of ${\bf u}-{\bf u}_{h}$ and $p-p_{h}$. Since (2) implies ${\rm div}\>{\bf u}=p/\lambda$ and since ${\rm div}\>{\bf u}_{h}=p_{h}/\lambda$ follows from (3), the energy norm is given by

[TABLE]

Local versions of the energy norm in (11), where the integration is limited to a subset $\omega\subset\Omega$, will be denoted by $|||(\>\cdot\>,\>\cdot\>)|||_{\omega}$. The definition of the stress directly leads to

[TABLE]

which implies

[TABLE]

and

[TABLE]

Note that (13) and (14) remain valid in the incompressibe limit $\lambda\rightarrow\infty$, where ${\cal A}$ tends to the projection onto the trace-free part, scaled by $1/(2\mu)$. Also note that our stress representation is purely two-dimensional while the true stress in plane-strain two-dimensional elasticity possesses a three-dimensional component $\sigma_{33}=p$. The following derivation may equivalently be done based on this three-dimensional stress using a slightly different mapping ${\cal A}_{3}$ instead of ${\cal A}$. At the end, however, the result is the same and we therefore stick to the unphysical but simpler two-dimensional stresses.

The inner product $(\>\cdot\>,\>\cdot\>)_{\cal A}:=(\>{\cal A}(\cdot)\>,\>\cdot\>)$ induces a norm $\|\>\cdot\>\|_{\cal A}$ on the divergence-free subspace of $H_{\Gamma_{N}}({\rm div},\Omega)^{2}$ due to

[TABLE]

which, combined with (10), gives

[TABLE]

again assuming ({\rm tr}\>\mbox{\boldmath\tau},1)=0 if $\Gamma_{N}=\emptyset$. For the derivation of our preliminary estimator in this section, however, we need to assume that $\Gamma_{N}\neq\emptyset$. This restriction will be overcome by the improved estimator in the next section. For the element-wise contribution to the energy norm in (11),

[TABLE]

holds, where $\|\>\cdot\>\|_{{\cal A},T}=({\cal A}(\cdot),\>\cdot\>)_{L^{2}(T)}^{1/2}$ denotes the corresponding element-wise norm. With the corresponding nonconforming version $\|\>\cdot\>\|_{{\cal A},h}=({\cal A}(\cdot),\>\cdot\>)_{h}^{1/2}$ , summing (17) over all elements leads to

[TABLE]

Our a posteriori error estimator will be based on \|\mbox{\boldmath\sigma}_{h}^{R}-2\mu\mbox{\boldmath\varepsilon}({\bf u}_{h})-p_{h}{\bf I}\|_{{\cal A},h}, the difference between post-processed and reconstructed stress and we start the derivation from this term. Let us assume that (2) holds with ${\bf f}={\cal P}_{h}{\bf f}$ and ${\bf g}={\cal P}_{h}^{\Gamma_{N}}{\bf g}$ since this approximation can be treated as an oscillation term (see the remarks at the end of this section). Inserting the relation \mbox{\boldmath\sigma}=2\mu\mbox{\boldmath\varepsilon}({\bf u})+p{\bf I} which holds for the exact solution, we obtain

[TABLE]

where (13) and (14) were used in the last equality. Our goal is to estimate the middle term on the right-hand side which by (18) coincides with the energy norm. The mixed term at the end of (19) can be rewritten as

[TABLE]

where $\llbracket\>\cdot\>\rrbracket_{E}$ denotes the jump across $E$ and the properties {\rm div}\>(\mbox{\boldmath\sigma}-\mbox{\boldmath\sigma}_{h}^{R})=0 in $\Omega$, (\mbox{\boldmath\sigma}-\mbox{\boldmath\sigma}_{h}^{R})\cdot{\bf n}=0 on $\Gamma_{N}$ were used. In (20), ${\bf u}_{h}^{C}\in H_{\Gamma_{D}}^{1}(\Omega)^{d}$ denotes any conforming approximation to ${\bf u}_{h}$. This leads to the upper bound

[TABLE]

Inserting this into (19) gives us

[TABLE]

which, combined with (18), implies

[TABLE]

Finally, the treatment of the right-hand side approximation as an oscillation term is described. To this end, let $({\bf u},p)$ and $(\widetilde{{\bf u}},\tilde{p})$ be the solutions of (2) with right-hand side data $({\bf f},{\bf g})$ and $({\cal P}_{h}{\bf f},{\cal P}_{h}^{\Gamma_{N}}{\bf g})$, respectively. Taking the difference and inserting $\widetilde{{\bf u}}-{\bf u}$ as test function leads to

[TABLE]

Using the local nature of the projections ${\cal P}_{h}$ and ${\cal P}_{h}^{\Gamma_{N}}$, (22) implies

[TABLE]

with a constant $C$. This proves that

[TABLE]

and therefore the right-hand side approximation may be treated as an oscillation term.

4 An improved and guaranteed a posteriori error estimator

We will now construct an improved a posteriori error estimator which avoids the unknown constants $C_{A}$ and $C_{K}$ present in (21). To this end, we perform two modifications, one associated with the stress reconstruction \mbox{\boldmath\sigma}_{h}^{R} and the other one with the conforming approximation ${\bf u}_{h}^{C}$.

The modified stress reconstruction \mbox{\boldmath\sigma}_{h}^{S} will be constructed in such a way that it satisfies

[TABLE]

How this can be achieved will be the topic of section 5. If (25) is satisfied, then the corresponding last term in (20) can be rewritten as

[TABLE]

for any $\alpha_{T}\in\mathds{R}$, $T\in{\cal T}_{h}$. Since for any \mbox{\boldmath\rho}_{T}\in RM(T), the space of rigid-body modes on $T$, we have \nabla\mbox{\boldmath\rho}_{T}=\alpha_{T}{\bf J} with $\alpha_{T}\in\mathds{R}$, (26) leads to

[TABLE]

Since \mbox{\boldmath\rho}_{T}\in RM(T) is arbitrary, this implies

[TABLE]

We may now use Korn’s inequality of the form

[TABLE]

with a constant $C_{K,T}^{\prime}$ which depends only on the shape of $T$, in particular, on the smallest interior angle (see [17, Sect. 3], [16] for detailed formulas). From (27) and (28) we are therefore led to

[TABLE]

with $C_{K}^{\prime}:=\max\{C_{K,T}^{\prime}:T\in{\cal T}_{h}\}$. The constant $C_{K}^{\prime}$ is therefore fully computable and, moreover, of moderate size for shape-regular triangulations.

The constant $C_{A}$ may be avoided by enforcing the constraint

[TABLE]

We will discuss possible approaches for the construction of such a conforming approximation in section 6. If (30) is satisfied, then the first term on the right-hand side in (20) can be rewritten as

[TABLE]

with $\alpha_{T}\in\mathds{R}$, $T\in{\cal T}_{h}$. Since $\alpha_{T}\in\mathds{R}$ can be chosen arbitrarily for each $T\in{\cal T}_{h}$, we obtain

[TABLE]

Since

[TABLE]

the term on the right-hand side of (32) is minimized, if ({\rm tr}(\mbox{\boldmath\sigma}-\mbox{\boldmath\sigma}_{h}^{S}-\alpha_{T}{\bf I}),1)_{L^{2}(T)}=0 holds. Using the version of the dev-div inequality from [7, Prop. 9.1.1] gives, on each element $T\in{\cal T}_{h}$,

[TABLE]

with a constant $C_{A,T}^{\prime}$ which again only depends on the shape-regularity of the triangulation ${\cal T}_{h}$. In fact, these constants coincide with those from (28) which is contained in the following result.

Proposition 4.1**.**

*For a two-dimensional triangulation ${\cal T}_{h}$, the constants $C_{K,T}^{\prime}$ in (28) and $C_{A,T}^{\prime}$ in (33) are identical. *

Proof 4.2**.**

We start from (33) and note that ({\rm tr}\>\mbox{\boldmath\tau},1)_{L^{2}(T)}=0 implies

[TABLE]

Moreover {\rm div}\>\mbox{\boldmath\tau}={\bf 0} implies the existence of \mbox{\boldmath\psi}\in H^{1}(T) with \mbox{\boldmath\tau}=\mbox{\boldmath\nabla}^{\perp}\mbox{\boldmath\psi} (cf. [14, Theorem I.3.1] and therefore (33) turns into

[TABLE]

Component-wise this has the form

[TABLE]

which may be rewritten as

[TABLE]

*This is equivalent to the inequality (28) with \mbox{\boldmath\psi}={\bf u}-{\bf u}_{h}^{C}. *

Combining (32) and (33) with $C_{A,T}^{\prime}$ replaced by $C_{K,T}^{\prime}$ leads to

[TABLE]

Combined with (16), this implies

[TABLE]

Inserting our improved estimates (29) and (34) into (19) and (20), we arrive at

[TABLE]

for any $\delta>0$. Since 2\mu\|\mbox{\boldmath\varepsilon}({\bf u}-{\bf u}_{h})\|_{h}^{2}\leq|||({\bf u}-{\bf u}_{h},p-p_{h})|||^{2} from the definition in (11), this finally leads to

[TABLE]

for any $\delta\in(0,1/2)$. The choice of $\delta$ may be optimized in dependence of the relative size of the individual error estimator terms

[TABLE]

Since these three contributions to the error estimator turn out to be of comparable size in our computations, choosing $\delta=1/4$ is sufficient for our purpose leading to the following guaranteed reliability result.

Theorem 4.3**.**

With the error estimator terms $\eta_{h}^{R}$, $\eta_{h}^{C}$ and $\eta_{h}^{S}$ defined in (36), the error $({\bf u}-{\bf u}_{h},p-p_{h})$, measured in the energy norm defined by (11), satisfies

[TABLE]

5 Construction of stresses with element-wise symmetry on average

This section provides a possible way for the construction of a modified stress reconstruction with the property ({\rm as}\>\mbox{\boldmath\sigma}_{h}^{S},{\bf J})_{L^{2}(T)}=0 for all $T\in{\cal T}_{h}$. To this end, we go back to the original stress reconstruction procedure at the end of section 2. In order to keep the equilibration property {\rm div}\>\mbox{\boldmath\sigma}_{h}^{R}+{\cal P}_{h}{\bf f}={\bf 0} unaffected, we compute a correction in

[TABLE]

where \mbox{\boldmath\Xi}_{h} is the standard conforming piecewise quadratic finite element space with zero boundary conditions on $\Gamma_{N}$. In order to retain the approximation properties of \mbox{\boldmath\sigma}_{h}^{R}, a first attempt would be to compute \mbox{\boldmath\sigma}_{h}^{R,\perp}\in\mbox{\boldmath\Sigma}_{h}^{R,\perp} such that

[TABLE]

Inserting \mbox{\boldmath\sigma}_{h}^{R,\perp}=\mbox{\boldmath\nabla}^{\perp}\mbox{\boldmath\chi}_{h} with \mbox{\boldmath\chi}_{h}\in\mbox{\boldmath\Xi}_{h}, (39) turns out to be equivalent to

[TABLE]

The solution \mbox{\boldmath\chi}_{h}^{\perp}\in\mbox{\boldmath\Xi}_{h} of (40) is uniquely determined by the associated KKT conditions

[TABLE]

for all \mbox{\boldmath\xi}_{h}\in\mbox{\boldmath\Xi}_{h} and $\rho_{h}\in Z_{h}$, where $Z_{h}$ denotes the space of scalar piecewise constant functions with respect to ${\cal T}_{h}$. In (41), $\nu_{h}\in Z_{h}$ plays the role of a Lagrange multiplier. The discrete inf-sup stability of the $P_{2}-P_{0}$ combination in two dimensions (cf. [7, Sect. 8.4.3]) and the coercivity in $H_{\Gamma_{N}}^{1}(\Omega)$ leads to the well-posedness of (41). The modified stress reconstruction is given by \mbox{\boldmath\sigma}_{h}^{S}=\mbox{\boldmath\sigma}_{h}^{R}+\mbox{\boldmath\nabla}^{\perp}\mbox{\boldmath\chi}_{h}.

It would be desirable to replace the global minimization problem (40) by a set of local ones like it will be done later in the next section for the estimation of the non-conformity error. The main incentive for such an approach would be the control of the correction \mbox{\boldmath\sigma}_{h}^{R,\perp} on an element $T$ by the right-hand side ({\rm as}\>\mbox{\boldmath\sigma}_{h}^{R},{\bf J}) in a neighborhood of $T$. Unfortunately, this is not possible, in general, since the following situation may occur in principle: On one element $T^{\ast}$ away from the boundary, ({\rm as}\>\mbox{\boldmath\sigma}_{h}^{R},{\bf J})_{L^{2}(T^{\ast})} does not vanish while ({\rm as}\>\mbox{\boldmath\sigma}_{h}^{R},{\bf J})_{L^{2}(T)}=0 for all $T\in{\cal T}_{h}\backslash\{T^{\ast}\}$. Then, it is not possible to find an admissible \mbox{\boldmath\chi}_{h} such that its support \mbox{supp}\>\mbox{\boldmath\chi}_{h} is contained in a subset $\omega\subset\Omega$ that does not touch the boundary, $\partial\omega\cap\partial\Omega=\emptyset$. This is due to the fact that, in this case,

[TABLE]

would hold. From the computational point of view, (40) constitutes a saddle point problem which requires much less effort to solve than the original one (3).

The following result gives an upper bound for the correction \mbox{\boldmath\nabla}^{\perp}\mbox{\boldmath\chi}_{h} which will later be used for showing global efficiency of the error estimator.

Proposition 5.1**.**

The correction \mbox{\boldmath\chi}_{h}^{\perp}\in\mbox{\boldmath\Xi}_{h} defined by (40) satisfies

[TABLE]

*where the constant $C$ depends only on the shape regularity of ${\cal T}_{h}$.

Proof 5.2**.**

Since \mbox{\boldmath\chi}_{h}^{\perp} is a solution of the KKT system (41), we obtain

[TABLE]

From the well-posedness of (41) we obtain

[TABLE]

with a constant $C$ (cf. [7, Theorem 5.2.1]). A combination of (44) and (45) implies

[TABLE]

Since {\rm as}\>\mbox{\boldmath\sigma}={\bf 0}, we obtain

[TABLE]

*where we used the fact that |{\rm as}\>(\mbox{\boldmath\sigma}-\mbox{\boldmath\sigma}_{h}^{R})|\leq|\mbox{\boldmath\sigma}-\mbox{\boldmath\sigma}_{h}^{R}| holds pointwise. *

6 Distance to divergence-constrained conformity

Our point of departure for the construction of a divergence-constrained conforming approximation is the minimization problem

[TABLE]

among all ${\bf u}_{h}^{C}\in{\bf V}_{h}^{C}$, the subspace of conforming piecewise quadratic functions. The solution of this global minimization problem can be replaced by local ones based on the partition of unity

[TABLE]

with respect to ${\cal V}_{h}^{\prime}=\{z\in{\cal V}_{h}:z\notin\Gamma_{D}\}$. Here, ${\cal V}_{h}$ denotes the set of vertices of the triangulation and $\phi_{z}$, $z\in{\cal V}_{h}^{\prime}$ are continuous piecewise linear functions with support restricted to

[TABLE]

For the partition of unity (49), the standard pyramid basis functions need to be extended for all vertices $z\in{\cal V}_{h}^{\prime}$ adjacent to a boundary vertex on $\Gamma_{D}$ such that it is constant along the connecting edge. This requires that the triangulation ${\cal T}_{h}$ is such that each vertex on $\Gamma_{D}$ is connected by an interior edge. From the decomposition

[TABLE]

we are led, for each $z\in{\cal V}_{h}^{\prime}$, to the problem

[TABLE]

among all ${\bf u}_{h,z}^{C}\in\widehat{{\bf V}}_{h,z}^{C}$, where $\widehat{{\bf V}}_{h,z}^{C}\subset H^{1}(\omega_{z})$ may be any space of conforming finite elements vanishing on all edges not adjacent to $z$. The compatibility condition for the constraint in (52) is satisfied since ${\bf u}_{h,z}$ and ${\bf u}_{h,z}^{C}$ both vanish on $\omega_{z}$. Since ${\bf u}_{h,z}={\bf u}_{h}\Phi_{z}$ is piecewise cubic, using conforming elements of polynomial degree 3 are used for $\widehat{{\bf V}}_{h,z}^{C}$ in order to secure the optimal approximation order. For each $z\in{\cal V}_{h}^{\prime}$, the solution ${\bf u}_{h,z}^{C}\in\widehat{{\bf V}}_{h,z}^{C}$ of the minimization problem (52) is obtained from the KKT system

[TABLE]

for all ${\bf v}_{h,z}^{C}\in\widehat{{\bf V}}_{h,z}^{C}$ and $\rho_{h,z}\in Z_{h,z}$. These local saddle-point problems for ${\bf u}_{h,z}^{C}\in\widehat{{\bf V}}_{h,z}^{C}$ and $\nu_{h,z}\in Z_{h,z}$ are again well-posed due to the inf-sup stability of these combinations of finite element spaces. For the conforming approximation

[TABLE]

$\widehat{{\bf V}}_{h}^{C}\subset H_{\Gamma_{D}}^{1}(\Omega)$ being the piecewise cubic finite element space, one obtains, for each $T\in{\cal T}_{h}$,

[TABLE]

i.e., it satisfies the constraint in (48).

Proposition 6.1**.**

The conforming approximation ${\bf u}_{h}^{C}\in\widehat{{\bf V}}_{h}^{C}$ defined by (53) and (54) satisfies

[TABLE]

*where the constant $C$ depends only on the shape regularity of ${\cal T}_{h}$.

Proof 6.2**.**

Since ${\bf u}_{h,z}^{C}\in\widehat{{\bf V}}_{h,z}^{C}$ solves (52), a simple scaling argument gives us

[TABLE]

for all $z\in{\cal V}_{h}^{\prime}$ (note that the right-hand side being zero implies ${\bf u}_{h,z}\in\widehat{{\bf V}}_{h,z}^{C}$ and therefore the left-hand side also vanishes). From (51) and (54) we get

[TABLE]

Since \langle\llbracket{\bf u}_{h}\rrbracket_{E},1\rangle_{E}{\color[rgb]{0,0,1}=0} is satisfied for all $E\in{\cal E}_{h}$, the same line of reasoning as in [1, Theorem 10] (cf. [15, Sect. 6]) implies that, for all $T\in{\cal T}_{h}$,

[TABLE]

*holds. Combining (57) and (58) finishes the proof. *

7 Global Efficiency of the Improved Estimator

Efficiency of the error estimator is shown if all three terms in (36), $\eta_{h}^{R}$, $\eta_{h}^{C}$ and $\eta_{h}^{S}$, can be bounded by the energy norm of the error multiplied with a constant which remains bounded in the incompressible limit.

Using the definition of ${\cal A}$ in (13) and (16), Proposition 5.1 implies

[TABLE]

The first estimator term $\eta_{h}^{R}$ can therefore be bounded in the form

[TABLE]

where (18) was used in the last equality.

The first term on the right-hand side in (60) can be treated by the following lemma. We will use the notation $a\lesssim b$ to indicate that $a$ is bounded by $b$ times a constant that is independent of the Lamé parameter $\lambda$.

Lemma 7.1**.**

The stress reconstruction computed by the algorithm at the end of Section 2 satisfies

[TABLE]

Proof 7.2**.**

The first step consists in the observation that for all functions \mbox{\boldmath\tau}_{h} which are element-wise of next-to-lowest order Raviart-Thomas type,

[TABLE]

holds. This can be shown in the usual way using a scaling argument and the finite dimension of the considered space. Applying (62) to \mbox{\boldmath\tau}_{h}=\mbox{\boldmath\sigma}_{h}^{R}-2\mu\mbox{\boldmath\varepsilon}({\bf u}_{h})-p_{h}{\bf I}, the right-hand side simplifies since, due to (6) and (7),

[TABLE]

holds, where ${\cal P}_{h}^{0}$ denotes the $L^{2}(\Omega)$-orthogonal projection onto piecewise constant functions. The identity {\rm div}(2\mu\mbox{\boldmath\varepsilon}({\bf u}_{h})+p_{h}{\bf I})=-{\cal P}_{h}^{0}{\bf f} follows from [13, Thm. 1]. From (62) we therefore get

[TABLE]

where the approximation estimate

[TABLE]

for the $L^{2}(T)$-orthogonal projections ${\cal P}_{h}^{0}$ and ${\cal P}_{h}$ onto polynomials of degree 0 and 1 was used. Arguing along the same lines as in [1, Theorem 6] one gets

[TABLE]

where $\omega_{E}$ denotes the union of the two elements adjacent to $E$. Summing over all $T$ leads to

[TABLE]

*which finishes the proof. *

For the second term in (36) we may use Proposition 6.1 to get

[TABLE]

Finally, the third term in (36) satisfies

[TABLE]

We summarize the global efficiency result in the following theorem.

Theorem 7.3**.**

The error estimator terms $\eta_{h}^{R}$, $\eta_{h}^{C}$ and $\eta_{h}^{S}$ defined in (36) satisfy

[TABLE]

The last term in (70) is of the same order as the approximation error is expected to decrease in the ideal case for the finite element spaces studied in this paper. It is, however, not an oscillation term of higher order. Nevertheless, 70 implies that the error estimator decreases proportionally to the approximation error.

8 Computational Results

This section presents computational results with the a posteriori error estimator studied in the previous sections. As an example, Cook’s membrane is considered which consists of the quadrilateral domain $\Omega\in\mathds{R}^{2}$ with corners $(0,0)$, $(0.48,0.44)$, $(0.48,0.6)$ and $(0,0.44)$, where $\Gamma_{D}$ coincides with the left boundary segment. The prescribed surface traction forces on $\Gamma_{N}$ are ${\bf g}={\bf 0}$ on the upper and lower boundary segments and ${\bf g}=(0,1)$ on the right. Starting from an initial triangulation with 44 elements, 17 adaptive refinement steps are performed based on the equilibration strategy, where a subset $\widetilde{{\cal T}}_{h}\subset{\cal T}_{h}$ of elements is refined such that

[TABLE]

holds with $\theta=0.5$ (cf. [23, Sect. 2.1]). Figures 1, 2 and 3 show the convergence behavior in terms of the error estimator for Poisson ratios $\nu=0.29$ (compressible case), $\nu=0.49$ (nearly incompressible case) and $\nu=0.5$ (incompressible case). Since the Poisson ratio is related to the Lamé parameters by $2\mu\nu=\lambda(1-2\nu)$ and since $\mu$ is set to 1 in our computations, this leads to the values $\lambda=1.381$, $\lambda=49$ and $\lambda=\infty$ in the three examples. The solid line (always in the middle) represents the estimator term $\eta_{h}^{R}$, the dashed line below stands for $\eta_{h}^{S}$ measuring the skew-symmetric part and the dotted line shows the values for $\eta_{h}^{C}$, the distance to the conforming space. In all cases, the optimal convergence behavior $\eta_{h}^{\Box}\sim N_{h}^{-1}$, if $N_{h}$ denotes the number of unknowns, is observed. For the investigation of the effectivity of error estimators of the type presented in this paper we refer to [2, 18] where the case of the Stokes equations with $\Gamma_{D}=\partial\Omega$ is treated. The fact that the estimator term $\eta_{h}^{S}$ measuring the symmetry is dominated by the other two contributions $\eta_{h}^{R}$ and $\eta_{h}^{C}$ suggests that the effectivity indices are comparable to those reported in these references.

For the incompressible case, the final triangulation after 17 adaptive refinement steps is shown in Figure 4. As expected, most of the refinement is happening in the vicinity of the strongest singularity at the upper left corner.

Final Remarks. We close our contribution with remarks on the generalization to three-dimensional elasticity computations. As already pointed out in the introduction the properties of the quadratic nonconforming element space were studied in [12] and its combination with piecewise linears again constitutes an inf-sup stable pair for incompressible linear elasticity. The stress reconstruction algorithm of [2] and [18] (see the end of Section 3) can also be generalized in a straightforward way to the three-dimensional case due to the fact that the corresponding conservation properties hold in a similar way on elements and faces as proven in [12]. For the improved and guaranteed estimator, the construction of stresses with element-wise symmetry on average becomes somewhat more complicated in the three-dimensional situation. This is due to the fact that the correction needs to be computed in the space of curls of Nédélec elements leading to a more involved local saddle point structure. For the computation of a divergence-constrained conforming approximation we see, however, no principal complications in three dimensions. Exploring the details of the associated analysis is the topic of ongoing work. The presentation of the results including three-dimensional computations are planned for a future paper.

Acknowledgement. We thank Martin Vohralík for elucidating discussions and for pointing out reference [1] to us. We are also grateful to two anonymous referees for the careful reading of our manuscript and for helpful suggestions. In particular, both of them found an error in an earlier version of the correction procedure in Section 5.