The corners the pressure pulls up

Moving a force found that two systems of forces are interchangeable on a rigid body when they have the same resultant and the same moment about every point. Equivalent in work, not in resultant found that a finite element is not a rigid body, and that the right substitute for a load spread along a beam element is the set of nodal actions that does the same work through every deflected shape the element can take. For a uniform load that set carries a couple at each end, $wL^2/12$, which no resultant can see, and with it a one-element cantilever is exact where the resultant split between the ends is a third wrong.

That essay closed on the case one dimension up. A plate or shell element carries load over an area, its shape functions are surfaces rather than curves, and the same principle applied to a uniform pressure produces nodal forces that “change sign and look like a mistake”. This essay is that case, and the mistake it looks like turns out to be three things at once: the correct load vector, the reason nodal-force output cannot be read as a load path, and — through the same integrals applied to mass — a negative mass that a whole family of dynamic calculations cannot run with.

One element, eight nodes, a downward pressure #

The eight-node serendipity element is the ordinary quadratic quadrilateral of a thousand finite-element programs: four corner nodes, four mid-side nodes, and a displacement inside it interpolated by eight shape functions, each equal to one at its own node and zero at the other seven. A uniform pressure $p$ on an element of area $A$ is replaced by nodal forces $f_i = \int N_i\,p\,dA$ — the work the pressure does through the shape of node $i$ — which is the whole of the work-equivalence principle for an area.

The mid-side nodes receive a third of the load each, and push down with it. The corners receive minus a twelfth each, and are pulled up. The eight shares add to one, as they must — four thirds less four twelfths — so the resultant is exact, and because the shape functions reproduce any linear field exactly, so is the moment about every axis. By the standard of moving a force the set is equivalent to the pressure. By the standard of every intuition about where a load goes, it is upside down at four of its eight nodes.

Why a corner has a negative share #

The share a node receives is the integral of its shape function over the element, so a negative share means a shape function that is mostly negative.

The corner function has to be one at its corner and zero at the seven other nodes, and with the eight terms the element allows, the only way to arrange that is a surface that dives below zero across the middle. It is zero along the straight line joining the two mid-side nodes beside its corner, positive in the small triangle between that line and the corner, and negative over the remaining seven eighths of the element. Integrated, the negative region wins: $0.032 - 0.115 = -1/12$.

So the corner’s negative load is not an artefact of a numerical rule or a sign convention. It is what the corner’s own deflected shape does under a uniform pressure: when the corner is pushed up by a unit amount with the other seven nodes held, most of the element moves down, and a downward pressure then does positive work. The pressure really does want to lift the corner, in the precise sense that lifting it releases potential energy.

Other elements, other signs #

The pattern belongs to the eight-node element, and two relatives show that it is a property of the shape functions and not of quadratic elements in general.

Add a ninth node at the centre and every share becomes positive: 1/36 at each corner, 1/9 at each mid-side, 4/9 at the centre. These are Simpson’s rule applied in both directions, because the nine-node element’s shape functions are products of one-dimensional quadratics, each of which integrates to Simpson’s weights of 1/6, 4/6, 1/6. The eight-node element has the same number of nodes on its boundary and no centre, and to cover the element without one, its corner functions have to be pulled negative across the middle.

The six-node triangle does something stranger still: its corners receive exactly nothing. Each corner function is positive near its corner and negative along the opposite side, and the two cancel exactly, so a uniform pressure on a mesh of quadratic triangles is carried entirely by the mid-side nodes. The four-node quadrilateral and the three-node triangle, whose shape functions never go negative, give the shares intuition expects: a quarter and a third to each node.

An assembled mesh, and what a printout shows #

In a mesh each node collects the shares of every element it belongs to. For the eight-node element that is uniformly unhelpful.

Every node at an element corner collects a negative load: −1/12 of an element’s pressure at the corners of the mesh, −1/6 along its edges where two elements meet, −1/3 inside where four meet. Every mid-side node collects a positive one, 1/3 on the boundary and 2/3 inside. The 15 negative loads and the positive ones add up to the whole pressure, which is all a load vector has to do.

This is the point at which the substitution stops being a numerical nicety. A printout of these forces is a correct load vector and not a picture of where the load goes. A designer reading a slab model’s nodal loads, or a mat foundation’s nodal reactions under a uniform bearing pressure, sees alternating signs and concludes that the corner nodes are in uplift — a pile under tension, a bearing lifting off — when the stress everywhere in the slab is compressive and smooth. The same substitution in reverse turns a reaction pattern into load shares that nobody should sum over a tributary area, which is the question the load a beam is given answers from the other direction: the load on a supporting member is an integral of a stress along its line, and a node is not a line.

Why intuition was trained on the wrong elements #

The intuition the corners offend — that a node carries the load on the area around it — is not wrong for the elements it was learned on. For the four-node quadrilateral and the three-node triangle, whose shape functions are planes or bilinear surfaces that never go negative, the work-equivalent share of a uniform pressure is exactly the node’s tributary area: a quarter of the element to each corner of a rectangle, a third to each corner of a triangle. Divide a slab into those elements, give each node the load on the area nearest it, and the result is the consistent vector to the last digit. Nothing distinguishes the two rules, so nothing teaches that they are different rules.

They part company as soon as the shape functions curve, because a curved function that is one at its node and zero at every other node must overshoot somewhere between them. On a quadratic element the overshoot is into negative values, and it is largest for the node whose neighbours are closest on every side — the corner of a serendipity element, flanked by two mid-side nodes and missing the centre node that would have taken some of the curvature. The tributary rule and the work rule were only ever the same rule for linear interpolation, and a quadratic element is chosen precisely because it is not linear.

The same thing happened one dimension down, where the beam element took couples at its ends that no tributary picture contains. There the missing term was invisible in any printout of forces, because a couple is not a force. Here it is visible, as a force with the wrong sign, which is why the area element is the one that alarms people and the beam element is the one that quietly gives the wrong deflection.

What the right vector is right about #

The work-equivalent vector is called consistent because it is the one the element’s own stiffness was derived with. The finite-element solution is then the displacement field, among all those the mesh can represent, that makes the total potential energy of the loaded structure smallest — the same principle the answer that depends on how it was divided described as convergence to a minimum over a restricted space. A load vector that is not work-equivalent is a different load, and the solution minimises the energy of that one instead.

A concrete wall 4 m long and 1 m deep, cantilevered and carrying its own weight, shows what that buys and what it does not. With one element along its length the work-equivalent vector gives a tip deflection of 0.799 of the fine-mesh answer. A split that gives an eighth of the weight to each node gives 0.996 — almost exact, with the wrong load. At two elements the work-equivalent vector gives 0.955 and the equal split 1.005, now past the answer; the corners-only split stays above it throughout, 1.113 falling to 1.001.

The consistent vector is not the most accurate on every mesh. It is the only one that promises which side of the answer it is on. A displacement element is too stiff, and with the load it was built for, the stiffness error has a sign: the strain energy is always underestimated, and a deflection that tracks it comes out low and rises monotonically as the mesh is refined, which it does here on every mesh. The equal split puts more weight at the corners — including the corners at the free end — and that extra lever happens to cancel most of the element’s stiffness on the coarsest mesh. The cancellation is luck, it changes sign at the next refinement, and a user with no fine-mesh answer to compare against cannot know which side of the truth a lumped load has landed on. The diagonal-scaled split happens to approach from below here too, but nothing promises it.

When one model’s reactions become another’s loads #

The alternating signs do their real damage at the boundary between two models. A slab is analysed with plate elements; its reactions along a supporting beam are read off node by node and applied to a separate model of the beam. With four-node elements the reactions are a sensible, smooth sequence of downward forces. With eight-node elements they alternate — large at the mid-side nodes, small or negative at the corners — and a beam model loaded by those point forces carries a sawtooth of shear and a moment diagram with a ripple in it, every value of which is a correct consequence of a load pattern that does not exist.

The total is right, so the check that every engineer makes — that the reactions sum to the load — passes. The moment at mid-span is nearly right, because the ripple averages out over the span. What goes wrong is local: the shear at a node, the bearing force at a connection that happens to sit at an element corner, the reaction at a pile that happens to sit at one. The repair is to hand the load across as what it is, a line load equal to the integral of the slab’s edge shear, rather than as the vector the slab’s own mesh needed; or to make the two models share a mesh, so that the vector is consistent for both. Either way the nodal forces of a quadratic mesh are an instruction to that mesh and not a description of the structure, and everything adds to nothing only in the sum.

The same integrals make a negative mass #

The integrals that give the consistent load, $\int N_i\,dA$, are the row sums of another matrix: the consistent mass matrix, $M_{ij} = \int \rho N_i N_j\,dA$, since the shape functions add to one everywhere. A dynamic calculation that marches through time explicitly needs a diagonal mass, one number per node, and the simplest way to make one is to sum each row of the consistent matrix onto its diagonal.

For the eight-node element that gives each corner a mass of minus a twelfth of the element’s, and a time-stepping scheme with a negative mass is unstable at any step. For the six-node triangle it gives the corners no mass at all, which makes the step infinitely short. So the lumping in general use scales the diagonal of the consistent matrix instead, in proportion, to the element’s total mass — Hinton, Rock and Zienkiewicz’s rule — which gives the eight-node corners 3/76 and the triangle’s 1/19, both small and both positive. The nine-node element needs no such device: its row sums are already positive, and the two rules agree.

The trade runs the opposite way from the one for loads. For a static load the consistent vector is the one to use, and a lumped one is an approximation with no promised sign. For an explicit dynamic calculation the consistent mass cannot be used at all with these elements, and the lumped one — chosen by a rule designed only to keep masses positive — is what makes the calculation possible, at a cost in accuracy that shows up in the structure’s higher periods.

One element, by hand #

For the eight-node element on the square $-1 \le \xi, \eta \le 1$, the shape function of the corner at $(-1, -1)$ is

$$N_1 = \tfrac14(1-\xi)(1-\eta)(-\xi-\eta-1).$$

Its integral over the square, of area 4, splits into three pieces: $\tfrac14\int\!\!\int (1-\xi)(1-\eta)(-1)\,d\xi\,d\eta = \tfrac14 \times 2 \times 2 \times (-1) = -1$, and two pieces of the form $\tfrac14\int\!\!\int (1-\xi)(1-\eta)(-\xi)\,d\xi\,d\eta = \tfrac14 \times \tfrac23 \times 2 = \tfrac13$ each, since $\int_{-1}^{1}(1-\xi)(-\xi)\,d\xi = \tfrac23$. The total is $-1 + \tfrac13 + \tfrac13 = -\tfrac13$, and as a share of the element’s area of 4 that is $-\tfrac1{12}$. The mid-side function $N_5 = \tfrac12(1-\xi^2)(1-\eta)$ integrates to $\tfrac12 \times \tfrac43 \times 2 = \tfrac43$, which is $\tfrac13$ of the area. Four of each: $4 \times (-\tfrac1{12}) + 4 \times \tfrac13 = 1$.

Straight sides, uniform pressure, and a static load #

The element is undistorted. On a square or rectangle the shares are exactly these fractions. A distorted element — a trapezium, a curved edge — has shares that depend on its shape, and the corners of a badly distorted eight-node element can receive shares more negative than a twelfth.

The pressure is uniform over each element. A pressure that varies linearly across an element shifts the shares toward the heavier side without changing their signs: for a pressure that falls from twice its mean to nothing across the element, the corners on the heavy side receive −1/18 of the element’s load and those on the light side −1/9.

The wall is a plane-stress continuum fixed along its whole end. Its fine mesh agrees with a deep cantilever’s bending-plus-shear deflection to a third of a per cent; the comparison between splits is a comparison of meshes of the same wall, and the ranking depends on the quantity compared. A split that is lucky for the tip deflection is not lucky for the stress at the root.

What the pictures cannot show #

What any of this does to stresses, which are what a designer usually reads. A consistent load gives stresses that converge smoothly with the mesh; a lumped one gives stresses that are wrong near every node where the load was moved, in the way a load spread out and the force that replaces it differ near the point of application. The error is local and dies away over an element or two, and on a coarse mesh an element or two is the whole structure.

Still open: the plate element whose rotations carry load #

Every element here has only translations at its nodes. A plate-bending element also has rotations there, and its work-equivalent load for a uniform pressure includes nodal moments, as the beam element’s did — couples that a resultant cannot see and that point, at some nodes, the way intuition would not. Whether those moments on a slab mesh read out as a pattern of fictitious hogging at every element corner, and whether the reactions of a slab modelled that way can be summed along a support line at all, is the question the plate that carries a third of its load by twisting asks of any mesh used to design one.