Locking and Hourglassing: Where the Privilege Ends

1. Why a separate engine?

The central result of FEM-01 was this: the cubic shape functions of the beam element are the exact solution of the Euler–Bernoulli equation when there is no distributed load inside the element. Hence a cantilever loaded at its tip is solved exactly with one element.

In a continuum — plane stress, plane strain, solid — that privilege does not exist. The shape functions cannot reproduce the true displacement field; they only approximate it. Mesh refinement is therefore genuinely necessary here.

But that is not the whole story. How bad the approximation can be is worse than most people expect.

2. How is an element's validity tested?

Whether a new element formulation is correct is tested not by comparison with commercial output but against the element's own mathematical properties. There are two basic criteria:

• Rigid body motion. Translating and rotating the element must produce no strain. If it does, the formulation is wrong and that element generates spurious internal forces.
• The patch test. A constant strain field must be reproducible exactly on an arbitrary mesh. This is the guarantee that the element will converge.

Our engine passes both to machine precision; the results are pinned in the test suite.

3. The patch test

The setup is as follows: all boundary nodes of a mesh are assigned displacements corresponding to a constant uniaxial stress state. Interior nodes are left free. If the element is correct, the interior stress comes out exactly constant.

4. Locking: why is the CST so stiff?

The constant strain triangle (CST), as its name says, takes the strain inside the element to be constant. The displacement field is linear, its derivative constant.

Now think about bending. In a bending beam the strain varies linearly through the depth: compression at the top, tension at the bottom, zero in the middle. An element with constant strain cannot represent that by definition.

5. Numerical comparison

Cantilever in plane stress. L = 1000, H = 100, t = 10 mm; tip load 1000 N. The reference, including shear deformation, is the Timoshenko solution: —.

Compare with the beam element in FEM-01: there one element gave the exact answer. Here four CST elements are 90% wrong. The difference between the two cases lies in whether the shape functions contain the solution space.

6. Convergence is always from below

Every value in the table is below the reference. This is not a coincidence but a variational result.

The practical value: displacement computed with compatible elements is a lower bound. If your mesh is coarse, the real deflection is larger than you computed — which means being on the unsafe side of a deflection limit check. The same simple relation does not hold for stress; there convergence may not be one-sided.

7. The cure: reduced integration

The source of locking is that the spurious shear energy is over-represented in the numerical integration. In a Q4 element the stiffness matrix is normally integrated with 2×2 Gauss points. Integrating the shear term at a single point (the centre) leaves most of that spurious energy out of account.

The third column of the table is the result: on the same mesh, reduced integration is markedly more accurate than full integration. On the 64×16 mesh the error is —, against — with full integration.

8. The side effect: the hourglass mode

Now look at the third cell in the first row of the table.

The reason is this: the stiffness matrix is built from the strains at the integration points. In a Q4 element integrated at a single point there exist displacement fields that produce no strain at all at the centre while being non-zero. Such fields cost the element no energy — they are spurious zero-energy modes. Every reduced-integration Q4 element has two of them.

On a sufficiently constrained mesh neighbouring elements block each other and these modes do not wake. On a coarse, single-row mesh there is nobody to block them; the system matrix becomes nearly singular and the solution blows up.

9. Measuring the convergence order

To know whether a model is converging properly you do not need the exact solution. It is enough to watch how the error shrinks across successive meshes — the same tool used in FEM-01.

An observed order that settles to a stable value is the sign of clean convergence. If it does not settle — if the order jumps from mesh to mesh — something is wrong with the model: a singularity, a distorted element, or a badly defined boundary condition.

10. What to do in practice

1. Do not use CST where bending dominates. Triangular elements are tempting because meshing is easy; in bending they are unacceptably stiff. At the very least use a second-order (six-node) triangle.
2. Do not use a single element through the thickness. The worst row in the table is 4×1. You need enough elements to represent the bending gradient.
3. If you use reduced integration, check for hourglassing. Look at the deformed shape: if neighbouring elements zigzag in opposite directions, the mode has woken.
4. Remember that displacement is a lower bound. On a coarse mesh the real deflection is larger than you computed.
5. Measure convergence, do not assume it. Solve at least three meshes and watch how the error shrinks.
6. Run the patch test when you move to new software. An element that cannot reproduce a constant stress field exactly is not trustworthy on any mesh.

11. Test yourself

1. Why is the CST element so stiff in bending? Where in the answer does "constant strain" appear?
2. From which principle does it follow that a displacement-based model is always stiffer than the real structure?
3. Why does reduced integration reduce locking? Which energy term is under-computed?
4. Why does the hourglass mode wake on a coarse mesh and not on a fine one?
5. You doubled a mesh and the error fell to a quarter. What is the observed order?
6. Does an element that passes the patch test give correct results on every problem?