Question 1 of 3: Ten short questions on finite element theory

Problem 1: Ten short questions on finite element theory (one of three equal-value problems)

Question text not reproduced: the examination questions are © Engineers and Geoscientists BC. Open the official past paper (linked at the top of this page) to read the question, then follow the worked solution below.

Given. No numerical data; the problem is a ten-part examination of the theory that underpins every displacement-based finite element code.

Find. A short, defensible answer to each of the ten parts, with the two “draw” and “identify” parts supported by a figure.

Part 1 — why the finite element solution does not satisfy equilibrium inside the domain

The displacement-based finite element method is not derived from the differential equations of equilibrium. It is derived from a variational statement — the principle of virtual work, equivalently the principle of minimum total potential energy — which is the weak (weighted-integral) form of those equations. The trial field is chosen to satisfy the kinematic requirements exactly: it is continuous across element boundaries, it satisfies the prescribed displacement boundary conditions, and the strains derived from it are compatible by construction. Equilibrium is then imposed only in the averaged sense

and only for the virtual displacements $\delta\{u\}$ that live in the same finite-dimensional space as the solution. Assembling that statement element by element produces $[K]\{d\}=\{F\}$, whose rows are nothing more than force balances at the nodes.

Consequently the stresses recovered as $\{\sigma\}=[D][B]\{d\}$ come from differentiating an approximate displacement field. Differentiation degrades the approximation, so in general

the computed traction on a loaded surface does not equal the applied traction, and the tractions transmitted across an inter-element boundary are discontinuous even though the displacements are continuous. Only three things are exact: compatibility, the essential boundary conditions, and nodal equilibrium of the assembled model. The equilibrium residual is driven to zero only in the limit of mesh refinement (or of increasing polynomial order), which is precisely what convergence means for this class of element.

Two useful corollaries follow. First, a constant-strain element satisfies the internal differential equations trivially — every derivative of a constant is zero — and yet still violates the traction boundary conditions; Problem 3.3 below measures exactly that violation. Second, the jump in recovered stress between adjacent elements is a free error estimator, and smoothing it is the basis of the Zienkiewicz–Zhu superconvergent patch recovery used in commercial adaptive codes.

Part 2 — the strain energy of a finite element solution

The admissible displacement fields available to the model form a subspace of the admissible fields available to the exact problem, so minimising the total potential energy $\Pi=U-W$ over the smaller space cannot do better than minimising it over the larger one: $\Pi_{h}\geq\Pi_{\text{exact}}$. For a structure loaded by prescribed forces the equilibrium value of the potential is $\Pi=-U$, and therefore

The physical reading is the familiar one: constraining the structure to deform only in the shapes the mesh can represent adds artificial constraint, so the finite element model is stiffer than the real structure, the displacements are under-predicted, and the stored energy is under-predicted with them. Convergence is monotone from below in the energy norm as the mesh is refined. The same argument run in reverse explains why a displacement-based model over-predicts natural frequencies: extra stiffness with the same mass raises every $\omega$.

Two carve-outs are worth remembering for the exam. The bound holds for the compatible displacement formulation with consistent nodal loads; incompatible-mode, reduced-integration and hybrid elements deliberately break the bound in order to soften an over-stiff element, and for those the strain energy may come out above the exact value.

Part 3 — shape functions for the two bar elements

A bar element carries one degree of freedom per node, the axial displacement, so the shape functions are the Lagrange polynomials that take the value one at their own node and zero at every other node. Figure 1-a has two nodes and therefore linear functions; Figure 1-b has three, so the functions are quadratic — and because the interior node sits at $L/3$ rather than at midspan, they are not the symmetric textbook quadratics.

For Figure 1-a, with $x$ measured from node 1,

two straight lines crossing at midspan. Their derivatives are constant, so this element can reproduce only a uniform axial strain — it is the bar analogue of the constant-strain triangle.

For Figure 1-b the nodes are at $x_{1}=0$, $x_{3}=L/3$ and $x_{2}=L$, and the Lagrange construction gives

Three features should appear in a hand sketch and are drawn above. Each curve is a parabola equal to one at its own node and zero at the other two; the three sum to unity at every point, which is what allows the element to represent a rigid-body translation; and because node 3 has been pushed towards node 1, $N_{1}$ swings negative over the right-hand part of the element, reaching its minimum

while $N_{3}$ overshoots to $1.125$ at midspan. A negative shape function is perfectly legitimate — only the partition-of-unity and Kronecker-delta properties are required — but it warns that an off-centre interior node degrades conditioning, and that placing it at $L/2$ is the numerically sensible default. The derivatives are now linear in $x$, so this element reproduces a linearly varying axial strain exactly.

Part 4 — what a hybrid finite element is

A hybrid element is one derived from a multi-field variational principle rather than from the single-field principle of minimum potential energy. Two or more fields are assumed independently — most commonly, in Pian's hybrid-stress element, an equilibrating stress field $\{\sigma\}=[P]\{\beta\}$ inside the element together with a compatible displacement field $\{u\}=[L]\{d\}$ on the element boundary only. The two are tied through the Hellinger–Reissner functional, the internal parameters $\{\beta\}$ are eliminated at element level by static condensation, and what emerges is an ordinary stiffness matrix in the usual nodal displacement degrees of freedom that can be assembled alongside conventional elements.

The label “hybrid” refers to that mix of an interior field with a separate boundary field; the closely related “mixed” elements interpolate both fields over the whole domain. The motivation is accuracy in stress rather than in displacement: because the interior field is chosen to satisfy the equilibrium equations exactly, hybrid elements give markedly better stresses on coarse meshes, are far less sensitive to element distortion, and can be made free of shear locking in plates and of volumetric locking as Poisson's ratio approaches 0.5. The price is a more elaborate formulation, the need to choose the stress basis $[P]$ carefully to avoid rank deficiency, and the loss of the guaranteed lower bound on strain energy noted in part 2.

Part 5 — what reduced integration means

The element stiffness matrix is evaluated numerically,

and “full” integration means using a Gauss rule of high enough order to evaluate that integrand exactly for an undistorted element — $2\times2$ for the bilinear quadrilateral, $3\times3$ for the eight-node quadrilateral. Reduced integration deliberately uses a lower-order rule: $1\times1$ for the Q4, $2\times2$ for the Q8.

There are three reasons to do it. It is cheaper, by a factor of four in two dimensions. It softens an element that low-order interpolation has made artificially stiff, and in doing so it cures the shear locking of thin bending-dominated elements and the volumetric locking of nearly incompressible material; the selective variant reduces only the shear or only the volumetric part of $[D]$ and integrates the rest fully. And the reduced Gauss points are the Barlow points, at which the sampled stresses are one order more accurate than anywhere else in the element, which is why codes extrapolate nodal stresses from them.

The danger is rank deficiency. Too few sampling points can leave deformation patterns that produce zero strain at every sampling point and therefore store no energy — the spurious zero-energy or hourglass modes. Part 7 quantifies this for the Q4. Production codes control it with hourglass stabilisation stiffness or by using at least one fully integrated element in each patch.

Part 6 — defects in the mesh of Figure 1-c

[Figure not reproduced: Figure 1.2 — Figure 1-c redrawn. Elements 1 and 2 are four-node bilinear; elements 3 and 4 carry mid-side nodes. The four flagged locations are the defects discussed below. See the official exam paper.]

Every defect in this mesh is a version of the same failure: the displacement field assumed along one side of an interface does not match the field assumed along the other side, so the mesh is non-conforming and the $C^{0}$ continuity on which the convergence proof rests is lost.

(A) Linear edges meeting a quadratic edge. Elements 1 and 2 are four-node bilinear, so displacement varies linearly along their right-hand edges. Element 3 carries a mid-side node on its left edge, so displacement there varies quadratically. The two fields agree at the three shared nodes and nowhere in between: a gap opens along one half of the interface and material overlaps along the other, in proportion to the quadratic bulge. Material is created and destroyed at the interface, the strain field has a spurious singularity there, and monotone convergence is no longer guaranteed.

(B) A hanging node. The top-left corner of element 4 lands part-way along the right edge of element 3, at that edge's mid-side node. Element 4 has no corresponding node there in its own connectivity sense — its edge runs only from that point downward — so the interface is a T-junction. The upper half of element 3's right edge is unmatched, and the displacement of element 4 along its left edge is linear against a quadratic in element 3.

(C) An unmatched mid-side node. The mid-side node on the bottom edge of element 3 has no partner in whatever lies below it. If nothing lies below, it is an unconstrained boundary node and harmless; if an element does adjoin, its edge is linear and the same incompatibility as (A) applies.

(D) An inconsistent node set. Element 4 carries mid-side nodes on two of its four edges only. That is not a member of the serendipity family: a valid isoparametric element must use one complete shape-function set, and mixing eight-node edges with four-node edges on the same element makes the interpolation neither $Q4$ nor $Q8$.

Three remedies are standard. Use a proper transition element — the five-, six- and seven-node serendipity variants have shape functions written specifically for an element with mid-side nodes on some edges and not others, and they restore compatibility exactly. Or impose multi-point constraints tying the mid-side node's displacement to the average of its two corner neighbours, which forces the quadratic edge to behave linearly. Or simply grade the mesh so that only elements of the same family ever share an edge, refining through a band of elements rather than at a single interface. The last is what a modern mesher does automatically, and it is the answer to give if only one is wanted.

Part 7 — zero eigenvalues of a square bilinear element in plane strain

A square four-node bilinear element has two degrees of freedom per node, so $[k]$ is $8\times8$. Integrated fully with the $2\times2$ Gauss rule its rank is 5, and therefore

The three modes are the rigid-body modes available to a body in a plane: translation in $x$, translation in $y$, and rotation about the out-of-plane axis. Each of them produces zero strain everywhere in the element, hence zero strain energy, hence $[k]\{d\}=\{0\}$. Their presence is not a defect — it is a necessary condition. An element whose stiffness matrix could not move rigidly without storing energy would fail the patch test and would lock; the three zeros disappear the moment the element is attached to supports that remove the rigid-body freedoms of the assembled structure. The count is unchanged by the choice of plane strain over plane stress, because the constitutive matrix $[D]$ is positive definite in both cases and only $[B]$ decides the null space.

The interesting comparison, and the reason the question is asked, is with reduced integration. Sampling at the single centre point makes the rank of $[k]$ only 3, so there are five zero eigenvalues: the three legitimate rigid-body modes plus two spurious hourglass modes, in which the element deforms into a bow-tie while producing zero strain at the centre. A mesh of such elements can develop a zero-energy, checkerboard deformation pattern under load, which is why hourglass control is mandatory whenever one-point integration is used.

Part 8 — consistent and diagonal mass matrices

In a dynamic analysis the semi-discrete equations are $[M]\{\ddot{d}\}+[C]\{\dot{d}\}+[K]\{d\}=\{F(t)\}$, and the question is how the inertia is distributed. The consistent mass matrix is obtained from the same interpolation used for the stiffness,

and it is symmetric, positive definite, banded and full within each element. The diagonal or lumped matrix concentrates the element mass at its nodes, by simple partitioning of the total mass, by row-summing the consistent matrix, or by the HRZ scheme that scales the consistent diagonal to preserve total mass.

The consistent matrix is the more accurate of the two and it is variationally correct: it is derived from the same kinetic-energy functional as the stiffness, it converges monotonically, and it gives an upper bound on the natural frequencies. Its disadvantages are cost and structure. It must be factorised, or at least solved with, at every time step, which makes it unattractive for explicit integration, and it costs more storage.

The lumped matrix is diagonal, so $[M]^{-1}$ is free and the central-difference explicit scheme becomes a sequence of scalar updates — the reason every explicit crash and blast code uses it. It is also more robust for wave propagation, where it damps the spurious high-frequency oscillation that the consistent matrix carries, and it tends to under-estimate frequencies, so the two together bracket the exact answer. Its drawbacks are lower accuracy for the same mesh, ambiguity over how much rotary inertia to assign to the rotational degrees of freedom of beams, plates and curved surface elements, and the fact that naive row-summing of a higher-order element such as the Q8 produces negative corner masses, which is why HRZ lumping exists. A common compromise in implicit codes is to average the two matrices.

Part 9 — analysing a plated reinforced concrete beam

The strengthened beam is a three-material composite — concrete, internal reinforcement, and a steel plate glued or bolted to the soffit — and the useful answer is a sequence of modelling decisions rather than a formula.

Idealisation and elements. A two-dimensional plane-stress model through the span captures flexure and shear economically; a three-dimensional solid model is needed only if the plate is narrower than the web, if the anchorage detail matters, or if torsion is present. Mesh the concrete with continuum elements, several through the depth so that the neutral axis and the compression block are resolved. Represent the internal bars either discretely, as two-node truss elements sharing nodes with the concrete mesh, or with the embedded-reinforcement option that most codes provide, which smears the bar stiffness into the parent element without requiring the mesh to follow the bar. Model the bonded plate with plane-stress or thin plate-bending elements of the true thickness.

The interface is the critical modelling choice. Plated beams almost never fail in flexure — they fail by debonding at the plate ends or by concrete-cover separation, so a perfectly bonded plate that shares nodes with the concrete will over-predict capacity and miss the governing mechanism entirely. Insert interface, cohesive-zone or spring elements between plate and concrete carrying a bond–slip law with a finite shear strength and a fracture energy, and check the interfacial shear against that law. Where the plate is bolted, model the bolts as connectors with their own slip.

Material models. Concrete needs a nonlinear law — smeared cracking or a concrete damaged-plasticity model — with a compressive stress–strain curve, a tensile strength, and tension softening defined through a fracture energy so that the answer is mesh-objective. Reinforcement and plate take an elastic–plastic law with strain hardening; an FRP plate would instead be linear elastic to brittle rupture.

Solution and verification. The analysis is incremental–iterative (Newton–Raphson, or arc-length once softening begins), under displacement control so that the post-peak branch can be traced. Refine the mesh at supports, under point loads and at the plate ends where the stress gradient is steepest, and run a mesh-convergence study. Validate the load–deflection curve and cracking moment against a hand section analysis and, where possible, against a test; then post-process deflections, crack pattern and width, steel and plate stresses, and interfacial shear. The design checks the model is feeding are the Canadian ones — CSA A23.3 for the reinforced concrete section, and the ISIS Canada design manuals (with CSA S806 for FRP) for the externally bonded strengthening.

Part 10 — Euler-Bernoulli and Timoshenko beam elements

The two elements differ in one kinematic assumption. Euler–Bernoulli theory requires plane sections to remain plane and normal to the deformed axis, so the transverse shear strain is identically zero and the section rotation is slaved to the slope, $\theta=dv/dx$. There is a single unknown field $v(x)$, its governing equation is fourth order, and the element must therefore provide $C^{1}$ continuity — hence the cubic Hermite shape functions and the two degrees of freedom per node, $v$ and $dv/dx$, of the classical beam element. It is exact for a prismatic beam under end loading, which is why the two-element model of Problem 2 below returns exact reactions.

Timoshenko theory keeps plane sections plane but drops the normality requirement, admitting a shear strain

with $v$ and $\theta$ now interpolated independently. Only $C^{0}$ continuity is needed, so simple Lagrange functions suffice, and the element requires a shear correction factor $k$ — $5/6$ for a rectangle, $9/10$ for a circle — to reconcile the assumed uniform shear stress with the true parabolic distribution. The extra flexibility lowers the stiffness and lowers the computed frequencies. Its own pathology is shear locking: with equal-order interpolation and full integration the element cannot represent pure bending without generating spurious shear, and becomes uselessly stiff as the beam gets thin. Reduced or selective integration of the shear term, or a linked interpolation, cures it.

The choice follows the slenderness. Use Euler–Bernoulli for slender members, roughly $L/h\geq10$ to 20 — ordinary building beams, columns and frames — where shear deformation contributes a per-cent or two at most. Use Timoshenko for deep beams and short spandrels, for sandwich and laminated composite members whose transverse shear modulus is low relative to the bending modulus, for thin-webbed plate girders, and for dynamic problems in which the higher modes matter, because shear and rotary inertia depress those frequencies substantially. Timoshenko degenerates to Euler–Bernoulli as $EI/(kGAL^{2})\to0$, so when in doubt the shear-flexible element is the safe default.