An Abstract Stabilization Method with Applications to Nonlinear Incompressible Elasticity

Many finite element schemes perform very well in terms of both accuracy and stability for linear elastic problems (see [1, 2, 3, 4, 5, 6, 7, 8, 9]), including in highly constrained situations such as incompressible cases. However, it is well-known that though such schemes can be extended to nonlinear elasticity cases (for instance, large deformation elasticity problems) the same level of stability is by no means guaranteed (see [10, 11, 12, 13]).

In [14, 15], a stabilized discontinuous Galerkin method for nonlinear compressible elasticity was introduced. Because the stabilization term adapts to the solution of the problem by locally changing the size of a penalty term on the appearance of discontinuities, it is called an adaptive stabilization strategy. This same adaptive stabilization strategy can also be found in [16]. Although all three papers discuss compressible and nearly incompressible nonlinear elasticity problems, they do not address the exactly incompressible problems.

In [17], an isogeometric “stream function” formulation was used to exactly enforce the linearized incompressibility constraint. In fact, the exact satisfaction of the linearized incompressibility constraint leads to a numerical approximation of the stability range that converges to the exact one (i.e., the one for the continuum problem). As the stream function is given by a fourth-order PDE, the high regularity of the NURBS shape functions [18] must be used there and the stream function in 3dimensions is not unique.

Brezzi, Fortin, and Marini [19] presented a modified mixed formulation for second-order elliptic equations and linear elasticity problems that automatically satisfied the coercivity condition on the discrete level. Using the same technique and applying stable Stokes elements to a modified formulation for Darcy-Stokes-Brinkman models, Xie, Xu and Xue [20] obtained a uniformly stable method.

In this paper, motivated by the modified formulation used in [19, 20] for second-order equations and the Darcy-Stokes-Brinkman models respectively, we devise a stabilization strategy for the classical mixed finite method designed for Stokes equations and obtain a modified method for nonlinear incompressible elasticity. As indicated in [10], a usual mixed finite element method that is stable for Stokes equations such as the MINI element is not stable for nonlinear incompressible elasticity even though the continuous problem is stable. We note that, for the usual mixed finite method, the unphysical instability is caused by the fact that the discrete solution is not exactly divergence-free. Hence, we rewrite the continuous problem. Based on that we design a stabilization strategy involving the divergence of the discrete solution. We prove that the modified mixed finite element method can remove the unphysical instability and lead to optimal convergence. Hence, all stable Stokes elements can also be made stable for discrete nonlinear elasticity problems.

The rest of the paper is organized as follows. In section 2, we present the finite strain incompressible elasticity problem and obtain the linearized formulation of the nonlinear elasticity problem. In section 3, we propose the stabilization strategy and derive the modified mixed finite element method in an abstract framework, and in section 4 we apply this modified method to the nonlinear incompressible elasticity problem. Furthermore, we obtain an optimal approximation of the modified finite element method. In section 5, we give two simple examples and the corresponding numerical experiments to verify the stability and accuracy of the modified finite element method. We present our conclusions in section 6.

To study the finite strain incompressible elasticity problem, we adopt what is known as the material description in this paper. Given a reference configuration $\Omega\subset R^{d}~{}(1\leq d\leq 3)$ for a d-dimensional bounded material body $\mathcal{B}$, the deformation of $\mathcal{B}$ can be described by the map $\hat{\varphi}:\Omega\rightarrow R^{d}$ defined by

$$\hat{\varphi}(X)=X+\hat{u}(X),$$

where $X=(X_{1},\cdots,X_{d})$ denotes the coordinates of a material point in the reference configuration and $\hat{u}(X)$ represents the corresponding displacement vector. Following the standard notation (see [17]), we introduce the deformation gradient $\hat{F}$ and the right Cauchy–Green deformation tensor $\hat{C}$ by setting

$$\hat{F}=I+\nabla\hat{u},~{}~{}~{}~{}\hat{C}=\hat{F}^{T}\hat{F},$$

(2.1)

where I is the second-order identity tensor and $\nabla$ is the gradient operator with respect to the coordinate $X$.

For a homogeneous neo-Hookean material, for example, latex and rubber, we define (see [21]) the potential energy function as

$$\psi(\hat{u})=\frac{1}{2}\mu[I:\hat{C}-d]-\mu\ln\hat{J}+\frac{\lambda}{2}(\ln\hat{J})^{2},$$

(2.2)

where $\lambda$ and $\mu$ are positive constants, “:”represents the usual inner product for second-order tensors, and $\hat{J}=\det\hat{F}$ is the deformation gradient Jacobian.

When we introduce he pressure-like variable (or simply pressure) $\hat{p}=\lambda\ln\hat{J}$, the potential energy (2.2) can be equivalently written as the following function of $\hat{u}$ and $\hat{p}$ (still denoted as $\psi$ with a little abuse of notation):

$$\psi(\hat{u},\hat{p})=\frac{1}{2}\mu[I:\hat{C}-d]-\mu\ln\hat{J}+\hat{p}\ln\hat{J}-\frac{1}{2\lambda}(\hat{p})^{2}.$$

When the body $\mathcal{B}$ is subject to a given load $b=b(X)$ per unit volume in the reference configuration, the total elastic energy functional reads as

$$\Pi(\hat{u},\hat{p})=\displaystyle\int_{\Omega}\psi(\hat{u},\hat{p})-\displaystyle\int_{\Omega}b\cdot\hat{u}.$$

(2.3)

According to the standard variational principles, the equilibrium is derived by searching for critical points of (2.3) in suitable admissible displacement and pressure spaces $\hat{V}$ and $\hat{Q}$. The corresponding Euler–Lagrange equations arising from (2.3) lead to this solution: Find $(\hat{u},\hat{p})\in\hat{V}\times\hat{Q}$ such that

$$\left\{\begin{array}[]{l}\mu\displaystyle\int_{\Omega}\hat{F}:\nabla v+\displaystyle\int_{\Omega}(\hat{p}-\mu)\hat{F}^{-T}:\nabla v=\displaystyle\int_{\Omega}b\cdot v~{}~{}\forall v\in V,\\ \displaystyle\displaystyle\int_{\Omega}(\ln\hat{J}-\frac{\hat{p}}{\lambda})q=0~{}~{}\forall q\in Q,\end{array}\right.$$

(2.4)

where $V$ and $Q$ are the admissible variation spaces for the displacements and pressures, respectively. Also note that in (2.4), the linearization of the deformation gradient Jacobian is

$$DJ(\hat{u})[v]=J(\hat{u})F(\hat{u})^{-T}:\nabla v=J(\hat{u})F(\hat{u}):\nabla v~{}~{}~{}\forall v\in V.$$

We now focus on the case of an incompressible material, which corresponds to taking the limit $\lambda\rightarrow+\infty$ in (2.4). Hence, our problem becomes: Find $(\hat{u},\hat{p})\in\hat{V}\times\hat{Q}$ such that

$$\left\{\begin{array}[]{l}\mu\displaystyle\int_{\Omega}\hat{F}:\nabla v+\displaystyle\int_{\Omega}(\hat{p}-\mu)\hat{F}^{-T}:\nabla v=\displaystyle\int_{\Omega}b\cdot v~{}~{}\forall v\in V,\\ \displaystyle\int_{\Omega}q\ln\hat{J}=0~{}~{}\forall q\in Q,\end{array}\right.$$

(2.5)

or, in residual form: Find $(\hat{u},\hat{p})\in\hat{V}\times\hat{Q}$ such that

$$\left\{\begin{array}[]{l}\mathcal{R}_{u}((\hat{u},\hat{p}),v)=0~{}~{}\forall v\in V,\\ \mathcal{R}_{p}((\hat{u},\hat{p}),q)=0~{}~{}\forall q\in Q,\end{array}\right.$$

(2.6)

where

$$\left\{\begin{array}[]{l}\mathcal{R}_{u}((\hat{u},\hat{p}),v):=\mu\displaystyle\int_{\Omega}\hat{F}:\nabla v+\displaystyle\int_{\Omega}(\hat{p}-\mu)\hat{F}^{-T}:\nabla v-\displaystyle\int_{\Omega}b\cdot v,\\ \mathcal{R}_{p}((\hat{u},\hat{p}),q):=\displaystyle\int_{\Omega}q\ln\hat{J}.\end{array}\right.$$

(2.7)

We now derive the linearization of problem (2.5) around a generic point $(\hat{u},\hat{p})$. Observing that

$$D\hat{F}(\hat{u})[u]=-\hat{F}^{-T}(\nabla u)^{T}\hat{F}^{-T}~{}~{}\forall u\in V,$$

we easily get the problem for the infinitesimal increment $(u,p)$: Find $(u,p)\in V\times Q$ such that

$$\left\{\begin{array}[]{l}\tilde{a}(u,v)+b(v,p)=-\mathcal{R}_{u}((\hat{u},\hat{p}),v)~{}\forall v\in V,\\ \tilde{b}(u,q)=-\mathcal{R}_{p}((\hat{u},\hat{p}),q)~{}~{}\forall q\in Q,\end{array}\right.$$

(2.8)

where

$$\left\{\begin{array}[]{l}\tilde{a}(u,v)=\mu\displaystyle\int_{\Omega}\nabla u:\nabla v+\displaystyle\int_{\Omega}(\mu-\hat{p})(\hat{F}^{-1}\nabla u)^{T}:\hat{F}^{-1}\nabla v,\\ \tilde{b}(u,q)=\displaystyle\int_{\Omega}q\hat{F}^{-T}:\nabla u.\end{array}\right.$$

(2.9)

We now note that the Piola identity $\hbox{div}(\hat{J}\hat{F}^{-T})=0$ and $\hat{J}=1$ give $\hbox{div}\hat{F}^{-T}=0$. Hence, we have

$$\hbox{div}(\hat{F}^{-T}u)=\hbox{div}\hat{F}^{-T}\cdot u+\hat{F}^{-T}:\nabla u=\hat{F}^{-T}:\nabla u.$$

Let $\hat{F}^{-T}u=w$, we then consider the following problem: Find $(w,p)\in V\times Q$ such that

$$\left\{\begin{array}[]{l}a(w,v)+b(v,p)=-\mathcal{R}_{u}((\hat{u},\hat{p}),v)~{}\forall v\in V,\\ b(w,q)=-\mathcal{R}_{p}((\hat{u},\hat{p}),q)~{}~{}\forall q\in Q,\end{array}\right.$$

(2.11)

where

$$\left\{\begin{array}[]{l}a(w,v)=\tilde{a}(\hat{F}w,\hat{F}v)=\mu\displaystyle\int_{\Omega}\nabla(\hat{F}w):\nabla(\hat{F}v)+\\ ~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\displaystyle\int_{\Omega}(\mu-\hat{p})(\hat{F}^{-1}\nabla(\hat{F}w))^{T}:\hat{F}^{-1}\nabla(\hat{F}v),\\ b(w,p)=\displaystyle\int_{\Omega}p\hbox{div}w.\end{array}\right.$$

(2.12)

Let $V_{h}\subset V~{}\hbox{and}~{}Q_{h}\subset Q$ be the compatible finite element spaces. The discrete problem of (2.11) reads: Find $(w_{h},p_{h})\in V_{h}\times Q_{h}$ such that

$$\left\{\begin{array}[]{l}a(w_{h},v_{h})+b(v_{h},p_{h})=-\mathcal{R}_{u}((\hat{u},\hat{p}),v_{h})~{}\forall v_{h}\in V_{h},\\ b(w_{h},q_{h})=-\mathcal{R}_{p}((\hat{u},\hat{p}),q_{h})~{}~{}\forall q_{h}\in Q_{h},\end{array}\right.$$

(2.13)

where $a(\cdot,\cdot)~{}\hbox{and}~{}b(\cdot,\cdot)$ are defined by (2.12).

In some cases, it is possible that the continuous problem (2.11) is stable (or well-posed) but that the discrete problem (2.13) is not (see [10]). In such cases, we say that there is an unphysical instability caused by the discretization.

As shown at the end of the previous section, it is necessary to establish stable discretization when the continuous problem (2.11) is stable. So in this section, we propose an abstract stability strategy framework for the mixed approximation.

Let $V$ and $Q$ each be a Hilbert space, and assume that

$$a:V\times V\rightarrow\mathbb{R},~{}b:V\times Q\rightarrow\mathbb{R}$$

are continuous bilinear forms. Let $f\in V^{{}^{\prime}}$ and $g\in Q^{{}^{\prime}}$. We denote both the dual pairing of $V$ with $V^{{}^{\prime}}$ and that of $Q$ with $Q^{{}^{\prime}}$ by $\langle\cdot,\cdot\rangle$. We consider the following problem: Find $(w,p)\in V\times Q$ such that

$$\left\{\begin{array}[]{l}a(w,v)+b(v,p)=\langle f,v\rangle~{}~{}~{}~{}\forall v\in V,\\ b(w,q)=\langle g,q\rangle~{}~{}~{}~{}\forall q\in Q.\end{array}\right.$$

(3.1)

We associate a mapping $B$ with the form $b$:

$$B:V\rightarrow Q^{{}^{\prime}}:\langle Bw,q\rangle=b(w,q)~{}~{}\forall q\in Q.$$

Define

$$K=\hbox{ker}(B):=\{v\in V:b(v,q)=0~{}~{}\forall q\in Q\}.$$

As is well-known, the continuous problem (3.1) is well-posed under the following assumptions:

• •

  $A_{1}:$ the inf-sup condition, i.e., it holds that

  $$\inf_{q\in Q}\sup_{v\in V}\frac{b(v,q)}{\|v\|_{V}\|q\|_{Q}}\geq\beta>0,$$

  (3.2)

  where $\|\cdot\|_{V}~{}\hbox{and}~{}\|\cdot\|_{Q}$ are norms on the spaces $V$ and $Q$, respectively;

• •

  $A_{2}:$ the coercivity on the kernel space, i.e., it holds that

  $$a(v,v)\geq\alpha\|v\|_{V}^{2}~{}~{}\forall~{}v\in K,\hbox{for some}~{}\alpha>0.$$

  (3.3)

Now we rewrite the problem (3.1) in an equivalent form as follows: Find $(w,p)\in V\times Q$ such that

$$\left\{\begin{array}[]{l}A(w,v)+b(v,p)=\langle f,v\rangle+M(g,Bv)_{Q^{{}^{\prime}}}~{}~{}~{}\forall v\in V,\\ b(w,q)=\langle g,q\rangle~{}~{}~{}\forall q\in Q,\end{array}\right.$$

(3.4)

where $A(w,v)=a(w,v)+M(Bw,Bv)_{Q^{{}^{\prime}}},~{}(\cdot,\cdot)_{Q^{{}^{\prime}}}$ denotes the inner product on $Q^{{}^{\prime}}$, and $M$ is a parameter to be chosen properly.

Let $V_{h}\subset V$ and $Q_{h}\subset Q$ be finite dimensional subspaces of $V$ and $Q$ such that the following assumption is satisfied:

• •

  $A_{3}:$ discrete inf-sup condition holds that

  $$\inf_{q_{h}\in Q_{h}}\sup_{v_{h}\in V_{h}}\frac{b(v_{h},q_{h})}{\|v_{h}\|_{V}\|q_{h}\|_{Q}}\geq\beta_{1}>0,$$

  (3.5)

  .

where $\beta_{1}$ is independent of $h$.

We consider the discretization of the problem (3.4) as follows: Find $(w_{h},p_{h})\in V_{h}\times Q_{h}$ such that

$$\left\{\begin{array}[]{l}A(w_{h},v_{h})+b(v_{h},p_{h})=\langle f,v_{h}\rangle+M(g,Bv_{h})_{Q^{{}^{\prime}}}~{}\forall v_{h}\in V_{h},\\ b(w_{h},q_{h})=\langle g,q_{h}\rangle~{}~{}\forall q_{h}\in Q_{h}.\end{array}\right.$$

(3.6)

• Proof

  For any $v_{h}\in V_{h}$, noting that $V_{h}\subset V$, we have $v_{h}=\varphi+\varphi^{\bot}$, where $\varphi\in K,\varphi^{\bot}\in K^{\bot}$ and $K^{\bot}$ denotes the orthogonal of $K$ in $V$. By Lemma 3.1, we have

  	$\displaystyle A(v_{h},v_{h})$	$\displaystyle=$	$\displaystyle a(v_{h},v_{h})+M(Bv_{h},Bv_{h})_{Q^{{}^{\prime}}}$
  		$\displaystyle=$	$\displaystyle a(\varphi+\varphi^{\bot},\varphi+\varphi^{\bot})+M(B\varphi^{\bot},B\varphi^{\bot})_{Q^{{}^{\prime}}}$
  		$\displaystyle\geq$	$\displaystyle a(\varphi,\varphi)+a(\varphi,\varphi^{\bot})+a(\varphi^{\bot},\varphi)+a(\varphi^{\bot},\varphi^{\bot})+M\beta^{2}\|\varphi^{\bot}\|_{V}.$

  By assumption $A_{2}$, the continuity of $a(\cdot,\cdot)$ and the Cauchy inequality, we obtain that

  	$\displaystyle A(v_{h},v_{h})$	$\displaystyle\geq$	$\displaystyle\alpha\|\varphi\|_{V}^{2}-C_{1}\|\varphi\|_{V}\|\varphi^{\bot}\|_{V}-C_{2}\|\varphi^{\bot}\|_{V}^{2}+M\beta^{2}\|\varphi^{\bot}\|_{V}^{2}$
  		$\displaystyle\geq$	$\displaystyle(\alpha-\varepsilon C_{1})\|\varphi\|_{V}^{2}-(C_{2}+\frac{C_{1}}{\varepsilon})\|\varphi^{\bot}\|^{2}_{V}+M\beta^{2}\|\varphi^{\bot}\|_{V}^{2}$
  		$\displaystyle\geq$	$\displaystyle(\alpha-\varepsilon C_{1})\|\varphi\|_{V}^{2}+(M\beta^{2}-(C_{2}+\frac{C_{1}}{\varepsilon}))\|\varphi^{\bot}\|_{V}^{2}.$

  Noting that $\|v_{h}\|^{2}_{V}=\|\varphi\|^{2}_{V}+\|\varphi^{\bot}\|^{2}_{V}$ and choosing $\varepsilon=\frac{\alpha}{2C_{1}},M_{0}=\displaystyle\frac{\frac{\alpha}{2}+C_{2}+\frac{2C_{1}^{2}}{\alpha}}{\beta^{2}}$, we have

  $\displaystyle A(v_{h},v_{h})\geq\alpha_{1}(\|\varphi\|^{2}_{V}+\|\varphi^{\bot}\|^{2}_{V})=\alpha_{1}\|v_{h}\|_{V}^{2},$

  which completes the proof.  

Furthermore, noting that $A(\cdot,\cdot)$ is continuous, by the classic Brezzi theory for the saddle point problem [22, 23, 24], we have

We now apply the abstract stabilized framework proposed in the previous section to nonlinear incompressible elasticity problems. Suppose that $\Omega$ is connected, we consider the d-dimensional Dirichlet boundary value problem (other boundary value problems are similar) of (2.11), which means that $V=(H_{0}^{1}(\Omega))^{d},Q=L^{2}_{0}(\Omega)$. We choose $V_{h}~{}\hbox{and}~{}Q_{h}$ as finite element spaces such that $A_{3}$ is satisfied. Furthermore, $K=\{v\in(H_{0}^{1}(\Omega))^{d}:\hbox{div}v=0\}$. We rewrite (2.11) as: Find $(w,p)\in V\times Q$ such that

$$\left\{\begin{array}[]{l}A(w,v)+b(v,p)=-\mathcal{R}_{u}((\hat{u},\hat{p}),v)-M\mathcal{R}_{p}((\hat{u},\hat{p}),\hbox{div}v)~{}~{}\forall v\in V,\\ b(w,q)=-\mathcal{R}_{p}((\hat{u},\hat{p}),q)~{}~{}\forall q\in Q,\end{array}\right.$$

(4.1)

where $A(w,v)=a(w,v)+M(\hbox{div}w,\hbox{div}v)$, $a(\cdot,\cdot)~{}\hbox{and}~{}b(\cdot,\cdot)$  are defined by (2.12), $(\cdot,\cdot)$ denotes the inner product on $L^{2}(\Omega)$, and $M$ is a parameter chosen properly.

The discretization of the problem (4.1) is as follows: Find $(w_{h},p_{h})\in V_{h}\times Q_{h}$ such that

$$\left\{\begin{array}[]{l}A(w_{h},v_{h})+b(v_{h},p_{h})=-\mathcal{R}_{u}((\hat{u},\hat{p}),v_{h})-M\mathcal{R}_{p}((\hat{u},\hat{p}),\hbox{div}v_{h})~{}~{}\forall v_{h}\in V_{h},\\ b(w_{h},q_{h})=-\mathcal{R}_{p}((\hat{u},\hat{p}),q_{h})~{}~{}\forall q_{h}\in Q_{h}.\end{array}\right.$$

(4.2)

Let $K^{0}=\{y\in(H^{-1}(\Omega))^{d}:<y,\phi>=0~{}~{}\forall~{}\phi\in K\}$, where $<\cdot,\cdot>$ denotes the duality paring between $(H^{-1}(\Omega))^{d}$ and $(H^{1}_{0}(\Omega))^{d}$.

• Proof

  By the abstract framework, the proof is obvious.  

Also from the abstract framework in the previous section, we can get an approximate accuracy result similar to (3.9).

In this section, we report our numerical experiments in regard to the performance of the stabilization strategy for two simple examples. We consider here simple problems using some mixed finite element formulations that are known to be optimal for Stokes equations. We first briefly present the finite element under consideration and then show the numerical results obtained using such element.

Within the framework of finite elasticity for incompressible materials, it is well known that the classical mixed finite element discretization can sometimes be unstable even though the continuous problem is stable. In this paper, we reformulated the continuous problem and proposed an abstract stabilization strategy based on the new continuous formulation and obtained a modified mixed finite element method. We proved theoretically in Section 3 that for a sufficiently large $M$, the modified mixed finite element method is stable whenever the continuous problem is stable, and the method maintains the optimal convergence of the classical one. We verified by numerical experiments in Section 5 that the modified mixed finite element method is much more stable than the classical one and is also locking-free. The $M_{0}$ provided in Section 3 always overestimates the parameter used in practical problems. However, we can choose the parameter heuristically by analyzing the stability and continuity of continuous problems in the numerical experiments presented in Section 5.