Weilin Deng | Miscellaneous notes on FEM

$\def\bs#1{\boldsymbol #1} \def\hbs#1{\hat{\boldsymbol {#1}}} \def\ul#1{\underline{\underline{#1}}} \newcommand{\pd}[2]{\frac{\partial #1}{\partial #2}}$

The matrix notations of Elasticity

Let the stress–strain ($\bs{\sigma}$–$\bs{\epsilon}$) relation ${\sigma} = [D]{\epsilon}$ be denoted as $\begin{align} \begin{Bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \sigma_{12} \\ \sigma_{23} \\ \sigma_{31} \end{Bmatrix} = \frac{E(1-\nu)}{(1+\nu)(1-2\nu)} \begin{bmatrix} 1 & & & & & \\ \frac{\nu}{1-\nu} & 1 & & & & \\ \frac{\nu}{1-\nu} & \frac{\nu}{1-\nu} & 1 & & & \\ 0 & 0 & 0 & \frac{1-2\nu}{2(1-\nu)} & & \\ 0 & 0 & 0 & 0 & \frac{1-2\nu}{2(1-\nu)} & \\ 0 & 0 & 0 & 0 & 0 & \frac{1-2\nu}{2(1-\nu)} \end{bmatrix} \begin{Bmatrix} \epsilon_{11} \\ \epsilon_{22} \\ \epsilon_{33} \\ \gamma_{12} \\ \gamma_{23} \\ \gamma_{31} \end{Bmatrix}, \label{eq:stress_strain} \end{align}$

where

$\begin{equation} \begin{Bmatrix} \epsilon_{11} \\ \epsilon_{22} \\ \epsilon_{33} \\ \gamma_{12} \\ \gamma_{23} \\ \gamma_{31} \end{Bmatrix} = \begin{Bmatrix} \pd{u}{x} \\ \pd{v}{y} \\ \pd{w}{z} \\ \pd{u}{y}+\pd{v}{x} \\ \pd{v}{z}+\pd{w}{y} \\ \pd{w}{x}+\pd{u}{z} \end{Bmatrix}. \label{eq:strain} \end{equation}$

For the plain strain problem, \eqref{eq:stress_strain} is reduced to

$\begin{equation} \begin{Bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{12} \end{Bmatrix} = \frac{E(1-\nu)}{(1+\nu)(1-2\nu)} \begin{bmatrix} 1 & & \\ \frac{\nu}{1-\nu} & 1 & \\ 0 & 0 & \frac{1-2\nu}{2(1-\nu)} \end{bmatrix} \begin{Bmatrix} \epsilon_{11} \\ \epsilon_{22} \\ \gamma_{12} \end{Bmatrix}. \label{eq:plain_strain} \end{equation}$

For the plain stress problem, we only need to replace $E$ with $E/(1-\nu^2)$ and $\nu$ with $\nu/(1-\nu)$ in \eqref{eq:plain_strain}.

In the case of axis-symmetric problem, we adopt the cylindrical coordinate and let $z$ be the symmetric axis. The non-zero components of stress are

$\begin{equation} {\sigma} = [\sigma_{rr}, \sigma_{\theta\theta}, \sigma_{zz}, \sigma_{rz}]^T \label{eq:axisym_sigma} \end{equation}$

and the non-zero components of strain are $\begin{align} {\epsilon} = [\epsilon_{rr}, \epsilon_{\theta\theta}, \epsilon_{zz}, \gamma_{rz}]^{T} = [\pd{u}{r}, \frac{u}{r}, \pd{w}{z}, \pd{u}{z}+\pd{w}{r}]^T. \label{eq:axisym_epsilon} \end{align}$

The elastic matrix $[D]$ is

$\begin{equation} [D] = \frac{E(1-\nu)}{(1+\nu)(1-2\nu)} \begin{bmatrix} 1 & & & \\ \frac{\nu}{1-\nu} & 1 & & \\ \frac{\nu}{1-\nu} & \frac{\nu}{1-\nu} & 1 & \\ 0 & 0 & 0 & \frac{1-2\nu}{1-\nu} \end{bmatrix}. \label{eq:axis_sysmmetric} \end{equation}$

Element stiffness matrix and force vector (2D)

Plain strain problem

Let $n$ be the number of node in an element. We adopt the isoparametric element, that is, the coordinate transformation and the displacement interpolation use the same shape functions. More precisely, the coordination transformation that transforms a global element to a standard, parent element through the isoparametric mapping is

$\begin{equation} x = \sum_{i=1}^{n}N_i(\xi,\eta) x_i, \quad y = \sum_{i=1}^{n}N_i(\xi,\eta) y_i, \label{eq:xy2xieta} \end{equation}$

where $(x_i,y_i)$ are the global coordinate of the $i_{th}$ node, $N_i$ is the shape function. The displacements are interpolated as

$\begin{equation} u(x,y) = \sum_{i=1}^{n}N_i(x,y) u_i, \quad v(x,y) = \sum_{i=1}^{n}N_i(x,y) v_i. \label{eq:uv_interpolate} \end{equation}$

where $(u_i,v_i)$ are the displacement of the $i_{th}$ node.

Now suppose $n = 4$ and consider the plain strain case. By substituting the interpolated displacements \eqref{eq:uv_interpolate} into \eqref{eq:strain} for the 2D case we obtain

$\begin{equation} \epsilon = [B] \delta^e = [B_1, B_2, B_3, B_4] \delta^e, \end{equation}$

where

$\begin{equation} [B_i] = \begin{bmatrix} N_{i,x} & 0 \\ 0 & N_{i,y} \\ N_{i,y} & N_{i,x} \end{bmatrix}, \end{equation}$

$\delta^e = [\delta_1^T, \delta_2^T, \delta_3^T, \delta_4^T]$, and $\delta_i = [u_i, v_i]^T$. Thus from \eqref{eq:plain_strain} we have

$\begin{equation} \sigma = [D][B]\delta^e. \end{equation}$

The element stiffness matrix is

$\begin{equation} [k] = \iint [B]^T [D] [B] \, dxdy = \int_{-1}^1\int_{-1}^1 [B]^T [D] [B] j \, d\xi d\eta, \label{eq:elem_stiff} \end{equation}$

where $j$ is the Jacobian determinant of

$\begin{equation} [J] = \begin{bmatrix} x_{,\xi} & y_{,\xi} \\ x_{,\eta} & y_{,\eta} \end{bmatrix}. \end{equation}$

Note that

$\begin{equation} \begin{Bmatrix} N_{i,x} \\ N_{i,y} \end{Bmatrix} = [J]^{-1} \begin{Bmatrix} N_{i,\xi} \\ N_{i,\eta} \end{Bmatrix}. \end{equation}$

We compute the element stiffness matrix \eqref{eq:elem_stiff} numerically using the Gauss–Legendre quadrature rule, i.e.,

$\begin{equation} [k] = \sum_{i=1}^{p} w_i [B]^T|_{(\xi_i,\eta_i)} [D] [B]|_{(\xi_i,\eta_i)} j|_{(\xi_i,\eta_i)}, \label{eq:elem_stiff_numerical} \end{equation}$

where $(\xi_i,\eta_i)$ and $w_i$ are the quadrature point and weight, respectively.

The force vector due to the body force is

$\begin{equation} [f] = \iint [N]^T [b] \, dxdy = \int_{-1}^1 \int_{-1}^1 [N]^T [b] j \, d\xi d\eta, \label{eq:elem_force} \end{equation}$

where $[N] = [N^1, N^2, N^3, N^4]$ and

$\begin{equation} [N^i] = \begin{bmatrix} N_i & 0 \\ 0 & N_i \end{bmatrix}, \quad [b] = [b_x, b_y]^T. \end{equation}$

Axis-symmetric problem

For the axis-symmetric problem, the coordinate transformation and displacement interpolation are

$\begin{equation} \begin{aligned} r &= \sum_{i=1}^{n}N_i(\xi,\eta) r_i, \quad z = \sum_{i=1}^{n}N_i(\xi,\eta) z_i, \\ u(r,z) &= \sum_{i=1}^{n}N_i(r,z) u_i, \quad w(r,z) = \sum_{i=1}^{n}N_i(r,z) w_i. \end{aligned} \end{equation}$

From \eqref{eq:axisym_epsilon} we have

$\begin{equation} [B_i] = \begin{bmatrix} N_{i,r} & 0 \\ \frac{N_i}{r} & 0 \\ 0 & N_{i,z} \\ N_{i,z} & N_{i,r} \end{bmatrix} \end{equation}$

and $\delta_i = [u_i, w_i]^T$. The element stiffness matrix becomes

$\begin{equation} \begin{aligned} [k] &= 2\pi \iint [B]^T [D] [B] r\, drdz = 2\pi \int_{-1}^1\int_{-1}^1 [B]^T [D] [B] r j \, d\xi d\eta, \\ &= 2\pi \sum_{i=1}^{p} w_i [B]^T|_{(\xi_i,\eta_i)} [D] [B]|_{(\xi_i,\eta_i)} r|_{(\xi_i,\eta_i)} j|_{(\xi_i,\eta_i)}. \label{eq:axisym_elem_stiff_numerical} \end{aligned} \end{equation}$

The force vector due to the body force becomes

$\begin{equation} [f] = 2\pi \iint [N]^T [b] r \, drdz = 2\pi \int_{-1}^1 \int_{-1}^1 [N]^T [b] r j \, d\xi d\eta. \label{eq:axisym_elem_force} \end{equation}$

Numerical integration

A good introduction to numerical introduction can be found in Sec. 5.5, p.455 of the book Finite Element Procedures (second edition) by Klaus-Jurgen Bathe. A polynomial up to the order $2N-1$ is integrated exactly when using a $N$-node Gauss quadrature. If a lower-order Gauss quadrature is used to compute the stiffness matrix, it can lead the stiffness matrix to be singular. In that case, a spurious zero energy mode or hourglass mode may exist. That is to say, there exists some particular displacement mode that has zero strain at the Gauss quadrature points. See more detailed discussion of spurious zero energy mode on p.239 of the book The Finite Element Method: Linear Static and Dynamic Finite Element Analysis by T. J. R. Hughes or here.

Nonlinear body force field

Consider a linear elastic, boundary value problem with the nonlinear body force $\bs{b}$ dependent on the displacement $\bs{u}$. the governing equation reads $\begin{align} \nabla \cdot \bs{\sigma}(\bs{u}) + \bs{b}(\bs{u}) = \bs{o}, \quad \text{in} \quad \Omega \label{eq:equilbrium} \end{align}$

with the boundary conditions

$\begin{equation} \begin{aligned} \bs{u} &= \bs{u}^*, \quad \text{on} \quad \Gamma_u, \\ \bs{\sigma} \cdot \bs{n} &= \bs{t}^*, \quad \text{on} \quad \Gamma_s. \label{eq:BC} \end{aligned} \end{equation}$

From the Principle of virtual work, we have

$\begin{equation} \begin{aligned} \int_{\Omega} \bs{\sigma}(\bs{u}):\delta\bs{\epsilon} \, dV &= \int_{\Gamma_s} \bs{t}^* \cdot \delta\bs{u} \, dS + \int_{\Omega} \bs{b}(\bs{u}) \cdot \delta\bs{u}\, dV. \end{aligned} \label{eq:virtualwork} \end{equation}$

where $\delta\bs{u}$ is a virtual displacement filed (test solution) which satisfies $\delta\bs{u} = \bs{o}$ on $\Gamma_u$, and $\delta\bs{\epsilon} = \text{sym}(\nabla\delta\bs{u})$. Since $\bs{\sigma} = \bs{C}:\nabla\bs{u}$, the \eqref{eq:virtualwork} becomes

$\begin{align} \int_{\Omega} \nabla \delta\bs{u}:\bs{C}:\nabla\bs{u} \, dV &= \int_{\Gamma_s} \bs{t}^* \cdot \delta\bs{u} \, dS + \int_{\Omega} \bs{b}(\bs{u}) \cdot \delta\bs{u}\, dV. \label{eq:virtualwork1} \end{align}$

We use the Newton-Raphson method to solve the nonlinear problem \eqref{eq:virtualwork1}. Suppose the displacement $\bs{u} + d\bs{u}$ satisfies \eqref{eq:virtualwork1}, then we have

$\begin{align} \int_{\Omega} \nabla \delta\bs{u}:\bs{C}:(\nabla\bs{u}+\nabla d\bs{u}) \, dV &= \int_{\Gamma_s} \bs{t}^* \cdot \delta\bs{u} \, dS + \int_{\Omega} \bs{b}(\bs{u}+d\bs{u}) \cdot \delta\bs{u}\, dV. \label{eq:virtualwork2} \end{align}$

To the first order approximation, in the index form, the \eqref{eq:virtualwork2} can be rewritten as

$\begin{align} \int_{\Omega} C_{ijkl} \delta u_{i,j} (u_{k,l} + du_{k,l}) \, dV &= \int_{\Gamma_s} t^*_i \delta u_i \, dS + \int_{\Omega} \delta u_i b_{i}(u_j) \, dV + \int_{\Omega} \delta u_i \pd{b_i}{u_j} du_j \, dV. \label{eq:virtualwork3} \end{align}$

Note that $C_{ijkl} = \lambda \delta_{ij} \delta_{kl} + \mu (\delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk})$. According to the Galerkin method in FEM, we choose

$\begin{equation} \delta u_i = \sum_{a=1}^{n} N^{(a)} \delta u_i^{(a)}, \quad u_j = \sum_{b=1}^{n} N^{(b)} u_j^{(b)}, \end{equation}$

and substitute $\delta u_i$ and $u_j$ into \eqref{eq:virtualwork3} to obtain

$\begin{equation} \begin{aligned} \int_{\Omega} C_{ijkl} \delta u^{(a)}_{i} (u^{(b)}_{k} + du^{(b)}_{k}) N^{(a)}_{,j} N^{(b)}_{,l} \, dV =& \int_{\Gamma_s} t^*_i \delta u^{(a)}_i N^{(a)} \, dS + \int_{\Omega} \delta u^{(a)}_i N^{(a)} b_{i}(N^{(b)} u_j^{(b)}) \, dV \\ &+ \int_{\Omega} \delta u^{(a)}_i \pd{b_i}{u_j} du^{(b)}_j N^{(a)} N^{(b)} \, dV. \end{aligned} \label{eq:virtualwork4} \end{equation}$