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 the first article in the theory series was called “Part 1 – Principle of Virtual Work“.
As I derived in the previous article, the principle of virtual work can be written: \[\displaystyle\sum_m \int_{V^{(m)}} \delta\boldsymbol{\varepsilon}^\top \boldsymbol{\sigma} dV = \sum_m\int_{V^{(m)}} \delta\mathbf{u}^\top \mathbf{b} dV + \sum_m \int_{A^{(m)}} \delta \mathbf{u}^\top \mathbf{t} dA + \sum \delta \mathbf{u}^\top \mathbf{f}.\]
The key step in solving this equation is to introduce isoparametric shape functions – which is the topic of this article. Note that the theory presented here is for 2D linear plane strain elements.
Isoparametric Shape Functions
The internal virtual work term in the Principle of Virtual Work has the form shown below. The challenge here is to calculate this integral since the element volume is not simple to integrate over.
One convenient way to fix this problem is to introduce a coordinate transformation (mapping). If we can create a transformation that maps the element shape into a square located between -1 and +1 then the integration would be much easier to perform. The shape function coordinates, which are also called Natural Coordinates, are \((\xi,\eta)\).
For the case of linear plane strain quadrilateral elements, one way to perform the mapping from the real coordinates to the shape function (natural) coordinates is to introduce the following shape functions:
With this definition the shape functions get the following values at the corner nodes.
In Python, the shape functions can be calculated from:
def shape(xi):
x,y = tuple(xi)
N = [(1.0-x)*(1.0-y), (1.0+x)*(1.0-y), (1.0+x)*(1.0+y), (1.0-x)*(1.0+y)]
return 0.25*np.array(N)
...
If we know the nodal coordinates ((\(x_1,y_1\)), (\(x_2,y_2\)), (\(x_3,y_3\)), (\(x_4,y_4\))), and a target location in the natural coordinates is (\(\xi,\eta\)), then the corresponding real location is:
Similarly, if the displacements of the nodes are known ((\(u_1,v_1\)), (\(u_2,v_2\)), (\(u_3,v_3\)), (\(u_4,v_4\))), and a target location in the natural coordinates is (\(\xi,\eta\)), then the displacement of the corresponding real location is:
\[u=N_1(\xi,\eta)u_1 + N_2(\xi,\eta)u_2 + N_3(\xi,\eta)u_3 + N_4(\xi,\eta)u_4\] \[v=N_1(\xi,\eta)v_1 + N_2(\xi,\eta)v_2 + N_3(\xi,\eta)v_3 + N_4(\xi,\eta)v_4\]
If you think about these equations some more it is easy to see that it is just a convenient mathematical mapping. The displacement field of a real physical specimen that is forced to undergo the same “nodal” displacement will not necessarily be the same as what is given by these equation. This is why it is important to use a refined FE mesh, and also why FE solvers support both linear and quadratic elements.
For numerical implementations it is convenient to use matrix notation. Define \(\mathbf{U}\) to be the displacements of the corner nodes: \[\mathbf{U} = [u_1, v_1, u_2, v_2, u_3, v_3, u_4, v_4]^\top ,\] and define \(\mathbf{H}\) from:
\[\mathbf{H}(\xi,\eta) = \begin{bmatrix} N_1 & 0 & N_2 & 0 & N_3 & 0 & N_4 & 0\\
0 & N_1 & 0 & N_2 & 0 & N_3 & 0 & N_4 \end{bmatrix}.\] The displacement of any point inside the element is then given by: \[\mathbf{u} = \begin{bmatrix}u\\v\end{bmatrix} = \mathbf{H} \mathbf{U}\]
In the next part in this series I will show how to calculate the strains from the nodal displacements.
