Finite Element: Q8 Quadrilaterals

Last updated August 30, 2023
By Ian Story

The following are notes that may be helpful when implementing the stiffness matrix for a Q8 quadrilateral. These notes supplement the following YouTube lectures by Dr. Clayton Pettit:

Finite Element Method: Quadrilateral Elements

Finite Element Method: Isoparametric Elements

Shape Functions

Throughout this derivation, we will use the following shape functions:

$N_1(\xi, \eta) = \frac{-(1-\eta)(1-\xi)(1+\xi+\eta)}{4}$

$N_2(\xi, \eta) = \frac{-(1-\eta)(1+\xi)(1-\xi+\eta)}{4}$

$N_3(\xi, \eta) = \frac{-(1+\eta)(1+\xi)(1-\xi-\eta)}{4}$

$N_4(\xi, \eta) = \frac{-(1+\eta)(1-\xi)(1+\xi-\eta)}{4}$

$N_5(\xi, \eta) = \frac{(1-\eta)(1-\xi)(1+\xi)}{2}$

$N_6(\xi, \eta) = \frac{(1-\eta)(1+\xi)(1+\eta)}{2}$

$N_7(\xi, \eta) = \frac{(1+\eta)(1-\xi)(1+\xi)}{2}$

$N_8(\xi, \eta) = \frac{(1-\eta)(1-\xi)(1+\eta)}{2}$

To simplify display in matrices, the function arguments $(\xi, \eta)$ are omitted in the following derivation. Just remember that everywhere $N_i$ appears, it actually means $N_i(\xi, \eta)$ , a function that must be evaluated at every point in the natural coordinate system.

Coordinate Mapping

$\setcounter{MaxMatrixCols}{20}x_{global}(\xi, \eta) =\begin{bmatrix}N_1 & N_2 & N_3 & N_4 & N_5 & N_6 & N_7 & N_8\end{bmatrix}\begin{bmatrix}x_1 \\x_2 \\x_3 \\x_4 \\x_5 \\x_6 \\x_7 \\x_8\end{bmatrix}$

$\setcounter{MaxMatrixCols}{20}y_{global}(\xi, \eta) =\begin{bmatrix}N_1 & N_2 & N_3 & N_4 & N_5 & N_6 & N_7 & N_8\end{bmatrix}\begin{bmatrix}y_1 \\y_2 \\y_3 \\y_4 \\y_5 \\y_6 \\y_7 \\y_8\end{bmatrix}$

Jacobian Matrix

$\setcounter{MaxMatrixCols}{20}[J](\xi, \eta) =\renewcommand*{\arraystretch}{2.5}\begin{bmatrix}\cfrac{\partial x}{\partial \xi} & \cfrac{\partial x}{\partial \eta} \\\cfrac{\partial y}{\partial \xi} & \cfrac{\partial y}{\partial \eta}\end{bmatrix}$

Where:

$\cfrac{\partial x}{\partial \xi} =\begin{bmatrix}\cfrac{\partial N_1}{\partial \xi} & \cfrac{\partial N_2}{\partial \xi} & \cfrac{\partial N_3}{\partial \xi} & \cfrac{\partial N_4}{\partial \xi} & \cfrac{\partial N_5}{\partial \xi} & \cfrac{\partial N_6}{\partial \xi} & \cfrac{\partial N_7}{\partial \xi} & \cfrac{\partial N_8}{\partial \xi}\end{bmatrix}\begin{bmatrix}x_1 \\x_2 \\x_3 \\x_4 \\x_5 \\x_6 \\x_7 \\x_8\end{bmatrix}$

The partial derivative for y and partial derivatives with respect to $\eta$ are identical to the above: just swap out the relevant symbols.

Making use of the $[\partial N]$ matrix, presented below, the Jacobian matrix can be expressed as follows:

$\setcounter{MaxMatrixCols}{20}[J](\xi, \eta) = [\partial N](\xi, \eta)\begin{bmatrix}x_1 & y_1 \\x_2 & y_2 \\x_3 & y_3 \\x_4 & y_4 \\x_5 & y_5 \\x_6 & y_6 \\x_7 & y_7 \\x_8 & y_8 \end{bmatrix}$

$\partial$ N Matrix

To simplify further calculations, we can define a $\partial$ N matrix as follows:

$[\partial N](\xi, \eta) =\renewcommand*{\arraystretch}{2.5}\begin{bmatrix}\cfrac{\partial N_1}{\partial \xi} & \cfrac{\partial N_2}{\partial \xi} & \cfrac{\partial N_3}{\partial \xi} & \cfrac{\partial N_4}{\partial \xi} & \cfrac{\partial N_5}{\partial \xi} & \cfrac{\partial N_6}{\partial \xi} & \cfrac{\partial N_7}{\partial \xi} & \cfrac{\partial N_8}{\partial \xi} \\\cfrac{\partial N_1}{\partial \eta} & \cfrac{\partial N_2}{\partial \eta} & \cfrac{\partial N_3}{\partial \eta} & \cfrac{\partial N_4}{\partial \eta} & \cfrac{\partial N_5}{\partial \eta} & \cfrac{\partial N_6}{\partial \eta} & \cfrac{\partial N_7}{\partial \eta} & \cfrac{\partial N_8}{\partial \eta}\end{bmatrix}$

Where:

$\cfrac{\partial N_1}{\partial \xi}(\xi, \eta) = \frac{-(\eta-1)(2\xi+\eta)}{4}$

$\cfrac{\partial N_2}{\partial \xi}(\xi, \eta) = \frac{-(\eta-1)(2\xi-\eta)}{4}$

$\cfrac{\partial N_3}{\partial \xi}(\xi, \eta) = \frac{(\eta+1)(2\xi+\eta)}{4}$

$\cfrac{\partial N_4}{\partial \xi}(\xi, \eta) = \frac{(\eta+1)(2\xi-\eta)}{4}$

$\cfrac{\partial N_5}{\partial \xi}(\xi, \eta) = \xi(\eta-1)$

$\cfrac{\partial N_6}{\partial \xi}(\xi, \eta) = \frac{-(\eta-1)(\eta+1)}{2}$

$\cfrac{\partial N_7}{\partial \xi}(\xi, \eta) = -\xi(\eta+1)$

$\cfrac{\partial N_8}{\partial \xi}(\xi, \eta) = \frac{(\eta-1)(\eta+1)}{2}$

$\cfrac{\partial N_1}{\partial \eta}(\xi, \eta) = \frac{-(\xi-1)(\xi+2\eta)}{4}$

$\cfrac{\partial N_2}{\partial \eta}(\xi, \eta) = \frac{-(\xi+1)(\xi-2\eta)}{4}$

$\cfrac{\partial N_3}{\partial \eta}(\xi, \eta) = \frac{(\xi+1)(\xi+2\eta)}{4}$

$\cfrac{\partial N_4}{\partial \eta}(\xi, \eta) = \frac{(\xi-1)(\xi-2\eta)}{4}$

$\cfrac{\partial N_5}{\partial \eta}(\xi, \eta) = \frac{(\xi-1)(\xi+1)}{2}$

$\cfrac{\partial N_6}{\partial \eta}(\xi, \eta) = -\eta(\xi+1)$

$\cfrac{\partial N_7}{\partial \eta}(\xi, \eta) = \frac{-(\xi-1)(\xi+1)}{2}$

$\cfrac{\partial N_8}{\partial \eta}(\xi, \eta) = \eta(\xi-1)$

B Matrix Derivation

$[\epsilon_{2D}](x,y) =\begin{bmatrix}\epsilon_{xx} \\\epsilon_{yy} \\\gamma_{xy}\end{bmatrix} =\begin{bmatrix}\epsilon_{xx} \\\epsilon_{yy} \\\epsilon_{xy} + \epsilon_{yx} \\\end{bmatrix} =\renewcommand*{\arraystretch}{2.5}\begin{bmatrix}\cfrac{du}{dx} \\\cfrac{dv}{dy} \\\cfrac{du}{dy} + \cfrac{dv}{dx} \\\end{bmatrix} =\renewcommand*{\arraystretch}{1.0}\begin{bmatrix}1 & 0 & 0 & 0 \\0 & 0 & 0 & 1 \\0 & 1 & 1 & 0\end{bmatrix}\renewcommand*{\arraystretch}{2.5}\begin{bmatrix}\cfrac{du}{dx} \\\cfrac{du}{dy} \\\cfrac{dv}{dx} \\\cfrac{dv}{dy}\end{bmatrix}$

Where (by multivariable chain rule):

$\cfrac{\partial u}{\partial x} = \cfrac{\partial u}{\partial \xi}\cfrac{\partial \xi}{\partial x} + \cfrac{\partial u}{\partial \eta}\cfrac{\partial \eta}{\partial x}$

$\cfrac{\partial u}{\partial y} = \cfrac{\partial u}{\partial \xi}\cfrac{\partial \xi}{\partial y} + \cfrac{\partial u}{\partial \eta}\cfrac{\partial \eta}{\partial y}$

$\cfrac{\partial v}{\partial x} = \cfrac{\partial v}{\partial \xi}\cfrac{\partial \xi}{\partial x} + \cfrac{\partial v}{\partial \eta}\cfrac{\partial \eta}{\partial x}$

$\cfrac{\partial v}{\partial y} = \cfrac{\partial v}{\partial \xi}\cfrac{\partial \xi}{\partial y} + \cfrac{\partial v}{\partial \eta}\cfrac{\partial \eta}{\partial y}$

Or, in matrix form:

$\renewcommand*{\arraystretch}{1.5}\begin{bmatrix}\cfrac{du}{dx} \\\cfrac{du}{dy} \\\cfrac{dv}{dx} \\\cfrac{dv}{dy}\end{bmatrix} =\renewcommand*{\arraystretch}{2.5}\begin{bmatrix}\cfrac{\partial \xi}{\partial x} & \cfrac{\partial \eta}{\partial x} & 0 & 0 \\\cfrac{\partial \xi}{\partial y} & \cfrac{\partial \eta}{\partial y} & 0 & 0 \\0 & 0 & \cfrac{\partial \xi}{\partial x} & \cfrac{\partial \eta}{\partial x} \\0 & 0 & \cfrac{\partial \xi}{\partial y} & \cfrac{\partial \eta}{\partial y}\end{bmatrix}\renewcommand*{\arraystretch}{2.5}\begin{bmatrix}\cfrac{du}{d\xi} \\\cfrac{du}{d\eta} \\\cfrac{dv}{d\xi} \\\cfrac{dv}{d\eta}\end{bmatrix}$

Conveniently, the first matrix can be simplified based on existing information, because:

$\renewcommand*{\arraystretch}{2.5}\begin{bmatrix}\cfrac{\partial \xi}{\partial x} & \cfrac{\partial \eta}{\partial x} \\\cfrac{\partial \xi}{\partial y} & \cfrac{\partial \eta}{\partial y}\end{bmatrix} =([J](\xi, \eta)^{-1})^T$

Because $u$ and $v$ are just the nodal displacements at $x$ and $y$ , the same coordinate mapping applies here:

$\setcounter{MaxMatrixCols}{20}(u\text{ or }v)_{global}(\xi, \eta) =\begin{bmatrix}N_1 & N_2 & N_3 & N_4 & N_5 & N_6 & N_7 & N_8\end{bmatrix}\begin{bmatrix}(u\text{ or }v)_1 \\(u\text{ or }v)_2 \\(u\text{ or }v)_3 \\(u\text{ or }v)_4 \\(u\text{ or }v)_5 \\(u\text{ or }v)_6 \\(u\text{ or }v)_7 \\(u\text{ or }v)_8\end{bmatrix}$

Therefore:

$\setcounter{MaxMatrixCols}{20}\cfrac{d(u\text{ or }v)}{d(\xi\text{ or }\eta)}(\xi, \eta) =\renewcommand*{\arraystretch}{2.5}\begin{bmatrix}\cfrac{\partial N_1}{\partial (\xi\text{ or }\eta)} & \cfrac{\partial N_2}{\partial (\xi\text{ or }\eta)} & \cfrac{\partial N_3}{\partial (\xi\text{ or }\eta)} & \cfrac{\partial N_4}{\partial (\xi\text{ or }\eta)} & \cfrac{\partial N_5}{\partial (\xi\text{ or }\eta)} & \cfrac{\partial N_6}{\partial (\xi\text{ or }\eta)} & \cfrac{\partial N_7}{\partial (\xi\text{ or }\eta)} & \cfrac{\partial N_8}{\partial (\xi\text{ or }\eta)}\end{bmatrix}\renewcommand*{\arraystretch}{1.0}\begin{bmatrix}(u\text{ or }v)_1 \\(u\text{ or }v)_2 \\(u\text{ or }v)_3 \\(u\text{ or }v)_4 \\(u\text{ or }v)_5 \\(u\text{ or }v)_6 \\(u\text{ or }v)_7 \\(u\text{ or }v)_8\end{bmatrix}$

Replacing the nodal displacements with a vector, $[u] =\begin{bmatrix}u_1 \\v_1 \\\vdots \\u_8 \\v_8 \\\end{bmatrix}$ ,

$\renewcommand*{\arraystretch}{2.5}\begin{bmatrix}\cfrac{du}{d\xi} \\\cfrac{du}{d\eta} \\\cfrac{dv}{d\xi} \\\cfrac{dv}{d\eta}\end{bmatrix} = \renewcommand*{\arraystretch}{2.5}\setcounter{MaxMatrixCols}{20}\begin{bmatrix}\cfrac{\partial N_1}{\partial \xi} & 0 & \cfrac{\partial N_2}{\partial \xi} & 0 & \cfrac{\partial N_3}{\partial \xi} & 0 & \cfrac{\partial N_4}{\partial \xi} & 0 & \cfrac{\partial N_5}{\partial \xi} & 0 & \cfrac{\partial N_6}{\partial \xi} & 0 & \cfrac{\partial N_7}{\partial \xi} & 0 & \cfrac{\partial N_8}{\partial \xi} & 0 \\\cfrac{\partial N_1}{\partial \eta} & 0 & \cfrac{\partial N_2}{\partial \eta} & 0 & \cfrac{\partial N_3}{\partial \eta} & 0 & \cfrac{\partial N_4}{\partial \eta} & 0 & \cfrac{\partial N_5}{\partial \eta} & 0 & \cfrac{\partial N_6}{\partial \eta} & 0 & \cfrac{\partial N_7}{\partial \eta} & 0 & \cfrac{\partial N_8}{\partial \eta} & 0 \\0 & \cfrac{\partial N_1}{\partial \xi} & 0 & \cfrac{\partial N_2}{\partial \xi} & 0 & \cfrac{\partial N_3}{\partial \xi} & 0 & \cfrac{\partial N_4}{\partial \xi} & 0 & \cfrac{\partial N_5}{\partial \xi} & 0 & \cfrac{\partial N_6}{\partial \xi} & 0 & \cfrac{\partial N_7}{\partial \xi} & 0 & \cfrac{\partial N_8}{\partial \xi} \\0 & \cfrac{\partial N_1}{\partial \eta} & 0 & \cfrac{\partial N_2}{\partial \eta} & 0 & \cfrac{\partial N_3}{\partial \eta} & 0 & \cfrac{\partial N_4}{\partial \eta} & 0 & \cfrac{\partial N_5}{\partial \eta} & 0 & \cfrac{\partial N_6}{\partial \eta} & 0 & \cfrac{\partial N_7}{\partial \eta} & 0 & \cfrac{\partial N_8}{\partial \eta}\end{bmatrix}[u]$

This 4×16 matrix is just the $[\partial N]$ matrix with zero columns inserted between each column, repeated and offset by one column. Therefore, let us define:

$[\partial N2] \equiv\renewcommand*{\arraystretch}{2.5}\setcounter{MaxMatrixCols}{20}\begin{bmatrix}\cfrac{\partial N_1}{\partial \xi} & 0 & \cfrac{\partial N_2}{\partial \xi} & 0 & \cfrac{\partial N_3}{\partial \xi} & 0 & \cfrac{\partial N_4}{\partial \xi} & 0 & \cfrac{\partial N_5}{\partial \xi} & 0 & \cfrac{\partial N_6}{\partial \xi} & 0 & \cfrac{\partial N_7}{\partial \xi} & 0 & \cfrac{\partial N_8}{\partial \xi} & 0 \\\cfrac{\partial N_1}{\partial \eta} & 0 & \cfrac{\partial N_2}{\partial \eta} & 0 & \cfrac{\partial N_3}{\partial \eta} & 0 & \cfrac{\partial N_4}{\partial \eta} & 0 & \cfrac{\partial N_5}{\partial \eta} & 0 & \cfrac{\partial N_6}{\partial \eta} & 0 & \cfrac{\partial N_7}{\partial \eta} & 0 & \cfrac{\partial N_8}{\partial \eta} & 0 \\0 & \cfrac{\partial N_1}{\partial \xi} & 0 & \cfrac{\partial N_2}{\partial \xi} & 0 & \cfrac{\partial N_3}{\partial \xi} & 0 & \cfrac{\partial N_4}{\partial \xi} & 0 & \cfrac{\partial N_5}{\partial \xi} & 0 & \cfrac{\partial N_6}{\partial \xi} & 0 & \cfrac{\partial N_7}{\partial \xi} & 0 & \cfrac{\partial N_8}{\partial \xi} \\0 & \cfrac{\partial N_1}{\partial \eta} & 0 & \cfrac{\partial N_2}{\partial \eta} & 0 & \cfrac{\partial N_3}{\partial \eta} & 0 & \cfrac{\partial N_4}{\partial \eta} & 0 & \cfrac{\partial N_5}{\partial \eta} & 0 & \cfrac{\partial N_6}{\partial \eta} & 0 & \cfrac{\partial N_7}{\partial \eta} & 0 & \cfrac{\partial N_8}{\partial \eta}\end{bmatrix}$

Putting all the pieces together gives

$[\epsilon_{2D}](x,y) =\begin{bmatrix}1 & 0 & 0 & 0 \\0 & 0 & 0 & 1 \\0 & 1 & 1 & 0\end{bmatrix}\begin{bmatrix}([J](\xi, \eta)^{-1})^T & \begin{bmatrix}0 & 0 \\0 & 0\end{bmatrix} \\\begin{bmatrix}0 & 0 \\0 & 0\end{bmatrix} & ([J](\xi, \eta)^{-1})^T\end{bmatrix}[\partial N2]$

By the definition of the kinematic matrix ( $[B]$ )

$[\epsilon_{2D}](x,y) = [B](\xi,\eta)[u]$

Which implies

$[B](\xi,\eta) =\begin{bmatrix}1 & 0 & 0 & 0 \\0 & 0 & 0 & 1 \\0 & 1 & 1 & 0\end{bmatrix}\begin{bmatrix}([J](\xi, \eta)^{-1})^T & \begin{bmatrix}0 & 0 \\0 & 0\end{bmatrix} \\\begin{bmatrix}0 & 0 \\0 & 0\end{bmatrix} & ([J](\xi, \eta)^{-1})^T\end{bmatrix}[\partial N2]$

C Matrix

Under plane stress,

$[\sigma_{2D}](x,y) =\begin{bmatrix}\sigma_{xx} \\\sigma_{yy} \\\tau_{xy}\end{bmatrix} =\cfrac{E}{1-\nu^2}\begin{bmatrix}1 & \nu & 0 \\\nu & 1 & 0 \\0 & 0 & 1-\nu\end{bmatrix}\begin{bmatrix}\epsilon_{xx} \\\epsilon_{yy} \\\epsilon_{xy}\end{bmatrix}$

By the definition of the constitutive matrix ( $[C]$ ):

$[\sigma_{2D}](x,y) = [C][\epsilon_{2D}](x,y)$

Therefore,

$[C] = \cfrac{E}{1-\nu^2}\begin{bmatrix}1 & \nu & 0 \\\nu & 1 & 0 \\0 & 0 & 1-\nu\end{bmatrix}$

Element Stiffness Matrix

$[k_e] = t\int_{-1}^{1}\int_{-1}^{1}[B](\xi,\eta)^T[C][B](\xi,\eta)\text{det}\bigl([J](\xi,\eta)\bigr) d\xi d\eta$

This integral can be simplified by using Gauss Integration (recommended n=3 for balance between performance and accuracy – this requires sampling the element at 9 points)

Tutorial: Dr. Clayton Pettit on Gaussian Integration

Post Processing

Once nodal deflections $[u]$ are known, the strains and stresses can be determined at any point on the element:

$[\epsilon_{2D}](x,y) =\begin{bmatrix}\epsilon_{xx} \\\epsilon_{yy} \\\gamma_{xy}\end{bmatrix} = [B](\xi,\eta)[u]$

$[\sigma_{2D}](x,y) = \begin{bmatrix}\sigma_{xx} \\\sigma_{yy} \\\tau_{xy}\end{bmatrix} =[C][\epsilon_{2D}](x,y)$

$\begin{bmatrix}\sigma_{xx} \\\sigma_{yy} \\\tau_{xy}\end{bmatrix} = [C][B](\xi,\eta)[u]$

Note that it is not easy to convert from the local/global coordinate systems to the natural coordinate system. Because all of the shape functions depend on inputs of $\xi$ and $\eta$ , working backwards requires a trial and error process until the $\xi, \eta$ guesses match the desired x, y outputs.

Therefore, it is more typical to calculate stresses in the natural coordinate system at regular intervals and place these at their corresponding places in the local coordinate system.

More Helpful Reference

https://engcourses-uofa.ca/books/introduction-to-solid-mechanics/finite-element-analysis/two-dimensional-solid-elements/quadrilateral-elements/