Hyperelasticity background

This example demonstrates the solution of a three-dimensional elasticity problem. In addition to illustrating how to use FunctionSpaces, Expressions and how to apply Dirichlet boundary conditions, it focuses on how to:

• Minimise a non-quadratic functional
• Use automatic computation of the directional derivative
• Solve a nonlinear variational problem
• Define compiled sub-domains
• Use specific form compiler optimization options

Equation and problem definition

By definition, boundary value problems for hyperelastic media can be expressed as minimisation problems, and the minimization approach is adopted in this example. For a domain \(\Omega \subset \mathbb{R}^{d}\), where \(d\) denotes the spatial dimension, the task is to find the displacement field \(u: \Omega \rightarrow \mathbb{R}^{d}\) that minimises the total potential energy \(\Pi\):

\[\min_{u \in V} \Pi,\]

where \(V\) is a suitable function space that satisfies boundary conditions on \(u\). The total potential energy is given by

\[\Pi = \int_{\Omega} \psi(u) \, {\rm d} x - \int_{\Omega} B \cdot u \, {\rm d} x - \int_{\partial\Omega} T \cdot u \, {\rm d} s,\]

where \(\psi\) is the elastic stored energy density, \(B\) is a body force (per unit reference volume) and \(T\) is a traction force (per unit reference area).

At minimum points of \(\Pi\), the directional derivative of \(\Pi\) with respect to change in \(u\)

(1)\[L(u; v) = D_{v} \Pi = \left. \frac{d \Pi(u + \epsilon v)}{d\epsilon} \right|_{\epsilon = 0}\]

is equal to zero for all \(v \in V\):

\[L(u; v) = 0 \quad \forall \ v \in V.\]

To minimise the potential energy, a solution to the variational equation above is sought. Depending on the potential energy \(\psi\), \(L(u; v)\) can be nonlinear in \(u\). In such a case, the Jacobian of \(L\) is required in order to solve this problem using Newton’s method. The Jacobian of \(L\) is defined as

(2)\[a(u; du, v) = D_{du} L = \left. \frac{d L(u + \epsilon du; v)}{d\epsilon} \right|_{\epsilon = 0} .\]

Elastic stored energy density

To define the elastic stored energy density, consider the deformation gradient \(F\)

\[F = I + \nabla u,\]

the right Cauchy-Green tensor \(C\)

\[C = F^{T} F,\]

and the scalars \(J\) and \(I_{c}\)

\[\begin{split}J &= \det(F), \\ I_{c} &= {\rm trace}(C).\end{split}\]

This demo considers a common neo-Hookean stored energy model of the form

\[\psi = \frac{\mu}{2} (I_{c} - 3) - \mu \ln(J) + \frac{\lambda}{2}\ln(J)^{2},\]

where \(\mu\) and \(\lambda\) are the Lame parameters. These can be expressed in terms of the more common Young’s modulus \(E\) and Poisson ratio \(\nu\) by:

\[\lambda = \frac{E \nu}{(1 + \nu)(1 - 2\nu)}, \quad \quad \mu = \frac{E}{2(1 + \nu)} .\]

Demo parameters

We consider a unit cube domain:

• \(\Omega = (0, 1) \times (0, 1) \times (0, 1)\) (unit cube)

We use the following definitions of the boundary and boundary conditions:

• \(\Gamma_{D_{0}} = 0 \times (0, 1) \times (0, 1)\) (Dirichlet boundary)

• \(\Gamma_{D_{1}} = 1 \times (0, 1) \times (0, 1)\) (Dirichlet boundary)

• \(\Gamma_{N} = \partial \Omega \backslash \Gamma_{D}\) (Neumann boundary)

• On \(\Gamma_{D_{0}}\): \(u = (0, 0, 0)\)

• On \(\Gamma_{D_{1}}\)

  \[\begin{split}u = (&0, \\ &(0.5 + (y - 0.5)\cos(\pi/3) - (z - 0.5)\sin(\pi/3) - y)/2, \\ &(0.5 + (y - 0.5)\sin(\pi/3) + (z - 0.5)\cos(\pi/3) - z))/2)\end{split}\]

• On \(\Gamma_{N}\): \(T = (0.1, 0, 0)\)

These are the body forces and material parameters used:

• \(B = (0, -0.5, 0)\)
• \(E = 10.0\)
• \(\nu = 0.3\)

With the above input the solution for \(u\) will look as follows: