Contents

Your beam deflects ten times less than the analytical solution. You refine the mesh. Still wrong. You add more elements. Still wrong. πŸ” At some point you start to wonder if you mistyped Young’s modulus.

You didn’t. The problem is shear locking β€” a pathology baked into bilinear quad elements that no amount of mesh refinement will fix. Welcome to one of FEM’s most reliably humbling gotchas.

The Problem in One Sentence

A bilinear quad element cannot represent pure bending without also generating fake shear strains that don’t exist. Those phantom strains steal energy from bending, making the element act artificially stiff.

Two Worlds: Bending vs Shear

Cross-section, bending stress, shear and moment diagrams

When your beam bends, two things happen in the cross-section:

Mode What it does Strain
Bending Fibres stretch/compress along the beam axis $\varepsilon_{xx} \neq 0$
Shear Cross-sections slide past each other $\gamma_{xy} \neq 0$

For a slender beam, shear deformation is negligible β€” nearly all strain energy goes into bending. A correct element must therefore represent pure bending with zero shear. The bilinear quad cannot. That’s the whole plot.

Why the Linear Quad Fails

Q4 bilinear element in reference coordinates

The Q4 element maps a unit square in reference coordinates $(\xi, \eta) \in [-1, 1]^2$ to the physical element. Each of its four nodes $a = 1,\ldots,4$ sits at a corner $(\xi_a, \eta_a) \in {-1,+1}^2$. The shape function for node $a$ is the product of two linear ramp functions β€” one in $\xi$, one in $\eta$ β€” chosen so that $N_a = 1$ at node $a$ and $N_a = 0$ at the other three [1, 2]:

\[N_a(\xi, \eta) = \frac{1}{4}(1 + \xi_a \xi)(1 + \eta_a \eta)\]

Here $\xi_a$ and $\eta_a$ are the reference coordinates of node $a$ (each $\pm 1$), so the product $(1+\xi_a\xi)$ is zero at the opposite edge and two at node $a$; the $\tfrac{1}{4}$ normalises the product to unity at the node.

Any displacement inside the element is then an interpolation of the four nodal displacements $(u_a, v_a)$:

\[u = \sum_{a=1}^{4} N_a\, u_a, \qquad v = \sum_{a=1}^{4} N_a\, v_a\]

Expanding those sums using the shape functions above β€” collecting the constant, $\xi$-, $\eta$-, and $\xi\eta$-terms β€” gives the polynomial form of the displacement field (mapped back to physical coordinates $x, y$):

\(u(x,y) = a_0 + a_1 x + a_2 y + a_3 xy\) \(v(x,y) = b_0 + b_1 x + b_2 y + b_3 xy\)

The coefficients $a_0,\ldots,a_3$ and $b_0,\ldots,b_3$ are linear combinations of the eight nodal displacements. The key point is the highest polynomial term: $xy$ β€” bilinear, not quadratic. There is no $x^2$ or $y^2$ term. (This polynomial form is exact for a rectangular element aligned with the axes; on a distorted quad the same conclusion holds via the isoparametric map, with extra algebra that doesn’t change the punchline.)

The shear strain. In small-strain continuum mechanics, engineering shear strain is the sum of the two off-diagonal displacement gradients:

\[\gamma_{xy} = \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}\]

Differentiating the bilinear fields above:

\[\frac{\partial u}{\partial y} = a_2 + a_3 x, \qquad \frac{\partial v}{\partial x} = b_1 + b_3 y\]

so

\[\gamma_{xy} = (a_2 + b_1) + a_3 x + b_3 y\]

The shear strain varies linearly across the element β€” it is not constant, and it cannot be made zero everywhere by any choice of nodal displacements unless $a_2 + b_1 = 0$, $a_3 = 0$, and $b_3 = 0$ simultaneously. Whether that is possible depends on the exact displacement field the physics demands.

The pure bending field. Now apply that test. In classical beam theory, pure bending of curvature $\kappa$ (where $\kappa = M/EI$, the ratio of bending moment to flexural rigidity) produces the following displacement field in the cross-section β€” derivable from the Euler–Bernoulli assumption that plane sections remain plane and perpendicular to the neutral axis:

\[u = \kappa\, x\, y, \qquad v = -\frac{\kappa}{2}(x^2 - \nu\, y^2)\]

where $\nu$ is Poisson’s ratio. The horizontal displacement $u$ is linear in both $x$ (distance along the beam) and $y$ (distance from the neutral axis) β€” fibres on the tension side stretch, fibres on the compression side shorten, proportional to both. The vertical displacement $v$ is quadratic in $x$ β€” the deflected shape of the beam β€” minus a Poisson term quadratic in $y$ (the cross-section expands laterally on the compression side and contracts on the tension side).

The exact shear strain in this field is:

\[\gamma_{xy}^{\text{exact}} = \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} = \kappa x + (-\kappa x) = 0\]

The two terms cancel exactly: $\partial u/\partial y = \kappa x$ (the rotation of the cross-section) and $\partial v/\partial x = -\kappa x$ (the slope of the deflected beam) are equal and opposite. This is the mathematical statement that in pure bending there is no shear β€” the cross-sections rotate, but they do not slide.

Why the quad cannot reproduce this. The exact $v$ field contains $x^2$ and $y^2$ β€” but the bilinear quad only has $x$, $y$, and $xy$. When the element fits the four nodal values of $v$, the interpolated $v^h$ is the best bilinear surface through those four points. That surface has $\partial v^h/\partial x = b_1 + b_3 y$ β€” a term that is at most linear in $y$ and constant in $x$. It cannot match $-\kappa x$, which varies linearly in $x$. The shortfall β€” call it the representation error in $\partial v^h/\partial x$ β€” appears directly in the shear strain:

\[\gamma_{xy}^{h} = \frac{\partial u^h}{\partial y} + \frac{\partial v^h}{\partial x} \neq 0 \quad \text{even though} \quad \gamma_{xy}^{\text{exact}} = 0\]

The element has no mechanism to turn it off. Think of it like trying to carry a tray perfectly level while your wrists insist on tilting β€” the geometry of your grip makes a tilt unavoidable no matter how carefully you try. That phantom shear absorbs energy, stiffens the element, and delivers deflections that are wrong in a hilariously consistent way. This is shear locking. πŸ”’

How Bad Is It?

Simply-supported beam with centre point load

Here’s how bad it gets. For a simply-supported beam of length $L$ and height $h$ with a centre point load $P$, the exact midspan deflection is:

\[\delta = \frac{PL^3}{48EI}\]

where $I = h^3/12$ per unit thickness (this is a 2D plane-stress problem; for a 3D beam of breadth $b$, use $I = bh^3/12$). With $P=1$, $E=1$, this reduces to $\delta = L^3/(4h^3)$, so a beam ten times thinner deflects a thousand times more β€” and your locked element won’t get within shouting distance.

With a single-layer mesh of $N$ bilinear quad elements, the ratio $\delta_\text{FEM}/\delta_\text{exact}$ scales roughly as:

\[\frac{\delta_\text{FEM}}{\delta_\text{exact}} \approx \frac{1}{1 + \alpha(h/L)^{-2}}\]

As your beam gets thinner ($h/L \to 0$), the error explodes. At $h/L = 0.1$ you might be off by 10Γ—; at $h/L = 0.03$ the beam barely twitches in the simulation. Adding more elements of the same type does not help β€” locking is a property of the element formulation, not mesh density. Yes, more elements. Still wrong. Refining helps a little β€” each element sees a smaller h/L, pushing the lock-up factor down β€” but the error never hits zero. You’re rearranging deck chairs.

Play with it below β€” but first, see exactly why in a single element.

Element Zoo: Q4 vs Q4-Reduced vs Q8 vs Q12

One element under pure bending. All four types side-by-side β€” drag ΞΊ and watch Q4 light up red with phantom shear while Q8 and Q12 sit there looking smug. Toggle overlays to see integration points and cross-section lines.

Pure bending field (plane stress):

\[u(\xi,\eta) = \kappa\,\xi\,\eta \qquad v(\xi,\eta) = -\frac{\kappa}{2}\!\left(\xi^2 - \nu\,\eta^2\right)\] \[\gamma_{xy}^{\text{exact}} = \frac{\partial u}{\partial \eta} + \frac{\partial v}{\partial \xi} = \kappa\xi - \kappa\xi = 0\]

The exact field has $\gamma_{xy}=0$ everywhere. Q4 with full 2Γ—2 Gauss integration cannot reproduce this β€” the bilinear shape functions produce parasitic shear away from the centroid, and all four integration points see non-zero shear. Q4 with 1Γ—1 integration samples only the centroid, where shear happens to be zero, so it unlocks β€” but opens the door to hourglass modes (a different headache). Q8 and Q12 have enough degrees of freedom to represent the bending field exactly, so their B-matrix gives zero shear throughout.

Q4 β€” full 2Γ—2 integration β€” LOCKS πŸ”’
Q4-Reduced β€” 1Γ—1 Gauss β€” unlocks, hourglass risk ⚠
Q8 β€” serendipity (8 nodes) β€” no locking βœ“
Q12 β€” serendipity (12 nodes) β€” no locking βœ“

Red fill = |Ξ³_xy| from B-matrix at each interior point. Q4 full: phantom shear throughout (locks). Q4 reduced: same strain field exists (dimmed) but integration only samples centroid where Ξ³=0 β€” that's why it unlocks, but misses hourglass modes. Q8/Q12: field is truly zero β€” higher-order DOFs eliminate the parasitic term. Circles = nodes (corner dark, midside amber). Crosses = Gauss points. Dashed = undeformed.

Interactive: Beam + Element View

Three beams stacked vertically: Q4 full 2Γ—2 (locked), Q8 serendipity (no locking), and the Euler-Bernoulli analytical solution. Same beam, same load β€” very different results. Each beam is colored by strain, and the zoom panel exposes the mid-span element where bending moment is highest and phantom shear is at its worst.

Beam geometry
Color map
low β†’ high
Results
Q4 full (locked) Q8 βœ“ Ξ΄_mid (FEM) Ξ΄_mid (exact) Ξ΄/Ξ΄_exact max Ξ³_xy verdict
⚠ Q4 full 2Γ—2 β€” LOCKED
βœ“ Q8 serendipity β€” no locking
β—‡ Analytical (Euler-Bernoulli exact)

Each column: full deformed beam (left) + zoom of one highlighted element (right). Phantom shear Ξ³_xy lights the locked beam red; selective and analytical stay blue.

Try this sequence β€” it tells the whole story in about 30 seconds:

  1. Set L/h = 20, elements = 6, layers = 1 β€” Q4 barely deflects; Q8 and analytical agree
  2. Color map Shear Ξ³_xy β€” Q4 mid-span element lights red with phantom shear; Q8 and analytical stay blue
  3. Push L/h to 30 β€” Q4 ratio collapses further; watch Ξ΄/Ξ΄_exact crater in the stats grid
  4. Color map Bending Ξ΅_xx β€” all three rows show the same bending gradient; only the shear field differs

Error vs Slenderness

Ξ΄_FEM/Ξ΄_exact vs L/h ratio, 6 elements. Full 2Γ—2 collapses for thin beams. Selective reduced + reduced 1Γ—1 stay near the analytical reference (purple). Reduced 1Γ—1 needs β‰₯2 layers to avoid hourglass singularity.

The Three Fixes

Fix 1: Reduced Integration (1Γ—1 Gauss)

Gauss integration schemes: full 2Γ—2, reduced 1Γ—1, selective

Use a single Gauss point at the element centroid. At the centroid, the offending $xy$ term is zero, so the phantom shear vanishes. Simple, elegant, slightly terrifying.

Scheme Points Catches locking? Problem
Full 2Γ—2 4 No β€” locks Accurate otherwise
Reduced 1Γ—1 1 Yes β€” unlocks Introduces hourglass modes

Hourglass modes are zero-energy deformation patterns that a single integration point cannot detect [3]. Adjacent elements can distort in a checkerboard pattern with no stiffness penalty β€” your mesh will look like modern art and your results will be garbage. In practice this requires explicit hourglass stabilisation, which adds complexity and tuning parameters.

Fix 2: Selective Reduced Integration

Integrate different strain components with different quadrature:

  • Normal strains ($\varepsilon_{xx}$, $\varepsilon_{yy}$): full 2Γ—2 Gauss
  • Shear strain ($\gamma_{xy}$): reduced 1Γ—1 Gauss

Underintegrating only the shear term kills the phantom shear without opening the door to hourglass modes. This is the standard in most commercial codes β€” ABAQUS CPE4R uses reduced integration with hourglass control; ANSYS PLANE182 with keyopt(1)=0 uses selective reduced. It’s the β€œjust enough surgery” option.

Fix 3: Higher-Order Elements

Q8 serendipity element in reference coordinates

The 8-node serendipity quad (Q8) has quadratic shape functions with $x^2$, $y^2$, $x^2y$, $xy^2$ terms β€” enough DOFs to represent pure bending exactly. No locking, no hourglass, no special treatment needed. Assembly costs more per element, but you typically need far fewer of them for the same accuracy in bending-dominated problems. When in doubt, reach for Q8.

When Does Locking Actually Matter?

Situation Locking risk What to do
Thick beam / block (L/h < 5) Low Q4 full integration is fine
Slender beam / shell (L/h > 10) High Use Q8 or selective reduced
Plane stress / plane strain Moderate Selective reduced is safe
Nearly-incompressible material (Ξ½ β†’ 0.5) Volumetric locking Different problem β€” use mixed or B-bar elements
3D hex elements, bending-dominated High ABAQUS C3D8R, ANSYS SOLID185 with B-bar

If your problem involves bending and your elements have aspect ratio above 5:1, suspect locking immediately. Always sanity-check at least one result against an analytical solution or a quadratic-element mesh β€” it takes five minutes and could save you days of debugging a phantom stiffness [4].

The same problem shows up in 3D. Your 8-node hex element (C3D8 in ABAQUS, SOLID185 in ANSYS) is the 3D analogue of Q4 β€” bilinear in all three directions β€” and it locks just as badly in bending-dominated problems. The fix is the same: C3D8R uses reduced integration with hourglass control, and it’s the default for structural problems in explicit dynamics for good reason. If you’re running an implicit 3D beam or plate model with C3D8, switch to C3D8R or the 20-node quadratic C3D20 (Q8’s 3D cousin) and watch your deflections correct themselves.

Detecting Locking in Your Own Model

Locking signature: Ξ΄/Ξ΄_exact vs L/h for Q4 and Q8

Before you trust any deflection number, run this three-step sanity check. It takes less time than a coffee break, and it’ll tell you whether you’re looking at physics or a pathology.

  1. Run with Q4 full 2Γ—2 integration. Note the midspan or tip deflection.
  2. Re-run with Q8 (or selective reduced). If deflection jumps by more than ~10%, you have locking. That jump isn’t mesh sensitivity β€” it’s the element formulation releasing the phantom stiffness it was holding onto.
  3. Plot Ξ΄_FEM/Ξ΄_exact versus L/h. If the ratio drops as the beam gets thinner, that’s the locking signature β€” not a coincidence, not a boundary condition error. The element is choking on slenderness.

An additional signal worth checking: strain energy. A locked mesh stores less strain energy than the exact solution β€” the β€œmissing” energy went into phantom shear that the physics doesn’t actually have. If your strain energy is suspiciously low compared to a reference or a coarser Q8 run, locking is the likely culprit.

No analytical reference to compare against? No problem. Run a coarse Q8 mesh and a fine Q4 mesh on the same problem. If the coarse Q8 is more accurate than the fine Q4, Q4 is locked β€” simple as that. The extra refinement bought you nothing because the problem was never mesh density in the first place.

And the one-liner to keep in your back pocket: if your beam gets stiffer as you make it thinner, something is wrong. Physics doesn’t work that way.

The Incompressible Cousin: Volumetric Locking

Volumetric locking: undeformed vs incompressible block

Shear locking has a sibling. For nearly-incompressible materials β€” rubber, biological tissue, metals in plastic flow β€” as $\nu \to 0.5$ the bulk modulus $K \to \infty$ and standard displacement elements lock volumetrically: your mesh cannot represent the divergence-free deformation field the physics demands. Different symptom, same root cause.

The fixes are structurally identical: the B-bar method (selective reduced integration on the volumetric term), mixed $u$-$p$ elements, or enhanced assumed strain (EAS). If you’ve fixed shear locking and your rubber seal still looks wrong, you’ve just met the cousin.

Which Element Should You Actually Use?

Β  Full 2Γ—2 Reduced 1Γ—1 Selective reduced Q8 serendipity
Shear locking Locks Fixed Fixed No locking
Hourglass modes None Possible None None
Accuracy (thick, L/h < 5) Good Good Good Good
Accuracy (thin, L/h > 10) Poor Good Good Excellent
Practical recommendation Avoid for bending Only with explicit hourglass stabilisation Default for Q4 problems Default for bending-dominated problems

In practice:

Q4 full 2Γ—2: safe only for chunky geometries (L/h below about 5) where bending is not the dominant mode. Never use it as your only element type for structural beams or plates without checking against a reference β€” this is the one that started your mesh-refinement spiral.

Q4 selective reduced: a solid default when you’re stuck with linear elements. Commercial codes like ABAQUS and ANSYS use it in their standard quad elements for a reason. Accurate and stable, though it’s still a linear element β€” so it needs more elements than Q8 to reach the same accuracy. Displacement error converges at O(hΒ²) β€” you need to halve element size to quarter the error.

Q4 reduced 1Γ—1: avoid unless your code has robust hourglass stabilisation already built in. The raw 1-point scheme is too fragile for production use. Don’t say you weren’t warned. ⚠️

Q8 serendipity: the best choice for bending-dominated problems. Exact for pure bending with 2Γ—2 Gauss integration, no hourglass risk, and a convergence rate that leaves Q4 far behind. Displacement error converges at O(hΒ³) with full 3Γ—3 Gauss β€” halve element size, error drops by 8Γ—. That’s why fewer elements get you there faster. The extra DOFs per element are almost always worth it when bending drives the response. Use 3Γ—3 Gauss for distorted meshes.

One bilinear term. One phantom shear. A factor of ten in your deflection. Now you know exactly where it comes from β€” and how to make it stop.

References

[1] Hughes, T. J. R. (2000). The Finite Element Method. Dover. Ch. 4 β€” reduced and selective integration, locking. [2] Zienkiewicz, O. C. & Taylor, R. L. (2000). The Finite Element Method Vol. 1 (5th ed.). Butterworth-Heinemann. Ch. 9. [3] Bathe, K. J. (1996). Finite Element Procedures. Prentice Hall. Ch. 4 β€” shear and volumetric locking. [4] MacNeal, R. H. & Harder, R. L. (1985). A proposed standard set of problems to test finite element accuracy. Finite Elements in Analysis and Design, 1(1), 3–20.