Engineering at Alberta Courses » FEA in Two and Three Dimensions

Finite Element Analysis: FEA in Two and Three Dimensions

The basic types of elements used in two dimensions are either the quadrilateral (four-sided) or the triangular (three-sided) elements. The quadrilateral elements can either have bilinear or quadratic displacement functions. The bilinear displacement function element possesses four nodes at the four corners of the element, while the quadratic displacement function element possesses eight nodes, four at the corners of the element and four at the midpoints of the sides. The triangular elements can have either bilinear or quadratic displacement functions. The bilinear displacement function element possesses three nodes at the corners of the triangle, and the quadratic displacement function element possesses six nodes, three nodes at the corners and three nodes at the midpoints of the sides. For three dimensions, the eight-node trilinear brick element is the three-dimensional extension of the bilinear quadrilateral and the 20-node brick element is the three-dimensional extension of the quadratic quadrilateral. The four-node tetrahedron is the three-dimensional extension of the bilinear triangular element, while the ten-node tetrahedron is the three-dimensional extension of the quadratic displacement function triangular element. The choice of the element type depends on many factors. For domains with irregular geometry, the triangle or the tetrahedron could be the favorable choice because of the ease of discretizing irregular domains with triangular meshes. However, linear triangular elements (as will be shown later) tend to be stiffer than other types of elements, and thus, finer meshes or the nonlinear triangular elements should be used, which in turn require greater computational resources. The quadrilateral element can provide better accuracy with fewer computational requirements but requires regular geometry to be discretized. In the following section, the general equations for two-dimensional problems in linear elasticity are formulated for general elements. Then, the derivation of the different shape functions for the bilinear and quadratic triangular and quadrilateral
two-dimensional elements is presented.

Stiffness Matrix and Nodal Forces Vector for a General 3D Linear Elastic Element

In three dimensions, the displacement vector $u_{3D}$ of an element has three components designated $u$ , $v$ , and $w$ such that:

$u_{3D} = \begin{pmatrix} u \\ v \\ w \\ \end{pmatrix}$

Assuming that the element has $n$ nodes, then, each node $i$ has 3 nodal degrees of freedom designated $u_i$ , $v_i$ , and $w_i$ . So, the nodal degrees of freedom vector $u_e$ can have the form:

$u_e = \begin{pmatrix} u_1 \\ v_1 \\ w_1 \\ u_2 \\ v_2 \\ w_2 \\ \vdots \\ u_n \\ v_n \\ w_n \\ \end{pmatrix}$

If $\forall i\leq n:N_i$ is the shape function associated with node $i$ , then the displacement vector function has the form:

$u_{3D} = \begin{pmatrix} u \\ v \\ w \\ \end{pmatrix} = \begin{pmatrix} u_1N_1 + u_2N_2 + \cdots + u_nN_n \\ v_1N_1 + v_2N_2 + \cdots + v_nN_n \\ w_1N_1 + w_2N_2 + \cdots + w_nN_n \\ \end{pmatrix}$

$= \begin{pmatrix} N_1 & 0 & 0 & N_2 & 0 & 0 & \cdots & N_n & 0 & 0 \\ 0 & N_1 & 0 & 0 & N_2 & 0 & \cdots & 0 & N_n & 0 \\ 0 & 0 & N_1 & 0 & 0 & N_2 & \cdots & 0 & 0 & N_n \\ \end{pmatrix} \begin{pmatrix} u_1 \\ v_1 \\ w_1 \\ u_2 \\ v_2 \\ w_2 \\ \vdots \\ u_n \\ v_n \\ w_n \\ \end{pmatrix}$

$= Nu_e$

For simplicity, the vector representation of the small strain matrix will be used. In this case, the small strain matrix can be written in vector form as a function of the nodal degrees of freedom vector as follows:

$\varepsilon = \begin{pmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ 2\varepsilon_{12} \\ 2\varepsilon_{13} \\ 2\varepsilon_{23} \\ \end{pmatrix} = \begin{pmatrix} \frac{\partial u}{\partial X_1} \\ \frac{\partial v}{\partial X_2} \\ \frac{\partial w}{\partial X_3} \\ \frac{\partial u}{\partial X_2} + \frac{\partial v}{\partial X_1} \\ \frac{\partial u}{\partial X_3} + \frac{\partial w}{\partial X_1} \\ \frac{\partial v}{\partial X_3} + \frac{\partial w}{\partial X_2} \\ \end{pmatrix}$

$= \begin{pmatrix} \frac{\partial N_1}{\partial X_1} & 0 & 0 & \frac{\partial N_2}{\partial X_1} & 0 & 0 & \cdots & \frac{\partial N_n}{\partial X_1} & 0 & 0 \\ 0 & \frac{\partial N_1}{\partial X_2} & 0 & 0 & \frac{\partial N_2}{\partial X_2} & 0 & \cdots & 0 & \frac{\partial N_n}{\partial X_2} & 0 \\ 0 & 0 & \frac{\partial N_1}{\partial X_3} & 0 & 0 & \frac{\partial N_2}{\partial X_3} & \cdots & 0 & 0 & \frac{\partial N_n}{\partial X_3} \\ \frac{\partial N_1}{\partial X_2} & \frac{\partial N_1}{\partial X_1} & 0 & \frac{\partial N_2}{\partial X_2} & \frac{\partial N_2}{\partial X_1} & 0 & \cdots & \frac{\partial N_n}{\partial X_2} & \frac{\partial N_n}{\partial X_1} & 0 \\ \frac{\partial N_1}{\partial X_3} & 0 & \frac{\partial N_1}{\partial X_1} & \frac{\partial N_2}{\partial X_3} & 0 & \frac{\partial N_2}{\partial X_1} & \cdots & \frac{\partial N_n}{\partial X_3} & 0 & \frac{\partial N_n}{\partial X_1} \\ 0 & \frac{\partial N_1}{\partial X_3} & \frac{\partial N_1}{\partial X_2} & 0 & \frac{\partial N_2}{\partial X_3} & \frac{\partial N_2}{\partial X_2} & \cdots & 0 & \frac{\partial N_n}{\partial X_3} & \frac{\partial N_n}{\partial X_2} \\ \end{pmatrix} \begin{pmatrix} u_1 \\ v_1 \\ w_1 \\ u_2 \\ v_2 \\ w_2 \\ \vdots \\ u_n \\ v_n \\ w_n \\ \end{pmatrix}$

$= Bu_e$

For a general linear elastic material, the vector representation of the stress matrix is related to the vector representation of the small strain matrix by a $6\times 6$ symmetric material coefficients matrix $C$ :

$\sigma = \begin{pmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \sigma_{12} \\ \sigma_{13} \\ \sigma_{23} \\ \end{pmatrix} = C \begin{pmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ 2\varepsilon_{12} \\ 2\varepsilon_{13} \\ 2\varepsilon_{23} \\ \end{pmatrix}$

Utilizing the principle of virtual work, the virtual displacement field $u_{3D}^*$ and the associated virtual strain field $\varepsilon^*$ can have the following forms with $u_e^*$ being an arbitrary “virtual” vector of degrees of freedom of the system:

$u_{3D}^* = \begin{pmatrix} u^* \\ v^* \\ w^* \\ \end{pmatrix} = u_e^* N^T$

$\varepsilon^{*T} = u_e^* B^T$

where

$u_e^* = \hspace{1mm} < \hspace{2mm} u_1^* \hspace{2mm} v_1^* \hspace{2mm} w_1^* \hspace{2mm} u_2^* \hspace{2mm} v_2^* \hspace{2mm} w_2^* \hspace{2mm} \cdots \hspace{2mm} u_n^* \hspace{2mm} v_n^* \hspace{2mm} w_n^* \hspace{2mm} >$

$\varepsilon^{*T} = \hspace{1mm} < \hspace{2mm} \varepsilon_{11}^* \hspace{2mm} \varepsilon_{22}^* \hspace{2mm} \varepsilon_{33}^* \hspace{2mm} 2\varepsilon_{12}^* \hspace{2mm} 2\varepsilon_{13}^* \hspace{2mm} 2\varepsilon_{23}^* \hspace{2mm} >$

The internal virtual work integral term for this particular element will have the form:

$IVW_e = \int_e \sum_{i,j=1}^3 \varepsilon_{ij}^* \sigma_{ij} dx = \int_e \varepsilon^* \cdot \sigma dx$

$=u_e^* \Big( \int_e B^T CBdx \Big) u_e$

where the integration is done on each term in the matrix $B^TCB$ and is performed over the volume (in 3D problems) of the element. The local stiffness matrix has dimensions $3n\times 3n$ and has the form:

$K^e = \int_e B^T CB dx$

The external virtual work integral term for this particular element will have the form:

$EVW_e = u_e^* \int_{\partial e} N^T t_n ds + u_e^* \int_e N^T \rho b dx$

Where $t_n\in\mathbb{R}^3$ is the traction vectors on the surface of the element, $\rho$ is the mass density, and $b\in\mathbb{R}^3$ is the body forces vector on the element. The integration in the first term is done over the element surface while in the second term is done over the element volume. The local nodal forces vector will $3n$ components and will have the form:

$f^e = \int_{\partial e} N^T t_n ds + \int_e N^T \rho b dx$

Stiffness Matrix and Nodal Forces Vector for a General 2D Linear Elastic Element

The equations in the previous section are repeated after reducing them to two dimensions. In this case, the displacement vector $u_{2D}$ of an element has two components designated $u$ and $v$ such that:

$u_{2D} = \begin{pmatrix} u \\ v \\ \end{pmatrix}$

Assuming that the element has n nodes, then, each node $i$ has 2 nodal degrees of freedom designated $u_i$ and $v_i$ . So, the nodal degrees of freedom vector $u_e$ can have the form:

$u_e = \begin{pmatrix} u_1 \\ v_1 \\ u_2 \\ v_2 \\ \vdots \\ u_n \\ v_n \\ \end{pmatrix}$

If $\forall i\leq n:N_i$ is the shape function associated with node $i$ , then the displacement vector function has the form:

$u_{2D} = \begin{pmatrix} u \\ v \\ \end{pmatrix} = \begin{pmatrix} u_1N_1 + u_2N_2 + \cdots + u_nN_n \\ v_1N_1 + v_2N_n + \cdots + v_nN_n \\ \end{pmatrix}$

$= \begin{pmatrix} N_1 & 0 & N_2 & 0 & \cdots & N_n & 0 \\ 0 & N_1 & 0 & N_2 & \cdots & 0 & N_n \\ \end{pmatrix} \begin{pmatrix} u_1 \\ v_1 \\ u_2 \\ v_2 \\ \vdots \\ u_n \\ v_n \\ \end{pmatrix}$

$= Nu_e$

$\varepsilon = \begin{pmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ 2\varepsilon_{12} \\ \end{pmatrix} = \begin{pmatrix} \frac{\partial u}{\partial X_1} \\ \frac{\partial v}{\partial X_2} \\ \frac{\partial u}{\partial X_2} + \frac{\partial v}{\partial X_1} \\ \end{pmatrix}$

$= \begin{pmatrix} \frac{\partial N_1}{\partial X_1} & 0 & \frac{\partial N_2}{\partial X_1} & 0 & \cdots & \frac{\partial N_n}{\partial X_1} & 0 \\ 0 & \frac{\partial N_1}{\partial X_2} & 0 & \frac{\partial N_2}{\partial X_2} & \cdots & 0 & \frac{\partial N_n}{\partial X_2} \\ \frac{\partial N_1}{\partial X_2} & \frac{\partial N_1}{\partial X_1} & \frac{\partial N_2}{\partial X_2} & \frac{\partial N_2}{\partial X_1} & \cdots & \frac{\partial N_n}{\partial X_2} & \frac{\partial N_n}{\partial X_1} \\ \end{pmatrix} \begin{pmatrix} u_1 \\ v_1 \\ u_2 \\ v_2 \\ \vdots \\ u_n \\ v_n \\ \end{pmatrix}$

$= Bu_e$

For a general plane strain or plain stress linear elastic material, the vector representation of the stress matrix is related to the vector representation of the small strain matrix by a $3\times 3$ symmetric material coefficients matrix $C$ :

$\sigma = \begin{pmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{12} \\ \end{pmatrix} = C \begin{pmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ 2\varepsilon_{12} \\ \end{pmatrix}$

Utilizing the principle of virtual work, the virtual displacement field $u_{2D}^*$ and the associated virtual strain field $\varepsilon^*$ can have the following forms with $u_e^*$ being an arbitrary “virtual” vector of degrees of freedom of the system:

$u_{2D}^* = \begin{pmatrix} u^* \\ v^* \\ \end{pmatrix} = u_e^* N^T$

$\varepsilon^{*T} = u_e^* B^T$

where

$u_e^* = \hspace{1mm} < \hspace{2mm} u_1^* \hspace{2mm} v_1^* \hspace{2mm} u_2^* \hspace{2mm} v_2^* \hspace{2mm} \cdots \hspace{2mm} u_n^* \hspace{2mm} v_n^* \hspace{2mm} > \hspace{1mm}$

$\varepsilon^{*T} = \hspace{1mm} < \hspace{2mm} \varepsilon_{11}^* \hspace{2mm} \varepsilon_{22}^* \hspace{2mm} 2\varepsilon_{12}^* \hspace{2mm} > \hspace{1mm}$

The internal virtual work integral term for this particular element will have the form:

$IVW_e = \int_e \sum_{i,j=1}^3 \varepsilon_{ij}^* \sigma_{ij} dx = \int_e \varepsilon^* \cdot \sigma dx$

$= u_e^*\Big(\int_e B^T CB dx \Big) u_e$

where the integration is done on each term in the matrix $B^TCB$ and is performed over the volume (which in the 2D case would be an integration over the area multiplied by the thickness) of the element. The local stiffness matrix has dimensions $2n\times 2n$ and has the form:

$K^e = \int_e B^T CB dx$

The external virtual work integral term for this particular element will have the form:

$EVW_e = u_e^* \int_{\partial e} N^T t_n ds + u_e^* \int_e N^T \rho b dx$

Where $t_n\in\mathbb{R}^2$ is the traction vectors on the surface of the element, $\rho$ is the mass density, and $b\in\mathbb{R}^2$ is the body forces vector on the element. The integration in the first term is done over the element surface while in the second term is done over the element volume. The local nodal forces vector will $2n$ components and will have the form: