Finite Element Method (FEM) - Corotated Linear 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 corotated linear elasticity model with triangular elements [SB12]:
  1. compute deformation gradient for each element
  2. extract rotation from deformation gradient 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. Extract rotation

To extract the rotational part from the deformation gradient, a polar decomposition is performed $$\mathbf F = \mathbf R \mathbf S$$ where $\mathbf R$ is an orthogonal matrix and $\mathbf S$ a symmetric matrix. If the element is not inverted, $\mathbf R$ is a rotation matrix. Otherwise $\mathbf R$ contains a reflection.

3. Compute Piola-Kirchhoff stress tensor

In this example we use the corotated linear elasticity constitutive model which defines the Piola-Kirchhoff stress tensor as $$\mathbf P (\mathbf F) = 2 \mu (\mathbf F - \mathbf R) + \lambda \text{tr}(\mathbf R^T \mathbf F - \mathbf I) \mathbf R,$$ 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