Finite Element Method (FEM)  Corotated Linear Elasticity Model


FEM algorithm:This example shows a finite element simulation using the corotated linear elasticity model with triangular elements [SB12]:
1. Compute deformation gradientThe 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 rotationTo 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 PiolaKirchhoff stress tensorIn this example we use the corotated linear elasticity constitutive model which defines the PiolaKirchhoff 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 forcesFinally, 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 integrationFinally, 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
