Finite Element Method (FEM) - Neohookean Elasticity Model

Your browser does not support the HTML5 canvas tag.
0.00 s
0.00 ms
0
0

FEM algorithm:

This example shows a finite element simulation using the Neohookean elasticity model with triangular elements [SB12]:
  1. compute deformation gradient for each element
  2. compute invariant for each element
  3. compute Piola-Kirchhoff stress tensor for each element
  4. compute elastic forces
  5. time integration to get new positions and velocities
Note that the simulation is only conditionally stable since a conditionally stable explicit time integration method is used and that inverted elements are not handled.

1. Compute deformation gradient

The deformation gradient is defined as $$\mathbf F = \frac{\partial \mathbf x(\mathbf X,t)}{\partial \mathbf X},$$ where $\mathbf X$ defines the reference configuration and $\mathbf x$ the deformed configuration of a body.

In this example we use linear triangular elements. The deformation gradient for such an element with the vertices $\mathbf x_1$, $\mathbf x_2$ and $\mathbf x_3$ is determined by $$\mathbf F = \mathbf D_s \mathbf D_m^{-1},$$ where $\mathbf D_s$ is the deformed shape matrix and $\mathbf D_m$ the constant reference shape matrix defined as $$\begin{align*} \mathbf D_s &= (\mathbf x_1 - \mathbf x_3 \quad \mathbf x_2 - \mathbf x_3) \\ \mathbf D_m &= (\mathbf X_1 - \mathbf X_3 \quad \mathbf X_2 - \mathbf X_3). \end{align*}$$ Note that $\mathbf D_m^{-1}$ can be precomputed.

2. Compute invariant

The three invariants are defined as $$I_1(\mathbf F) = \text{tr}(\mathbf F^T \mathbf F),\quad I_2(\mathbf F) = \text{tr}\left( (\mathbf F^T \mathbf F)^2 \right ), \quad I_3(\mathbf F) = \text{det}(\mathbf F^T \mathbf F.)$$

3. Compute Piola-Kirchhoff stress tensor

In this example we use the Neohookean elasticity constitutive model which defines the Piola-Kirchhoff stress tensor as $$\begin{align*} \mathbf P (\mathbf F) &= \mu \mathbf F - \mu \mathbf F^{-T} + \frac{\lambda \text{log}(I_3)}{2} \mathbf F^{-T} \\ &= \mu \mathbf F + \left (\frac{\lambda \text{log}(I_3)}{2} - \mu \right ) \mathbf F^{-T}, \end{align*}$$ where $\mu$ and $\lambda$ are the Lamé parameters. These parameters are related to the Young’s modulus $k$ (measure of stretch resistance) and the Poisson’s ratio $\nu$ (measure of incompressibility) as: $$\mu = \frac{k}{2 (1 + \nu)}, \quad\quad \lambda = \frac{k \nu}{(1+\nu)(1 - 2\nu)}.$$

4. Compute elastic forces

Finally, the elastic forces for the three nodes of our linear triangular elements can be computed as the negative gradients of the strain energy which yields $$\begin{align*} \left [ \mathbf f_1 \quad \mathbf f_2 \right ] &= -A \mathbf P(\mathbf F) \mathbf D_m^{-T} \\ \mathbf f_3 &= -\sum_{i=1}^2 \mathbf f_i, \end{align*}$$ where $A$ is the area of the triangular element in rest configuration and $\mathbf D_m^{-T}$ is the inverted and transposed the reference shape matrix.

2. Time integration

Finally, the particles are advected by numerical time integration. In our case we use a symplectic Euler method: $$\begin{align*} \mathbf v(t + \Delta t) &= \mathbf v(t) + \frac{\Delta t}{m} \left (\mathbf f(t) + \mathbf f^{\text{ext}}(t) \right ) \\ \mathbf x(t + \Delta t) &= \mathbf x(t) + \Delta t \mathbf v(t + \Delta t), \end{align*}$$ where $\mathbf f^{\text{ext}}$ are the external forces.

References

  • [SB12] Eftychios Sifakis, Jernej Barbic. FEM Simulation of 3D Deformable Solids. ACM SIGGRAPH Courses, 2012