# What exactly causes mesh locking in thin plate bending problems?

In thin plate bending problems, it would be very nice if we could model a thin rectangular plate with a mesh of many elements in the transverse directions, but only a single element in the longitudinal direction. However, when try to model a plate this way and use piecewise linear polynomials, the resulting numerically obtained displacements are much smaller when compared to known analytical solutions. This phenomenon is called "mesh locking". It can easily be resolved by either using higher order polynomials basis functions, mixed formulations, specialized plate/shell elements, or by utilizing more elements in the longitudinal direction.

I know that "mesh locking" happens in thin plate bending problems when you use low order elements and few elements in the longitudinal direction, but I'm still confused as to why it happens. What's the fundamental principle that prevents the low order basis functions from approximating the analytical solution? How do specialized formulations/elements and/or more elements in the longitudinal direction remedy this problem?

Let me try to expound a very simple argument.

A 3D element with piecewise linear shape functions can exactly represent a constant strain solution. (Consider e.g. the patch test.) Loosely speaking this means that when the exact solution exhibits a strain gradient, mesh size has to be adapted so that there is not too much change of strain across one single element. This holds true also for hexahedral elements, where some linear terms arise in the strain expression: however linear hex element cannot represent exactly an arbitrary linear strain field.

If we admit that at the limit of very thin slabs the 3D solution approaches the Kirchhoff–Love plate solution, then we have a very peculiar situation: mild strain gradient in the in-plane directions $x_1, x_2$, but a very steep gradient in the transverse (thickness) direction: remember that strain (stress) changes sign across $x_3$.

Now it should be clear that this strain/stress state cannot be sensibly approximated by a single linear FE in the $x_3$ direction.

In other terms, for thin plates (Kirchhoff-Love) the most simple straining solution is constant curvature $\kappa$: no gradient in $x_1, x_2$, but a gradient equal to $\kappa$ along $x_3$ of the in plane longitudinal strains $\varepsilon_{11}$, $\varepsilon_{22}$ and vanishing transverse shear strain $\gamma_{13}=\gamma_{23}=0$. A coarse mesh should be able to represent the "simple straining" state, but this is not possible with linear elements: in fact while, as pointed out by Bill Green, it is possible to have $\varepsilon_{11}$ and $\varepsilon_{22}$ linear across $x_3$, it is impossible to meet at the same time the vanishing transverse shear strain condition. These spurious shear stains are sometimes called "parasitic."

As a result the approximated elastic strain energy is wrong, and the finite element patch (under the assumption of complete integration) exhibits a stiffness which is bigger than expected. This point could be easily developed by computing the elastic strain energy of the exact solution, and comparing with the elastic strain energy of a single-element-across-the-thickness patch.

By the way this is a common situation, with any FE approximation. Actually what makes the bending problem pathological is the fact that for a single linear element in the thickness direction $x_3$, the strain energy is dominated by the parasitic transverse shear deformation, so that the displacement components are heavily underestimated.

Remedy for this problem are

• structural elements (plate/shell) in which Kirchhoff-Love or Reissner-Mindlin theory is incorporated as an assumed displacement field. (Gradient in $x_3$ direction is exactly modelled by this assumption.)

• more linear element in the thickness (more than 5) or use of higher order elements in order to better catch the linear strain/stress gradient along $x_3$, without the introduction of parasitic shear strains.

• selective reduced integration, which amounts to integrating the transverse shear field with a single Gauss point at the centroid of the element (where the parasitic shear fields vanishes) while still integrating other terms with more gauss points.

## Final notes

Thanks to Bill Green who criticised my first (simplistic) answer and forced me to write a longer (and I hope also better) answer. It was not my intention to give a complete account of mesh locking phenomena, which are indeed well understood and described in the literature.

I'm still convinced that, without going into the full details of shear locking, the starting point for explaining this phenomenon should be the observation that for pure bending the limiting solution for thin plates has very high strain gradients in the thickness directions, compared to the in plane ones. Linear elements are "constantish" in strain. Linear strains for hex elements arise from $x_1\cdot x_3$ (and similar mixed terms): the lack of squared terms (like $x_1^2$) in fact makes it impossible to properly decouple shear and longitudinal strains...

• I hope that my answer is not too simplistic, but the question is really very broad, and a short treatise could be written on this argument. – Stefano M Jun 15 '15 at 22:13
• This is wrong. The normal strains, e.g. $\sigma_{11}$ are well-approximated by the basic shape functions. The problem is that large, erroneous transverse shear stresses are also introduced. This phenomenon is often referred to as "shear locking" or "parasitic shear" in the literature. Many improvements to the formulation have been proposed including a simple one of integrating the transverse shear terms in the weak form with only a single integration point. A definitive discussion of this problem can be found in the book by MacNeal, Finite Elements: Their Design and Performance. – Bill Greene Jun 16 '15 at 0:37
• @BillGreene you are right, and your comment has proven that my argument can be misleading. But let me add that I know that when you compute the strain energy for a single element subjected to bending, you have linear $\varepsilon_{11}$ along $x_3$ and shear strain $\gamma_{13}$, and that the contribution of $G\gamma_{13}^2$ is dominant in the strain energy (hence the "shear locking"). Nevertheless I would hold the argument that "linear elements" cannot represent strain gradients exactly, in the very simple sense that linear normal strain but vanishing shear strain is not possible. – Stefano M Jun 16 '15 at 6:20

You just need a finite element that does not manifest shear locking. There is software that implements such element, based on the Mindlin model:

http://members.ozemail.com.au/~comecau/plate_analytical.htm

This software was developed by our company; but it is impossible to provide here full mathematical details of the method implemented in the program. As an alternative, you may consider using an approach used in the following book (it is devoted to shell elements, i.e with curved surfaces, and is based on the Kirchhoff-Love geometric hypothesis; still, that element does not manifest shear locking either, and you may consider adaptation of that element for your purposes):

http://www.amazon.com/Computational-Geometry-Surfaces-Application-Analysis/dp/0646930818

There is one more thing to be considered in terms of shear locking. In the theory of plates and shells, the displacements (deflections) field is the principal function that defines all other parameters of the mechanical behavior. Rotations of the nodes represent the 1st derivatives of the displacements; bending moments are the second derivatives; shear forces correspond to the 3rd derivatives. Quite obvious that if the deflections are not interpolated by polynomials of sufficiently high order, shear forces would just not be evaluated with any acceptable degree of precision at all.

• Welcome to Scicomp.SE. I appreciate that you're (also) trying to be helpful, but the answer-to-advertisement ratio of your posts dangerously low; this is not what this site is for. (And at the very least, some disclosure that this is your software would be in order.) – Christian Clason Apr 17 '16 at 10:36
• Thank you, Christian; I have expanded the answer. Hopefully is now corresponds to the criteria of acceptability. – SparseSolverCodes Apr 17 '16 at 11:13