Introduction
This article is part of my series on basic FEA theory. The series started with two Python code examples (Full Finite Element Solver in 100 Lines of Python and Full Finite Element Solver in 200 Lines of Python), and articles on (1) virtual work; and (2) isoparametric shape functions.
The goal of this article is to determine a matrix equation of the form: \( \boldsymbol{\varepsilon} = \mathbf{B}^{(m)} \mathbf{U}\). The vector \( \mathbf{U}\) contains all nodal displacement components (of an element) and is a 8×1 vector that was introduced in my previous article. The strain vector has shape 3×1, therefore the matrix \(\mathbf{B}^{(m)}\) has to have a 3×8 shape. In this article I will derive the components of this \( \mathbf{B}^{(m)}\) matrix.
Calculate Strains from Nodal Displacements
The strains for a 2D element at small displacements are given by the partial derivatives of the displacement with respect to the original coordinates: \[ \varepsilon_{ij} = \begin{bmatrix}\varepsilon_{11}\\ \varepsilon_{22} \\ \varepsilon_{12} \end{bmatrix} = \begin{bmatrix} \partial u/ \partial x \\ \partial v / \partial y \\ \frac{1}{2}(\partial u / \partial y + \partial v / \partial x) \end{bmatrix}.\]
Since we are using a parametric shape function to map the nodal coordinates (see this article), I will start by looking a the partial derivatives of the displacements with respect to the natural coordinates (and not the real coordinates): \[ \begin{cases} \displaystyle\frac{\partial u}{\partial \xi} = \frac{\partial u}{\partial x} \frac{\partial x}{\partial \xi} + \frac{\partial u}{\partial y} \frac{\partial y}{\partial \xi} \\ \displaystyle\frac{\partial u}{\partial \eta} = \frac{\partial u}{\partial x} \frac{\partial x}{\partial \eta} + \frac{\partial u}{\partial y} \frac{\partial y}{\partial \eta} \end{cases} \]
This equation can be written in matrix form as: \[ \begin{bmatrix} \partial u / \partial \xi \\ \partial u / \partial \eta \end{bmatrix} = \begin{bmatrix} \partial x / \partial \xi & \partial y / \partial \xi \\ \partial x / \partial \eta & \partial y / \partial \eta \end{bmatrix} \begin{bmatrix} \partial u / \partial x \\ \partial u / \partial y\end{bmatrix} \equiv \begin{bmatrix} J_{11} & J_{12} \\ J_{21} & J_{22} \end{bmatrix} \begin{bmatrix} \partial u / \partial x \\ \partial u / \partial y \end{bmatrix}. \] Note that the Jacobian matrix \( \mathbf{J} \) is determined by the known nodal coordinates and the point \( (\xi,\eta) \).
From my previous article we know that the position of any point inside the element that is specified by \((\xi,\eta)\) is given by \[ \begin{cases} x&=N_1(\xi,\eta)x_1 + N_2(\xi,\eta)x_2 + N_3(\xi,\eta)x_3 + N_4(\xi,\eta)x_4\\
y&=N_1(\xi,\eta)y_1 + N_2(\xi,\eta)y_2 + N_3(\xi,\eta)y_3 + N_4(\xi,\eta)y_4\end{cases} \]
If I let \(\mathbf{x}\) be the nodal coordinates \[ \mathbf{x} = \begin{bmatrix} x_1 & y_1 \\ x_2 & y_2 \\ x_3 & y_3 \\ x_4 & y_4\end{bmatrix},\] and if I define the gradient of the shape functions by \[ d\mathbf{N}(\xi,\eta) = \begin{bmatrix} \partial N_1/\partial \xi & \partial N_2/\partial \xi & \partial N_3/\partial \xi & \partial N_4/\partial \xi \\ \partial N_1/\partial \eta & \partial N_2/\partial \eta & \partial N_3/\partial \eta & \partial N_4/\partial \eta \end{bmatrix}, \] then the \(\mathbf{J}\) matrix becomes \[ \mathbf{J} = d\mathbf{N}\, \mathbf{x} \]
Here’s a Python function that calculates the gradient of the shape functions (\(d\mathbf{N}\)):
def gradshape(xi):
x,y = tuple(xi)
dN = [[-(1.0-y), (1.0-y), (1.0+y), -(1.0+y)],
[-(1.0-x), -(1.0+x), (1.0+x), (1.0-x)]]
return 0.25*np.array(dN)
...
Hence, the partial derivatives with respect to the displacements can be written: \[\begin{bmatrix} \partial u/\partial x \\ \partial u/\partial y\end{bmatrix} = \mathbf{J}^{-1} \begin{bmatrix} \partial u/\partial \xi \\ \partial u/\partial \eta \end{bmatrix} = \mathbf{J}^{-1} \, d\mathbf{N} \, \begin{bmatrix}u_1 \\ u_2 \\ u_3 \\ u_4 \end{bmatrix},\] \[\begin{bmatrix} \partial v/\partial x \\ \partial v/\partial y\end{bmatrix} = \mathbf{J}^{-1} \begin{bmatrix} \partial v/\partial \xi \\ \partial v/\partial \eta \end{bmatrix} = \mathbf{J}^{-1} \, d\mathbf{N} \, \begin{bmatrix}v_1 \\ v_2 \\ v_3 \\ v_4 \end{bmatrix}.\]
If I then define \( \mathbf{M} \equiv \mathbf{J}^{-1} d\mathbf{N}\), then the strain can be calculated from: \[ \begin{bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{12} \end{bmatrix} = \begin{bmatrix} M_{11} & 0 & M_{12} & 0 & M_{13} & 0 & M_{14} & 0 \\ 0 & M_{21} & 0 & M_{22} & 0 & M_{23} & 0 & M_{24} \\ M_{11} & M_{21} & M_{12} & M_{22} & M_{13} & M_{23} & M_{14} & M_{24} \end{bmatrix} \begin{bmatrix} u_1 \\ v_1 \\ u_2 \\ v_2 \\ u_3 \\ v_3 \\ u_4 \\ v_4 \end{bmatrix} = \mathbf{B}^{(m)} \,\mathbf{U} \]
This equation show how we can calculate the strain at any point (specified by \((\xi,\eta)\)) as long as we know the initial nodal locations, and the nodal displacements. The next part in this series covers how to calculate the stresses and element stiffness matrix.
1 thought on “Basic FEA Theory – Part 3 – Strain Calculation”
Please correct me, if I’m wrong, but based on the conventional definition of the Jabobian (https://mathworld.wolfram.com/Jacobian.html), the matrix in the third equation should be called the transpose of the Jacobian. Are you using a different definition? Thank you in advance.