Theory

Mixed formulation for incompressible elasticity

The governing equations and boundary conditions for the linear elastic problem are:

\[ \begin{align} -\textrm{div} \boldsymbol{\sigma(u)} = \boldsymbol{f} & \quad \textrm{in}\; \Omega \label{eq:momentum} \\ \boldsymbol{u} = \boldsymbol{\bar{u}} & \quad \textrm{on}\; \partial\Omega_u, \label{eq:dirichlet}\\ \boldsymbol{\sigma \cdot n} = \boldsymbol{t} & \quad \textrm{on}\; \partial\Omega_n, \label{eq:neumann} \end{align} \]

With Hooke's constitutive equation defined as:

\begin{align} \label{eq:constitutive}\boldsymbol{\sigma} = \lambda tr(\boldsymbol{\varepsilon(u)})\boldsymbol{I} + 2 \mu\boldsymbol{\varepsilon(u)} \end{align}

Where \(\boldsymbol{\sigma}\) is the Cauchy stress, \(\boldsymbol{\varepsilon}\) is the strain, \(\boldsymbol{\varepsilon(u)}= \frac{1}{2}(\nabla\boldsymbol{u}+[\nabla\boldsymbol{u}]^T)\). The Lame constants \(\lambda\) and \(\mu\) are related through Young's modulus \(E\) and the Poisson ratio \(\nu\):

\[ \begin{align} \label{eq:lamelambda}\lambda = \frac{\nu E}{(1 + \nu)(1 - 2\nu)} \end{align} \]

\[ \begin{align} \label{eq:lamemu}\mu = \frac{E}{2(1 + \nu)} \end{align} \]

The problem above of compressible elasticity is solved using finite elements. However, for the incompressible state (where \(\nu\) = 0.5), a mixed formulation of the linear elastic problem can be considered through the introduction of a new unknown for pressure, \(p\). With \(p = \lambda \textrm{div}\boldsymbol{u}\), the constitutive relation is now defined as:

\[ \begin{align} \label{eq:incompressible}\boldsymbol{\sigma} = 2 \mu\boldsymbol{\varepsilon} + \lambda \textrm{div}\boldsymbol{u}\boldsymbol{I} = 2 \mu\boldsymbol{\varepsilon} + p\boldsymbol{I} \end{align} \]

Therefore, we add the following equation to the strong form of the problem:

\[ \begin{align} \label{eq:incompressiblestrong}-\textrm{div} \boldsymbol{u} + \frac{1}{\lambda} p = 0 & \quad \textrm{in}\; \Omega \end{align} \]

As standard for Galerkin finite elements, we multiply Eqs. \eqref{eq:momentum} and \eqref{eq:incompressiblestrong} by the appropriate test functions \(\delta\mathbf{u}\) and \(\delta p\) and integrate over the domain \(\Omega\):

\[ \begin{align} \label{eq:momentum_int} 2\mu\int_{\Omega}\mathbf{\varepsilon (u)}:\mathbf{\varepsilon(\delta u)} \,\textrm{d}\Omega + \int_{\Omega}p \,\textrm{div}\,\delta\mathbf{u} \,\textrm{d}\Omega &= \int_{\Omega}\delta\mathbf{u}\cdot\mathbf{f}\,\textrm{d}\Omega + \int_{\partial\Omega}\delta\mathbf{u}\cdot\mathbf{t}\,\textrm{d}\partial\Omega \\ \int_{\Omega}\textrm{div}\,\mathbf{u}\,\delta p \,\textrm{d}\Omega -\int_{\Omega}\frac{1}{\lambda} p\, \delta p\,\textrm{d}\Omega &= 0 \end{align} \]

As a result, we arrive at the following global linear system of equations:

\[ \begin{align} \begin{bmatrix} {\bf{A}} & {\bf{C}} \\ {\bf{C^T}} & {\bf{B}} \end{bmatrix} \begin{bmatrix} \bf{u} \\ \bf{p} \\ \end{bmatrix} = - \begin{bmatrix} {\bf{F}} \\ {\bf{0}} \\ \end{bmatrix} \end{align} \]

where

\[ \begin{align} \mathbf{A} = 2\mu\int_{\Omega}\mathbf{\varepsilon (u)}:\mathbf{\varepsilon(\delta u)}\,\textrm{d}\Omega \end{align} \]

\[ \begin{align} \mathbf{C^T} = \int_{\Omega}\textrm{div}\,\mathbf{ u}\,\delta p\,\textrm{d}\Omega \end{align} \]

\[ \begin{align} \mathbf{B} = -\int_{\Omega}\frac{1}{\lambda} p\, \delta p\,\textrm{d}\Omega \end{align} \]

\[ \begin{align} \mathbf{F} = \int_{\Omega}\delta\mathbf{u}\cdot\mathbf{f}\,\textrm{d}\Omega + \int_{\partial\Omega}\delta\mathbf{u}\cdot\mathbf{t}\,\textrm{d}\partial\Omega \end{align} \]

In the case of incompressible material, the lower right diagonal block of our matrix, \({\bf{B}} = 0\), resulting in the well-known saddle point problem. Careful choice has to be made for the finite element approximations of the displacement and pressure fields to ensure that the discrete problem has a valid and stable solution, i.e. the inf-sup conditions have to be satisfied.

One element in particular that can be adopted is the well-known Taylor-Hood element, whereby both displacement and pressure fields are continuous and are approximated in Sobolev \(H^1\) functional space. One must note that the polynomial order for the approximation of the displacement field must be greater by one than that of the pressure field to achieve inf-sup stability. Such choice of the element is often denoted as \(P_k - P_{k-1}\) for \(k \geq 2\). See [Boffi, 2013] for more details on Taylor-Hood elements.

Another element that can be utilised is the Crouzeix-Raviart element. With this element the displacement is approximated in Sobolev \(H^1\) functional space and the pressure field is discontinuous in \(L^2\) space. Two order difference is required between displacement and pressure fields for stability, i.e. such element is denoted as \(P_k - P_{k-2}\). For the 2D case, \(k \geq 2\) is sufficient but for 3D where we require degrees of freedom on the faces, \(k \geq 3\).

References:

[Boffi, 2013] Boffi D, Brezzi F, Fortin M. Mixed finite element methods and applications, Springer, 2013