Compatible-Strain Mixed Finite Element Methods for 3D Compressible and Incompressible Nonlinear Elasticity

TL;DR

This paper introduces compatible-strain mixed finite element methods (CSFEMs) for 3D nonlinear elasticity, effectively modeling large strains and complex geometries in both compressible and incompressible regimes with high accuracy and stability.

Contribution

The paper develops a new family of mixed finite element methods using compatible strains, stabilizing terms, and conforming interpolation for nonlinear elasticity in three dimensions, ensuring stability and efficiency.

Findings

Accurate stress and pressure approximation in large strain regimes

No observed numerical artifacts such as locking or hourglass instability

Effective modeling of heterogeneous solids and complex geometries

Abstract

A new family of mixed finite element methods$-$compatible-strain mixed finite element methods (CSFEMs)$-$are introduced for three-dimensional compressible and incompressible nonlinear elasticity. A Hu-Washizu-type functional is extremized in order to obtain a mixed formulation for nonlinear elasticity. The independent fields of the mixed formulations are the displacement, the displacement gradient, and the first Piola-Kirchhoff stress. A pressure-like field is also introduced in the case of incompressible elasticity. We define the displacement in $H^1$, the displacement gradient in $H(curl)$, the stress in $H(div)$, and the pressure-like field in $L^2$. In this setting, for improving the stability of the proposed finite element methods without compromising their consistency, we consider some stabilizing terms in the Hu-Washizu-type functional that vanish at its critical points. Using a conforming interpolation, the solution and the test spaces are approximated with some piecewise polynomial subspaces of them. In three dimensions, this requires using the N\'ed\'elec edge elements for the displacement gradient and the N\'ed\'elec face elements for the stress. This approach results in mixed finite element methods that satisfy the Hadamard jump condition and the continuity of traction on all internal faces of the mesh. This, in particular, makes CSFEMs quite efficient for modeling heterogeneous solids. We assess the performance of CSFEMs by solving several numerical examples, and demonstrate their good performance for bending problems, for bodies with complex geometries, and in the near-incompressible and the incompressible regimes. Using CSFEMs, one can capture very large strains and accurately approximate stresses and the pressure field. Moreover, in our numerical examples, we do not observe any numerical artifacts such as checkerboarding of pressure, hourglass instability, or locking.

Figures

Equations

$$(\mathbf{grad}\mathbf{U}{)}^{IJ}=\partial U^I/\partial {\mathrm{X}}^J,(\mathbf{curl}\mathbf{T}{)}^{IJ}={\epsilon }_{JKL}\partial T^{IL}/\partial {\mathrm{X}}^K,(\mathbf{div}\mathbf{T}{)}^I=\partial T^{IJ}/\partial {\mathrm{X}}^J,$$

$$H^1(T\mathcal{B})$$

$$H^{\mathbf{c}}(\mathcal{B})$$

$$H^{\mathbf{d}}(\mathcal{B})$$

$$\begin{gathered}\scalebox{0.9}{\lx@xy@svg{\hbox{\raise 0.0pt\hbox{\kern 5.25pt\hbox{\ignorespaces\ignorespaces\ignorespaces\hbox{\vtop{\kern 0.0pt\offinterlineskip\halign{\entry@#!@&&\entry@@#!@\cr&&&&&\\&&&&&\\&&&&&\crcr}}}\ignorespaces{\hbox{\kern-3.0pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{}$}}}}}}}{\hbox{\kern 16.87495pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{\text{displacements}\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces{\hbox{\kern 90.39996pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces{\hbox{\kern 47.04999pt\raise-4.75pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\kern 47.04999pt\raise-28.30101pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\lx@xy@drawline@}}{\hbox{\kern 90.39996pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{\text{disp.~{}gradients}\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces{\hbox{\kern 173.03137pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces{\hbox{\kern 122.08598pt\raise-4.75pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\kern 122.08598pt\raise-28.84001pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\lx@xy@drawline@}}{\hbox{\kern 173.03137pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{\text{compatibility}\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces{\hbox{\kern 202.03142pt\raise-4.75pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\kern 202.03142pt\raise-28.175pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\lx@xy@drawline@}}{\hbox{\kern 271.12837pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{}$}}}}}}}{\hbox{\kern 313.39085pt\raise 0.0pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{}$}}}}}}}{\hbox{\kern-5.25pt\raise-39.07498pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{\mathbf{0}\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces{\hbox{\kern 28.74155pt\raise-39.07498pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}{\hbox{\kern 28.74155pt\raise-39.07498pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{H^{1}(T\mathcal{B})\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces{\hbox{\kern 77.13158pt\raise-33.27498pt\hbox{{}\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\hbox{\hbox{\kern 0.0pt\raise-1.57501pt\hbox{$\scriptstyle{\mathbf{grad}}$}}}\kern 3.0pt}}}}}}\ignorespaces{\hbox{\kern 107.17224pt\raise-39.07498pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}{\hbox{\kern 107.17224pt\raise-39.07498pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{H^{\mathbf{c}}(\mathcal{B})\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces{\hbox{\kern 153.65996pt\raise-33.88748pt\hbox{{}\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\hbox{\hbox{\kern 0.0pt\raise-2.1875pt\hbox{$\scriptstyle{\mathbf{curl}}$}}}\kern 3.0pt}}}}}}\ignorespaces{\hbox{\kern 186.83768pt\raise-39.07498pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces{\hbox{\kern 122.08598pt\raise-44.32498pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\kern 122.08598pt\raise-68.89996pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\lx@xy@drawline@}}{\hbox{\kern 186.83768pt\raise-39.07498pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{H^{\mathbf{d}}(\mathcal{B})\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces{\hbox{\kern 229.75691pt\raise-33.88748pt\hbox{{}\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\hbox{\hbox{\kern 0.0pt\raise-2.1875pt\hbox{$\scriptstyle{\mathbf{div}}$}}}\kern 3.0pt}}}}}}\ignorespaces{\hbox{\kern 256.86366pt\raise-39.07498pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces{\hbox{\kern 202.03142pt\raise-44.32498pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\kern 202.03142pt\raise-68.89996pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\lx@xy@drawline@}}{\hbox{\kern 256.86366pt\raise-39.07498pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{L^{2}(T\mathcal{B})\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces{\hbox{\kern 311.14085pt\raise-39.07498pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}\ignorespaces\ignorespaces\ignorespaces{\hbox{\kern 274.12837pt\raise-44.32498pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\kern 274.12837pt\raise-68.89996pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}\ignorespaces\ignorespaces{\hbox{\lx@xy@drawline@}}\ignorespaces{\hbox{\lx@xy@drawline@}}{\hbox{\kern 311.14085pt\raise-39.07498pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{\mathbf{0}}$}}}}}}}{\hbox{\kern-3.0pt\raise-78.14996pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{}$}}}}}}}{\hbox{\kern 44.04999pt\raise-78.14996pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{}$}}}}}}}{\hbox{\kern 88.84998pt\raise-78.14996pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{\text{stress~{}functions}\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces{\hbox{\kern 166.94695pt\raise-78.14996pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}{\hbox{\kern 166.94695pt\raise-78.14996pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{\text{first~{}PK~{}stresses}\ignorespaces\ignorespaces\ignorespaces\ignorespaces}$}}}}}}}\ignorespaces\ignorespaces\ignorespaces\ignorespaces{}{\hbox{\lx@xy@droprule}}\ignorespaces{\hbox{\kern 248.74084pt\raise-78.14996pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\lx@xy@tip{1}\lx@xy@tip{-1}}}}}}{\hbox{\lx@xy@droprule}}{\hbox{\lx@xy@droprule}}{\hbox{\kern 248.74084pt\raise-78.14996pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{\text{equilibrium}}$}}}}}}}{\hbox{\kern 313.39085pt\raise-78.14996pt\hbox{\hbox{\kern 0.0pt\raise 0.0pt\hbox{\hbox{\kern 3.0pt\raise 0.0pt\hbox{$\textstyle{}$}}}}}}}\ignorespaces}}}}\ignorespaces}\end{gathered}$$

$$\mathcal{I}(\mathbf{U},\mathbf{K},\mathbf{P})={\int }_{\mathcal{B}}W(\mathbf{X},\mathbf{K})dV-\text{llangle}\mathbf{P},\mathbf{K}-\mathbf{grad}\mathbf{U}\text{rrangle}-\text{llangle}{\rho }_0\mathbf{B},\mathbf{U}\text{rrangle}-{\int }_{{\Gamma }_t}⟨\overline{\mathbf{T}},\mathbf{U}⟩dA,$$

$$\overline{\mathcal{I}}(\boldsymbol{U},\boldsymbol{K},\boldsymbol{P},p)=\mathcal{I}(\boldsymbol{U},\boldsymbol{K},\boldsymbol{P})\Big{|}_{J(\boldsymbol{K})=1}+\int_{\mathcal{B}}p\,C\big{(}J(\boldsymbol{K})\big{)}dV,$$

$$\mathcal{J}(\mathbf{U},\mathbf{K},\mathbf{P},p)=\overline{\mathcal{I}}(\mathbf{U},\mathbf{K},\mathbf{P},p)+\frac{\alpha }{2}\text{llangle}\mathbf{K}-\mathbf{grad}\mathbf{U},\mathbf{K}-\mathbf{grad}\mathbf{U}\text{rrangle},$$

$$\text{llangle}\mathbf{P},\mathbf{grad}\mathbf{{\rm Y}}\text{rrangle}+\alpha s_1\left(\mathbf{U},\mathbf{K},\mathbf{{\rm Y}}\right)$$

$$\text{llangle}\tilde{\mathbf{P}}(\mathbf{K}),\mathbf{\kappa }\text{rrangle}-\text{llangle}\mathbf{P},\mathbf{\kappa }\text{rrangle}+\text{llangle}p\mathbf{Q}(\mathbf{K}),\mathbf{\kappa }\text{rrangle}+\alpha s_2\left(\mathbf{U},\mathbf{K},\mathbf{\kappa }\right)$$

$$\text{llangle}\mathbf{grad}\mathbf{U},\mathbf{\pi }\text{rrangle}-\text{llangle}\mathbf{K},\mathbf{\pi }\text{rrangle}$$

$$\text{llangle}C(J),q\text{rrangle}$$

$$f(\mathbf{{\rm Y}})=\text{llangle}{\rho }_0\mathbf{B},\mathbf{{\rm Y}}\text{rrangle}+{\int }_{{\Gamma }_t}⟨\overline{\mathbf{T}},\mathbf{{\rm Y}}⟩dA,$$

$$s_1\left(\mathbf{U},\mathbf{K},\mathbf{{\rm Y}}\right)$$

$$s_2\left(\mathbf{U},\mathbf{K},\mathbf{\kappa }\right)$$

$$\text{llangle}\mathbf{P},\mathbf{grad}\mathbf{{\rm Y}}\text{rrangle}+\alpha s_1\left(\mathbf{U},\mathbf{K},\mathbf{{\rm Y}}\right)$$

$$\text{llangle}\hat{\mathbf{P}}(\mathbf{K}),\mathbf{\kappa }\text{rrangle}-\text{llangle}\mathbf{P},\mathbf{\kappa }\text{rrangle}+\alpha s_2\left(\mathbf{U},\mathbf{K},\mathbf{\kappa }\right)$$

$$\text{llangle}\mathbf{grad}\mathbf{U},\mathbf{\pi }\text{rrangle}-\text{llangle}\mathbf{K},\mathbf{\pi }\text{rrangle}$$

$${\mathcal{P}}_r^{-}(T{\mathbb{R}}^3)$$

$${\mathcal{P}}_r^{\ominus }(T{\mathbb{R}}^3)$$

$${\mathcal{P}}_r^{-}({\text{⊗}}^{2}T{\mathbb{R}}^3)$$

$${\mathcal{P}}_r^{\ominus }({\text{⊗}}^{2}T{\mathbb{R}}^3)$$

$$\dim {\mathcal{P}}_r({\text{⊗}}^{2}T{\mathbb{R}}^3)=3\dim {\mathcal{P}}_r(T{\mathbb{R}}^3)=9\dim {\mathcal{P}}_r({\mathbb{R}}^3)$$

$$\dim {\mathcal{P}}_r^{-}({\text{⊗}}^{2}T{\mathbb{R}}^3)=3\dim {\mathcal{P}}_r^{-}(T{\mathbb{R}}^3)$$

$$\dim {\mathcal{P}}_r^{\ominus }({\text{⊗}}^{2}T{\mathbb{R}}^3)=3\dim {\mathcal{P}}_r^{\ominus }(T{\mathbb{R}}^3)$$

$$\left(\hat{\mathcal{T}},{\mathcal{P}}_r(T\hat{\mathcal{T}}),{\text{textSigma}}^{\hat{\mathcal{T}},1}\right),$$

$$\left(\hat{\mathcal{T}},{\mathcal{P}}_r^{-}({\text{⊗}}^{2}T\hat{\mathcal{T}}),{\text{textSigma}}^{\hat{\mathcal{T}},\mathbf{c}-}\right),\left(\hat{\mathcal{T}},{\mathcal{P}}_r({\text{⊗}}^{2}T\hat{\mathcal{T}}),{\text{textSigma}}^{\hat{\mathcal{T}},\mathbf{c}}\right),$$

$$\left(\hat{\mathcal{T}},{\mathcal{P}}_r^{\ominus }({\text{⊗}}^{2}T\hat{\mathcal{T}}),{\text{textSigma}}^{\hat{\mathcal{T}},\mathbf{d}-}\right),\left(\hat{\mathcal{T}},{\mathcal{P}}_r({\text{⊗}}^{2}T\hat{\mathcal{T}}),{\text{textSigma}}^{\hat{\mathcal{T}},\mathbf{d}}\right),$$

$$\left(\hat{\mathcal{T}},{\mathcal{P}}_r(\hat{\mathcal{T}}),{\text{textSigma}}^{\hat{\mathcal{T}},\ell }\right),$$

$${\mathbf{h}}_{3i-2}^{\hat{\mathcal{T}}}=\left[\begin{array}{c}{l}_i^r\\ 0\\ 0\end{array}\right],{\mathbf{h}}_{3i-1}^{\hat{\mathcal{T}}}=\left[\begin{array}{c}0\\ l_i^r\\ 0\end{array}\right],{\mathbf{h}}_{3i}^{\hat{\mathcal{T}}}=\left[\begin{array}{c}0\\ 0\\ l_i^r\end{array}\right],i=1,2,...,m.$$

$$l_1^1=1-{\xi }^1-{\xi }^2-{\xi }^3,l_2^1={\xi }^1,l_3^1={\xi }^2,l_4^1={\xi }^3.$$

$${\varphi }_{I,J}^{\hat{\mathcal{T}},{\hat{\mathcal{E}}}_k}(\mathbf{T})$$

$${\varphi }_{I,J}^{\hat{\mathcal{T}},{\hat{\mathcal{F}}}_l}(\mathbf{T})$$

$${\varphi }_{I,J}^{\hat{\mathcal{T}},\hat{\mathcal{T}}}(\mathbf{T})$$

$${\mathbf{r}}_{1,J}^{\hat{\mathcal{T}},{\hat{\mathcal{E}}}_k}={\left[\begin{array}{c}{\mathbf{v}}_J^{\hat{\mathcal{T}},{\hat{\mathcal{E}}}_k}\\ 0\\ 0\end{array}\right]}_{3\times 3},{\mathbf{r}}_{2,J}^{\hat{\mathcal{T}},{\hat{\mathcal{E}}}_k}={\left[\begin{array}{c}0\\ {\mathbf{v}}_J^{\hat{\mathcal{T}},{\hat{\mathcal{E}}}_k}\\ 0\end{array}\right]}_{3\times 3},{\mathbf{r}}_{3,J}^{\hat{\mathcal{T}},{\hat{\mathcal{E}}}_k}={\left[\begin{array}{c}0\\ 0\\ {\mathbf{v}}_J^{\hat{\mathcal{T}},{\hat{\mathcal{E}}}_k}\end{array}\right]}_{3\times 3}.$$

$$\displaystyle\phi^{\widehat{\mathcal{T}},\widehat{\mathcal{E}}_{p}}_{M,N}\Big{(}\boldsymbol{r}^{\widehat{\mathcal{T}},\widehat{\mathcal{E}}_{k}}_{I,J}\Big{)}$$

$$\displaystyle\phi^{\widehat{\mathcal{T}},\widehat{\mathcal{F}}_{q}}_{M,N}\Big{(}\boldsymbol{r}^{\widehat{\mathcal{T}},\widehat{\mathcal{F}}_{l}}_{I,J}\Big{)}$$

$$\displaystyle\phi^{\widehat{\mathcal{T}},\widehat{\mathcal{T}}}_{M,N}\Big{(}\boldsymbol{r}^{\widehat{\mathcal{T}},\widehat{\mathcal{T}}}_{I,J}\Big{)}$$

$${\psi }_{I,J}^{\hat{\mathcal{T}},{\hat{\mathcal{F}}}_l}(\mathbf{T})$$

$${\psi }_{I,J}^{\hat{\mathcal{T}},\hat{\mathcal{T}}}(\mathbf{T})$$

Full text

Keywords:

Mixed finite element methods; finite element exterior calculus; nonlinear elasticity; incompressible elasticity; Hilbert complex.

1. 1 Introduction

2. 2 A Mixed Formulation for Nonlinear Elasticity

3. 2.1 Preliminaries

4. 2.2 Mixed Formulations

5. 3 Finite Element Approximations

6. 3.1 Finite Elements

7. 3.2 Finite Element Spaces

8. 3.3 The Compatible-Strain Mixed Finite Element Methods

9. 3.4 Matrix Formulation of CSFEMs

10. 3.5 Solvability and Stability

11. 4 Numerical Examples

12. 5 Concluding Remarks

1 Introduction

It is known that the standard finite elements formulated in terms of the displacement field are not effective for various problems in nonlinear elasticity such as nearly incompressible or incompressible solids, bending analyses, capturing very large strains, and accurate calculation of stress. Developing finite element methods using the mixed formulations of elasticity is one path to overcome these limitations. However, it is a challenge to develop a robust and efficient mixed finite element method for nonlinear elasticity free from numerical instabilities and artifacts. This is more pronounced for problems in 3D as there is a much wider range of deformations in dimension three and 3D problems require more expensive computations. We recently proposed a new family of mixed finite element methods — compatible-strain mixed finite element methods — for 2D compressible [1] and incompressible [2] nonlinear elasticity. Our observations in several numerical examples indicated that these mixed methods have excellent performance in solving various 2D problems and do not suffer from numerical instabilities and artifacts including the difficulties mentioned earlier. In this paper, we extend these mixed methods to 3D compressible and incompressible nonlinear elasticity.

Over the years different approaches have been proposed in the finite element literature to capture large deformations of solids. Here we focus on some well-know works that are based on a mixed formulation and have proved promising for 3D nonlinear problems (see also our literature review in [1] and [2]). Mixed formulations are based on a saddle-point variational principle such as the two-field Hellinger-Reissner principle or the three-field Hu-Washizu principle, see [3] and [4, §1.5]. One of the most commonly used schemes in the literature has been the enhanced strain method originally introduced by Simo and Rifai [5] for infinitesimal strains and later extended to 2D and 3D nonlinear elasticity by Simo and Armero [6] and Simo et al. [7]. In these methods, strain is assumed to be additively decomposed into a compatible part associated with the displacement field, and an enhanced part. The problem is then written as a two-field mixed formulation in terms of the displacement and the enhanced strain, which is derived from a three-field Hu-Washizu-type mixed formulation after eliminating the stress assuming that the enhanced strain and the stress are $L^{2}$-orthogonal. See [8] for a detailed discussion of early developments of enhanced strain methods and their locking and stability. Using an interpolation of strain and stress different from those in the original enhanced strain method, Kasper and Taylor [9] proposed a new mixed method for 2D and 3D linear elasticity called mixed-enhanced strain method and extended it to nonlinear elasticity in [10]. Lamichhane et al. [11] proposed a parameter-dependent modification of the standard Hu-Washizu mixed formulation for 2D linear elasticity and studied its uniform convergence in the incompressible limit for different interpolations. Their study also incorporates the enhanced strain methods proposed in [5] and [9]. Chavan et al. [12] extended the approach introduced in [11] to 3D nonlinear elasticity considering a Mooney-Rivlin material model. By solving several 2D and 3D problems, they demonstrated the good performance of their method in bending problems and for nearly incompressible solids. Reese et al. [13] introduced a new reduced-integration stabilized brick element method for 3D finite elasticity whose stabilization is based on the enhanced strain method. Their scheme is numerically efficient and shows robust performance in bending of thin shells, compression tests, and the near incompressible regime.

The well-posedness of a mixed formulation requires that certain pairs of independent variables are defined in compatible spaces. This is commonly written as an inf-sup condition also known as the LBB condition named after the works of Ladyzhenskaya [14], Babuška [15], and Brezzi [16]. At the discrete level, the satisfaction of LBB condition is a necessary condition for the stability of mixed finite element methods. Thus, only particular combinations of finite element spaces for the independent variables result in convergent methods. Because of theses difficulties, the mixed formulations of 2D problems cannot simply be extended to 3D problems. In other words, one cannot use a mixed formulation with finite element spaces that converge in 2D and only switch the 2D elements with the counterpart 3D elements to obtain a convergent method. In [2] we formulated $96$ different four-field mixed finite element methods for 2D incompressible nonlinear elasticity by considering different combinations of first and second-order finite element spaces to independently approximate displacement, displacement gradient, stress, and pressure. By examining the linearized discrete systems of those mixed methods, we showed that $75$ out of $96$ of them result in singular tangent stiffness (Jacobian) matrices for any mesh and that only the remaining $21$ cases may result in convergent schemes. In this paper, using the same approach, we show that in 3D all the $96$ possible choices of the first and second-order four-field mixed methods lead to singular tangent stiffness matrices for any mesh and regardless of its size. To overcome this difficulty, we add some stabilization terms to the mixed formulations without compromising the consistency of their discretization schemes. This can also help to introduce a convergent mixed method with a fewer degrees of freedom, which is greatly beneficial for computationally expensive 3D problems. An example of such modification is the work of Hughes et al. [17] on the Stokes problem, where they introduced a stabilized mixed finite element method using an equal-order $C^{0}$ interpolation of both velocity and pressure. Furthermore, inspired by the work of Hughes et al. [17], Franca et al. [18] developed a mixed finite element method for nearly incompressible linear elastic solids by adding stabilization terms to the weak formulation associated with the critical point of the Hellinger-Reissner principle. Klaas et al. [19] developed a stabilized displacement-pressure mixed finite element method for 3D finite elasticity by using linear shape functions for both displacement and pressure. In these works, the combinations of the finite element spaces are unstable according to the LBB condition and result in unphysical solutions. However, adding the stabilization terms resulted in convergent mixed methods.

This paper is organized as follows. In §2, we discuss the mixed formulations that are used in the 3D CSFEMs. In §2.1, we review some preliminaries and definitions. In §2.2, by defining suitable Hu-Washizu-type energy functionals we derive a three-field mixed formulation for compressible elastostatics and a four-field mixed formulation for incompressible elastostatics. In §3, we discuss the finite element approximations for the proposed mixed formulations. In §3.1, we define the finite elements (shape functions and degrees of freedom) for the displacement, displacement gradient, stress, and pressure. In §3.2, we define the finite element approximation spaces and use them to introduce the mixed finite element methods in §3.3. Next, the matrix formulation of the mixed finite elements are discussed in §3.4. In §3.5, we investigate singularities of the mixed methods for some combinations of finite element spaces and explain how the stabilizing terms remove those singularities. To study the performance of the 3D CSFEMs, we present several numerical examples in §4 for both compressible and incompressible solids in dimension three. The paper ends by some concluding remarks in §5.

2 A Mixed Formulation for Nonlinear Elasticity

In this section, following [2], we present two mixed formulations one for 3D compressible nonlinear elasticity and one for 3D incompressible nonlinear elasticity.

2.1 Preliminaries

Suppose $\mathbf{X}=(\mathrm{X}^{1},\mathrm{X}^{2},\mathrm{X}^{3})\in\mathbb{R}^{3}$ is the position of a material point in the reference configuration $\mathcal{B}\subset\mathbb{R}^{3}$ with boundary $\partial\mathcal{B}$. For any vector field $\boldsymbol{U}$ and any $\binom{2}{0}$-tensor field $\boldsymbol{T}$, one can define $\binom{2}{0}$-tensors $\mathbf{grad}\,\boldsymbol{U}$ and $\mathbf{curl}\,\boldsymbol{T}$ and a vector field $\mathbf{div}\,\boldsymbol{T}$ with components

[TABLE]

where $\varepsilon_{JKL}$ is the standard permutation symbol, and summation convention for repeated indices is assumed. Suppose $L^{2}(\mathcal{B})$, $L^{2}(T\mathcal{B})$, and $L^{2}(\otimes^{2}T\mathcal{B})$ are the spaces of square integrable scalar fields, vector fields, and $\binom{2}{0}$-tensor fields, respectively. Define the following spaces:

[TABLE]

In the above spaces, $\mathbf{grad}$, $\mathbf{curl}\,$, and $\mathbf{div}$ are defined in the distributional sense. Recalling that $\mathbf{curl}\,(\mathbf{grad}\,\boldsymbol{Y})=\boldsymbol{0}$ and $\mathbf{div}(\mathbf{curl}\,\boldsymbol{T})=\boldsymbol{0}$, one writes the following differential complex [20, 21]:

[TABLE]

where the first arrow is a trivial operator sending zero to zero, and the last arrow indicates the zero operator mapping the $L^{2}$-space to zero. The physical interpretation of this differential complex is as follows: Let $\boldsymbol{U}(\mathbf{X}):=\varphi(\mathbf{X})-\mathbf{X}$, $\mathbf{X}\in\mathcal{B}$, be the displacement field associated with a motion $\varphi:\mathcal{B}\rightarrow\mathbb{R}^{3}$. Then, $\boldsymbol{K}:=\mathbf{grad}\,\boldsymbol{U}$ is the displacement gradient and $\mathbf{curl}\,\boldsymbol{K}=\boldsymbol{0}$ is the necessary condition for the compatibility of $\boldsymbol{K}$. Moreover, given a first Piola-Kirchhoff stress tensor $\boldsymbol{P}$, the equilibrium equation $\mathbf{div}\,\boldsymbol{P}=\boldsymbol{0}$ is the necessary condition for the existence of a stress function $\boldsymbol{\Psi}$ such that $\boldsymbol{P}=\mathbf{curl}\,\boldsymbol{\Psi}$. This holds whenever $\boldsymbol{U}\in H^{1}(T\mathcal{B})$, $\boldsymbol{K}\in\operatorname{ker}(\mathbf{curl}\,)\subset H^{\mathbf{c}}(\mathcal{B})$, and $\boldsymbol{P}\in\operatorname{ker}(\mathbf{div})\subset H^{\mathbf{d}}(\mathcal{B})$. The deformation gradient is defined as $\boldsymbol{F}:=\boldsymbol{I}+\boldsymbol{K}$, where $\boldsymbol{I}$ is the identity tensor, and $J:=\operatorname{det}\boldsymbol{F}$ (in Cartesian coordinates for both the reference and current configurations). One can show that $dv=JdV$, where $dV$ and $dv$ are the volume elements of the undeformed and deformed configurations, respectively. For incompressible solids, $J=1$. To weakly impose $J-1=0$, one considers a Lagrange multiplier $p$ as an independent field variable, which physically is realized as a pressure-like variable. At the discrete level, the restriction of $J$ to an element is a scaler describing the change of volume of that element [22]. Hence, one can assume that discrete pressure $p$ is also defined on each element, and in general, is not continuous across the element interfaces. Therefore, as a discontinuous scalar-valued field, $p\in L^{2}(\mathcal{B})$.

2.2 Mixed Formulations

Let $\rho_{0}$ be the mass density of the body $\mathcal{B}$ and $\boldsymbol{B}$ be the body force per unit mass. Assume that the boundary of the body is a disjoint union of two subsets $\partial\mathcal{B}=\Gamma_{d}\sqcup\Gamma_{t}$ and is subjected to the displacement boundary condition \boldsymbol{U}\big{|}_{\Gamma_{d}}=\overline{\boldsymbol{U}} and the traction boundary condition (\boldsymbol{P}\boldsymbol{N})\big{|}_{\Gamma_{t}}=\overline{\boldsymbol{T}}, where $\boldsymbol{N}$ is the unit outward normal vector field of $\partial\mathcal{B}$ in the reference configuration. Also, define $H^{1}(T\mathcal{B},\Gamma_{d},\overline{\boldsymbol{U}}):=\left\{\boldsymbol{U}\in H^{1}(T\mathcal{B}):\boldsymbol{U}|_{\Gamma_{d}}=\overline{\boldsymbol{U}}\right\}$ and $H^{1}(T\mathcal{B},\Gamma_{d}):=H^{1}(T\mathcal{B},\Gamma_{d},\boldsymbol{0})$, where $\overline{\boldsymbol{U}}$ is of $H^{1/2}$-class. Suppose $\langle,\rangle$ is the standard inner product of $\mathbb{R}^{3}$ and let $\llangle,\rrangle$ stand for the $L^{2}$-inner products of scalar, vector, and tensor fields, that is, $\llangle f,g\rrangle:=\int_{\mathcal{B}}fg\,dV$, $\llangle\boldsymbol{Y},\boldsymbol{Z}\rrangle:=\int_{\mathcal{B}}Y^{I}Z^{I}dV$, and $\llangle\boldsymbol{S},\boldsymbol{T}\rrangle:=\int_{\mathcal{B}}S^{IJ}T^{IJ}dV$. Then, one can define a Hu-Washizu-type functional $\mathcal{I}:H^{1}(T\mathcal{B},\Gamma_{d},\overline{\boldsymbol{U}})\times H^{\mathbf{c}}(\mathcal{B})\times H^{\mathbf{d}}(\mathcal{B})=:\mathcal{D}\rightarrow\mathbb{R}$ as

[TABLE]

where $W(\mathbf{X},\boldsymbol{K})$ is the stored energy function of a hyperelastic material. For an isotropic solid, the energy function can be written as $W=\widehat{W}(\mathbf{X},I_{1},I_{2},I_{3})$, where $I_{1}=\operatorname{tr}\boldsymbol{C}$, $I_{2}=\frac{1}{2}[(\operatorname{tr}\boldsymbol{C})^{2}-\operatorname{tr}\boldsymbol{C}^{2}]$, and $I_{3}=\operatorname{det}\boldsymbol{C}$ are the invariants of the right Cauchy-Green deformation tensor $\boldsymbol{C}=\boldsymbol{F}^{\mathsf{T}}\boldsymbol{F}$. For incompressible solids, $J=\sqrt{I_{3}}=1$, and one modifies (2.1) by defining

[TABLE]

where $C:\mathbb{R}^{+}\rightarrow\mathbb{R}$ is a smooth function such that $C(J)=0$ if and only if $J=1$ and $p\in L^{2}(\mathcal{B})$ is a pressure-like scalar field, to which we may refer simply as pressure. For 3D computations, in order to improve the stability of the mixed finite element methods, a stabilizing term is added to (2.2) as

[TABLE]

where $\alpha\geq 0$ is a penalty constant for enforcing $\boldsymbol{K}=\mathbf{grad}\,\boldsymbol{U}$. Extremizing (2.3), as discussed in [2, §2.2], results in the following weak formulation of the boundary-value problem for incompressible nonlinear elastostatics:

Given a body force $\boldsymbol{B}$ of $L^{2}$-class, a boundary displacement $\overline{\boldsymbol{U}}$ on $\Gamma_{d}$ of $H^{1/2}$-class, a boundary traction $\overline{\boldsymbol{T}}$ on $\Gamma_{t}$ of $L^{2}$-class, and a stability constant $\alpha\geq 0$, find $(\boldsymbol{U},\boldsymbol{K},\boldsymbol{P},p)\in H^{1}(T\mathcal{B},\Gamma_{d},\overline{\boldsymbol{U}})\times H^{\mathbf{c}}(\mathcal{B})\times H^{\mathbf{d}}(\mathcal{B})\times L^{2}(\mathcal{B})$ such that

[TABLE]

where

[TABLE]

and

[TABLE]

In (2.4)2, $\widetilde{\boldsymbol{P}}(\boldsymbol{K})={\partial\widetilde{W}}/{\partial\boldsymbol{K}}$ with \widetilde{W}=\widehat{W}(\mathbf{X},I_{1},I_{2},I_{3})\big{|}_{I_{3}=1} is the constitutive part of the stress, and $\boldsymbol{Q}(\boldsymbol{K})={\partial C}/{\partial\boldsymbol{K}}=C^{\prime}(J)(\boldsymbol{F}^{-1})^{\mathsf{T}}$ is the contribution of the incompressibility constraint $J=1$. Note that setting $\alpha=0$ results in the standard weak formulation of incompressible nonlinear elastostatics [2, §2.2]. The solutions of the above weak formulation are the critical points of the functional (2.3). Using Green’s formula $\llangle\mathbf{div}\boldsymbol{P},\boldsymbol{\Upsilon}\rrangle=-\llangle\boldsymbol{P},\mathbf{grad}\boldsymbol{\Upsilon}\rrangle+\int_{\partial\mathcal{B}}\langle\boldsymbol{P}\boldsymbol{N},\boldsymbol{\Upsilon}\rangle\,dA,~{}\forall\boldsymbol{\Upsilon}\in H^{1}(T\mathcal{B},\Gamma_{d})$ and assuming that $\int_{\partial\mathcal{B}}\langle\boldsymbol{P}\boldsymbol{N},\boldsymbol{\Upsilon}\rangle\,dA=\int_{\Gamma_{t}}\langle\overline{\boldsymbol{T}},\boldsymbol{\Upsilon}\rangle\,dA$ holds, $\forall\boldsymbol{\Upsilon}\in H^{1}(T\mathcal{B},\Gamma_{d})$, one can show that (2.4) results in the following set of governing equations for incompressible nonlinear elastostatics:

[TABLE]

Conversely, one can obtain (2.4) from (2.7), see [1, §2.2]. Note that (2.7) is the constitutive relation of an incompressible solid, which in terms of the Cauchy stress reads $\boldsymbol{\sigma}=\widetilde{\boldsymbol{P}}(\boldsymbol{K})\boldsymbol{F}^{\mathsf{T}}+\bar{p}\boldsymbol{I}$, where $\bar{p}=p\,C^{\prime}(J)$. Note that adding the stabilizing terms (2.6) to the weak formulation (2.4) does not change the set of governing equations (2.7). In other words, these terms will vanish for the exact solutions of (2.4). Hence, with proper discretization, the extra terms (2.6) may improve the stability of the resulting mixed finite element methods without compromising their consistency (see [23] for consistency and stability). We discuss this further in §3.5.

By setting $p=q=0$ in (2.4) and replacing $\widetilde{\boldsymbol{P}}(\boldsymbol{K})$ with $\widehat{\boldsymbol{P}}(\boldsymbol{K})={\partial\widehat{W}}/{\partial\boldsymbol{K}}$, one can readily arrive at the following weak formulation of the boundary-value problem of compressible nonlinear elastostatics:

Given $\boldsymbol{B}$, $\overline{\boldsymbol{U}}$, $\overline{\boldsymbol{T}}$, and $\alpha\geq 0$, find $(\boldsymbol{U},\boldsymbol{K},\boldsymbol{P})\in H^{1}(T\mathcal{B},\Gamma_{d},\overline{\boldsymbol{U}})\times H^{\mathbf{c}}(\mathcal{B})\times H^{\mathbf{d}}(\mathcal{B})$ such that

[TABLE]

Similarly, one can show that (2.8) results in the following set of governing equations for compressible nonlinear elastostatics:

[TABLE]

3 Finite Element Approximations

3.1 Finite Elements

Suppose $\mathcal{P}_{r}(\mathbb{R}^{3})$ is the space of $\mathbb{R}$-valued polynomials in three variables $\{\mathrm{X}^{1},\mathrm{X}^{2},\mathrm{X}^{3}\}$ of degree at most $r\geq 0$ and suppose $\mathcal{H}_{r}(\mathbb{R}^{3})\subset\mathcal{P}_{r}(\mathbb{R}^{3})$ is the space of homogeneous polynomials of degree $r$, that is, all the terms of the members of $\mathcal{H}_{r}(\mathbb{R}^{3})$ are of degree $r$. For $r<0$, these spaces are assumed to be empty. By $\mathcal{P}_{r}(T\mathbb{R}^{3})$ and \mathcal{P}_{r}(\text{\large\otimes}^{2}T\mathbb{R}^{3}) we denote the spaces of polynomial vector and $\binom{2}{0}$-tensor fields in $\mathbb{R}^{3}$ with Cartesian components in $\mathcal{P}_{r}(\mathbb{R}^{3})$. The spaces $\mathcal{H}_{r}(T\mathbb{R}^{3})$ and \mathcal{H}_{r}(\text{\large\otimes}^{2}T\mathbb{R}^{3}) are defined similarly. Next define the following subspaces of $\mathcal{P}_{r}(T\mathbb{R}^{3})$:

[TABLE]

where $(\mathrm{L}_{1}(\boldsymbol{Y}))^{I}=\varepsilon_{IJL}\mathrm{X}^{L}Y^{J}$ for any vector field $\boldsymbol{Y}$, and $(\mathrm{L}_{2}(f))^{I}=\mathrm{X}^{I}f$ for any scalar field $f$. Similarly, one defines the following subspaces of \mathcal{P}_{r}(\text{\large\otimes}^{2}T\mathbb{R}^{3}):

[TABLE]

where $(\mathbf{L}_{1}(\boldsymbol{T}))^{IJ}=\varepsilon_{JLK}\mathrm{X}^{K}T^{IL}$ for any $\binom{2}{0}$-tensor field $\boldsymbol{T}$, and $(\mathbf{L}_{2}(\boldsymbol{Y}))^{IJ}=\mathrm{X}^{J}Y^{I}$ for any vector field $\boldsymbol{Y}$. One can show that

[TABLE]

Let $\widehat{\mathcal{T}}$ be a reference tetrahedral element with coordinates $\boldsymbol{\xi}=(\xi^{1},\xi^{2},\xi^{3})$ shown in Figure 1. We denote the edges of $\widehat{\mathcal{T}}$ by $\widehat{\mathcal{E}}_{i},i=1,2,\dots,6$ and their corresponding lengths by $\hat{\ell}_{i},i=1,2,\dots,6$, and the faces of $\widehat{\mathcal{T}}$ by $\widehat{\mathcal{F}}_{i},i=1,2,3,4$ and their corresponding areas by $\hat{A}_{i},i=1,2,3,4$. For an edge joining two vertices $i$ and $j$, one defines a unique orientation as $i\rightarrow j$, where $i<j$. We also define a unit tangent vector $\hat{\boldsymbol{\mathsf{t}}}_{i}$ on each edge such that it agrees with the edge orientation. Moreover, on each face containing three edges $\widehat{\mathcal{E}}_{i}$, $\widehat{\mathcal{E}}_{j}$, and $\widehat{\mathcal{E}}_{k}$, we define a unit normal vector $\hat{\boldsymbol{\mathsf{n}}}_{l}=\hat{\boldsymbol{\mathsf{t}}}_{i}\times\hat{\boldsymbol{\mathsf{t}}}_{j}$, where $i<j<k$.

Following [24, 25], we define a finite element as a triplet $(\mathcal{T},\mathcal{P}(\mathcal{T}),\text{{\textSigma}})$, where $\mathcal{T}$ is a tetrahedron in $\mathbb{R}^{3}$, $\mathcal{P}(\mathcal{T})$ is a space of polynomials on $\mathcal{T}$, and \textSigma is a set of $\mathbb{R}$-valued linear functionals acting on the members of $\mathcal{P}(\mathcal{T})$. The members of \textSigma are called the local degrees of freedom (DOF) and the local shape functions form a basis for $\mathcal{P}(\mathcal{T})$ (see [2, §3.1]). We consider the following reference finite elements for the four field variables:

[TABLE]

Note that \mathcal{P}_{r}(\widehat{\mathcal{T}})=\mathcal{P}_{r}(\mathbb{R}^{3})\big{|}_{\widehat{\mathcal{T}}} and $\mathcal{P}_{r}(T\widehat{\mathcal{T}})$, \mathcal{P}_{r}(\text{\large\otimes}^{2}T\widehat{\mathcal{T}}), \mathcal{P}^{-}_{r}(\text{\large\otimes}^{2}T\widehat{\mathcal{T}}), and \mathcal{P}^{\ominus}_{r}(\text{\large\otimes}^{2}T\widehat{\mathcal{T}}) are defined similarly. The finite element for $\boldsymbol{U}$ is based on the standard Lagrange finite elements. For a vector field $\boldsymbol{V}:\widehat{\mathcal{T}}\rightarrow\mathbb{R}^{3}$, the set of local degrees of freedom is $\text{{\textSigma}}^{\widehat{\mathcal{T}},1}=\{V^{1}(\boldsymbol{\xi}_{1}),V^{2}(\boldsymbol{\xi}_{1}),V^{3}(\boldsymbol{\xi}_{1}),\dots,V^{1}(\boldsymbol{\xi}_{m}),V^{2}(\boldsymbol{\xi}_{m}),V^{3}(\boldsymbol{\xi}_{m})\}$, where $\boldsymbol{\xi}_{i}$ contains the coordinates of the $i$-th node of the $m$-node $\widehat{\mathcal{T}}$, where $m=4$ ($m=10$) for $r=1$ ($r=2$). For $r=1,2$, a basis of the polynomial space $\mathcal{P}_{r}(T\widehat{\mathcal{T}})$ includes

[TABLE]

The Lagrange polynomials $l^{r}_{i}$ for the four-node reference tetrahedron $\widehat{\mathcal{T}}$ are

[TABLE]

For the ten-node $\widehat{\mathcal{T}}$, the Lagrange polynomials are $l^{2}_{i}=l^{1}_{i}(2l^{1}_{i}-1)$ for the nodes at the vertices $i=1,2,3,4$ and $l^{2}_{k}=4l^{1}_{i}l^{1}_{j}$ for the middle node of each edge joining vertices $i$ and $j$ as shown in Figure 1. We will use $\mathcal{P}_{r}(T\widehat{\mathcal{T}}),r=1,2$ spanned by $\boldsymbol{h}^{\widehat{\mathcal{T}}}_{l},l=1,2,...,3m$ to construct the approximation space of $\boldsymbol{U}$. To interpolate $\boldsymbol{K}\in H^{\mathbf{c}}$, we define two finite elements given in (3.1)2 based on the Nédélec $1^{\mathrm{st}}$-kind edge elements in $\mathbb{R}^{3}$ (NE1) [26] and the Nédélec $2^{\mathrm{nd}}$-kind edge elements in $\mathbb{R}^{3}$ (NE2) [27], respectively. Let $\overrightarrow{\boldsymbol{T}}_{\!\!I}:=\begin{bmatrix}T^{I1}&T^{I2}&T^{I3}\end{bmatrix}^{\mathsf{T}}$ be a vector containing the elements of the $I$-th row of a $\binom{2}{0}$-tensor $\boldsymbol{T}$. The set of the local degrees of freedom $\text{{\textSigma}}^{\widehat{\mathcal{T}},\mathbf{c}-}$ ($\text{{\textSigma}}^{\widehat{\mathcal{T}},\mathbf{c}}$) in (3.1)2 is defined as $\left\{\phi^{\widehat{\mathcal{T}},\widehat{\mathcal{E}}_{k}}_{I,J},\phi^{\widehat{\mathcal{T}},\widehat{\mathcal{F}}_{l}}_{I,J},\phi^{\widehat{\mathcal{T}},\widehat{\mathcal{T}}}_{I,J}\right\}$, where

[TABLE]

We next discuss the corresponding local shape functions of the two finite elements given in (3.1)2. Let $\boldsymbol{v}_{J}^{\widehat{\mathcal{T}},\widehat{\mathcal{E}}_{k}}$, $\boldsymbol{v}_{J}^{\widehat{\mathcal{T}},\widehat{\mathcal{F}}_{l}}$, and $\boldsymbol{v}_{J}^{\widehat{\mathcal{T}},\widehat{\mathcal{T}}}$ denote the shape functions of those Nédélec edge elements that are associated to the $k-$th edge of $\widehat{\mathcal{T}}$, the $l-$th face of $\widehat{\mathcal{T}}$, and the entire $\widehat{\mathcal{T}}$, respectively. We consider these vector-valued polynomials in $\mathbb{R}^{3}$ as row vectors and define the following tensorial shape functions:

[TABLE]

Similarly, we define $\boldsymbol{r}_{I,J}^{\widehat{\mathcal{T}},\widehat{\mathcal{F}}_{l}}$ and $\boldsymbol{r}_{I,J}^{\widehat{\mathcal{T}},\widehat{\mathcal{T}}}$ for $I=1,2,3$ using $\boldsymbol{v}_{J}^{\widehat{\mathcal{T}},\widehat{\mathcal{F}}_{l}}$ and $\boldsymbol{v}_{J}^{\widehat{\mathcal{T}},\widehat{\mathcal{T}}}$, respectively. The polynomial spaces \mathcal{P}^{-}_{r}(\text{\large\otimes}^{2}T\widehat{\mathcal{T}}) and \mathcal{P}_{r}(\text{\large\otimes}^{2}T\widehat{\mathcal{T}}) in (3.1)2 are spanned by a basis $\left\{\boldsymbol{r}_{I,J}^{\widehat{\mathcal{T}},\widehat{\mathcal{E}}_{k}},\boldsymbol{r}_{I,J}^{\widehat{\mathcal{T}},\widehat{\mathcal{F}}_{l}},\boldsymbol{r}_{I,J}^{\widehat{\mathcal{T}},\widehat{\mathcal{T}}}\right\}$ that are respectively based on the shape functions of NE1 and NE2. Moreover, the following relations hold:

[TABLE]

Later in this section, we will provide explicit expressions for some of the above shape functions that are used in our numerical examples. To interpolate $\boldsymbol{P}\in H^{\mathbf{d}}$, we define the two finite elements given in (3.1)3 based on respectively the Nédélec $1^{\mathrm{st}}$-kind face elements in $\mathbb{R}^{3}$ (NF1) [26], and the Nédélec $2^{\mathrm{nd}}$-kind face elements in $\mathbb{R}^{3}$ (NF2) [27]. We denote the set of the local degrees of freedom $\text{{\textSigma}}^{\widehat{\mathcal{T}},\mathbf{d}-}$ ($\text{{\textSigma}}^{\widehat{\mathcal{T}},\mathbf{d}}$) by $\left\{\psi^{\widehat{\mathcal{T}},\widehat{\mathcal{F}}_{l}}_{I,J},\psi^{\widehat{\mathcal{T}},\widehat{\mathcal{T}}}_{I,J}\right\}$, where

[TABLE]

We denote the set of the local shape functions of (3.1)3 by $\left\{\boldsymbol{s}_{I,J}^{\widehat{\mathcal{T}},\widehat{\mathcal{F}}_{l}},\boldsymbol{s}_{I,J}^{\widehat{\mathcal{T}},\widehat{\mathcal{T}}}\right\}$. Both $\boldsymbol{s}_{I,J}^{\widehat{\mathcal{T}},\widehat{\mathcal{F}}_{l}}$ and $\boldsymbol{s}_{I,J}^{\widehat{\mathcal{T}},\widehat{\mathcal{T}}}$ are defined similar to (3.5) but using the vector-valued shape functions of Nédélec face elements, which we denote by $\boldsymbol{u}_{J}^{\widehat{\mathcal{T}},\widehat{\mathcal{F}}_{l}}$ and $\boldsymbol{u}_{J}^{\widehat{\mathcal{T}},\widehat{\mathcal{T}}}$. Also, one has

[TABLE]

For the reference finite element of pressure (3.1)4, we have $\text{{\textSigma}}^{\widehat{\mathcal{T}},\ell}=\{\omega^{\widehat{\mathcal{T}}}_{i}\}$, where $\omega^{\widehat{\mathcal{T}}}_{i}(f)=\frac{1}{\hat{A}}\int_{\widehat{\mathcal{T}}}p_{i}f\,d\hat{V}$ for all the polynomials $p_{i}$ that form a basis for \mathcal{P}_{r}(\mathbb{R}^{3})\big{|}_{\widehat{\mathcal{T}}}. Also, the set of local shape functions \big{\{}t^{\widehat{\mathcal{T}}}_{i}\big{\}}, which spans $\mathcal{P}_{r}(\widehat{\mathcal{T}})$, is $\{1\}$ for $r=0$, $\{1,\xi^{1},\xi^{2},\xi^{3}\}$ for $r=1$, and $\{1,\xi^{1},\xi^{2},(\xi^{1})^{2},(\xi^{2})^{2},(\xi^{3})^{2},\xi^{1}\xi^{2},\xi^{1}\xi^{3},\xi^{2}\xi^{3}\}$ for $r=2$. Choosing $p_{i}$ properly, one can show that $\omega^{\widehat{\mathcal{T}}}_{i}(t_{j}^{\widehat{\mathcal{T}}})=\delta_{ij}$.

To extend 2D CSFEMs [2] to 3D, we first followed the same approach we had proposed in [2] and considered $r=1,2$ for the finite elements of $\boldsymbol{U}$, $\boldsymbol{K}$, and $\boldsymbol{P}$ and $r=0,1,2$ for the finite elements of $p$ in (3.1). This provides $96$ combinations of elements for discretizing the boundary-value problem (2.4). Using the matrix formulation of the linearization of (2.4) for $\alpha=0$ using the approach discussed in [2, §3.5], we concluded that all the $96$ choices lead to strictly singular or unstable methods in 3D. We will discuss this further in §3.5. To overcome this singularity issue, we modify a suitable combination of elements among the aforementioned unstable $96$ choices and propose the following convergent finite elements for $\boldsymbol{U}$, $\boldsymbol{K}$,$\boldsymbol{P}$, and $p$:

[TABLE]

where \overline{\mathcal{P}}_{3}(\text{\large\otimes}^{2}T\widehat{\mathcal{T}}):=\mathcal{P}_{1}(\text{\large\otimes}^{2}T\widehat{\mathcal{T}})\oplus\text{span}\left\{\boldsymbol{r}^{\widehat{\mathcal{T}},\widehat{\mathcal{T}}}_{I,J}\right\}_{I,J=1,2,3} for \boldsymbol{r}^{\widehat{\mathcal{T}},\widehat{\mathcal{T}}}_{I,J}\in\mathcal{P}^{-}_{3}(\text{\large\otimes}^{2}T\widehat{\mathcal{T}}), and $\overline{\text{{\textSigma}}}^{\widehat{\mathcal{T}},\mathbf{c}}$ is the union of $\text{{\textSigma}}^{\widehat{\mathcal{T}},\mathbf{c}}$ for $r=1$ and $\left\{\phi^{\widehat{\mathcal{T}},\widehat{\mathcal{T}}}_{I,J}\right\}_{I,J=1,2,3}\subset\text{{\textSigma}}^{\widehat{\mathcal{T}},\mathbf{c}-}$ for $r=3$. The schematic diagram of (3.9) is illustrated in Figure 2. Moreover, the shape functions for the finite element of $\boldsymbol{U}$ in (3.9) are given by (3.2) for $r=2$ and $m=10$, and the shape function for the finite element of $p$ in (3.9) is simply $t^{\widehat{\mathcal{T}}}=1$. We use the results of Arnold et al. [28] to provide the explicit expression for the shape functions for the finite elements of $\boldsymbol{K}$ and $\boldsymbol{P}$. The finite element of $\boldsymbol{K}$ in (3.9) has $6$ shape functions associated to each edge $\widehat{\mathcal{E}}_{k}$ of $\widehat{\mathcal{T}}$ and $9$ shape functions associated to $\widehat{\mathcal{T}}$. Let us ignore the superscript of $l^{1}_{i},i=1,2,3,4$ in (3.3) and consider $\nabla l_{i}=\begin{bmatrix}\partial l_{i}/\partial\xi^{1}&\partial l_{i}/\partial\xi^{2}&\partial l_{i}/\partial\xi^{3}\end{bmatrix}$ as a row vector. Then, for an edge $\widehat{\mathcal{E}}_{k}$ joining two vertices $i$ and $j$ as shown in Figure 1, the $6$ shape functions $\boldsymbol{r}_{I,2}^{\widehat{\mathcal{T}},\widehat{\mathcal{E}}_{k}},\boldsymbol{r}_{I,3}^{\widehat{\mathcal{T}},\widehat{\mathcal{E}}_{k}},I=1,2,3$ are obtained using (3.5) and the following vector-valued shape functions for NE2 of order $1$ [28]:

[TABLE]

The $9$ remaining shape functions $\boldsymbol{r}^{\widehat{\mathcal{T}},\widehat{\mathcal{T}}}_{I,1},\boldsymbol{r}^{\widehat{\mathcal{T}},\widehat{\mathcal{T}}}_{I,2},\boldsymbol{r}^{\widehat{\mathcal{T}},\widehat{\mathcal{T}}}_{I,3},I=1,2,3$ are obtained similar to (3.5) and using the following vector-valued shape functions for NE1 of order $3$ [28]:

[TABLE]

where $\boldsymbol{w}_{ij}=l_{i}\nabla l_{j}-l_{j}\nabla l_{i}$. The finite element of $\boldsymbol{P}$ in (3.9) has $3$ shape functions $\boldsymbol{s}_{I}^{\widehat{\mathcal{T}},\widehat{\mathcal{F}}_{l}},I=1,2,3$ associated to each face $\widehat{\mathcal{F}}_{l}$ of $\widehat{\mathcal{T}}$ that contains the three vertices $i$, $j$, and $k$ according to Figure 1. These shape functions are obtained similar to (3.5) and using the following vector-valued shape function for NF1 of order $1$[28]:

[TABLE]

Next we explain how to calculate the shape functions in an arbitrary element in a mesh using the shape functions of the reference finite elements. Let $\mathcal{B}_{h}$ be a triangulation of the reference configuration $\mathcal{B}$ consisting of arbitrary tetrahedra $\mathcal{T}$ such that the intersection of any two distinct tetrahedra is either empty or a common face/edge/vertex of each. The discretization parameter $h$ is defined as $h:={\max}\,\mathrm{diam}\,\mathcal{T},\,\forall{\mathcal{T}\in\mathcal{B}_{h}}$. We define a local ordering for vertices of each $\mathcal{T}\in\mathcal{B}_{h}$ by assigning the numbers $1,2,3,4$ to them. We then denote the Cartesian coordinates of the $i$-th vertex of $\mathcal{T}$ by a column vector $\boldsymbol{\mathsf{X}}^{\mathcal{T}}_{i}=[\mathsf{X}^{1,\mathcal{T}}_{i}~{}\mathsf{X}^{2,\mathcal{T}}_{i}~{}\mathsf{X}^{3,\mathcal{T}}_{i}]^{\mathsf{T}}$ and define the following affine transformation:

[TABLE]

where $\boldsymbol{\mathsf{J}}_{\mathcal{T}}=[\boldsymbol{\mathsf{X}}^{\mathcal{T}}_{2}-\boldsymbol{\mathsf{X}}^{\mathcal{T}}_{1}~{}\boldsymbol{\mathsf{X}}^{\mathcal{T}}_{3}-\boldsymbol{\mathsf{X}}^{\mathcal{T}}_{1}~{}\boldsymbol{\mathsf{X}}^{\mathcal{T}}_{4}-\boldsymbol{\mathsf{X}}^{\mathcal{T}}_{1}]_{3\times 3}$. Note that $\boldsymbol{\mathsf{T}}_{\mathcal{T}}$ is bijective and $\boldsymbol{\mathsf{J}}_{\mathcal{T}}$ is invertible. For an element $\mathcal{T}\in\mathcal{B}_{h}$, we denote its edges by $\mathcal{E}^{\mathcal{T}}_{i}=\boldsymbol{\mathsf{T}}_{\mathcal{T}}(\widehat{\mathcal{E}}_{i})$, $i=1,2,\dots,6$, and its faces by $\mathcal{F}^{\mathcal{T}}_{i}=\boldsymbol{\mathsf{T}}_{\mathcal{T}}(\widehat{\mathcal{F}}_{i})$, $i=1,2,3,4$. We assume that $\mathcal{E}^{\mathcal{T}}_{i}$ inherits the orientation of its reference counterpart $\widehat{\mathcal{E}}_{i}$. Moreover, the tangent vector $\boldsymbol{\mathsf{t}}_{i}$ defined on $\mathcal{E}^{\mathcal{T}}_{i}$ inherits the orientation of $\mathcal{E}^{\mathcal{T}}_{i}$, and the normal vector on $\mathcal{F}^{\mathcal{T}}_{l}$ containing three edges $\mathcal{E}^{\mathcal{T}}_{i}$, $\mathcal{E}^{\mathcal{T}}_{j}$, and $\mathcal{E}^{\mathcal{T}}_{k}$ such that $i<j<k$ is defined as $\boldsymbol{\mathsf{n}}_{l}=\boldsymbol{\mathsf{t}}_{i}\times\boldsymbol{\mathsf{t}}_{j}$. One can show that $\boldsymbol{\mathsf{t}}_{i}=\boldsymbol{\mathsf{J}}_{\mathcal{T}}\hat{\boldsymbol{\mathsf{t}}}_{i}$ and $\boldsymbol{\mathsf{n}}_{i}=\mathrm{det}\,\boldsymbol{\mathsf{J}}_{\mathcal{T}}\boldsymbol{\mathsf{J}}_{\mathcal{T}}^{\mathsf{-T}}\hat{\boldsymbol{\mathsf{n}}}_{i}$. For efficient assembly of the finite elements of $\boldsymbol{K}\in H^{\mathbf{c}}$ and $\boldsymbol{P}\in H^{\mathbf{d}}$, we use the numbering scheme discussed in [29]. Using this scheme, one first assumes that every vertex in a mesh $\mathcal{B}_{h}$ has a distinct global number and then the local ordering of four vertices of every tetrahedron in that mesh $\mathcal{T}\in\mathcal{B}_{h}$ agree with the ascending order of the global numbers of its four vertices. Considering the edge orientations of the reference element shown in Figure 1, the orientation of every edge in the mesh is from a vertex with a smaller global number to a vertex with a larger global number. The advantage of this scheme is that the orientation of a common edge between elements in a mesh is uniquely defined and is identical to that of the edge in any of those elements. It follows that some elements in a mesh sharing a common edge have an identical tangent vector on that edge, and any two elements with a common face have an identical normal vector on that face. For an illustration of this, see [29, Figue 5.2]. Note that using this scheme, the normal vectors of some of the exterior faces of the mesh are not pointed outward, and $o_{\mathcal{T}}=\mathrm{sign}\left(\mathrm{det}\,\boldsymbol{\mathsf{J}}_{\mathcal{T}}\right)$ can be either $1$ or $-1$.

Consider the following mappings:

[TABLE]

where $\boldsymbol{\mathsf{T}}^{\mathbf{c}}_{\mathcal{T}}$ and $\boldsymbol{\mathsf{T}}^{\mathbf{d}}_{\mathcal{T}}$ are known as the Piola transforms. For a $\binom{2}{0}$-tensor $\boldsymbol{T}$, one calculates the Piola transforms separately for each row:

[TABLE]

and $\boldsymbol{\mathsf{T}}^{\mathbf{d}}_{\mathcal{T}}(\boldsymbol{T})$ is calculated similarly. Using [2, Proposition 8], (3.11), and the local shape functions in the reference element $\widehat{\mathcal{T}}$, one can obtain the local shape functions in any element $\mathcal{T}\in\mathcal{B}_{h}$ enabling one to locally interpolate the four field variables $(\boldsymbol{U},\boldsymbol{K},\boldsymbol{P},p)$ over that element. In particular, the local shape functions for $\boldsymbol{U}$ are obtained as $\boldsymbol{h}^{\mathcal{T}}_{k}=\boldsymbol{\mathsf{T}}^{1}_{\mathcal{T}}\left(\boldsymbol{h}^{\widehat{\mathcal{T}}}_{k}\right)$; the local shape functions for $\boldsymbol{K}$ are $\boldsymbol{r}_{I,J}^{\mathcal{T},\mathcal{E}_{k}}=\boldsymbol{\mathsf{T}}^{\mathbf{c}}_{\mathcal{T}}\left(\boldsymbol{r}_{I,J}^{\widehat{\mathcal{T}},\widehat{\mathcal{E}}_{k}}\right)$, $\boldsymbol{r}_{I,J}^{\mathcal{T},\mathcal{F}_{l}}=\boldsymbol{\mathsf{T}}^{\mathbf{c}}_{\mathcal{T}}\left(\boldsymbol{r}_{I,J}^{\widehat{\mathcal{T}},\widehat{\mathcal{F}}_{l}}\right)$, and $\boldsymbol{r}_{I,J}^{\mathcal{T},\mathcal{T}}=\boldsymbol{\mathsf{T}}^{\mathbf{c}}_{\mathcal{T}}\left(\boldsymbol{r}_{I,J}^{\widehat{\mathcal{T}},\widehat{\mathcal{T}}}\right)$; the local shape functions for $\boldsymbol{P}$ are $\boldsymbol{s}_{I,J}^{\mathcal{T},\mathcal{F}_{l}}=\boldsymbol{\mathsf{T}}^{\mathbf{d}}_{\mathcal{T}}\left(\boldsymbol{s}_{I,J}^{\widehat{\mathcal{T}},\widehat{\mathcal{F}}_{l}}\right)$ and $\boldsymbol{s}_{I,J}^{\mathcal{T},{\mathcal{T}}}=\boldsymbol{\mathsf{T}}^{\mathbf{d}}_{\mathcal{T}}\left(\boldsymbol{s}_{I,J}^{\widehat{\mathcal{T}},\widehat{\mathcal{T}}}\right)$; and $t_{i}^{\mathcal{T}}=\boldsymbol{\mathsf{T}}^{\ell}_{\mathcal{T}}\left(t_{i}^{\widehat{\mathcal{T}}}\right)$ gives the local shape functions for $p$. Using [2, Proposition 8] and (3.11), one can also obtain the local degrees of freedom for the finite elements of any element $\mathcal{T}\in\mathcal{B}_{h}$. For example, considering (3.7), $\psi^{{\mathcal{T}},{\mathcal{F}}_{l}}_{I,J}(\boldsymbol{T})=\left(\psi^{\widehat{\mathcal{T}},\widehat{\mathcal{F}}_{l}}_{I,J}\circ{\boldsymbol{\mathsf{T}}^{\mathbf{d}}_{\mathcal{T}}}^{-1}\right)(\boldsymbol{T})$ and $\psi^{{\mathcal{T}},\mathcal{T}}_{I,J}(\boldsymbol{T})=\left(\psi^{\widehat{\mathcal{T}},\widehat{\mathcal{T}}}_{I,J}\circ{\boldsymbol{\mathsf{T}}^{\mathbf{d}}_{\mathcal{T}}}^{-1}\right)(\boldsymbol{T})$ are the degrees of freedom for the finite element of $\mathcal{T}$ that we use for $\boldsymbol{P}$. The other degrees of freedom can be written similarly using their reference counterparts. In this work, the traction boundary conditions are imposed weakly through (2.5); one does not need to impose them directly by calculating the related degrees of freedom on the boundary of the mesh. Thus, in practice, all degrees of freedom, even those on the boundary of the mesh, are obtained by solving the final discrete system; calculating their explicit expressions are not required.

3.2 Finite Element Spaces

Next, some conforming finite element spaces are introduced in order to discretize (2.4) and (2.8). Let $\mathcal{F}_{h}^{i}$ be the set of all interior faces of a 3D mesh $\mathcal{B}_{h}$. Given a face $\mathcal{F}\in\mathcal{F}_{h}^{i}$, there are two elements $\mathcal{T},\mathcal{T}^{\prime}\in\mathcal{B}_{h}$ such that $\mathcal{F}=\mathcal{T}\cap\mathcal{T}^{\prime}$. Suppose $\boldsymbol{V}$ is a vector-valued function and $\boldsymbol{T}$ is a tensor-valued function both defined on $\mathcal{B}_{h}$ with limits on both sides of $\forall\mathcal{F}\in\mathcal{F}_{h}^{i}$. We define the following notions of jump across a face $\mathcal{F}\in\mathcal{F}_{h}^{i}$:

[TABLE]

where $\boldsymbol{V}_{\mathcal{T}}:=\boldsymbol{V}|_{\mathcal{T}}$, $\boldsymbol{T}_{\mathcal{T}}:=\boldsymbol{T}|_{\mathcal{T}}$, and $\boldsymbol{V}_{\mathcal{T^{\prime}}}$ and $\boldsymbol{T}_{\mathcal{T}^{\prime}}$ are defined similarly. $\boldsymbol{\mathsf{t}}$ ($\boldsymbol{\mathsf{n}}$) is a unit vector tangent (normal) to $\mathcal{F}$. We write $\llbracket\mathsf{t}\boldsymbol{T}\rrbracket_{\mathcal{F}}=\boldsymbol{0}$ ($\llbracket\mathsf{n}\boldsymbol{T}\rrbracket_{\mathcal{F}}=\boldsymbol{0}$), if the jump is zero for any unit vector $\boldsymbol{\mathsf{t}}$ ($\boldsymbol{\mathsf{n}}$) on $\mathcal{F}$. Note that all the above jumps are vector-valued functions in 3D. Consider the following finite element spaces:

[TABLE]

Note that the above sapces are conforming, i.e., $V_{h,r}^{1}\subset H^{1}(T\mathcal{B}_{h})$, $V_{h,r}^{\mathbf{c}{-}}\subset V_{h,r}^{\mathbf{c}}\subset H^{\mathbf{c}}(\mathcal{B}_{h})$, $V_{h,r}^{\mathbf{d}{-}}\subset V_{h,r}^{\mathbf{d}}\subset H^{\mathbf{d}}(\mathcal{B}_{h})$, and $V_{h,r}^{\ell}\subset L^{2}(\mathcal{B}_{h})$. Recalling the definition of \overline{\mathcal{P}}_{3}(\text{\large\otimes}^{2}T\mathcal{T}) in (3.9), we define

[TABLE]

Note that \big{(}\boldsymbol{r}^{\mathcal{T},\mathcal{T}}_{I,J}\boldsymbol{\mathsf{t}}\big{)}\big{|}_{\mathcal{F}}=\boldsymbol{0} for $I,J=1,2,3$, and for every $\mathcal{F}$ in the mesh and $\overline{V}_{h,3}^{\mathbf{c}}\subset H^{\mathbf{c}}(\mathcal{B}_{h})$.

3.3 The Compatible-Strain Mixed Finite Element Methods

We write the following mixed finite element method for the boundary-value problem of incompressible nonlinear elastostatics (2.4) based on the reference elements (3.9) and the corresponding approximation spaces ($V_{h,2}^{1}$, $\overline{V}_{h,3}^{\mathbf{c}}$, $V_{h,1}^{\mathbf{d}-}$, $V_{h,0}^{\ell}$) defined in the previous section:

Given a body force $\boldsymbol{B}$ of $L^{2}$-class, a boundary displacement $\overline{\boldsymbol{U}}$ on $\Gamma_{d}$ of $H^{1/2}$-class, a boundary traction $\overline{\boldsymbol{T}}$ on $\Gamma_{t}$ of $L^{2}$-class, and a stability constant $\alpha\geq 0$, find $(\boldsymbol{U}_{h},\boldsymbol{K}_{h},\boldsymbol{P}_{h},p_{h})\in V_{h,2}^{1}(\Gamma_{d},\overline{\boldsymbol{U}})\times\overline{V}_{h,3}^{\mathbf{c}}\times V_{h,1}^{\mathbf{d}-}\times V_{h,0}^{\ell}$ such that

[TABLE]

where

[TABLE]

and

[TABLE]

Similarly, one can define the following mixed finite element method for the boundary-value problem of compressible nonlinear elastostatics (2.8):

Given ($\boldsymbol{B},\overline{\boldsymbol{U}},{\boldsymbol{T}}$) and $\alpha\geq 0$, find $(\boldsymbol{U}_{h},\boldsymbol{K}_{h},\boldsymbol{P}_{h})\in V_{h,2}^{1}(\Gamma_{d},\overline{\boldsymbol{U}})\times\overline{V}_{h,3}^{\mathbf{c}}\times V_{h,1}^{\mathbf{d}-}$ such that

[TABLE]

The above mixed finite element methods are extensions of the compatible-strain mixed finite element methods (CSFEMs) introduced in [2] and [1] to three dimensional problems.

Remark 1** (Compatibility of Strain and Continuity of Traction).**

1. $(i)$

Recalling (3.12) and considering a displacement gradient field $\boldsymbol{K}_{h}$ on a 3D mesh $\mathcal{B}_{h}$, the Hadamard jump condition is defined as the zero jump $\llbracket\mathsf{t}\boldsymbol{K}_{h}\rrbracket_{\mathcal{F}}=\boldsymbol{0}$ for any tangent vector on ${\mathcal{F}}$ and the three edges enclosing ${\mathcal{F}}$. A necessary condition for the existence of $\boldsymbol{U}_{h}\in H^{1}(T\mathcal{B}_{h})$ such that $\boldsymbol{K}_{h}=\mathbf{grad}\,\boldsymbol{U}_{h}$ is that the Hadamard jump condition holds $\forall\mathcal{F}\in\mathcal{F}_{h}^{i}$ [30]. By construction, the mixed finite element methods (3.14) and (3.15) satisfy this necessary condition as $\boldsymbol{K}_{h}\in\overline{V}_{h,3}^{\mathbf{c}}\subset H^{\mathbf{c}}(\mathcal{B}_{h})$. 2. $(ii)$

Let $\boldsymbol{P}_{h}$ be a stress field on a 3D mesh $\mathcal{B}_{h}$. The localization of the balance of linear momentum requires that $\llbracket\mathsf{n}\boldsymbol{P}_{h}\rrbracket_{\mathcal{F}}=\boldsymbol{0}$, $\forall\mathcal{F}\in\mathcal{F}_{h}^{i}$, that is, the traction vector associated with $\boldsymbol{P}_{h}$ is continuous across all the internal faces of $\mathcal{B}_{h}$. By construction, (3.14) and (3.15) satisfy this requirement as $\boldsymbol{P}_{h}\in V_{h,1}^{\mathbf{d}-}\subset H^{\mathbf{d}}(\mathcal{B}_{h})$.

3.4 Matrix Formulation of CSFEMs

The procedure of obtaining the matrix formulation of (3.14) or (3.15) is similar to that of 2D CSFEMs, which we discussed in detail in [2, §3.4]. In this section, we only write the final formulations needed for the implementation and studying the stability of the 3D CSFEMs. One can write (3.14) in the following matrix form

[TABLE]

where

[TABLE]

The column vectors $\boldsymbol{\mathsf{q}}^{1}_{h}$, $\boldsymbol{\mathsf{q}}^{\mathbf{c}}_{h}$, $\boldsymbol{\mathsf{q}}^{\mathbf{d}}_{h}$, $\boldsymbol{\mathsf{q}}^{\ell}_{h}$ contain all the unknown global degrees of freedom for $\boldsymbol{U}$, $\boldsymbol{K}$, $\boldsymbol{P}$, and $p$, respectively. Let $n$ be the total number of nodes in $\mathcal{B}_{h}$ except those lying on the displacement boundary $\Gamma_{d}$, and let $n_{\mathcal{E}}$, $n_{\mathcal{F}}$, and $n_{\mathcal{T}}$ be the total numbers of edges, faces, and elements in $\mathcal{B}_{h}$, respectively. The number of degrees of freedom in $\boldsymbol{\mathsf{q}}^{1}_{h}$, $\boldsymbol{\mathsf{q}}^{\mathbf{c}}_{h}$, $\boldsymbol{\mathsf{q}}^{\mathbf{d}}_{h}$, and $\boldsymbol{\mathsf{q}}^{\ell}_{h}$ is $3n$, $6n_{\mathcal{E}}+9n_{\mathcal{T}}$, $3n_{\mathcal{F}}$, and $n_{\mathcal{T}}$, respectively, see Figure 2. The total number of degrees of freedom is $N=3n+6n_{\mathcal{E}}+3n_{\mathcal{F}}+10n_{\mathcal{T}}$. The size of $\mathbb{K}_{h}$ is $N\times N$, and the size of $\mathbb{Q}_{h}$, $\mathbb{N}_{h}$, and $\mathbb{F}_{h}$ is $N\times 1$. Let us define $\boldsymbol{V}_{\mathcal{T}}:=\boldsymbol{V}_{h}|_{\mathcal{T}}$ and $\boldsymbol{T}_{\mathcal{T}}:=\boldsymbol{T}_{h}|_{\mathcal{T}}$ for any discrete vector field $\boldsymbol{V}_{h}$ and any discrete tensor field $\boldsymbol{T}_{h}$. The global sparse matrices $\boldsymbol{\mathsf{M}}^{11}_{h}$, $\boldsymbol{\mathsf{M}}^{1\mathbf{c}}_{h}$, $\boldsymbol{\mathsf{M}}^{\mathbf{c}\mathbf{c}}_{h}$, $\boldsymbol{\mathsf{K}}^{1\mathbf{d}}_{h}$, and $\boldsymbol{\mathsf{K}}^{\mathbf{c}\mathbf{d}}_{h}$ in $\mathbb{K}_{h}$ are the result of assembling a set of $n_{\mathcal{T}}$ local matrices that are obtained from calculating respectively $\llangle\mathbf{grad}\,\boldsymbol{U}_{\mathcal{T}},\mathbf{grad}\,\boldsymbol{\Upsilon}_{\mathcal{T}}\rrangle$, $-\llangle\boldsymbol{K}_{\mathcal{T}},\mathbf{grad}\,\boldsymbol{\Upsilon}_{\mathcal{T}}\rrangle$, $\llangle\boldsymbol{K}_{\mathcal{T}},\boldsymbol{\kappa}_{\mathcal{T}}\rrangle$, $\llangle\boldsymbol{P}_{\mathcal{T}},\mathbf{grad}\,\boldsymbol{\Upsilon}_{\mathcal{T}}\rrangle$, and $-\llangle\boldsymbol{P}_{\mathcal{T}},\boldsymbol{\kappa}_{\mathcal{T}}\rrangle$, $\forall\mathcal{T}\in\mathcal{B}_{h}$. Moreover, $\mathbb{K}_{h}$ is a symmetric matrix and $\boldsymbol{\mathsf{M}}^{\mathbf{c}1}_{h}=(\boldsymbol{\mathsf{M}}^{1\mathbf{c}}_{h})^{\mathsf{T}}$, $\boldsymbol{\mathsf{K}}^{\mathbf{d}1}_{h}=(\boldsymbol{\mathsf{K}}^{1\mathbf{d}}_{h})^{\mathsf{T}}$, and $\boldsymbol{\mathsf{K}}^{\mathbf{d}\mathbf{c}}_{h}=(\boldsymbol{\mathsf{K}}^{\mathbf{c}\mathbf{d}}_{h})^{\mathsf{T}}$. For given $\boldsymbol{\mathsf{q}}^{\mathbf{c}}_{h}$ and $\boldsymbol{\mathsf{q}}^{\ell}_{h}$, one obtains the global vectors $\boldsymbol{\mathsf{N}}^{\mathbf{c}}_{h}(\boldsymbol{\mathsf{q}}^{\mathbf{c}}_{h},\boldsymbol{\mathsf{q}}^{\ell}_{h})$ and $\boldsymbol{\mathsf{N}}^{\ell}_{h}(\boldsymbol{\mathsf{q}}^{\mathbf{c}}_{h})$ in $\mathbb{N}_{h}$ by assembling a set of $n_{\mathcal{T}}$ local vectors that are obtained form calculating the nonlinear terms $\llangle\widetilde{\boldsymbol{P}}(\boldsymbol{K}_{\mathcal{T}})+p_{\mathcal{T}}\boldsymbol{Q}(\boldsymbol{K}_{\mathcal{T}}),\boldsymbol{\kappa}_{\mathcal{T}}\rrangle$ and $\llangle C(J_{\mathcal{T}}),q_{\mathcal{T}}\rrangle$, $\forall\mathcal{T}\in\mathcal{B}_{h}$, respectively. Similarly, for a given body force $\boldsymbol{B}$, one obtains $\boldsymbol{\mathsf{F}}^{1}_{h}$ in $\mathbb{F}_{h}$ by calculating $\llangle\rho_{0}\boldsymbol{B},\boldsymbol{\Upsilon}_{\mathcal{T}}\rrangle$, $\forall\mathcal{T}\in\mathcal{B}_{h}$. Finally, for a given traction $\overline{\boldsymbol{T}}$ on $\Gamma_{t}$, one obtains $\boldsymbol{\mathsf{F}}^{1}_{\Gamma_{t}}$ in $\mathbb{F}_{h}$ through assembling all the local vectors obtained from $\int_{\mathcal{F}^{\mathcal{T}}_{t}}\langle\overline{\boldsymbol{T}},\boldsymbol{\Upsilon}|_{\mathcal{F}^{\mathcal{T}}_{t}}\rangle\,dA$ for every face $\mathcal{F}^{\mathcal{T}}_{t}$ lying on $\Gamma_{t}$. See [2, (3.21)-(3.33)] for details of calculating the local matrices and vectors in each element. To obtain the matrix form of (3.15) for compressible solids, we modify (3.16) by setting $p_{h}=0$ ($\boldsymbol{\mathsf{q}}^{\ell}_{h}=\boldsymbol{0}$) and removing the fourth row and the fourth column of $\mathbb{K}_{h}$ and the forth entries of $\mathbb{Q}_{h}$, $\mathbb{N}_{h}$, and $\mathbb{F}_{h}$. We also use $\widehat{\boldsymbol{P}}(\boldsymbol{K})$ instead of $\widetilde{\boldsymbol{P}}(\boldsymbol{K})$ in our calculations.

Using Newton’s method, one can approximate the solution of the nonlinear equation (3.16) iteratively using $\mathbb{Q}^{(i+1)}_{h}=\mathbb{Q}^{(i)}_{h}-\mathbb{K}^{-1}_{t_{h}}\left(\mathbb{Q}^{(i)}_{h}\right)\mathbb{R}_{h}\left(\mathbb{Q}^{(i)}_{h}\right)$, where $\mathbb{R}_{h}(\mathbb{Q}_{h})=\mathbb{K}_{h}\mathbb{Q}_{h}+\mathbb{N}_{h}(\mathbb{Q}_{h})-\mathbb{F}_{h}$ is the residual vector and $\mathbb{K}_{{t}_{h}}$ is the tangent stiffness matrix (Jacobian matrix) given by

[TABLE]

The matrix $\boldsymbol{\mathsf{H}}^{\mathbf{c}\mathbf{c}}_{h}(\boldsymbol{\mathsf{q}}^{\mathbf{c}}_{h},\boldsymbol{\mathsf{q}}^{\ell}_{h})$ \big{(}\boldsymbol{\mathsf{H}}^{\mathbf{c}\ell}_{h}(\boldsymbol{\mathsf{q}}^{\mathbf{c}}_{h})\big{)} contains the derivative of components of $\boldsymbol{\mathsf{N}}^{\mathbf{c}}_{h}(\boldsymbol{\mathsf{q}}^{\mathbf{c}}_{h},\boldsymbol{\mathsf{q}}^{\ell}_{h})$ in $\mathbb{N}_{h}$ with respect to components of $\boldsymbol{\mathsf{q}}^{\mathbf{c}}_{h}$ ($\boldsymbol{\mathsf{q}}^{\ell}_{h}$). Also, $\boldsymbol{\mathsf{H}}^{\ell\mathbf{c}}_{h}(\boldsymbol{\mathsf{q}}^{\mathbf{c}}_{h})$ contains the derivative of components of $\boldsymbol{\mathsf{N}}^{\ell}_{h}(\boldsymbol{\mathsf{q}}^{\mathbf{c}}_{h})$ in $\mathbb{N}_{h}$ with respect to components of $\boldsymbol{\mathsf{q}}^{\mathbf{c}}_{h}$. Linearizing $\llangle\widetilde{\boldsymbol{P}}(\boldsymbol{K})+p\boldsymbol{Q}(\boldsymbol{K}),\boldsymbol{\kappa}\rrangle$ in (2.4)2 at a given displacement gradient $\boldsymbol{K}^{0}\in H^{\mathbf{c}}(\mathcal{B})$ and a given pressure $p^{0}\in L^{2}(\mathcal{B})$ gives $\llangle\boldsymbol{\widetilde{A}}(\boldsymbol{K}^{0},p^{0})\!:\!\boldsymbol{K},\boldsymbol{\kappa}\rrangle+\llangle p\boldsymbol{Q}(\boldsymbol{K}^{0}),\boldsymbol{\kappa}\rrangle$, where $\boldsymbol{\widetilde{A}}$ is the elasticity tensor and $(\boldsymbol{\widetilde{A}}\!:\!\boldsymbol{K})^{IJ}=\widetilde{A}^{IJMN}K^{MN}$. Also, linearizing $\llangle C(J(\boldsymbol{K})),q\rrangle$ in (2.4)4 at $\boldsymbol{K}^{0}$ results in $\llangle\boldsymbol{Q}(\boldsymbol{K}^{0})\!:\!\boldsymbol{K},q\rrangle=\llangle\boldsymbol{K},q\boldsymbol{Q}(\boldsymbol{K}^{0})\rrangle$, where $\boldsymbol{Q}\!:\!\boldsymbol{K}={Q}^{IJ}{K}^{IJ}$. After discretization, for given $\boldsymbol{K}^{0}_{h}$ and $p^{0}_{h}$ (or $\boldsymbol{\mathsf{q}}^{\mathbf{c}}_{h}$ and $\boldsymbol{\mathsf{q}}^{\ell}_{h}$), one can calculate the local matrices for $\llangle\boldsymbol{\widetilde{A}}_{\mathcal{T}}(\boldsymbol{K}^{0}_{\mathcal{T}},p^{0}_{\mathcal{T}})\!:\!\boldsymbol{K}_{\mathcal{T}},\boldsymbol{\kappa}_{\mathcal{T}}\rrangle$, $\llangle p_{\mathcal{T}}\boldsymbol{Q}(\boldsymbol{K}^{0}_{\mathcal{T}}),\boldsymbol{\kappa}_{\mathcal{T}}\rrangle$, and $\llangle\boldsymbol{K}_{\mathcal{T}},q_{\mathcal{T}}\boldsymbol{Q}(\boldsymbol{K}^{0}_{\mathcal{T}})\rrangle$ $\forall\mathcal{T}\in\mathcal{B}_{h}$ to assemble $\boldsymbol{\mathsf{H}}^{\mathbf{c}\mathbf{c}}_{h}$, $\boldsymbol{\mathsf{H}}^{\mathbf{c}\ell}_{h}$, and $\boldsymbol{\mathsf{H}}^{\ell\mathbf{c}}_{h}$, respectively. For more details, see [2, (3.36)], and note that $\boldsymbol{\mathsf{H}}^{\mathbf{c}\ell}_{h}=(\boldsymbol{\mathsf{H}}^{\ell\mathbf{c}}_{h})^{\mathsf{T}}$. For compressible solids, (3.17) simplifies to

[TABLE]

where $\widehat{\boldsymbol{\mathsf{H}}}_{h}^{\mathbf{c}\mathbf{c}}(\boldsymbol{\mathsf{q}}^{\mathbf{c}}_{h})$ is obtained by linearizing $\llangle\widehat{\boldsymbol{P}}(\boldsymbol{K}_{h}),\boldsymbol{\kappa}_{h}\rrangle$ in (3.15)2.

We use neo-Hookean materials for our numerical examples. However, note that our formulation can use any elastic constitutive equation. For compressible solids, we use the energy function

[TABLE]

where $\mu$ and $\kappa$ are the shear and bulk moduli at the ground state, respectively. Recalling that $\boldsymbol{F}=\boldsymbol{I}+\boldsymbol{K}$, the constitutive relation reads

[TABLE]

To calculate $\mathbb{K}_{t_{h}}$ defined in (3.18), one needs to obtain the elasticity tensor $\boldsymbol{\widehat{A}}(\boldsymbol{K})$ by taking the derivative of components of $\widehat{\boldsymbol{P}}(\boldsymbol{K})$ with respect to components of $\boldsymbol{K}$. In the implementation (see [2, (3.36)]), it is more convenient to represent the elasticity tensor as a matrix $\widehat{\boldsymbol{\mathsf{A}}}$, whose size is $9\times 9$ in 3D. Let $\lceil\boldsymbol{T}\rceil$ be a vector representation of a tensor $\boldsymbol{T}$, and let $\left[\boldsymbol{V}\right]_{\!\times}$ be a skew-symmetric matrix representing a vector $\boldsymbol{V}$ that are defined as

[TABLE]

Considering (3.19), one obtains

[TABLE]

where $\boldsymbol{\mathsf{I}}$ is the $9\times 9$ identity matrix and

[TABLE]

For incompressible solids with $I_{3}=1$, the following neo-Hookean energy function is used.

[TABLE]

The constitutive part of stress is $\widetilde{P}(\boldsymbol{K})=\mu(\boldsymbol{I}+\boldsymbol{K})$. To impose the incompressibility constraint $J=1$, we use the constraint functions $C_{1}(J)=J-1$ or $C_{2}(J)=\ln J$; we choose the function that results in a better numerical performance of the method in a given example. To obtain $\mathbb{N}_{h}$ in (3.16) and $\mathbb{K}_{t_{h}}$ in (3.17), one needs the following matrices:

[TABLE]

3.5 Solvability and Stability

Theorem 2**.**

Let $N^{1}$, $N^{c}$, $N^{d}$, and $N^{\ell}$ be the numbers of degrees of freedom in $\boldsymbol{\mathsf{q}}^{1}_{h}$, $\boldsymbol{\mathsf{q}}^{\mathbf{c}}_{h}$, $\boldsymbol{\mathsf{q}}^{\mathbf{d}}_{h}$, and $\boldsymbol{\mathsf{q}}^{\ell}_{h}$, respectively. For $\alpha>0$, the tangent stiffness matrix $\mathbb{K}_{t_{h}}$ of incompressible solids (3.17) is non-singular if and only if the following conditions hold.

1. $(i)$

$\operatorname{ker}(\boldsymbol{\mathsf{H}}^{\mathbf{c}\ell}_{h})=\{\boldsymbol{0}_{N^{\ell}\times 1}\}$, 2. $(ii)$

$\operatorname{ker}(\boldsymbol{\mathsf{K}}^{1\mathbf{d}}_{h})\cap\operatorname{ker}(\boldsymbol{\mathsf{B}}^{\ell\mathbf{c}}_{0}\boldsymbol{\mathsf{K}}^{\mathbf{c}\mathbf{d}}_{h})=\{\boldsymbol{0}_{N^{d}\times 1}\}$,  $\left(\operatorname{ker}(\boldsymbol{\mathsf{K}}^{\mathbf{c}\mathbf{d}}_{h})\subseteq\operatorname{ker}(\boldsymbol{\mathsf{B}}^{\ell\mathbf{c}}_{0}\boldsymbol{\mathsf{K}}^{\mathbf{c}\mathbf{d}}_{h})\right)$, 3. $(iii)$

$\operatorname{ker}\left(\boldsymbol{\mathsf{H}}^{\mathbf{c}\mathbf{c}}_{h}+\alpha\boldsymbol{\mathsf{M}}^{\mathbf{c}\mathbf{c}}_{h}-\alpha\boldsymbol{\mathsf{M}}^{\mathbf{c}1}_{h}\left(\boldsymbol{\mathsf{M}}^{11}_{h}\right)^{-1}\boldsymbol{\mathsf{M}}^{1\mathbf{c}}_{h}\right)\cap\operatorname{ker}\left(\boldsymbol{\mathsf{K}}^{\mathbf{d}\mathbf{c}}_{h}-\boldsymbol{\mathsf{K}}^{\mathbf{d}1}_{h}\left(\boldsymbol{\mathsf{M}}^{11}_{h}\right)^{-1}\boldsymbol{\mathsf{M}}^{1\mathbf{c}}_{h}\right)\cap\operatorname{ker}\left(\boldsymbol{\mathsf{H}}^{\ell\mathbf{c}}_{h}\right)=\{\boldsymbol{0}_{N^{c}\times 1}\}$,

where $\boldsymbol{\mathsf{B}}^{\ell\mathbf{c}}_{0}$ is a matrix whose rows form a basis for $\operatorname{ker}(\boldsymbol{\mathsf{H}}^{\ell\mathbf{c}}_{h})$. For $\alpha=0$, $\mathbb{K}_{t_{h}}$ is non-singular if and only if $(i)$ and $(ii)$ and the following conditions hold.

1. $(iii)^{\prime}$

$\operatorname{ker}(\boldsymbol{\mathsf{K}}^{\mathbf{d}1}_{h})=\{\boldsymbol{0}_{N^{1}\times 1}\}$, 2. $(iv)^{\prime}$

$\operatorname{ker}(\boldsymbol{\mathsf{H}}^{\mathbf{c}\mathbf{c}}_{h})\cap\operatorname{ker}\left(\boldsymbol{\mathsf{B}}^{1\mathbf{d}}_{0}\boldsymbol{\mathsf{K}}^{\mathbf{d}\mathbf{c}}_{h}\right)\cap\operatorname{ker}\left(\boldsymbol{\mathsf{H}}^{\ell\mathbf{c}}_{h}\right)=\{\boldsymbol{0}_{N^{c}\times 1}\}$,

where $\boldsymbol{\mathsf{B}}^{1\mathbf{d}}_{0}$ is a matrix whose rows form a basis for $\operatorname{ker}(\boldsymbol{\mathsf{K}}^{1\mathbf{d}}_{h})$.

Proof.

Rearrange the rows and the columns of $\mathbb{K}_{t_{h}}$ to obtain

[TABLE]

Then, according to [31, Theorem 3.2.1], the matrix $\mathbb{K}_{t_{h}}$ is non-singular if and only if the following holds:

1. $(1)$

The restriction of $\mathbb{A}_{h}$ to $\operatorname{ker}(\mathbb{B}_{h})$ is surjective (or equivalently injective), 2. $(2)$

$\mathbb{B}_{h}$ is surjective (or equivalently $\mathbb{B}^{\mathsf{T}}_{h}$ is injective or $\operatorname{ker}(\mathbb{B}^{\mathsf{T}}_{h})=\{\boldsymbol{0}\}$).

Consider the following sets:

[TABLE]

One can show that

[TABLE]

Therefore, the requirement $(2)$ ($\operatorname{ker}(\mathbb{B}^{\mathsf{T}}_{h})=\{\boldsymbol{0}\}$) is equivalent to $S^{\prime}_{1}=S^{\prime}_{2}=\emptyset$. It is straightforward to show that $S^{\prime}_{1}=\emptyset$ is equivalent to $(i)$. We write $S^{\prime}_{2}=\emptyset$ as $\operatorname{ker}(\boldsymbol{\mathsf{K}}^{1\mathbf{d}}_{h})\cap s^{\prime}_{2}=\{\boldsymbol{0}\}$, where $s^{\prime}_{2}=\left\{\boldsymbol{\mathsf{X}}:\boldsymbol{\mathsf{K}}^{\mathbf{c}\mathbf{d}}_{h}\boldsymbol{\mathsf{X}}=-\boldsymbol{\mathsf{H}}^{\mathbf{c}\ell}_{h}\boldsymbol{\mathsf{Y}}\text{ for some }\boldsymbol{\mathsf{Y}}\in\mathbb{R}^{N^{\ell}}\right\}$. Using $\operatorname{Im}(\boldsymbol{\mathsf{H}}^{\mathbf{c}\ell}_{h})=\operatorname{ker}(\boldsymbol{\mathsf{H}}^{\ell\mathbf{c}}_{h})^{\perp}$, where the superscript $\perp$ indicates the orthogonal complement, one can readily show that $s^{\prime}_{2}=\left\{\boldsymbol{\mathsf{X}}:\boldsymbol{\mathsf{b}}_{0}^{\mathsf{T}}\boldsymbol{\mathsf{K}}^{\mathbf{c}\mathbf{d}}_{h}\boldsymbol{\mathsf{X}}=0,\forall\boldsymbol{\mathsf{b}}_{0}\in\operatorname{ker}(\boldsymbol{\mathsf{H}}^{\ell\mathbf{c}}_{h})\right\}$. Let $\boldsymbol{\mathsf{B}}^{\ell\mathbf{c}}_{0}$ be a matrix whose rows form a basis for $\operatorname{ker}(\boldsymbol{\mathsf{H}}^{\ell\mathbf{c}}_{h})$, then $s^{\prime}_{2}=\operatorname{ker}(\boldsymbol{\mathsf{B}}^{\ell\mathbf{c}}_{0}\boldsymbol{\mathsf{K}}^{\mathbf{c}\mathbf{d}}_{h})$. Therefore, $S^{\prime}_{2}=\emptyset$ is equivalent to $(ii)$. Next, we assume that $\alpha>0$ and show that $(1)$ is equivalent to $(iii)$. For $\alpha>0$, one can write

[TABLE]

where use was made of the fact that $\boldsymbol{\mathsf{M}}^{11}_{h}$ is a Gram matrix and positive-definite by construction (and hence injective). Note that $(1)$ is equivalent to $\operatorname{ker}(\mathbb{A}_{h})\cap\operatorname{ker}(\mathbb{B}_{h})=(\operatorname{ker}(\mathbb{A}_{h})\cap\{\boldsymbol{0}\})\sqcup(\operatorname{ker}(\mathbb{A}_{h})\cap S_{1})\sqcup(\operatorname{ker}(\mathbb{A}_{h})\cap S_{2})=\{\boldsymbol{0}\}$, which is equivalent to $\operatorname{ker}(\mathbb{A}_{h})\cap S_{1}=\operatorname{ker}(\mathbb{A}_{h})\cap S_{2}=\emptyset$. We know that $\mathbb{A}_{h}\boldsymbol{\mathsf{Q}}\neq\boldsymbol{0}$, $\forall\boldsymbol{\mathsf{Q}}\in S_{1}$, due to injectivity of $\boldsymbol{\mathsf{M}}^{11}_{h}$, so $\operatorname{ker}(\mathbb{A}_{h})\cap S_{1}=\emptyset$ is trivial. The remaining condition $\operatorname{ker}(\mathbb{A}_{h})\cap S_{2}=\emptyset$ simplifies to $(iii)$. For $\alpha=0$, one can write

[TABLE]

Now, $\operatorname{ker}(\mathbb{A}_{h})\cap S_{1}=\emptyset$ simplifies to $(iii)^{\prime}$, and $\operatorname{ker}(\mathbb{A}_{h})\cap S_{2}=\emptyset$ simplifies to $(iv)^{\prime}$. ∎

Corollary 3**.**

For $\alpha>0$, the tangent stiffness matrix $\mathbb{K}_{t_{h}}$ of compressible solids (3.18) is non-singular if and only if the following conditions hold:

1. $(i)$

$\operatorname{ker}(\boldsymbol{\mathsf{K}}^{1\mathbf{d}}_{h})\cap\operatorname{ker}(\boldsymbol{\mathsf{K}}^{\mathbf{c}\mathbf{d}}_{h})=\{\boldsymbol{0}_{N^{d}\times 1}\}$, 2. $(ii)$

$\operatorname{ker}\left(\boldsymbol{\mathsf{H}}^{\mathbf{c}\mathbf{c}}_{h}+\alpha\boldsymbol{\mathsf{M}}^{\mathbf{c}\mathbf{c}}_{h}-\alpha\boldsymbol{\mathsf{M}}^{\mathbf{c}1}_{h}\left(\boldsymbol{\mathsf{M}}^{11}_{h}\right)^{-1}\boldsymbol{\mathsf{M}}^{1\mathbf{c}}_{h}\right)\cap\operatorname{ker}\left(\boldsymbol{\mathsf{K}}^{\mathbf{d}\mathbf{c}}_{h}-\boldsymbol{\mathsf{K}}^{\mathbf{d}1}_{h}\left(\boldsymbol{\mathsf{M}}^{11}_{h}\right)^{-1}\boldsymbol{\mathsf{M}}^{1\mathbf{c}}_{h}\right)=\{\boldsymbol{0}_{N^{c}\times 1}\}$.

For $\alpha=0$, $\mathbb{K}_{t_{h}}$ is non-singular if and only if $(i)$ and the following conditions hold:

1. $(ii)^{\prime}$

$\operatorname{ker}(\boldsymbol{\mathsf{K}}^{\mathbf{d}1}_{h})=\{\boldsymbol{0}_{N^{1}\times 1}\}$, 2. $(iii)^{\prime}$

$\operatorname{ker}(\boldsymbol{\mathsf{H}}^{\mathbf{c}\mathbf{c}}_{h})\cap\operatorname{ker}\left(\boldsymbol{\mathsf{B}}^{1\mathbf{d}}_{0}\boldsymbol{\mathsf{K}}^{\mathbf{d}\mathbf{c}}_{h}\right)=\{\boldsymbol{0}_{N^{c}\times 1}\}$.

Corollary 4**.**

If the tangent stiffness matrix $\mathbb{K}_{t_{h}}$ is non-singular, then

1. $(1)$

$N^{d}\leq N^{c}+N^{1}$* for $\alpha\geq 0$ and for both compressible and incompressible solids,* 2. $(2)$

$N^{\ell}\leq N^{c}$* only for incompressible solids,* 3. $(3)$

$N^{1}\leq N^{d}$* only for $\alpha=0$.*

Proof.

Noting that $\operatorname{ker}(\boldsymbol{\mathsf{K}}^{\mathbf{c}\mathbf{d}}_{h})\subseteq\operatorname{ker}(\boldsymbol{\mathsf{B}}^{\ell\mathbf{c}}_{0}\boldsymbol{\mathsf{K}}^{\mathbf{c}\mathbf{d}}_{h})$, both Theorem 2 $(ii)$ and Corollary 3 $(i)$, imply $(1)$. Theorem 2 (i) implies $(2)$. Both Theorem 2 $(iii)^{\prime}$ and Corollary 3 $(ii)^{\prime}$, imply $(3)$. ∎

In view of Theorem 2, one can see how adding (2.6) to the weak formulation (2.4) may improve the stability and the performance of the resulting finite element methods. Without the stabilization terms ($\alpha=0$), the violation of $\operatorname{ker}(\boldsymbol{\mathsf{K}}^{\mathbf{d}1}_{h})=\{\boldsymbol{0}_{N^{1}\times 1}\}$ or more strongly having $N^{1}>N^{d}$ leads to a singular tangent stiffness matrix $\mathbb{K}_{t_{h}}$. This restricts the choices of finite elements for the displacement and stress in both 2D and 3D. In particular, considering $(\boldsymbol{U}_{h},\boldsymbol{P}_{h})$ in $V_{h,m}^{1}\times V_{h,n}^{\mathbf{d}-}$ or $V_{h,m}^{1}\times V_{h,n}^{\mathbf{d}}$ such that $m>n$ results in a singular $\mathbb{K}_{t_{h}}$ in both 2D and 3D independent of the size of the mesh. Adding the stabilization terms (2.6) ($\alpha>0$) overcomes this limitation and enables one to improve the convergence of the displacement field by discretizing it using second-order shape functions without the need for modifying the finite elements of other fields. For instance, for $\alpha=0$, the finite elements (3.9) result in a singular system, but they converge to correct solutions for large values of $\alpha$. To avoid the singularity of (3.9) for $\alpha=0$, we have no choice but to approximate the displacement field using first-order polynomials and to compromise the rate of convergence of the method. Also, note that approximating the displacement $\boldsymbol{U}$ in a second-order polynomial space leads to a more accurate discretization of $\boldsymbol{K}=\mathbf{grad}\,\boldsymbol{U}$ as the intersection of image of $\mathbf{grad}$ and the approximation space of $\boldsymbol{K}$ becomes larger at the discrete level.

We next discuss how modifying the finite element of the displacement gradient $\boldsymbol{K}$ in (3.9) and its resulting finite element space (3.13) lead to solvability of the mixed finite element methods (3.14) and (3.15). Let $(\boldsymbol{U}_{h},\boldsymbol{K}_{h},\boldsymbol{P}_{h})\in V_{h,m}^{1}\times V_{h,n}^{\mathbf{c}-}(V_{h,n}^{\mathbf{c}})\times V_{h,k}^{\mathbf{d}-}(V_{h,k}^{\mathbf{d}})$ for $m,n,k=1,2$, which results in $32$ different combinations (note that pressure is not relevant here). In 3D, all the $32$ combinations except $V_{h,m}^{1}\times V_{h,2}^{\mathbf{c}}\times V_{h,1}^{\mathbf{d}-},m=1,2$ result in a singular $\mathbb{K}_{t_{h}}$. These combinations either give $N^{d}>N^{c}+N^{1}$ for any mesh or their smallest singular value of $\begin{bmatrix}\boldsymbol{\mathsf{K}}^{\mathbf{d}\mathbf{c}}_{h}&\boldsymbol{\mathsf{K}}^{\mathbf{d}1}_{h}\end{bmatrix}^{\mathsf{T}}$ goes to zero as one refines the mesh (see [2, Remark 16]). Any of these two cases is a violation of Theorem 2 $(ii)$ (or Corollary 3 $(i)$). The two remaining choices $V_{h,m}^{1}\times V_{h,2}^{\mathbf{c}}\times V_{h,1}^{\mathbf{d}-},m=1,2$ are not practical as they have poor performances considering their expensive computational cost. $V_{h,2}^{\mathbf{c}}$ of displacement gradient has $90$ degrees of freedom per element, which significantly increases the computational cost, but paired with the lowest-order space of stress $V_{h,1}^{\mathbf{d}-}$, it cannot improve the overall convergence of the method. To resolve this issue, we proposed $\overline{V}_{h,3}^{\mathbf{c}}$ in (3.13) and considered $(\boldsymbol{U}_{h},\boldsymbol{K}_{h},\boldsymbol{P}_{h})\in V_{h,2}^{1}\times\overline{V}_{h,3}^{\mathbf{c}}\times V_{h,1}^{\mathbf{d}-}$. Note that, in each element, $\overline{V}_{h,3}^{\mathbf{c}}$ has only $9$ degrees of freedom more than the first-order space $V_{h,1}^{\mathbf{c}}$ with $36$ degrees of freedom (see Figure 2 ). Hence, it does not increase the computational cost of the method significantly. Moreover, we observe that the smallest singular value of $[\boldsymbol{\mathsf{K}}^{\mathbf{d}\mathbf{c}}_{h}~{}\boldsymbol{\mathsf{K}}^{\mathbf{d}1}_{h}]^{\mathsf{T}}$ for $V_{h,2}^{1}\times\overline{V}_{h,3}^{\mathbf{c}}\times V_{h,1}^{\mathbf{d}-}$ remains positive as we refine different arbitrary meshes.

So far, we have discussed that, for $\alpha>0$, (3.14) and (3.15) do not result in a singular $\mathbb{K}_{t_{h}}$ even if $\operatorname{ker}(\boldsymbol{\mathsf{K}}^{\mathbf{d}1}_{h})\neq\{\boldsymbol{0}_{N^{1}\times 1}\}$, and they result in $N^{d}\leq N^{c}+N^{1}$ and $\operatorname{ker}(\boldsymbol{\mathsf{K}}^{1\mathbf{d}}_{h})\cap\operatorname{ker}(\boldsymbol{\mathsf{K}}^{\mathbf{c}\mathbf{d}}_{h})=\{\boldsymbol{0}_{N^{d}\times 1}\}$, which are required for satisfying Theorem 2 $(ii)$ or Corollary 3 $(i)$. These have been made possible through studying the linear operators $\boldsymbol{\mathsf{K}}^{\mathbf{c}\mathbf{d}}_{h}$ and $\boldsymbol{\mathsf{K}}^{1\mathbf{d}}_{h}$ in $\mathbb{K}_{t_{h}}$, which are independent of the physics of the problem. The stability of (3.14) requires that all the conditions of Theorem 2 hold as one refines the mesh. However, given the nonlinear nature of the problems of interest here, this is difficult to check. In particular, the nonlinear operators $\boldsymbol{\mathsf{H}}^{\mathbf{c}\mathbf{c}}_{h}(\boldsymbol{\mathsf{q}}^{\mathbf{c}}_{h},\boldsymbol{\mathsf{q}}^{\ell}_{h})$ and $\boldsymbol{\mathsf{H}}^{\mathbf{c}\ell}_{h}(\boldsymbol{\mathsf{q}}^{\mathbf{c}}_{h})$ in $\mathbb{K}_{t_{h}}$ depend on the material properties of the body and its state of deformation. Therefore, one cannot draw a general conclusion for stability or convergence of the mixed methods only by studying the formulations and without considering the physics of the problem. Based on the various numerical examples presented in the next section, we have concluded that (3.14) and (3.15) have an overall good performance in capturing the large deformations of incompressible and compressible solids in 3D.

4 Numerical Examples

In this section, we consider several examples to assess the performance of the mixed finite elements (3.14) and (3.15) in modeling compressible and incompressible solids in 3D. We use the Frobenius norm $\|\boldsymbol{T}\|:=(\sum\nolimits_{I,J}T^{IJ}T^{IJ})^{\frac{1}{2}}$ for $\boldsymbol{K}_{h}$ and $\boldsymbol{P}_{h}$ in the deformed configurations. We use the $L^{2}$-norm for $\boldsymbol{U}_{h}$, $\boldsymbol{K}_{h}$, $\boldsymbol{P}_{h}$, and $p_{h}$ over the entire mesh in convergence analyses. We use $\alpha=1\times 10^{6}$ in all the examples (the solutions actually converge for smaller values of $\alpha$ in each example; assuming larger values does not change the solutions).

Example 1: Inflation of a Hollow Spherical Ball.

Let us consider an incompressible hollow spherical ball shown in Figure 3. We assume that the inner boundary of the ball is subjected to the displacement boundary condition $\boldsymbol{U}_{\text{in}}=(\lambda-1)\mathbf{X}$, the outer boundary is traction free, and there are no body forces. This is an example of a universal deformation [32] and the exact solution reads

[TABLE]

where $R=\|\mathbf{X}\|$, $r(R)=\left(R^{3}+(\lambda^{3}-1)R_{\text{in}}^{3}\right)^{\frac{1}{3}}$, and $g(R)=R\left(3r^{3}(R)+(\lambda^{3}-1)R^{3}_{\text{in}}\right)/r^{4}(R)$. It follows that $\boldsymbol{K}_{e}=\mathbf{grad}\,\boldsymbol{U}_{e}$, and $\boldsymbol{P}_{e}=\widetilde{\boldsymbol{P}}(\boldsymbol{K}_{e})+p_{e}\boldsymbol{Q}(\boldsymbol{K}_{e})$. Having the exact solution, we assess the accuracy and convergence of CSFEM given in (3.14). For our computations, we consider the neo-Hookean energy function (3.20) with $\mu=1~{}\mathrm{N}/\mathrm{mm}^{2}$, the constraint function $C(J)=J-1$, $R_{\text{in}}=0.5\,\mathrm{mm}$, $R_{\text{out}}=1\,\mathrm{mm}$, and $\lambda=3$. Using symmetry, we model only $1/24$ of a hemisphere as shown in Figure 3. To study the convergence order of (3.14), we plot the relative errors of the field variables versus the maximum diameter $h$ of some unstructured meshes in a log-log graph in Figure 4. The convergence order of the displacement $\boldsymbol{U}_{h}$ is close to $2$, and those of the displacement gradient $\boldsymbol{K}_{h}$, the stress $\boldsymbol{P}_{h}$, and the pressure-like variable $p_{h}$ are almost $1$. Figure 5 shows the reference and the deformed configurations of the four unstructured meshes given in Figure 3 obtained using CSFEM in (3.14) for $\lambda=3$. Colors show the values of $\|\boldsymbol{K}_{h}\|$ in the first row and the values of $p_{h}$ in the second row with lighter colors associated with the larger values.

Example 2: $\boldsymbol{3}$D Cook’s Membrane.

In this example, the 3D Cook’s membrane problem depicted in Figure 6 is analyzed in order to study the performance of CSFEMs in bending analysis. We consider two cases of tractions imposed on the right side of the membrane (on $16\,\mathrm{mm}\times 10\,\mathrm{mm}$ face): $\overline{\boldsymbol{T}}_{1}=\left(0,f,0\right)$ and $\overline{\boldsymbol{T}}_{2}=\left(0,2f,f\right)$. We use the energy function (3.20) with $\mu=1\,\mathrm{N}/\mathrm{mm}^{2}$ and $C(J)=\ln J$ to impose the incompressibility constraint. Figure 7 shows the convergence of the vertical displacement of point $A$ indicated in Figure 6 for different values of traction $\overline{\boldsymbol{T}}_{1}=\left(0,f,0\right)$ using the mixed method (3.14). Since the membrane deforms in two dimensions, the results of the 3D analysis using (3.14) are compared to those obtained by a 2D analysis using $\mathsf{H}2\mathsf{c}2\mathsf{d}\bar{2}\mathsf{L}1$ in [2]. The comparison shows a good agreement between the two analyses. Considering $\overline{\boldsymbol{T}}_{2}=\left(0,2f,f\right)$, the membrane deforms in three dimensions, for which a similar convergence graph for point $A$ is presented in Figure 8. The convergence of the independent field variables $(\boldsymbol{U}_{h},\boldsymbol{K}_{h},\boldsymbol{P}_{h},p_{h})$ obtained using the mixed method (3.14) is illustrated in Figure 9 for different values of $\overline{\boldsymbol{T}}_{2}=\left(0,2f,f\right)$. One observes that $\boldsymbol{U}_{h}$ and $\boldsymbol{K}_{h}$ have a faster convergence in comparison with $\boldsymbol{P}_{h}$ or $p_{h}$. The deformed configurations of the four meshes in Figure 6 using the mixed method (3.14) are given in Figure 10 and Figure 11 for $\overline{\boldsymbol{T}}_{1}=\left(0,0.3,0\right)\,\mathrm{N}/\mathrm{mm}^{2}$ and $\overline{\boldsymbol{T}}_{2}=\left(0,0.2,0.1\right)\,\mathrm{N}/\mathrm{mm}^{2}$, respectively. In both figures, colors indicate the values of $\|\boldsymbol{p}_{h}\|$ in the first row and the values of $p_{h}$ in the second row with lighter colors corresponding to larger values. It is well-known that the standard displacement-pressure mixed methods for incompressible materials approximate displacement accurately but suffer form numerical artifacts in approximating pressure (they are unable to provide an approximation of stress either). By contrast, Figures 10 and 11 clearly show that the mixed method (3.14) does not suffer from any numerical artifacts in approximating the stresses and the pressure in a large deformation of an incompressible solid even for relatively coarse meshes.

Example 3. Compression of a Near-Incompressible Block.

Let us consider a block under compression as shown in Figure 12. The length and the width of the block are $2\,\mathrm{mm}$ and its hight is $1\,\mathrm{mm}$. The loading square surface on the upper face of the block has an edge of $1\,\mathrm{mm}$ and is subjected to a traction $\overline{\boldsymbol{T}}=(0,0,f)$. The vertical (horizontal) displacement at the bottom (top) of the block is zero. As shown in Figure 12, using symmetry the meshes are generated for only a quarter of the block.

In this example we test the performance of the mixed method (3.15) in the near-incompressible regime. Note that many of the existing finite element methods are unable to solve this problem or suffer from numerical artifacts. Reese et al. [13] developed a reduced-integration stabilized brick element and used it to solve this problem. To compare our numerical results to those of [13], we consider the energy function (3.19) with $\lambda=400889.806\,\mathrm{N}/\mathrm{mm}^{2}$ and $\mu=80.194\,\mathrm{N}/\mathrm{mm}^{2}$. Figure 13 illustrates the convergence of the vertical displacement of point $A$ (see Figure 12) for different values of $\overline{\boldsymbol{T}}=(0,0,f)$. The results obtained using (3.15) agree with those reported by Reese et al. [13]. Figure 14 depicts the deformed configuration of the block for $\overline{\boldsymbol{T}}=(0,0,320)~{}\mathrm{N}/\mathrm{mm}^{2}$. Colors show the values of $\|\boldsymbol{K}_{h}\|$, where lighter colors are assigned to larger values.

Example 4. Stretching a Heterogeneous Block.

As was mentioned in Remark 1, CSFEMs (3.14) and (3.15), by construction, satisfy the Hadamard jump condition and the continuity of traction on all the internal faces in a given mesh. This provides an efficient framework to model heterogeneous solids provided that the constituent materials do not slide at their interfaces, i.e., the displacement field is continuous at the material interfaces. One can generate a 3D mesh such that some of the internal faces of the mesh closely approximate the given material interfaces and assign the material model of each inhomogeneity to its corresponding region in the mesh. Using CSFEMs guarantees that the necessary kinematic and kinetic conditions are automatically satisfied at the material interfaces.

We consider an incompressible cubic block of edge $1\,\mathrm{mm}$ with a spherical inhomogeneity of diameter $0.5\,\mathrm{mm}$ at its center as shown in 15. The bottom of the block at $Z=-0.5\,\mathrm{mm}$ and the top of the block at $Z=0.5\,\mathrm{mm}$ are subjected to displacement boundaries $(0,0,-0.5)\,\mathrm{mm}$ and $(0,0,0.5)\,\mathrm{mm}$, respectively (stretch = $2$), and the other four faces are traction free. Using symmetry, we model only $1/8$ of the block as shown in Figure 15. The energy function (3.20) is considered for the block with $\mu=1\,\mathrm{N}/\mathrm{mm}^{2}$ for the matrix, and $\mu=\bar{\mu}$ for the spherical inhomogeneity. $C(J)=J-1$ is used for imposing the incompressibility constraint. We study four different cases: (i) a homogeneous block with $\bar{\mu}=1\,\mathrm{N}/\mathrm{mm}^{2}$, (ii) a very soft inhomogeneity with $\bar{\mu}=1e-5\,\mathrm{N}/\mathrm{mm}^{2}$, (iii) a reinforced block with $\bar{\mu}=4\,\mathrm{N}/\mathrm{mm}^{2}$, and (iv) a rigid inhomogeneity with $\bar{\mu}=1e5\,\mathrm{N}/\mathrm{mm}^{2}$. Figure 16 illustrates the convergence of the $L^{2}$-norm of the field variables $(\boldsymbol{U}_{h},\boldsymbol{K}_{h},\boldsymbol{P}_{h},p_{h})$ calculated in the matrix for all the four cases (the values of $p_{h}$ become disproportionately large in the inhomogeneity for case (iv)). One can see that a significant change in the material properties of the inhomogeneity only slightly changes the convergence of the method.

Figure 17 shows the deformed configurations of $1/8$ of the block for all the four cases for a mesh consisting of $5450$ elements. This corresponds to the last points on the convergence graphs given in Figure 16. Colors indicate the values of $\|\boldsymbol{K}_{h}\|$, $\|\boldsymbol{p}_{h}\|$, and $p_{h}$ in the first, second, and third row, respectively, where the lighter colors are associated with larger values. As expected, the values of $\|\boldsymbol{K}_{h}\|$ ($\|\boldsymbol{P}_{h}\|$) in the inhomogeneity decrease (increase) as the inhomogeneity becomes stiffer. In contrast to case (i), one can see a discontinuous change of color from the matrix to the inhomogeneity in cases (ii)-(iv). As expected, the values of $\boldsymbol{K}_{h}$, $\boldsymbol{P}_{h}$, and $p_{h}$ are continuous at the interface of the two regions in case (i) (homogeneous block) but they are discontinuous in cases (ii)-(iv) (heterogeneous blocks). Nevertheless, in all the four cases, the interface conditions are satisfied, i.e., $\boldsymbol{K}_{h}\boldsymbol{T}$ and $\boldsymbol{P}_{h}\boldsymbol{N}$ are continuous at the interface of the two regions, where $\boldsymbol{T}$ and $\boldsymbol{N}$ are respectively a tangent vector field and a normal vector field on the interface. For case (ii), one observes that $\|\boldsymbol{p}_{h}\|$ is almost uniformly zero in the spherical inhomogeneity. Hence, the traction field on the interface of the two regions is zero as well, which must be the case as a very soft inhomogeneity behaves like a hole. We solved another example by considering a block with the same geometry and the same boundary conditions but with an actual hole. It was observed that the $L^{2}$-norm of all the four field variables are equal to those calculated in the matrix for the case (ii).

Example 5: Stretching a Block with Randomly Distributed Holes.

Next, we assess the performance of CSFEM for very large strains in a complex geometry. Let us consider an incompressible cubic block of edge $1\,\mathrm{mm}$ with $6$ spherical holes as shown in Figure 18. The coordinates of the centers of the holes are $(0.25,0.6,0.6)$, $(0.7,0.5,0.3)$, $(0.6,0.2,0.7)$, $(0.2,0.2,0.2)$, $(0.3,0.8,0.2)$, $(0.8,0.75,0.7)$ and their diameters are respectively $0.4$, $0.4$, $0.3$, $0.3$, $0.3$, $0.3$. The left face of the block is fixed, the right face is subjected to a uniform displacement boundary $(u,0,0)$, and the other four faces are traction free. We use the energy function (3.20) with $\mu=1\,\mathrm{N}/\mathrm{mm}^{2}$, and $C(J)=J-1$ to impose the incompressibility constraint. The reference and the deformed configurations of the block obtained using (3.14) for $u=2\,\mathrm{mm}$ are shown in Figure 18. The mesh consists of $11756$ elements and colors indicate the values of $\|\boldsymbol{K}_{h}\|$ with lighter colors corresponding to larger values. Note that this result corresponds to the last points on the convergence graphs given in Figure 19. One can see that all the holes are stretched severely along the $x$-axis. Hence, relative to the $x$-axis, the beginning and the end portions of the boundary of each hole have the lower values of $\|\boldsymbol{K}_{h}\|$ while the middle portion has the larger values of $\|\boldsymbol{K}_{h}\|$. Figure 19 illustrates the convergence of (3.14) for different values of the displacement boundary condition $(u,0,0)$ imposed on the right face of the block. For all values of $u$, one observes that CSFEM given in (3.14) has good convergence considering all the four independent variables $(\boldsymbol{U}_{h},\boldsymbol{K}_{h},\boldsymbol{P}_{h},p_{h})$.

5 Concluding Remarks

A new mixed finite element method for 3D compressible and incompressible nonlinear elasticity was introduced. This work is an extension of [1] and [2] to three-dimensional nonlinear elasticity problems. We proposed a new four-field mixed formulation for incompressible nonlinear elasticity in terms of the displacement $\boldsymbol{U}$, the displacement gradient $\boldsymbol{K}$, the first Piola-Kirchhoff stress $\boldsymbol{P}$, and a pressure-like field $p$. By setting $p=0$ in this formulation, one can readily obtain a three-filed mixed formulation for compressible solids. In the present formulation it is assumed that $(\boldsymbol{U},\boldsymbol{K},\boldsymbol{P},p)\in H^{1}(T\mathcal{B})\times H^{\mathbf{c}}(\mathcal{B})\times H^{\mathbf{d}}(\mathcal{B})\times L^{2}(\mathcal{B})$. The new formulation has some additional terms compared with those used for 2D finite elements in [2] that vanish for the exact solutions. Provided with a proper discretization, the extra terms improve the stability of the resulting mixed finite element methods without compromising consistency. To obtain the mixed finite element methods, first four conforming finite element spaces were defined and then were used for approximating the four field variables. The discrete fields of the CSDEMs are: $\boldsymbol{U}_{h}\in V_{h,2}^{1}\subset H^{1}(T\mathcal{B}_{h})$, $\boldsymbol{K}_{h}\in\overline{V}_{h,3}^{\mathbf{c}}\subset H^{\mathbf{c}}(\mathcal{B}_{h})$, $\boldsymbol{P}_{h}\in V_{h,1}^{\mathbf{d}-}\subset H^{\mathbf{d}}(\mathcal{B}_{h})$, and $p_{h}\in V_{h,0}^{\ell}\subset L^{2}(\mathcal{B}_{h})$. The discrete spaces $V_{h,2}^{1}$, $V_{h,1}^{\mathbf{d}-}$, and $V_{h,0}^{\ell}$ are constructed using the second-order Lagrange elements, the first-order Nédélec $1^{\mathrm{st}}$-kind face elements, and the piecewise constant elements, respectively. The discrete space $\overline{V}_{h,3}^{\mathbf{c}}$ is constructed using the first-order Nédélec $2^{\mathrm{nd}}$-kind edge elements and is enriched by volume-based third-order shape functions of Nédélec $1^{\mathrm{st}}$-kind edge elements. Due to interelement continuities of these conforming spaces, our proposed mixed methods by construction provide a continuous approximation of the displacement field and satisfy both the Hadamard jump condition and the continuity of traction at the discrete level. We solved several 3D numerical examples using CSFEMs. Our observations indicate that CSFEMs have a robust performance for bending, tension, and compression problems, and in the near-incompressible and the incompressible regimes. They are also capable of modeling problems with very large strains and accurately approximating stresses. Moreover, they seem to be free from numerical artifacts such as checkerboarding of pressure, hourglass instability, and locking.

Acknowledgments.

This research was supported by AFOSR – Grant No. FA9550-12-1-0290.