Where Shear Locking Comes From

In my case study on shear locking, a finite element (FEM) simulation of a bent flat steel bar came out about 60% too stiff with linear bricks on a coarse mesh, and still about 10% off after refining it to four elements through the thickness. Shear locking is an effect in simulations where the modelled rod resists bending more than a real rod does. The case study showed it experimentally, and here I want to show where it comes from mathematically.

Imagine taking a cylindrical rod and bending it into an arc of a circle by applying a moment, a pure turning effect with no net force. There is no shear force, no twist and no stretching of the axis. That is pure bending. In practice it rarely appears on its own, and anything worth simulating is more complicated. But any slender part that bends contains this deformation at its core. If an element cannot represent pure bending, it cannot represent the bending part of a more complicated problem either, so the defect appears there too. How big it is will differ from problem to problem.

FEM cuts the body into small elements and replaces the unknown displacement inside each one by a simple function fixed by its values at the nodes. The solver finds those nodal values. For bending, brick (hexahedral) elements are standard engineering practice, and there are two common classes, linear and quadratic. Measure the coordinates from the centre of the element, $\tilde x=x-X$, $\tilde y=y-Y$, $\tilde z=z-Z$. In the linear brick every displacement component is a sum of eight terms,

$$u_x=\sum_{i,j,k\in\{0,1\}}\alpha_{ijk}\,\tilde x^{i}\tilde y^{j}\tilde z^{k},$$

and the same for $u_y$ with coefficients $\beta_{ijk}$ and $u_z$ with $\gamma_{ijk}$. Eight coefficients need eight known values, so the element has 8 nodes, the corners, and the nodal values fix the coefficients one to one. In the quadratic brick the powers run over $\{0,1,2\}$, which gives 27 coefficients and 27 nodes. The question I am asking is whether the exact solution of the bending problem can be written in the form an element can build. If it can, the solver is able to find it. If it cannot, no choice of nodal values will ever reproduce it.

So first, the exact solution. In §17 of Landau and Lifshitz, Theory of Elasticity (Chapter II, Bending of rods), a rod bent with radius of curvature $R$ and made of a material with Poisson's ratio $\sigma$ has the displacement field

$$u_x=\frac{xy}{R},\qquad u_y=-\frac{x^2+\sigma(y^2-z^2)}{2R},\qquad u_z=-\frac{\sigma yz}{R}.$$

Here $u_x$ and $u_z$ contain only products of different coordinates, but $u_y$ contains $x^2$, $y^2$ and $z^2$.

The linear brick cannot build these squares. Along one edge it has two nodes, at $\tilde x=\pm a$, and at both of them $\tilde x^2=a^2$. The element cannot tell the square from a constant. Now let's see what this does to the strain. The shear strain is

$$u_{xy}=\frac12\left(\frac{\partial u_x}{\partial y}+\frac{\partial u_y}{\partial x}\right),$$

and in pure bending it must be zero everywhere. Differentiating the element's own sums,

$$\frac{\partial u_x}{\partial y}=\sum_{i,k}\alpha_{i1k}\,\tilde x^{i}\tilde z^{k},\qquad \frac{\partial u_y}{\partial x}=\sum_{j,k}\beta_{1jk}\,\tilde y^{j}\tilde z^{k}.$$

In the second sum $\tilde x$ has disappeared, because $u_y$ is only linear in $\tilde x$. So only $u_x$ can produce an $\tilde x$ in the shear:

$$2u_{xy}=\alpha_{110}\,\tilde x+\alpha_{111}\,\tilde x\tilde z+\ (\text{terms in }\tilde y,\tilde z\text{ only}).$$

No shear means $u_{xy}=0$ for all $x$ in the element, which requires $\alpha_{110}=0$. But $\alpha_{110}$ is the coefficient of $\tilde x\tilde y$ in $u_{xx}$:

$$u_{xx}=\frac{\partial u_x}{\partial x}=\sum_{j,k}\alpha_{1jk}\,\tilde y^{j}\tilde z^{k} .$$

Pure bending needs $u_{xx}=y/R=(Y+\tilde y)/R$, which requires $\alpha_{110}=1/R$. The two conditions contradict each other: $\alpha_{110}=1/R$ for the right bending strain, $\alpha_{110}=0$ for no shear. This holds whatever the nodal values are. The exact solution does not belong to the set of functions a linear brick can build. In the exact field the $x^2$ term in $u_y$ cancels the shear that the $xy$ term in $u_x$ creates, and the linear element has no such term. The element therefore reports a shear that does not exist, that fake shear stores energy, and the element comes out too stiff. That is shear locking.

In the quadratic brick every term in the three components has each coordinate at power at most two, so the whole field is in the family, $x^2$, $y^2$ and $z^2$ included. The solution belongs to the set of functions a quadratic brick can build.

What remains is to be sure that the solver actually finds it. A real elastic body settles into the state of minimum total potential energy, which is the elastic strain energy minus the work done by the external loads. FEM minimizes the same quantity, but only over the functions the mesh can build. If the exact solution is among them, it is also the minimum among them, so the solver finds it. This holds for a linear elastic material with proper boundary conditions and nothing too exotic in the loads.

So the linear brick does not merely give a less accurate answer in bending. Bending is a state it cannot represent, so a stiffness error is built in, and refining the mesh only shrinks it, as the bar showed, from about 60% on the coarse mesh to about 10% on the fine one. The quadratic brick contains that state, which is why it is the better element for this kind of problem. And as the case study showed, a coarse mesh of quadratic bricks landed within 1% of the exact answer and ran two to seven times faster than the refined linear mesh.