Weakly symmetric stress equilibration for hyperelastic materialmodels

TL;DR

This paper introduces a novel weakly symmetric stress equilibration method for hyperelastic materials, improving surface traction force calculations through a vertex-patch-wise $H(div)$-conforming stress approximation.

Contribution

It develops a new stress equilibration procedure that ensures weak symmetry and enhances the accuracy of stress and traction force computations in hyperelastic models.

Findings

Improved surface traction force results in computational experiments.

The method guarantees weak symmetry in the stress approximation.

Analytical and numerical validation of momentum balance properties.

Abstract

A stress equilibration procedure for hyperelastic material models is proposed andanalyzed in this paper. Based on the displacement-pressure approximation computed with a stable finite element pair, it constructs, in a vertex-patch-wise manner, an $H(div)$-conforming approximation to the first Piola-Kirchhoff stress. This is done in such a way that its associated Cauchy stress is weakly symmetric in the sense that its anti-symmetric part is zero tested against continuous piecewise linear functions. Our main result is the identification of the subspace of test functions perpendicular to the range of the local equilibration system on each patch which turn out to be rigid body modes associated with the current configuration. Momentum balance properties are investigated analytically and numerically and the resulting stress reconstruction is shown to provide improved results for surface traction forces by computational experiments.

Tables

Table 1. Table 1: Approximated resultant normal traction force for 𝐏 ^ h ​ ( 𝐮 h ) subscript ^ 𝐏 ℎ subscript 𝐮 ℎ \widehat{{\bf P}}_{h}({\bf u}_{h})

${ℐ}_{D,n}\left({\widehat{𝐏}}_h\left({𝐮}_h\right)\right)$	$\gamma =0.05$	$\gamma =0.2$	$\gamma =0.5$
${𝒯}_3$	$1.69\cdot {10}^{-3}$	$8.29\cdot {10}^{-3}$	$2.31\cdot {10}^{-2}$
${𝒯}_5$	$9.59\cdot {10}^{-4}$	$5.40\cdot {10}^{-3}$	$2.59\cdot {10}^{-3}$

Table 2. Table 2: Approximated resultant normal traction force for 𝐏 h R superscript subscript 𝐏 ℎ 𝑅 {\bf P}_{h}^{R}

${ℐ}_{D,n}\left({𝐏}_h^R\right)$	$\gamma =0.05$	$\gamma =0.2$	$\gamma =0.5$
${𝒯}_3$	$1.74\cdot {10}^{-11}$	-$3.80\cdot {10}^{-9}$	-$2.87\cdot {10}^{-10}$
${𝒯}_5$	$8.11\cdot {10}^{-11}$	$5.96\cdot {10}^{-10}$	$2.92\cdot {10}^{-9}$

Table 3. Table 3: Approximated resultant normal traction force for 𝐏 h R superscript subscript 𝐏 ℎ 𝑅 {\bf P}_{h}^{R}

$\gamma =0.2$	${𝒯}_3$	${𝒯}_4$	${𝒯}_5$
${\Vert \left({𝐏}_h^R-{𝐏}_{2h}^R\right)\cdot 𝐧\Vert }_{-1/2,{\Gamma }_D}$	$4.4075\cdot {10}^{-3}$	$2.3382\cdot {10}^{-3}$	$1.2470\cdot {10}^{-3}$
rate $\alpha$		0.915	0.907
$\parallel (𝐏\left({𝐮}_h\right)-𝐏\left({𝐮}_{2h}\right)\cdot 𝐧{\parallel }_{-1/2,{\Gamma }_D}$	$2.7801\cdot {10}^{-3}$	$2.4555\cdot {10}^{-3}$	$2.1781\cdot {10}^{-3}$
rate $\alpha$		0.179	0.173

Equations

$$({\bf P}({\bf u}),\mbox{\boldmath$\nabla$}{\bf v})=({\bf f},{\bf v})+\langle{\bf g},{\bf v}\rangle_{0,\Gamma_{N}}$$

$${\psi }_{NH}(\mathbf{B})=\frac{1}{2}\left(\mu \mathrm{tr}\mathbf{B}+\frac{\lambda }{2}\det (\mathbf{B})-\left(\mu +\frac{\lambda }{2}\right)\ln (\det (\mathbf{B}))\right).$$

$$\mathbf{P}(\mathbf{u})={\partial }_{\mathbf{F}}{\psi }_{NH}(\mathbf{B}(\mathbf{u}))=\mu \mathbf{F}(\mathbf{u})+\left(\frac{\lambda }{2}(\det (\mathbf{B}(\mathbf{u}))-1)-\mu \right)\mathbf{F}(\mathbf{u}{)}^{-T}$$

$$\mathbf{P}(\mathbf{u},p)=\mu \mathbf{F}(\mathbf{u})+\left(p\left(1+\frac{p}{2\lambda }\right)-\mu \right)\mathbf{F}(\mathbf{u}{)}^{-T}$$

$$\begin{split}({\bf P}({\bf u},p),\mbox{\boldmath$\nabla$}{\bf v})&=({\bf f},{\bf v})+\langle{\bf g},{\bf v}\rangle_{0,\Gamma_{N}}\mbox{ for all }{\bf v}\in{\bf V}\>,\\ (\det({\bf F}({\bf u}))-1,q)-\frac{1}{\lambda}(p,q)&=0\hskip 65.44142pt\mbox{ for all }q\in Q^{\prime}\end{split}$$

$$(\mathbf{P}({\mathbf{u}}_h),\nabla {\mathbf{v}}_h)=(\mathbf{f},{\mathbf{v}}_h)+⟨\mathbf{g},{\mathbf{v}}_h{⟩}_{0,{\Gamma }_N}$$

$$\begin{split}({\bf P}({\bf u}_{h},p_{h}),\mbox{\boldmath$\nabla$}{\bf v}_{h})&=({\bf f},{\bf v}_{h})+\langle{\bf g},{\bf v}_{h}\rangle_{0,\Gamma_{N}}\mbox{ for all }{\bf v}_{h}\in{\bf V}_{h}\>,\\ (\det({\bf F}({\bf u}_{h}))-1,q_{h})-\frac{1}{\lambda}(p_{h},q_{h})&=0\hskip 64.01869pt\mbox{ for all }q_{h}\in Q_{h}\end{split}$$

$$\displaystyle\mbox{\boldmath$\Pi$}_{h}^{\Delta}=\{{\bf P}_{h}:\Omega\rightarrow{\rm I\kern-2.5ptR}^{d\times d}\mbox{ with }\left.{\bf P}_{h}\right|_{T}\in P_{k}(T)^{d\times d}+P_{k}(T)^{d}{\bf x}^{T}\}\>,$$

$${\mathcal{S}}_h^{\ast }:=\{S\in {\mathcal{S}}_h:S⊈{\Gamma }_D\}.$$

$$[[{\mathbf{P}}_h\cdot \mathbf{n}{]]}_S={{\mathbf{P}}_h\cdot \mathbf{n}|}_{T_{-}}-{{\mathbf{P}}_h\cdot \mathbf{n}|}_{T_{+}}.$$

$$[[{\mathbf{P}}_h\cdot \mathbf{n}{]]}_S={{\mathbf{P}}_h\cdot \mathbf{n}|}_{T_{-}}.$$

$$\llbracket{\bf P}_{h}\cdot{\bf n}\rrbracket_{S}^{\ast}=\left\{\begin{array}[]{lcr}\left.{\bf P}_{h}\cdot{\bf n}\right|_{T_{-}}-{\bf g}&,\mbox{ if }&S\subset\Gamma_{N}\>,\\ \llbracket{\bf P}_{h}\cdot{\bf n}\rrbracket_{S}&,\mbox{ if }&S\nsubseteq\Gamma_{N}\>.\end{array}\right.$$

$$\begin{split}({\rm div}\>{\bf P}_{h}^{\Delta},{\bf z}_{h})_{h}&=-({\bf f}+{\rm div}\>\widehat{\bf P}_{h}({\bf u}_{h}),{\bf z}_{h})_{h}\hskip 11.38092pt\mbox{ for all }{\bf z}_{h}\in{\bf Z}_{h}\>,\\ \langle\llbracket{\bf P}_{h}^{\Delta}\cdot{\bf n}\rrbracket_{S},\mbox{\boldmath$\zeta$}\rangle_{S}&=-\langle\llbracket\widehat{\bf P}_{h}({\bf u}_{h})\cdot{\bf n}\rrbracket_{S}^{\ast},\mbox{\boldmath$\zeta$}\rangle_{S}\hskip 17.07182pt\mbox{ for all }\mbox{\boldmath$\zeta$}\in P_{k}(S)^{d}\>,\>S\in{\cal S}_{h}^{\ast}\>,\\ ({\bf P}_{h}^{\Delta}{\bf F}({\bf u}_{h})^{T},{\bf J}(\mbox{\boldmath$\gamma$}_{h}))&=-(\widehat{{\bf P}}_{h}({\bf u}_{h}){\bf F}({\bf u}_{h})^{T},{\bf J}(\mbox{\boldmath$\gamma$}_{h}))\mbox{ for all }\mbox{\boldmath$\gamma$}_{h}\in{\bf X}_{h}\>.\end{split}$$

$${\bf J}(\theta):=\begin{pmatrix}\;\;0&\theta\\ -\theta&0\end{pmatrix}\mbox{ for }d=2\mbox{ and }{\bf J}(\mbox{\boldmath$\theta$}):=\begin{pmatrix}\;\;0&\;\;\theta_{3}&-\theta_{2}\\ -\theta_{3}&\;\;0&\;\;\theta_{1}\\ \;\;\theta_{2}&-\theta_{1}&\;\;0\end{pmatrix}\mbox{ for }d=3$$

$$1\equiv \sum _{z\in {\mathcal{V}}_h}{\tilde{\varphi }}_z\text{mbox}on\Omega ,$$

$${\tilde{\omega }}_z:=\bigcup \{T\in {\mathcal{T}}_h:z\text{mbox}isavertexofT\}.$$

$$1\equiv \sum _{z\in {\mathcal{V}}_h^{\prime }}{\varphi }_z\text{mbox}on\Omega .$$

$${\omega }_z:=\bigcup \{T\in {\mathcal{T}}_h:{\varphi }_z=1\text{mbox}foratleastonevertexz\text{mbox}ofT\}.$$

$$\mbox{\boldmath$\Pi$}_{h,z}^{\Delta}=\{{\bf q}_{h}\in\mbox{\boldmath$\Pi$}_{h}^{\Delta}:{\bf q}_{h}\cdot{\bf n}={\bf 0}\mbox{ on }\partial\omega_{z}\backslash\partial\Omega\>,\>{\bf q}_{h}\equiv{\bf 0}\mbox{ on }\Omega\backslash\overline{\omega}_{z}\}$$

$${\mathbf{P}}_h^{\Delta }=\sum _{z\in {\mathcal{V}}_h^{\prime }}{\mathbf{P}}_{h,z}^{\Delta },$$

$$\begin{split}({\rm div}\>{\bf P}_{h,z}^{\Delta},{\bf z}_{h,z})_{\omega_{z},h}&=-(({\bf f}+{\rm div}\>\widehat{{\bf P}}_{h}({\bf u}_{h}))\phi_{z},{\bf z}_{h,z})_{\omega_{z},h}\mbox{ for all }{\bf z}_{h,z}\in{\bf Z}_{h,z}\>,\\ \langle\llbracket{\bf P}_{h,z}^{\Delta}\cdot{\bf n}\rrbracket_{S},\mbox{\boldmath$\zeta$}\rangle_{S}&=-\langle\llbracket\widehat{{\bf P}}_{h}({\bf u}_{h})\cdot{\bf n}\rrbracket_{S}\,\phi_{z},\mbox{\boldmath$\zeta$}\rangle_{S}\hskip 24.18501pt\mbox{ for all }\mbox{\boldmath$\zeta$}\in P_{k}(S)^{d}\>,\>S\in{\cal S}_{h,z}\>,\\ ({\bf P}_{h,z}^{\Delta}{\bf F}({\bf u}_{h})^{T},{\bf J}(\mbox{\boldmath$\gamma$}_{h,z}))_{\omega_{z}}&=-(\widehat{{\bf P}}_{h}({\bf u}_{h}){\bf F}({\bf u}_{h})^{T}\phi_{z},{\bf J}(\mbox{\boldmath$\gamma$}_{h,z}))_{\omega_{z}}\hskip 11.38092pt\mbox{ for all }\mbox{\boldmath$\gamma$}_{h,z}\in{\bf X}_{h,z}\>.\end{split}$$

$$\mathrm{div}{\mathbf{P}}_{h,z}^{\Delta }$$

$$[[{\mathbf{P}}_{h,z}^{\Delta }\cdot \mathbf{n}{]]}_S$$

$$({\bf P}({\bf z}),\mbox{\boldmath$\nabla$}{\bf v})_{\omega_{z}}=({\bf f},{\bf v})_{\omega_{z}}+\langle{\bf g},{\bf v}\rangle_{\partial\omega\cap\Gamma_{N}}\mbox{ for all }{\bf v}\in{\bf V}_{z}$$

$$({\bf P}({\bf z}),\nabla{\bf v})_{\omega_{z}}=({\bf P}({\bf z}){\bf F}({\bf z})^{T},\mbox{\boldmath$\nabla$}{\bf v}{\bf F}({\bf z})^{-1})_{\omega_{z}}$$

$$\text{mbox}span\{\left(\begin{array}{c}1\\ 0\end{array}\right),\left(\begin{array}{c}0\\ 1\end{array}\right),\left(\begin{array}{c}{x}_2+z_2\\ -(x_1+z_1)\end{array}\right)\}$$

$$\text{mbox}span\{\left(\begin{array}{c}1\\ 0\\ 0\end{array}\right),\left(\begin{array}{c}0\\ 1\\ 0\end{array}\right),\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right),\left(\begin{array}{c}0\\ x_3+z_3\\ -(x_2+z_2)\end{array}\right),\left(\begin{array}{c}-(x_3+z_3)\\ 0\\ x_1+z_1\end{array}\right),\left(\begin{array}{c}{x}_2+z_2\\ -(x_1+z_1)\\ 0\end{array}\right)\}.$$

$$\begin{split}{\bf R}_{h,z}^{\perp}=\{&({\bf z}_{h,z},{\bf s}_{h,z},\mbox{\boldmath$\gamma$}_{h,z})\in{\bf Z}_{h,z}\times{\bf S}_{h,z}\times{\bf X}_{h,z}:\\ &({\rm div}\>{\bf P}_{h,z}^{\Delta},{\bf z}_{h,z})_{\omega_{z},h}-\sum_{S\in{\cal S}_{z,h}}\langle\llbracket{\bf P}_{h,z}^{\Delta}\cdot{\bf n}\rrbracket_{S},{\bf s}_{h,z}\rangle_{S}+({\bf P}_{h,z}^{\Delta},{\bf J}(\mbox{\boldmath$\gamma$}_{h,z}){\bf F}({\bf u}_{h}))_{\omega_{z}}=0\mbox{ for all }{\bf P}_{h,z}^{\Delta}\in\mbox{\boldmath$\Pi$}_{h,z}^{\Delta}\}\>,\end{split}$$

$$\begin{split}{\bf R}_{h,z}^{\perp}&=\{({\cal P}_{h,z}^{k}\mbox{\boldmath$\rho$},\{{\cal P}_{h,S}^{k}\mbox{\boldmath$\rho$}\}_{S\in{\cal S}_{h,z}},\mbox{\boldmath$\theta$}):(\mbox{\boldmath$\rho$},\mbox{\boldmath$\theta$})\in{\bf R}{\bf M}({\bf u}_{h})\times{\rm I\kern-2.5ptR}^{d(d-1)/2}\mbox{ such that }{\bf J}(\mbox{\boldmath$\theta$}){\bf F}({\bf u}_{h})=\mbox{\boldmath$\nabla$}\mbox{\boldmath$\rho$}\}\mbox{ if }|\partial\omega_{z}\cap\Gamma_{D}|=0\>,\\ {\bf R}_{h,z}^{\perp}&=\{({\bf 0},{\bf 0},{\bf 0})\}\mbox{ if }|\partial\omega_{z}\cap\Gamma_{D}|>0\>.\end{split}$$

$${\bf J}(\mbox{\boldmath$\gamma$}_{h,z}){\bf F}({\bf u}_{h})=\nabla(\gamma_{1}\mbox{\boldmath$\rho$}_{1}+\gamma_{2}\mbox{\boldmath$\rho$}_{2}+\gamma_{3}\mbox{\boldmath$\rho$}_{3})\>.$$

$$({\rm div}\>{\bf P}_{h,z}^{\Delta},{\bf z}_{h,z})_{\omega_{z}}+({\bf P}_{h,z}^{\Delta},{\bf J}(\mbox{\boldmath$\gamma$}_{h,z}){\bf F}({\bf u}_{h}))_{\omega_{z}}=0\>.$$

Full text

††articletype: Original Paper

1 Introduction

This paper is concerned with a stress equilibration procedure for hyperelastic material models in nonlinear solid mechanics. It extends the approach proposed and studied in our earlier work 1 to the case of geometrically and materially nonlinear elasticity in the form of a hyperelastic material law. Due to the fact that the symmetry condition does not hold for the first Piola-Kirchhoff stress (which is the result from the reconstruction process) but for the Cauchy stress, the use of symmetric stress elements is not feasible anymore in the hyperelastic case. The weak symmetry condition from linear elasticity can, however, be generalized to a suitable constraint for the Piola-Kirchhoff stress as is done in this contribution. To the best of our knowledge, our contribution is the first attempt to develop a stress equilibration procedure for the hyperelastic situation. Our hope is that this will be of use for the development of an a posteriori error estimator for hyperelastic problems in the future. The issue of a posteriori error estimation and adaptive refinement is, however, beyond the scope of this contribution.

Expressing the internal forces of a material, the components of the stress-tensor are crucial for the prediction of the weakening of a material, including plastic behavior or damage. A specific application area where this is an issue is associated with implant shape design which constitutes an optimal control problem, see 2. Therefore, the accurate approximation of the stress-tensor is of strong importance in numerous applications and in particular in the hyperelastic material model this paper is concerned with. The mathematical foundations of hyperelastic material models in solid mechanics are covered, e.g., in the books by Marsden and Hughes 3 and Ciarlet 4. The numerical treatment of the associated variational problems are investigated in detail by Le Tallec 5. Specifically for incompressible hyperelasticity, issues connected to the use of displacement-pressure formulations are discussed in 6. A priori analysis of numerical methods are available under restrictive assumptions, see Carstensen and Dolzmann 7 and, for a least-squares finite element approach, Müller et.al. 8.

Common displacement-based approaches or, in the incompressible regime, mixed displacement-pressure formulations for this model lead to approximations of the stresses that are not $H({\rm div})$-conforming, i.e., have discontinuities of the normal components on the interface between two elements. In particular, this means on the one hand that they do not control momentum conservation and on the other hand that the normal component of the boundary traces are not well-defined implying that the approximation of the surface traction forces can also not be guaranteed. In contrast to variational principles involving a direct approximation of the stress in an $H({\rm div})$-conforming space (see Chapter 9 of the monograph 9 for an overview), this paper proposes an algorithm to obtain an $H({\rm div})$-conforming approximation of the stress-tensor by post-processing the displacement-based approximation.

The idea of reconstructing the matrix-valued stress and vector-valued flux goes back to the hypercircle theorem by Prager and Synge 10 (see also Section III.9 in Braess’ book 11 for a presentation in modern mathematical language). Besides the accurate approximation in an $H({\rm div})$-conforming space, the stress or flux reconstruction builds the basis of an a posteriori error estimator, which was actually already one of the motivations of Prager and Synge 10. Over the years, a posteriori error estimators based on flux reconstruction were explored in detail in many contributions 12, 13, 14, 15, 16. An important algorithmic innovation was given by Braess and Schöberl 17 by the equilibration procedure which is completely local and provides the link to residual error estimation. An important aspect of the use of reconstruction-based error estimation of the above type is that it provides guaranteed upper bounds for the error with accessible constants. Another important aspect is that these a posteriori error estimators are valid for any approximation that is inserted into the procedure. In particular, it does not assume that the underlying finite-dimensional variational problems are solved to high precision. The extension of reconstruction strategies to linear elasticity was the subject of a number of contributions in the last two decades 18, 19, 20, 21, 22, stress reconstruction in the context of Stokes flow was also studied recently 23. More recently, a posteriori error estimation based on the reconstruction of weakly symmetric stresses was investigated in our earlier work 24 and 1. In particular, the stress equilibration procedure considered in our recent contribution 1 serves as a point of departure for our treatment of hyperelastic material models in the present paper. The recent paper by Botti and Riedlbeck 25 should also be mentioned here. It treats nonlinear elasticity restricted to a geometrically linear situation. In that case, the (Piola-Kirchhoff) stress is still symmetric which allows the use of symmetric stress elements as it is done in the approach by Botti and Riedlbeck 25.

We emphasize once more that our paper does not discuss the issue of a posteriori error estimation. The development of an a posteriori error estimator based on the stress equilibration for hyperelastic material models and, in particular, its analysis are expected to be rather involved and to require rather restrictive assumptions. After all, it is well-known that the solution of the variational problem may not be unique (see the examples in Chapter 5 in 4). We nevertheless hope that our stress reconstruction procedure will be of use for the future study of such an a posteriori error estimator. For the time being, we concentrate on other motivations for the use of equilibrated stresses like the enhanced accuracy of surface force approximations which will be studied in detail. Other approaches to the direct finite element approximation of stresses in geometrically nonlinear elasticity can be found e.g. in 26 and 8.

Besides the fact that the symmetry condition for the stresses becomes more complicated in the geometrically and materially nonlinear situation associated with hyperelastic models which was already mentioned above, other challenging issues arise if one wants to extend the stress equilibration procedure from our recent work 1 to that case. The stresses computed directly from displacement and, possibly, pressure approximations are no longer piecewise polynomial due to the nonlinearity of the model. Therefore, in order to get a stress reconstruction in an appropriate $H({\rm div})$-conforming finite element space, a suitable projection to piecewise polynomial stresses need to be carried out first. Another problem is concerned with the subspace of test functions which are perpendicular to the range of the local equilibration systems for vertex patches not connected to the Dirichlet boundary. The main result of this contribution is the identification of these subspaces as associated with rigid body modes in the current configuration, i.e., involving the displacement approximations. The right-hand sides arising from straightforward piecewise polynomial projections of the stresses are shown to have components outside of these ranges which means that the local equilibration systems possess an additional compatibility error. We propose a remedy involving a more complicated test space to overcome this problem. This leads to compatible local problems and thus to a truly equilibrated stress reconstruction.

The outline of this paper is as follows. We start with the variational formulation of elastic deformations governed by hyperelastic material models and the weakly symmetric stress reconstruction in Section 2. Section 3 presents the local equilibration algorithm. The solvability of the local problems on vertex patches is analysed in 4. In particular, the subspace of test functions orthogonal to the range of the local operators associated with equilibration is identified and this result is used for the investigation of the compatibility of the right-hand side. Section 5 proposes our remedy to deal with this problem and derives a more complicated test space which leads to compatible local equilibration systems for which an inf-sup condition holds. The improved accuracy of the surface forces associated with the equilibrated stresses will be the topic of Section 6. Finally, computational results illustrating the properties of the equilibrated stresses are collected in Section 7.

2 Hyperelasticity and weakly symmetric stress reconstruction

The hyperelastic problems under our consideration are based on an open, bounded and connected domain $\Omega\subset{\rm I\kern-2.5ptR}^{d}$ ($d=2,3$) with Lipschitz-continuous boundary which constitutes the reference configuration of the undeformed state. The boundary is divided into two disjoint non-empty subsets $\Gamma_{D}$ and $\Gamma_{N}$. On $\Gamma_{D}$, homogeneous displacement boundary conditions ${\bf u}={\bf 0}$ are imposed, while surface traction forces ${\bf P}\cdot{\bf n}={\bf g}$ are prescribed on $\Gamma_{N}$. For an appropriate subspace ${\bf V}\subset H_{\Gamma_{D}}^{1}(\Omega)^{d}$, the boundary value problem of hyperelasticity then consists in the variational problem of finding ${\bf u}\in{\bf V}$ such that

[TABLE]

holds for all ${\bf v}\in{\bf V}$. Here, ${\bf P}({\bf u})=\partial_{\bf F}\psi({\bf B})$ denotes the first Piola-Kirchhoff stress tensor with respect to the stored energy function $\psi:{\rm I\kern-2.5ptR}^{d\times d}_{\text{sym}}\to{\rm I\kern-2.5ptR}$, where the deformation gradient is given by {\bf F}({\bf u})={\bf I}+\mbox{\boldmath\nabla}{\bf u} and the left Cauchy-Green strain tensor is defined as ${\bf B}({\bf u})={\bf F}({\bf u}){\bf F}({\bf u})^{T}$. Simple brackets $(\>\cdot\>,\>\cdot\>)$ as in (1) will from now on always abbreviate the inner product in $L^{2}(\Omega)$ with respect to the reference configuration; ${\bf f}$ and ${\bf g}$ stand for volume and surface loads, transformed back to the reference configuration. An example of a stored energy function which we will also use later in our computations in Section 7 is associated with the Neo-Hookean model

[TABLE]

In this case, the Piola-Kirchhoff stress tensor is given by

[TABLE]

and ${\bf V}=W_{\Gamma_{D}}^{1,4}(\Omega)^{d}$ would be sufficient for the variational problem (1) to be properly defined. In order to deal with materials in the incompressible parameter regime ($\lambda\gg\mu$), a pressure-like variable may be introduced, e.g. by setting $p=\lambda(\det({\bf F}({\bf u}))-1)$. Note that with this choice, $p$ does not really stand for the physical pressure but that it is possible to obtain the pressure from $p$ even in the incompressible limit. The above choice is motivated from the fact that it turns into the constraint $p={\rm div}\>{\bf u}$, familiar from linear elasticity, in the small strain limit. Other options for the definition of $p$ are possible and may have advantages. The Piola-Kirchhoff stress is now given in terms of ${\bf u}$ and $p$ which, in the Neo-Hookean example, reads

[TABLE]

due to the fact that $\det({\bf B}({\bf u}))-1=\det({\bf F}({\bf u}))^{2}-1=(\det({\bf F}({\bf u}))-1)(\det({\bf F}({\bf u}))+1)$ holds. With a pressure space $Q$ ($Q=L^{4/3}(\Omega)$ would be appropriate in the Neo-Hooke case), the variational problem turns into one of saddle point type which consists in finding ${\bf u}\in{\bf V}$ and $p\in Q$ such that

[TABLE]

with $Q^{\prime}$ denoting the dual space of $Q$ ($Q^{\prime}=L^{4}(\Omega)$ with the above choices for the Neo-Hooke case).

For $k\geq 1$, let ${\bf V}_{h}\subset{\bf V}$ be the subspace of continuous piecewise polynomials of degree $k+1$ with respect to a triangulation ${\cal T}_{h}$ for each component of ${\bf V}_{h}$. Our finite-dimensional variational problem for hyperelasticity consists in finding ${\bf u}_{h}\in{\bf V}_{h}$ such that

[TABLE]

holds for all ${\bf v}_{h}\in{\bf V}_{h}$. In the incompressible regime, a discrete pressure space $Q_{h}$ consisting of continuous piecewise polynomials of degree $k$ may be used to define a corresponding discrete saddle point problem. It consists in finding $({\bf u}_{h},p_{h})\in{\bf V}_{h}\times Q_{h}$ such that

[TABLE]

is satisfied. The direct use of ${\bf P}({\bf u}_{h})$ or, in the incompressible regime, ${\bf P}({\bf u}_{h},p_{h})$ as an approximation for the Piola-Kirchhoff stress, has, however, certain deficiencies which are already known from the linear elasticity situation. Most importantly, ${\bf P}({\bf u}_{h})\cdot{\bf n}$ is not continuous at interfaces between elements of the underlying triangulation implying that traction forces are not well-defined. It also means that ${\bf P}({\bf u}_{h})$ is not $H({\rm div})$-conforming and that the conservation of momentum is not controlled. This motivates the need to construct an $H({\rm div})$-conforming stress reconstruction ${\bf P}_{h}^{R}$ with all these desired properties.

The idea of equilbration is to compute the reconstructed stress ${\bf P}_{h}^{R}$ in the $H({\rm div})$-conforming Raviart-Thomas space of degree $k$ as an additive correction to ${\bf P}({\bf u}_{h})$. This is done using the broken Raviart-Thomas space of degree $k$ for each row leading to

[TABLE]

where $P_{k}(T)$ denotes the space of polynomials of degree $k$ on the triangle ($d=2$) or tetrahedron ($d=3$) $T$. In other words, each row of the stress tensor {\bf P}_{h}\in\mbox{\boldmath\Pi}_{h}^{\Delta} is element-wise given by a function in the Raviart-Thomas space. Unfortunately, in contrast to the linear elasticity situation, {\bf P}({\bf u}_{h})\in\mbox{\boldmath\Pi}^{\Delta} does not hold, in general, due to the nonlinearity of the stress-strain relation. Obviously, for the Neo-Hookean model in (3), ${\bf P}({\bf u}_{h})$ is not even piecewise polynomial. Therefore, ${\bf P}({\bf u}_{h})$ needs to be projected first to an element \widehat{{\bf P}}_{h}({\bf u}_{h})\in\mbox{\boldmath\Pi}_{h}^{\Delta}. An obvious candidate would be to set $\widehat{{\bf P}}_{h}({\bf u}_{h})={\cal P}_{h}^{k}{\bf P}({\bf u}_{h})$, where ${\cal P}_{h}^{k}$ denotes the component-wise and element-wise $L^{2}$-orthogonal projection onto $P_{k}(T)$. We will stick with this choice of $\widehat{{\bf P}}_{h}({\bf u}_{h})$ for the moment until we present an alternative one in Section 5 as a remedy for certain deficiencies associated with it.

Following the weakly symmetric equilibration procedure from 1, we perform the construction for the difference ${{\bf P}}_{h}^{\Delta}:={\bf P}_{h}^{R}-\widehat{{\bf P}}_{h}({\bf u}_{h})$ between the reconstructed and the projected original stress. Recall that the extension of the hypercircle theorem to linear elasticity requires a symmetric reconstruction satisfying the equilibration condition ${\rm div}\>{\bf P}_{h}^{\Delta}=-{\bf f}-{\rm div}\>\widehat{\bf P}_{h}({\bf u}_{h})$ in each triangle and the jump condition allowing ${\bf P}_{h}^{R}$ to be $H({\rm div})$-conforming. In order to write this jump condition in a precise way, let ${\cal S}_{h}$ denote the set of all sides (edges in 2D and faces in 3D) of the triangulation ${\cal T}_{h}$ and ${\cal S}_{h}^{\ast}$ the set of sides not contained in $\Gamma_{D}$

[TABLE]

Further, for all sides $S\in{\cal S}_{h}$, let ${\bf n}$ be the normal direction associated with $S$ (depending on its orientation), $T_{+}$ and $T_{-}$ the elements adjacent to $S$ (such that ${\bf n}$ points into $T_{+}$) and the jump of ${\bf P}_{h}$ over $S$ defined by

[TABLE]

For sides $S\subset\Gamma_{N}$ located on the Neumann boundary we assume that ${\bf n}$ points outside of $\Omega$ and define the jump by

[TABLE]

In order to use the same formulas also for patches adjacent to the Neumann boundary $\Gamma_{N}$ we define the auxiliary jump by

[TABLE]

With this, the jump condition for the correction reads $\llbracket{\bf P}_{h}^{\Delta}\cdot{\bf n}\rrbracket_{S}=-\llbracket\widehat{\bf P}_{h}({\bf u}_{h})\cdot{\bf n}\rrbracket_{S}^{\ast}$ for all sides $S\in{\cal S}_{h}^{\ast}$.

Similarly as in 1, the symmetry condition will be imposed weakly in order to obtain a reconstructed stress with reasonable symmetry properties. In the hyperelastic setting, symmetry does not hold for ${\bf P}({\bf u})$ but instead for the related Cauchy stress tensor \mbox{\boldmath\sigma}({\bf u})={\bf P}({\bf u}){\bf F}({\bf u})^{T}/\det({\bf F}({\bf u})) which adequately describes stresses in the deformed configuration. Rewritting the equilibration and jump conditions in a weak form and applying the weak symmetry condition to ${\bf P}_{h}^{R}{\bf F}({\bf u})^{T}$ leads to the following conditions for ${\bf P}_{h}^{\Delta}$:

[TABLE]

${\bf Z}_{h}$ may be chosen to be the space of discontinuous $d$-dimensional vector functions which are piecewise polynomial of degree $k$, and ${\bf X}_{h}$ may stand for the continuous $d(d-1)/2$-dimensional vector functions which are piecewise polynomial of degree $k$ with {\bf J}(\mbox{\boldmath\theta}) being defined by

[TABLE]

for every $d(d-1)/2$-dimensional vector $\theta$. This choice is motivated by the inf-sup stability of the corresponding combination with the use of Raviart-Thomas element of degree $k\geq 1$ as stress approximation space in the Hellinger-Reissner formulation (see Boffi, Brezzi and Fortin 27).

3 Local stress equilibration algorithm

For the sake of the efficient computation of the stress reconstruction, we localize the problem using a partition of unity. The commonly used partition of unity with respect to the set ${\cal V}_{h}$ of all vertices of ${\cal T}_{h}$,

[TABLE]

consists of continuous piecewise linear functions $\tilde{\phi}_{z}$. In this case, the support of $\tilde{\phi}_{z}$ is restricted to

[TABLE]

In analogy to the stress equilibration procedure described in 1 for the linear elasticity case, we modify this classical partition of unity in order to exclude patches formed by vertices $z\in\Gamma_{N}$, where the local problems may possess to few degrees of freedom to be solvable. To this end, let ${\cal V}_{h}^{\prime}=\{z\in{\cal V}_{h}:z\notin\Gamma_{N}\}$ denote the subset of vertices which are not located on a side (edge/face) of $\Gamma_{N}$. The modified partition of unity is defined by

[TABLE]

For $z\in{\cal V}_{h}^{\prime}$ not connected by an edge to $\Gamma_{N}$ the function $\phi_{z}$ is equal to $\tilde{\phi}_{z}$. Otherwise, the function $\phi_{z}$ has to be modified in order to account for unity at the connected vertices on $\Gamma_{N}$. For each $z_{N}\in\Gamma_{N}$ one vertex $z_{I}\notin\Gamma_{N}$ connected by an edge with $z_{N}$ is chosen and $\tilde{\phi}_{z_{I}}$ is extended by the value $1$ along the edge from $z_{I}$ to $z_{N}$ to obtain the modified function $\phi_{z_{I}}$. The support of $\phi_{z}$ is denoted by

[TABLE]

For the partition of unity (14) to hold, we require the triangulation ${\cal T}_{h}$ to be such that each vertex on $\Gamma_{N}$ is connected to an interior edge. For the localized equilibration algorithm, we will also need the local subspaces

[TABLE]

for all $z\in{\cal V}_{h}^{\prime}$. Moreover, we need to work with the local sets of sides ${\cal S}_{h,z}:=\{S\in{\cal S}_{h}:S\subset\overline{\omega}_{z}\}$ and the restrictions ${\bf Z}_{h,z}$ and ${\bf X}_{h,z}$ to $\omega_{z}$ of the test spaces ${\bf Z}_{h}$ and ${\bf X}_{h}$, respectively. The conditions in (10) can be restated for a sum of patch-wise contributions

[TABLE]

where, for each $z\in{\cal V}_{h}^{\prime}$, {\bf P}_{h,z}^{\Delta}\in\mbox{\boldmath\Pi}_{h,z}^{\Delta} is computed such that $\|{\bf P}_{h,z}^{\Delta}\|_{\omega_{z}}^{2}$ is minimized subject to the following constraints:

[TABLE]

The minimization in (18) is necessary since solutions to (18) are not expected to be unique, in general, similarly to the linear elasticity case treated in our earlier work 1. At this point, we may introduce the local orthogonal projections ${\cal P}_{h,z}^{k}:L^{2}(\omega_{z})\rightarrow{\bf Z}_{h,z}$ and ${\cal P}_{h,S}^{k}:L^{2}(S)\rightarrow P_{k}(S)^{d}$ which means that the first two conditions in (18) can be written shortly as

[TABLE]

For each $z\in{\cal V}_{h}^{\prime}$, (18) constitutes a low-dimensional quadratic minimization problem with linear constraints for which standard methods are available for the efficient solution. Note that it is not guaranteed at this point that (18) has a solution at all. In fact, it does not, in general, as will become clear from the results of the next section. This is the reason why we will modify the test space in Section 5 in order to have well-posed local patch problems.

To get an idea about the structure of the system (18) and as a motivation for the result in the next section, we consider its underlying continuous problem. On the continuous level, the system (18) constitutes the stress-based dual formulation of the variational problem (1) restricted to $\omega_{z}$. With a suitable subspace ${\bf V}_{z}\subset H_{\Gamma_{D}\cap\partial\omega_{z}}^{1}(\omega_{z})^{d}$ this means that ${\bf z}\in{\bf V}_{z}$ is sought such that

[TABLE]

holds. On vertex patches with $\Gamma_{D}\cap\partial\omega_{z}=\emptyset$, there is a non-trivial subspace ${\bf V}_{z}^{\circ}\subset{\bf V}_{z}$ of test functions such that ({\bf P}({\bf z}),\mbox{\boldmath\nabla}{\bf v})=0 for ${\bf v}\in{\bf V}_{z}^{\circ}$. Obviously, all constants are contained in ${\bf V}_{z}^{\circ}$. Moreover, since

[TABLE]

holds and since ${\bf P}({\bf z}){\bf F}({\bf z})^{T}$ is a symmetric matrix, also all ${\bf v}$ with \mbox{\boldmath\nabla}{\bf v}{\bf F}({\bf z})^{-1} being skew-symmetric will be contained in ${\bf V}_{z}^{\circ}$. In two dimensions, we arrive at

[TABLE]

being contained in ${\bf V}^{\circ}$, and for $d=3$ this is true for

[TABLE]

These are exactly the rigid body modes associated with the current configuration deformed by \mbox{\boldmath\varphi}({\bf x})={\bf x}+{\bf z} which we would like to denote by ${\bf R}{\bf M}({\bf z})$ from now on. From the above derivation, it should not be surprising that the corresponding rigid body mode spaces ${\bf R}{\bf M}({\bf u}_{h})$ will appear in the investigation of the well-posedness of the discrete local problems (18) in the following section.

4 Solvability of the local problems on vertex patches

We turn our attention to the solvability of the local minimization problem subject to the constraints (18). To this end, we need to guarantee that for every right hand side, a function {\bf P}_{h,z}^{\Delta}\in\mbox{\boldmath\Pi}_{h,z}^{\Delta} exists such that the constraints (18) are satisfied. The left-hand side in (18) defines a linear operator ${\cal L}_{h,z}:\Pi_{h,z}^{\Delta}\rightarrow{\bf Z}_{h,z}^{\prime}\times{\bf S}_{h,z}^{\prime}\times{\bf X}_{h,z}^{\prime}$, where {\bf S}_{h,z}=\{\mbox{\boldmath\zeta}\in P_{k}(S)^{d}:S\in{\cal S}_{h,z}\} denotes the trace space on the interior sides and $(\>\cdot\>)^{\prime}$ stands for the dual space. The subspace ${\bf R}_{h,z}^{\perp}\subseteq{\bf Z}_{h,z}\times{\bf S}_{h,z}\times{\bf X}_{h,z}$ orthogonal to the range of ${\cal L}_{h,z}$, i.e., the null space of its adjoint ${\cal L}_{h,z}^{\ast}$, is obviously of interest for the solvability since the linear functionals on the right-hand side in (18) need to vanish on ${\bf R}_{h,z}^{\perp}$. This subspace can be characterized as follows.

Proposition 4.1**.**

The subspace

[TABLE]

i.e., the null space of the adjoint operator ${\cal L}_{h,z}^{\ast}$ associated with the constraints (18), can be characterized as follows:

[TABLE]

Here, $|\>\cdot\>|$ denotes the $d-1$-dimensional measure of boundary curves or surfaces, respectively.

Proof 4.2**.**

The proof is carried out for $d=3$; the two-dimensional case is much easier and can be derived from the three-dimensional one in the usual way by setting $u_{3}\equiv 0$ and all other functions to be independent of $x_{3}$ (with appropriate modifications of operators such as ${\rm div}$, $\nabla$, ${\bf curl}$, etc.).

1st Step.* We start by showing that the component \mbox{\boldmath\gamma}_{h,z} of ${\bf R}_{h,z}^{\perp}$ in (22) needs to satisfy*

[TABLE]

Let us restrict ourselves to the $H({\rm div})$-conforming subspace of \mbox{\boldmath\Pi}_{h,z}^{\Delta}, i.e., with the property that $\llbracket{\bf P}_{h,z}^{\Delta}\cdot{\bf n}\rrbracket_{S}={\bf 0}$ for all $S\in{\cal S}_{h,z}$. Then, the condition in (22) for the definition of ${\bf R}_{h,z}^{\perp}$ turns into

[TABLE]

By definition, we can write

[TABLE]

We may restrict ourselves further to divergence-free ${\bf P}_{h,z}^{\Delta}$ with ${\bf P}_{h,z}\cdot{\bf n}={\bf 0}$ on the entire boundary $\partial\omega_{z}$. These stress approximations can be written as {\bf P}_{h,z}^{\Delta}={\bf curl}\>\mbox{\boldmath\psi}_{h,z} with \mbox{\boldmath\psi}_{h,z} in the Nédélec space $N_{k}({\cal T}_{h})^{d}$ (cf. 9 Corollary 2.3.2) with boundary conditions {\bf n}\times\mbox{\boldmath\psi}_{h,z}={\bf 0} on $\partial\omega_{z}$. Inserting this into (25) and integrating by parts leads to

[TABLE]

where we used the fact that {\bf curl}\>\mbox{\boldmath\nabla}\mbox{\boldmath\rho}_{1}={\bf curl}\>\mbox{\boldmath\nabla}\mbox{\boldmath\rho}_{2}={\bf curl}\>\mbox{\boldmath\nabla}\mbox{\boldmath\rho}_{3}={\bf 0}. It can be shown that (27) can only hold for all \mbox{\boldmath\psi}_{h,z} if $\nabla\gamma_{1}=\nabla\gamma_{2}=\nabla\gamma_{3}={\bf 0}$ in the following way: In the lowest-order case $k=1$, one may insert as test functions \mbox{\boldmath\psi}_{h,z} with tangential component \mbox{\boldmath\psi}_{h,z}\cdot{\bf t}_{E}\equiv{\bf e}_{i} for $i=1,2,3$ on an interior edge $E\subset\omega_{z}\backslash\partial\omega_{z}$ and \mbox{\boldmath\psi}_{h,z}\cdot{\bf t}_{E^{\prime}}\equiv{\bf 0} on all the other interior edges $E^{\prime}$. If $(\#E)_{z}$ denotes the number of interior edges in $\omega_{z}$, this gives $3(\#E)_{z}$ linearly independent conditions for the $3(\#E)_{z}$ constant values $(\nabla\gamma_{i})\cdot{\bf t}_{E}$ for $i=1,2,3$ on all interior edges $E$. Therefore, (27) implies that the tangential derivatives of $\gamma_{1}$, $\gamma_{2}$ and $\gamma_{3}$ vanish along all interior edges $E$ which implies that $\gamma_{1}$, $\gamma_{2}$ and $\gamma_{3}$ are themselves constant. For the higher-order case, each increase of the polynomial degree from $k-1$ to $k$ gives additional degrees of freedom to be controlled: For each of the three components, one per edge (including edges on $\partial\omega_{z}$), additionally $k-2$ per face (including faces on $\partial\omega_{z}$) and additionally $(k-2)(k-3)/2$ per tetraeder. This is more than compensated for by the additional test functions available in the Nédélec space $N_{k}({\cal T}_{h})$, see 9 Proposition 2.3.5 so that $\gamma_{1}$, $\gamma_{2}$ and $\gamma_{3}$ are still forced by (27) to remain constant. Finally, the fact that $\gamma_{1}$, $\gamma_{2}$ and $\gamma_{3}$ need to be constant implies that (26) can be written as (24).

2nd Step.* Inserting (24) into (25) and, restricting ourselves to {\bf P}_{h,z}\in\mbox{\boldmath\Pi}_{h,z}^{\Delta} with, in addition to $\llbracket{\bf P}_{h,z}\cdot{\bf n}\rrbracket_{S}={\bf 0}$ for all $S\in{\cal S}_{h,z}$, ${\bf P}_{h,z}\cdot{\bf n}={\bf 0}$ on all of $\partial\omega_{z}$ (which is automatically satisfied if $\partial\omega_{z}\cap\Gamma_{D}=\emptyset$), integration by parts leads to*

[TABLE]

with an arbitrary constant ${\bf a}\in{\rm I\kern-2.5ptR}^{3}$. The range of the divergence operator satisfies

[TABLE]

If ${\cal P}_{h,z}^{k,0}$ denotes the $L^{2}(\omega_{z})$-orthogonal projection to ${\bf Z}_{h,z}^{0}$, then (28) implies that {\bf z}_{h,z}={\cal P}_{h,z}^{k,0}(\gamma_{1}\mbox{\boldmath\rho}_{1}+\gamma_{2}\mbox{\boldmath\rho}_{2}+\gamma_{3}\mbox{\boldmath\rho}_{3})+{\bf a} which means that {\bf z}_{h,z}={\cal P}_{h,z}^{k}(\gamma_{1}\mbox{\boldmath\rho}_{1}+\gamma_{2}\mbox{\boldmath\rho}_{2}+\gamma_{3}\mbox{\boldmath\rho}_{3}+\widetilde{{\bf a}}) with some $\widetilde{{\bf a}}\in{\rm I\kern-2.5ptR}^{3}$. Since all rigid body modes \mbox{\boldmath\rho}\in{\bf R}{\bf M}({\bf u}_{h}) which can be written as \mbox{\boldmath\rho}=\gamma_{1}\mbox{\boldmath\rho}_{1}+\gamma_{2}\mbox{\boldmath\rho}_{2}+\gamma_{3}\mbox{\boldmath\rho}_{3}+\widetilde{{\bf a}}, we have the corresponding representation of ${\bf z}_{h,z}$ in (23).

3rd Step.* Now we need to consider the two cases in (23) separately. If $|\partial\omega_{z}\cap\Gamma_{D}|=0$, we have indeed that every pair (\mbox{\boldmath\rho},\mbox{\boldmath\theta})\in{\bf R}{\bf M}({\bf u}_{h})\times{\rm I\kern-2.5ptR}^{d(d-1)/2} with {\bf J}(\mbox{\boldmath\theta}){\bf F}({\bf u}_{h})=\mbox{\boldmath\nabla}\mbox{\boldmath\rho} gives rise to a solution of (22) in the form ({\cal P}_{h,z}^{k}\mbox{\boldmath\rho},\{{\cal P}_{h,S}^{k}\mbox{\boldmath\rho}\}_{S\in{\cal S}_{h,z}},\mbox{\boldmath\theta}). This is due to the fact that, for all {\bf P}_{h,z}^{\Delta}\in\mbox{\boldmath\Pi}_{h,z}^{\Delta},*

[TABLE]

holds. On the other hand, in the case $|\partial\omega_{z}\cap\Gamma_{D}|>0$,

[TABLE]

holds for all {\bf P}_{h,z}^{\Delta}\in\mbox{\boldmath\Pi}_{h,z}^{\Delta}. Choosing {\bf P}_{h,z}^{\Delta}\in\mbox{\boldmath\Pi}_{h,z}^{\Delta} appropriately, this implies that {\cal P}_{h,S}^{k}\mbox{\boldmath\rho}={\bf 0} must hold on all $S\subset\Gamma_{D}$. Since there is at least one side $S\subset\Gamma_{D}$ and due to the special structure of the space ${\bf R}{\bf M}({\bf u}_{h})$, the only possibility is \mbox{\boldmath\rho}={\bf 0}.

Remark 4.3**.**

In the linear elasticity case, Proposition 1 turns into the corresponding result from our earlier work 1, where ({\bf z}_{h,z},{\bf s}_{h,z},\mbox{\boldmath\gamma}_{h,z})=(\mbox{\boldmath\rho},\{\left.\mbox{\boldmath\rho}\right|_{S}\}_{S\in{\cal S}_{h}},\mbox{\boldmath\theta}) for (\mbox{\boldmath\rho},\mbox{\boldmath\theta})\in{\bf R}{\bf M}({\bf u}_{h})\times{\rm I\kern-2.5ptR}^{d(d-1)/2} with {\bf J}(\theta)=\mbox{\boldmath\nabla}\mbox{\boldmath\rho}.

Basic linear algebra tells us that the right-hand side of the linear system (18) is in the range of the operator ${\cal L}_{h,z}$ if it is orthogonal to ${\bf R}_{h,z}^{\perp}$, the null space of ${\cal L}_{h,z}^{\ast}$. Using Proposition 4.1 this is obviously the case for patches $\omega_{z}$ with $|\partial\omega_{z}\cap\Gamma_{D}|>0$ since ${\bf R}_{h,z}^{\perp}$ only contains zero in that case. In the case of interior patches $\omega_{z}$ in the sense that $|\partial\omega_{z}\cap\Gamma_{D}|=0$, we may insert the representation of ${\bf R}_{h,z}^{\perp}$ into the right-hand side of (18). This leads to

[TABLE]

for all ({\bf z}_{h,z},{\bf s}_{h,z},\mbox{\boldmath\gamma}_{h,z})\in{\bf R}_{h,z}^{\perp}. The first two terms vanish since

[TABLE]

holds due to the definition of $\widehat{{\bf P}}_{h}({\bf u}_{h})$ as projection onto piecewise polynomials of degree $k$ and the Galerkin condition (6) which holds for piecewise polynomials of degree $k+1$ (of which \phi_{z}{\cal P}_{h,z}^{k}\mbox{\boldmath\rho} is a fine specimen). Therefore, for each ({\bf z}_{h,z},{\bf s}_{h,z},\mbox{\boldmath\gamma}_{h,z})\in{\bf R}_{h,z}^{\perp}, we end up with the expression

[TABLE]

for the inconsistency of the right-hand side in (18). This motivates the choice of a modified test space such that the term in (30) actually vanishes.

5 A Modification Leading to Equilibrated Stresses

Our construction so far is based on using the simple component-wise $L^{2}(\Omega)$-projection $\widehat{{\bf P}}_{h}({\bf u}_{h})={\cal P}_{h}^{k}{\bf P}({\bf u}_{h})$ onto the space of piecewise polynomials of degree $k$. Due to the incompatibility of the right-hand sides in the local equilibration systems (18) on interior vertex patches, these problems do not possess a solution, in general. It is certainly possible to solve these systems in a least-squares sense but that would mean that we do not get equilibrated stresses from this procedure. In particular, this means that momentum conservation would not be satisfied locally on each element. We will therefore take up our findings from Section 4 and derive a modification of $\widehat{{\bf P}}_{h}({\bf u}_{h})$ such that the right-hand side in (18) becomes compatible. In view of (29), it is reasonable to choose the test space ${\bf Z}_{h,z}$ as well as the test functions $\zeta$ on the sides $S\in{\cal S}_{z,h}$ in such a way that they contain the rigid body modes \mbox{\boldmath\rho}\in{\bf R}{\bf M}({\bf u}_{h}) of the deformed configuration. With this choice, {\cal P}_{h,z}^{k}\mbox{\boldmath\rho}={\cal P}_{h,S}^{k}\mbox{\boldmath\rho}=\mbox{\boldmath\rho} and the sum over the sides in (29) vanishes. A straightforward way to do this consists in building the test spaces on the basis of piecewise polynomials in the deformed variables \mbox{\boldmath\varphi}({\bf x})={\bf x}+{\bf u}_{h}({\bf x}) instead of ${\bf x}$. This choice also makes sense in view of the fact that the quantities ${\rm div}\>{\bf P}_{h}^{R}$ and $\llbracket{\bf P}_{h}^{R}\cdot{\bf n}\rrbracket$ which are actually tested are mappings from the reference configuration to (forces in) the current configuration. In fact, this modification of the test spaces is only needed for the subspace of polynomials of degree 1 and one can use a hierarchical construction where the enrichment to polynomials of higher degree is again based on the reference coordinates from ${\bf x}$.

Let us assume, for the moment, that also the test space in the Galerkin formulation (6) would contain the rigid body modes of the deformed configuration. Then, the compatibility condition in (29) would turn into

[TABLE]

if $\widehat{{\bf P}}_{h}({\bf u}_{h})$ is defined as the $L^{2}(\omega_{z})$-orthogonal projection with respect to piecewise polynomials in the deformed coordinates. The compatibility term in (31) does indeed miraculously cancel out, if \phi_{z}\mbox{\boldmath\rho} is assumed to be in the test space of the Galerkin formulation (6). Using such a test space is not as far-fetched as one might think. It would ensure invariance with respect to the rigid body modes in the deformed configuration which is not fulfilled for the use of standard polynomial-based finite elements. However, such an approach is expected to be too complicated for practical use and therefore we need to come up with a suitable choice for $\widehat{{\bf P}}({\bf u}_{h})$ leading to a compatible right-hand side in the absence of this ideal situation.

We restrict our construction to the lowest-order case $k=1$ and consider the following slightly more general formulation of (18):

[TABLE]

where we are still free to construct $\widehat{{\bf f}}$ in an appropriate way from ${\bf f}$. As test space in the first equation of (32),

[TABLE]

could be chosen, where $\varphi$ again denotes the mapping from the reference to the (approximated) deformed configuration given by \mbox{\boldmath\varphi}({\bf x})={\bf x}+{\bf u}_{h}({\bf x}). The test space for the second equation in (32) would then be given component-wise by transformed polynomials of the form

[TABLE]

However, in order to make sure that the rigid body modes associated with the deformed configuration ${\bf R}{\bf M}({\bf u}_{h})$ are contained in the test space, it is sufficient to replace the original undeformed rigid body modes ${\bf R}{\bf M}({\bf 0})$ in the piecewise polynomial test space by ${\bf R}{\bf M}({\bf u}_{h})$. The test space ${\bf X}_{h,z}$ for the third equation in (32), the weak symmetry condition, may remain unchanged since only constant rotations appear in the compatibility conditions resulting from Proposition 4.1. For these spaces, the compatibility condition

[TABLE]

for all ({\bf z}_{h,z},{\bf s}_{h,z},\mbox{\boldmath\gamma}_{h,z})\in{\bf R}_{h,z}^{\perp} is therefore equivalent to

[TABLE]

due to (29). Making use of the Galerkin condition (6), the compatibility condition (36) is certainly fulfilled if, on all elements $T\in{\cal T}_{h}$,

[TABLE]

holds. Note that \mbox{\boldmath\rho}\circ\mbox{\boldmath\varphi}^{-1}\in{\bf R}{\bf M}({\bf 0}) holds with the original undeformed rigid body modes. The first relation in (37) constitutes $d(d+1)^{2}/2$ conditions (9 in two dimensions, 24 in three dimensions) and can thus be fulfilled by choosing $\widehat{{\bf f}}\in P_{2}(T)^{d}$. The spare degrees of freedoms may be used to minimize $\|\widehat{{\bf f}}-{\bf f}\|_{T}$ among all $\widehat{{\bf f}}\in P_{2}(T)^{d}$ satisfying the constraints. The second relation in (37) constitutes $d(d+1)^{2}/2-d$ conditions (7 in two dimensions, 21 in three dimensions) since the constant rigid body modes gives zero on both sides. These conditions can be fulfilled by $\widehat{{\bf P}}_{h}\in P_{1}(T)^{d\times d}$. Again, a reasonable elimination of the spare degrees of freedoms consists in minimizing $\|\widehat{{\bf P}}_{h}-{\bf P}_{h}({\bf u}_{h})\|_{T}$ among all $\widehat{{\bf P}}_{h}\in P_{1}(T)^{d\times d}$ satisfying the constraints.

We end this section with a remark on the inf-sup stability of the system (32) which follows along the same lines as in 27 for the linear elasticity formulation. It is easy to see that the null space associated with the first and second equation in (32)

[TABLE]

remains unchanged by the modification of the test spaces, i.e.,

[TABLE]

where \mbox{\boldmath\Xi}_{h,z} is the subspace of Nédélec elements (of the first kind) on $\omega_{z}$ with vanishing tangential trace on $\partial\omega_{z}$. All that is left to show for the inf-sup stability of (32) is therefore that

[TABLE]

holds with a constant $\beta>0$. If we define \mbox{\boldmath\xi}_{h,z}^{\varphi}:\mbox{\boldmath\varphi}(\omega_{z})\rightarrow{\rm I\kern-2.5ptR}^{d\times d} by \mbox{\boldmath\xi}_{h,z}^{\varphi}\circ\mbox{\boldmath\varphi}=\mbox{\boldmath\xi}_{h,z}{\bf F}({\bf u}_{h})^{-1}, then, according to the transformation rule of the curl operator (cf. 9 Sect. 2.1.3), \mbox{\boldmath\xi}_{h,z}^{\varphi}\in H({\bf curl}^{\varphi},\mbox{\boldmath\varphi}(\omega_{z})) and

[TABLE]

where ${\bf curl}^{\varphi}$ denotes the curl with respect to the mapped coordinates. The inf-sup condition (40) is therefore equivalent to the existence of a constant $\beta>0$ such that

[TABLE]

with the mapped Nédélec space \mbox{\boldmath\Xi}_{h,z}^{\varphi} holds. This is exactly the inf-sup condition for the original spaces from 9 in mapped coordinates using parametric Raviart-Thomas elements 28 for the stress approximation.

The combination of the inf-sup stability of the system ((32) with the fact that our right-hand side is guaranteed to be in its range ensures that there is a correction ${\bf P}_{z,h}^{\Delta}$ in the broken Raviart-Thomas space leading to an equilibrated stress ${\bf P}_{h}^{R}$ in the end.

6 Improved Approximation of Surface Traction Forces

One of the motivations for the construction of equilibrated stresses is that this leads to approximations of the surface traction forces with an ensured convergence rate. The divergence theorem implies that

[TABLE]

holds for all ${\bf v}\in H^{1}(\Omega)^{d}$. If we assume that $({\bf P}-{\bf P}_{h}^{R})\cdot{\bf n}={\bf 0}$ on $\Gamma_{N}$ and ${\rm div}({\bf P}-{\bf P}_{h}^{R})={\bf 0}$ in $\Omega$ holds (for example, since and ${\bf g}$ are piecewise constant), then (43) turns into

[TABLE]

This implies that

[TABLE]

is satisfied which means that the approximation of the surface traction forces, measured in the $H^{-1/2}(\Gamma)$ norm, converges at least as fast as the stress approximation in the $L^{2}(\Omega)$ norm. Since, by construction, $\|{\bf P}_{h}^{R}-{\bf P}_{h}({\bf u}_{h})\|$ is expected to be locally an $O(h^{2})$-approximation, the term on the right-hand side in (45) will converge at the same order as $\|{\bf P}-{\bf P}({\bf u}_{h})\|$, in general.

If we insert ${\bf v}\in{\bf R}{\bf M}({\bf u}_{h})$, the rigid body modes in the deformed configuration, into the numerators in the middle of (45), then

[TABLE]

since \mbox{\boldmath\nabla}{\bf v}{\bf F}({\bf u}_{h})^{-1}={\bf J}(\mbox{\boldmath\theta}), which constitutes a global version of (24), and $({\bf P}-{\bf P}_{h}^{R}){\bf F}({\bf u}_{h})^{T}$ is weakly symmetric in the sense of (32).

7 Computational Results

We tested our stress equilibration procedure for the well-known Cook’s membrane example with a quadrilateral geometry. The corners of the domain are located at $(0,0)$, $(0.48,0.44)$, $(0.48,0.6)$ and $(0,0.44)$ and the boundary is divided into the left line segment $\Gamma_{D}$ and the lower, right, and upper segments which together form $\Gamma_{N}$. Figure 1 shows this geometry and the triangulation ${\cal T}_{3}$ which is the result of three levels of uniform refinement. The surface traction force on the right boundary segment is ${\bf g}=(0,\gamma)^{T}$ with different values $\gamma>0$, while the upper and lower boundary parts are traction-free; the volume forces ${\bf f}$ are set to zero. In order to test the robustness of our approach with respect to the incompressibility, we set $\mu=1$ and $\lambda=\infty$ in the Neo-Hookean law (4) and use the displacement-pressure approximation from (7) as starting point for our stress equilibration procedure. All our computations are for the lowest-order case $k=1$ using the Taylor-Hood combination of finite element spaces.

Of particular interest is the distribution of the traction forces on the left boundary including the singularity with infinite stress components at the upper left corner. The distribution of the normal traction force along the left boundary is shown in Figure 2, for the load value $\gamma=0.2$, on the triangulation ${\cal T}_{5}$ which results from two further uniform refinements of ${\cal T}_{3}$. The left graph shows the values for ${\bf n}\cdot(\widehat{{\bf P}}_{h}({\bf u}_{h})\cdot{\bf n})$, corresponding to the projected Piola-Kirchhoff stress from the Galerkin approximation. The right graph shows ${\bf n}\cdot({\bf P}_{h}^{R}\cdot{\bf n})$ for the reconstructed stress. Both pictures represent piecewise affine traction force distributions along the vertical axis. At a first glance, one may get the impression that the left distribution “looks better than” the right one. However, at closer inspection it becomes obvious that the reconstructed stress in the right graph is better able to represent the singular behavior at the upper end. More importantly, the surface forces obtained from the reconstructed Piola-Kirchhoff stress ${\bf P}_{h}^{R}$ recover the correct resultant force

[TABLE]

This is a consequence of the divergence theorem which implies

[TABLE]

The approximations ${\cal I}_{D,n}(\widehat{{\bf P}}_{h}({\bf u}_{h}))$ and ${\cal I}_{D,n}({\bf P}_{h}^{R})$ are shown for the two triangulations ${\cal T}_{3}$ and ${\cal T}_{5}$ and several values of $\gamma$ in Tables 1 and 2. Apparently, the values produced by $\widehat{{\bf P}}_{h}$ are not exact while those coming from the stress reconstruction differ from zero only in the range of machine precision.

The reference and the deformed configuration are shown in Figure 3 for $\gamma=0.2$. The picture clearly indicates that this example is well inside the geometrically nonlinear regime. Table 3 compares the convergence of $\|({\bf P}-{\bf P}_{h}^{R})\cdot{\bf n}\|_{-1/2,\Gamma_{D}}$ versus $\|({\bf P}-{\bf P}({\bf u}_{h}))\cdot{\bf n}\|_{-1/2,\Gamma_{D}}$ on a sequence of meshes. Since we do not know the exact values of ${\bf P}\cdot{\bf n}$ on $\Gamma_{D}$, we access the convergence behavior by the computation of $\|({\bf P}_{h}^{R}-{\bf P}_{2h}^{R})\cdot{\bf n}\|_{-1/2,\Gamma_{D}}$ and $\|({\bf P}({\bf u}_{h})-{\bf P}({\bf u}_{2h}))\cdot{\bf n}\|_{-1/2,\Gamma_{D}}$, respectively. The norm is evaluated approximately by

[TABLE]

where $V_{h}$ denotes the space of continuous piecewise linear functions on ${\cal T}_{h}$. The values in Table 3 indicate that the convergence for the equilibrated stresses is quite a bit faster than the $O(h^{\alpha})$-behavior with $\alpha\approx 0.544$ expected from the regularity of the problem. It can also be seen that the convergence rate is much higher than the one obtained for the original stresses.

8 Conclusions

In this paper, a stress equilibration procedure for hyperelastic material models was proposed and investigated. It is necessarily based on a weakly symmetric stress formulation and treats geometrically and materially nonlinear elasticity problems. Our main contribution is the identification of the subspace of test functions perpendicular to the range of the equilibration system on vertex patches not attached to the Dirichlet boundary. This result is then used to propose an appropriate projection for the Piola-Kirchhoff stress in order to get compatible patch problems. For the moment, this stress equilibration procedure is used for its own sake, for example, in order to obtain better approximations of traction forces. Our future goal will be to develop an a posteriori error estimator on the basis of stress equilibration for hyperelastic material models. Clearly, this will only be possible under restrictive assumptions excluding all the known situations where uniqueness of the solution does not hold.