Finite Element Method (FEM) - St. Venant-Kirchhoff 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 St. Venant-Kirchhoff model with triangular elements [SB12]:
  1. compute deformation gradient for each element
  2. compute Green strain tensor 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.

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 Green strain tensor

The strain is measured using the Green strain tensor which is defined as $$\mathbf E = \frac12 (\mathbf F^T \mathbf F - \mathbf I)$$ where $\mathbf I$ is the identity matrix.

3. Compute Piola-Kirchhoff stress tensor

In this example we use the St. Venant-Kirchhoff constitutive model which defines the Piola-Kirchhoff stress tensor as $$\mathbf P (\mathbf F) = \mathbf F \left ( 2 \mu \mathbf E + \lambda \text{tr}(\mathbf E) \mathbf I \right ),$$ 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