Question 4 of 4: Four-node isoparametric element — strains at the centre and the zero-energy mode

Problem 4: Four-node isoparametric element — strains at the centre and the zero-energy mode (equal value)

Question text not reproduced: the examination questions are © Engineers and Geoscientists BC. Open the official past paper (linked at the top of this page) to read the question, then follow the worked solution below.

Given. A four-node isoparametric quadrilateral that was a square of side two units before loading, with the four nodal displacement pairs listed in Table 2 and the bilinear shape functions printed with the question.

Find. The three strain components at the element centre, $\xi=\eta=0$, and an interpretation of what those values mean for the element.

Approach. Interpolate the displacement field with the printed shape functions, differentiate it through the Jacobian to get the strains as functions of $\xi$ and $\eta$, then evaluate at the centre and interpret the result in terms of numerical integration.

1. Establish the Jacobian. The parent square runs from $-1$ to $+1$ in each natural coordinate while the physical element runs over two units, so with the half-side $a = 1$ the mapping is simply $x = x_0+a\xi = x_0+\xi$ and $y = y_0+\eta$. The Jacobian is therefore the identity, $$[J]=\begin{bmatrix}\partial x/\partial\xi & \partial y/\partial\xi\\ \partial x/\partial\eta & \partial y/\partial\eta\end{bmatrix}=\begin{bmatrix}1 & 0\\ 0 & 1\end{bmatrix},$$ so $\partial/\partial x = \partial/\partial\xi$ and $\partial/\partial y = \partial/\partial\eta$ throughout, and no inverse mapping is needed.
2. Interpolate the displacement field. With $u=\sum N_i u_i$ and only $u_2=-c$ and $u_3=+c$ non-zero, $$u(\xi,\eta)=c\left(N_3-N_2\right)=\frac{c}{4}\left[(1+\xi)(1+\eta)-(1+\xi)(1-\eta)\right]=\frac{c}{2}(1+\xi)\eta,$$ and with only $v_3=-c$ and $v_4=+c$ non-zero, $$v(\xi,\eta)=c\left(N_4-N_3\right)=\frac{c}{4}\left[(1-\xi)(1+\eta)-(1+\xi)(1+\eta)\right]=-\frac{c}{2}\xi(1+\eta).$$ Substituting the four nodal stations reproduces Table 2 exactly, which is the check that the interpolation has been set up correctly.
3. Differentiate to obtain the strain field. Using the identity Jacobian from Step 1, $$\varepsilon_x=\frac{\partial u}{\partial x}=\frac{\partial u}{\partial\xi}=\frac{c\,\eta}{2},\ \ \varepsilon_y=\frac{\partial v}{\partial y}=\frac{\partial v}{\partial\eta}=-\frac{c\,\xi}{2},$$ $$\gamma_{xy}=\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}=\frac{c(1+\xi)}{2}-\frac{c(1+\eta)}{2}=\frac{c}{2}\left(\xi-\eta\right).$$ Every component is linear in the natural coordinates, which is the familiar result that the bilinear quadrilateral carries linearly varying strain even though it uses a bilinear displacement field.
4. Evaluate at the centre. Setting $\xi=\eta=0$ in Step 3, $$\boxed{\varepsilon_x=0,\ \ \varepsilon_y=0,\ \ \gamma_{xy}=0\ \ \text{at the element centre}}$$ All three strains vanish simultaneously, and with them the stresses $\sigma_x=\sigma_y=\tau_{xy}=0$ at that point, for any value of $c$.
5. Confirm that the element really is deformed. The result is not a rigid-body motion in disguise. A rigid-body field would have $u$ and $v$ at most linear in position with a single antisymmetric rotation term, whereas Step 2 gives products such as $\xi\eta$; and the strain field of Step 3 is non-zero everywhere except along the lines $\eta=0$, $\xi=0$ and $\xi=\eta$. At the four Gauss stations of a two-by-two rule, $\xi,\eta=\pm1/\sqrt3$, the strains are $$\varepsilon_x=\pm\frac{c}{2\sqrt3}=\pm0.2887c,\ \ \varepsilon_y=\mp0.2887c,$$ so a two-by-two integration sees real straining while a single central point sees none.
6. Interpret the result (4.2). The displacement pattern of Table 2 is a combination of the two hourglass or zero-energy modes of the four-node quadrilateral — deformation patterns whose strain field is proportional to $\xi$ or $\eta$ and therefore vanishes identically at the centroid. If the element stiffness matrix is formed by one-point (reduced) Gauss integration, $[k]=\left|J\right|\,w\,[B]^{T}[D][B]$ evaluated at $\xi=\eta=0$, then this pattern produces no strain at the sampling point, no strain energy, and hence no restoring force: the stiffness matrix is rank-deficient and admits a spurious mechanism in addition to the three genuine rigid-body modes. In a mesh these modes chain together into the characteristic zig-zag hourglass pattern and the computed displacements can be wildly too large.
7. State the remedies. Three are standard. Integrate the element fully with a two-by-two Gauss rule, which samples the strains of Step 5 and restores the correct rank at the cost of some shear locking in bending-dominated problems. Keep one-point integration for its speed but add an hourglass-control stiffness or viscosity that penalises exactly this deformation pattern, which is what explicit codes do. Or replace the element with a formulation that does not suffer the defect, such as the eight-node quadrilateral or an assumed-strain or incompatible-modes quadrilateral.