1. Introduction

2. Governing Equations

Let the body occupied by $\Omega$ be under some external forces. Let $\Omega’$ denote the current deformed configuration. Since we wish to solve the system of equations to obtain the displacement field, we consider the system of equations in the reference configuration, and focus on the finite element formulation. The reader is referred to well-known references ^{1 2 3} for… theory, and the list is by no means (even vaguely) exhaustive.

2.1. Kinematics

We denote the position of a particle in the reference state $\Omega$ by $X \in \Omega$. In the deformed state, we denote the position by $x = \phi (X) \in \Omega’$, where $\phi : \Omega \rightarrow \mathbb{R}^3$ denotes the deformation, and $\Omega’ = \phi(\Omega)$. Further, let $F = \nabla \phi \in \mathbb{R}^{3 \times 3}$ denote the deformation gradient, and let $u(X) = \phi(X) - I$ denote the displacement.

We denote the right Cauchy-Green strain tensor by $C = F^T F$, and the Green-St-Venant strain tensor by $E = \frac{1}{2}\left( C - I \right)$, where $I \in \mathbb{R}^{3 \times 3}$ is the identity matrix. This gives us

Finally, the linearized strain is given by $\epsilon = \frac{1}{2}\left(\nabla u + \nabla u^T \right)$, and it can be seen that $\epsilon \approx E$ for small displacements $u$.

2.2. Equilibrium Equations

For a pure displacement problem, under an external body force $f : \Omega \rightarrow \mathbb{R}^3$, the equation of equilibrium in the reference configuration is given by [Ciarlet’ 1988, Theorem 2.6-1] ¹

where $T : \Omega \rightarrow \mathbb{R}^3$ is the second order symmetric stress tensor, and in particular the first Piola-Kirchhoff stress tensor. We complete the system \ref{eq:stress_eq} with homogeneous Dirichlet boundary conditions

We now list the constitutive relations for the first Piola-Kirchhoff stress tensor.

2.2.1. Hyperelasticity: St-Venant Kirchhoff Material

In this post, we consider the simplest constitutive relation for a hyperelastic material, the St-Venant Kirchhoff relation. For the St-Venant Kirchhoff material, the first Piola-Kirchhoff stress tensor is given by

where $\lambda \in \mathbb{R}$ and $\mu \in \mathbb{R}$ are Lamé parameters related to the Youngs modulus $E_Y \in \mathbb{R}$ and Poisson ratio $\nu \in \mathbb{R}$ by

Note on the choice of the St-Venant Kirchhoff material. It can be shown that the St-Venant Kirchhoff material is not the most sound way of approximating physically observed material behaviour. Indeed, it can be shown that such a material, i.e., where the the stress relation is given by \ref{eq:st_venant_fpks}, can undergo extreme deformation to arbitrarily small volumes in a finite by expending finite energy. This, however, is physically inconsistent with natural materials. Better choices include Neo-Hookean materials, Ogden… For more information, see… Here we make use of \ref{eq:st_venant_fpks} due to its simple form.

2.2.2. Linear Elasticity

In the linear elastic framework, the first Piola-Kirchhoff stress tensor is given by

2.3. Weak Formulation and Existence of Solution

We now present the weak formulation for \ref{eq:stress_eq}-\ref{eq:Dirichlet_boundary}. In this post, we consider the problem in 1D, and seek the solution $u \in H_0^1(\Omega)$. The continuous weak formulation is as follows: find $u \in H_0^1(\Omega)$ such that ¹, [Oden’ 1979, Eq. 7.1] ⁴,

The existence of a solution to \ref{eq:continuous_weak_form} is established using the implicit function theorem in [Ciarlet’ 1988, Theorem 6.4-1] for the St-Venant Kirchhoff material \ref{eq:st_venant_fpks}, and using Gårding operators in [Oden’ 1979, Theorem 7.1] under more regularity assumptions.

We now present the fully discrete formulation. In this post, we focus on 1D, and in $P_1$ elements. Let $\Omega = (0, 1)$, and let $\Omega$ be divided into $M$ cells $(x_{j}, x_{j+1})$, $0 \leq j \leq M-1$, of uniform width denoted by $h = \frac{1}{M}$. Let $V_h$ is the subspace of piecewise-linear functions with basis functions given by

The fully discrete form of \ref{eq:continuous_weak_form} reads as follows: find $u_h \in V_h$ such that

Existence of solution. Before we set up a numerical experiment, we first prove the existence of a solution to \ref{eq:discrete_weak_form} for appropriate forces $f$. We can rewrite \ref{eq:discrete_weak_form} as

where $U \in \mathbb{R}^{M-1}$ collects the entries of $u_h = \sum_{j} U_j \phi_j$, i.e., $U = [U_1 \; U_2 \; \dots \; U_{M-1}]^T$, and $\mathcal{T} : \mathbb{R}^{M-1} \rightarrow \mathbb{R}^{M-1}$ is a nonlinear map with entries

and the vector $L_f \in \mathbb{R}^{M-1}$ collects the entries

We now present a little existence result.

Theorem 2.3.1. Let $f \in C^0(\Omega)$. Then, for $\lVert f \rVert_\infty$ small enough, there exists a solution to \ref{eq:nonlinear_map_eq}.

Proof We make use of the inverse function theorem. First note that, for $U = 0$ (here we mean $0 \in \mathbb{R}^{M-1}$), we have $\mathcal{T}(U) = 0$. Now, consider the Jacobian $\mathcal{J}$ of $\mathcal{T}$

By definition, since $E(u_h) = \frac{1}{2}(F(u_h)^2 - 1)$, we have

Since $u_h = \sum_j U_j \phi_j$, by the chain rule, we have

For $U = 0$, we have $F(u_h) = 1$, and thus we have from \ref{eq:proof_jacobian} that $J_T(0) = \frac{\left(\lambda + 2\mu \right)}{h}\text{tri}(1, 2, 1)$ is a tri-diagonal matrix such that $\mathcal{J}_T(0)$ is symmetric positive definite. Hence $\mathcal{J}(0)$ is invertible. Thus, by the inverse function theorem ⁵ $\exists$ open neighborhoods $R_1, R_2 \subset \mathbb{R}^{M-1}$ $0 \in R_1$, $0 \in R_2$, $\mathcal{T}$ is one-one on $O_1$, and

That is, for any $f \in R_2$, i.e., if $\lVert f \rVert_\infty$ is small enough, $\exists U_f \in R_1$ such that $\mathcal{T}(U_f) = f$. This completes the proof.

□

Note that the above theorem does not prove the non-existence of solutions for any arbitrary $f$. In fact, in our numerical experiments we have obtained solutions for large $\lVert f \rVert_\infty$, however, another issue lurks with the nonlinear system above. A crucial point to consider now is how we still have not mentioned uniqueness for our hyperelastic system. For linear elasticity, both uniqueness and existence is well-established and follows from Korn’s inequality ¹, but for hyperelastic system this is not the case. As we shall explore below, uniqueness for hyperelastic systems indeed isn’t guaranteed and leads to spurious oscillations.

3. Numerical Implementation and Experiments

We now describe the details of our numerical implementation. We make use of Newton’s method to solve the system \ref{eq:nonlinear_map_eq}. Starting with an initial guess $U^{(0)}$, we perform the step

We perform the Newton step till we obtain convergence of the residuals $\lVert \mathcal{T}(U^{(m-1)}) - L_f \rVert_\infty < \epsilon_{tol}$ or $\lvert \delta U^{(m)} \rVert_\infty < \epsilon_{tol}$, for some prescribed tolerance $\epsilon_{tol} > 0$.

Since we are using $P_1$ elements for the displacement $u$, this means that $F$ is a piecewise-constant on each grid cell. This makes numerical integration straightforward when computing the Jacobians of $\mathcal{T}$ since from \ref{eq:proof_jacobian}

References