Prev
Unit 10 · Boundary Value Problems and Partial Differential Equations
The Finite Element Method in Two Dimensions
70m100 XP
Next
🌡️Boundary Value Problems and Partial Differential Equations

The Finite Element Method in Two Dimensions

The weak form in the plane, linear triangles, assembly, Prandtl torsion, the patch test, corner singularities and mesh convergence

Gandsfjord, 23 August 1991
Sleipner A is to be a gas platform in the North Sea. Its gravity base — a cluster of twenty-four concrete cells, each about 24 m across, standing on the sea bed — is built in a dry dock and then floated to a fjord near Stavanger, where it is ballasted down step by step in a controlled test, to check its strength before the deck is installed. Partway through the test a loud noise is heard. A wall of one of the 'tricells' — the triangular spaces where three cylinders meet — has failed. Water floods in faster than the pumps can remove it. Within minutes the base sinks to the bottom of the fjord, 220 m down, and hits the sea floor with the force of a small earthquake. Nobody is hurt. The economic loss is estimated at around 700 million dollars. The investigation by SINTEF traced the failure to the design calculations. A finite element analysis of the cells, using elements that were too coarse and too distorted for the steep stress gradients in the tricell walls, had underestimated the shear stresses by 47%. The walls were designed too thin. The replacement platform was designed with a refined analysis and installed in 1993. Sleipner became the standard case study on the most important habit in computational engineering: check that your mesh has converged.
Lesson 10-5 built the method on a line. In the plane the recipe is identical — energy or weak form, shape functions on elements, element matrices, assembly — and only the geometry changes. The model problem is Poisson's equation, which covers steady heat conduction, seepage (lesson 10-4) and the torsion of shafts:
−∇2u=f  in a region Ω,u=0  on its boundary.-\nabla^2u = f \ \text{ in a region } \Omega, \qquad u = 0 \ \text{ on its boundary}.
Its energy is Π(u)=∫Ω12∣∇u∣2 dA−∫Ωf u dA\Pi(u) = \int_\Omega \tfrac12|\nabla u|^2\,dA - \int_\Omega f\,u\,dA, and the region is now covered by a mesh of triangles that can follow any outline: a curved wall, a hole, a fillet.
Derivation
The weak form in two dimensions
Multiply −∇2u=f-\nabla^2u = f by a test function vv that vanishes on the boundary and integrate over Ω\Omega. The two-dimensional version of integration by parts is Green's identity (the divergence theorem applied to v∇uv\nabla u):
∫Ω−∇2u  v dA=∫Ω∇u⋅∇v dA−∮∂Ωv ∂u∂n ds.\int_\Omega -\nabla^2u\;v\,dA = \int_\Omega \nabla u\cdot\nabla v\,dA - \oint_{\partial\Omega} v\,\frac{\partial u}{\partial n}\,ds.
The boundary integral vanishes because v=0v = 0 there, leaving the weak form
∫Ω∇u⋅∇v dA=∫Ωf v dAfor every admissible v,\int_\Omega \nabla u\cdot\nabla v\,dA = \int_\Omega f\,v\,dA \qquad\text{for every admissible } v,
exactly the bar's ∫EA u′v′=∫fv\int EA\,u'v' = \int fv with the derivative replaced by the gradient. Writing u=∑jujφju = \sum_j u_j\varphi_j and taking v=φiv = \varphi_i gives Ku=FK\mathbf u = \mathbf F with Kij=∫Ω∇φi⋅∇φj dAK_{ij} = \int_\Omega\nabla\varphi_i\cdot\nabla\varphi_j\,dA. Again only first derivatives appear, so functions that are merely continuous — planes glued along edges — are allowed.
💡The linear triangle
Courant's 1941 choice: on each triangle, uu is a plane, u=α+βx+γyu = \alpha + \beta x + \gamma y, fixed by its values at the three corners. Its shape functions N1,N2,N3N_1, N_2, N_3 are the planes that equal 1 at one corner and 0 at the other two, so
  • the gradient of uu is constant on each triangle — flux, stress or shear is one value per element, the 2-D version of the bar's constant stress;
  • neighbouring planes agree along their shared edge (a straight line is fixed by its two end values), so uu is continuous across the mesh;
  • the element stiffness is ke=Ae GTGk_e = A_e\,G^{\mathsf T}G, where AeA_e is the triangle's area and GG is the 2×32 \times 3 matrix whose columns are the three constant gradients.
Derivation
The triangle's gradients and stiffness
The gradients. Write each shape function as Ni=ai+bix+ciyN_i = a_i + b_ix + c_iy. The conditions '1 at corner ii, 0 at the others' for all three functions at once say that the matrix
M=[1x1y11x2y21x3y3]M = \begin{bmatrix} 1 & x_1 & y_1 \\ 1 & x_2 & y_2 \\ 1 & x_3 & y_3 \end{bmatrix}
times the coefficient matrix [a1a2a3b1b2b3c1c2c3]\begin{bmatrix} a_1 & a_2 & a_3 \\ b_1 & b_2 & b_3 \\ c_1 & c_2 & c_3 \end{bmatrix} is the identity. So the coefficients are M−1M^{-1}, and its last two rows are the gradients ∇Ni=(bi,ci)\nabla N_i = (b_i, c_i) — the matrix GG. (The lab computes it exactly this way. ∣det⁡M∣=2Ae|\det M| = 2A_e, so a nearly flat triangle makes MM nearly singular and the gradients huge: distorted elements are badly conditioned.)
The stiffness. The integrand ∇Ni⋅∇Nj\nabla N_i\cdot\nabla N_j is constant over the triangle, so the integral is that constant times the area: kij=Ae (bibj+cicj)k_{ij} = A_e\,(b_ib_j + c_ic_j), or ke=Ae GTGk_e = A_e\,G^{\mathsf T}G.
Two built-in checks. The three shape functions add up to 1 everywhere, so their gradients add up to zero: every row of kek_e sums to zero — a constant temperature stores no energy and drives no flux. And because a linear field has constant gradient, a mesh of linear triangles can represent any linear field exactly. That is the patch test: impose a linear field on the boundary of any mesh, however distorted, and a correct element must reproduce it everywhere to rounding error.
💡Torsion of a shaft: Prandtl's stress function
In 1903 Ludwig Prandtl showed that twisting a prismatic shaft of any cross-section reduces to a Poisson problem for a stress function ϕ(x,y)\phi(x, y) over the section: −∇2ϕ=2Gθ  in the section,ϕ=0  on its outline,-\nabla^2\phi = 2G\theta \ \text{ in the section}, \qquad \phi = 0 \ \text{ on its outline}, with GG the shear modulus and θ\theta the twist per unit length. The two shear stresses are the derivatives of ϕ\phi (their size is ∣∇ϕ∣|\nabla\phi|), and the torque the shaft carries is T=2∫ϕ dAT = 2\int\phi\,dA, which defines the torsion constant JJ through T=GJθT = GJ\theta. Prandtl's membrane analogy makes it visible: ϕ\phi has the shape of a soap film stretched over a hole of the section's outline and inflated by a uniform pressure — the steepest slopes of the film are where the shear stress is highest. For a square of side aa the exact values (from series) are J=0.14058 a4J = 0.14058\,a^4 and a peak shear of 0.6753 Gθa0.6753\,G\theta a at the middle of each side.
📐Sleipner's lesson in miniature
Torsion of a square shaft, solved with linear triangles for Prandtl's stress function. The torsion constant JJ — a global, energy-like quantity — is within 5% on a 8×88 \times 8 mesh and within 0.1% on 64×6464 \times 64. The peak shear stress at the middle of each side — the local quantity a designer checks — is underestimated by 28% on a 4×44 \times 4 mesh, 18% on 8×88 \times 8, 9% on 16×1616 \times 16, 4.6% on 32×3232 \times 32 and still 2.3% on 64×6464 \times 64. Every coarse answer looked plausible. Only the sequence of refinements shows how far from the truth each one was.
A map is not the territory it represents, but, if correct, it has a similar structure to the territory, which accounts for its usefulness.
– Alfred Korzybski, Science and Sanity (1933)
💡How fast does it converge — and where does it not?
For a smooth solution and linear triangles of size hh, the error in uu falls like h2h^2, but the error in its gradient — the stress, the flux, the shear — only like hh: halving the elements halves the stress error and quarters the displacement error, the same pattern as the bar. Energy-type quantities such as the torsion constant converge like h2h^2, because errors in the gradient enter the energy squared. Corners spoil this. At a re-entrant corner — one that points into the material, such as the inside corner of an L-shaped plate or the root of an unfilleted bracket — the exact solution behaves like r2/3r^{2/3} in the distance rr from the corner, so its gradient behaves like r−1/3r^{-1/3} and is infinite at the corner itself. The computed peak stress there never converges: each halving of hh raises it by about 21/3=1.262^{1/3} = 1.26. And the singularity slows convergence everywhere else as well. Real designs round such corners with fillets; analysts refine the mesh locally around them, or report stresses at a distance from the corner rather than at it.
🕰Courant's 1941 address
On 3 May 1941 Richard Courant gave an address to the American Mathematical Society on variational methods, published in 1943. To show what Ritz's method could do with simple local functions, he computed the torsion of a square shaft with a square hole using functions that were linear on each triangle of a mesh, on four successively finer meshes: the torsional stiffness came out 0.344, 0.352, 0.353, 0.353. It is the earliest finite element computation that historians of the method have found — done by hand, fifteen years before engineers at Boeing rediscovered the idea with computers.
⚠Checking a finite element model
Before believing any FE result: (1) check reactions sum to the applied loads; (2) compare with a hand calculation or a simpler model; (3) run the patch test on a new element or code; (4) refine the mesh, especially at corners, holes and load points, until the quantity of interest converges — stresses at re-entrant corners may never converge, because the exact stress is infinite there; (5) check element shapes — badly distorted elements degrade accuracy, as at Sleipner; (6) check units. Most famous FE failures were not software bugs but unchecked models.
💡Where you will use this
Structures and Mechanical design: every modern stress analysis, from brackets to bridges, is finite elements — triangles and quadrilaterals in 2-D, tetrahedra and bricks in 3-D, with element integrals evaluated at the Gauss points of lesson 8-4. Geotechnical Engineering: foundations, tunnels, slopes and dams with nonlinear soil models. Aircraft and Spacecraft structures: stress, buckling and modal analysis (Unit 6) of airframes and launch vehicles. Heat transfer, electromagnetics, CFD: the same machinery with different weak forms.
📚Further reading
  • Gilbert Strang and George Fix, An Analysis of the Finite Element Method (2nd ed., Wellesley–Cambridge, 2008) — convergence rates, the patch test and singularities, with proofs.
  • Thomas J. R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis (Dover, 2000), chapters 2–3 — two-dimensional elements, isoparametric mapping and Gauss quadrature.
  • Richard Courant, 'Variational methods for the solution of problems of equilibrium and vibrations', Bulletin of the American Mathematical Society 49, 1–23 (1943) — the address with the first triangle computation.
  • Bernt Jakobsen and Finn Rosendahl, 'The Sleipner platform accident', Structural Engineering International 4(3), 190–193 (1994).
↓But there is a problem
Everything in the course so far has answered the question given this model, what happens? — how the beam deflects, how the orbit evolves, how the water seeps. Engineering asks the opposite question: which design is best? The lightest bracket that does not yield, the fuel-optimal trajectory to orbit, the gain that minimises tracking error, the parameters that best fit a test. These are optimisation problems — and in 1847 Augustin-Louis Cauchy, studying the orbits of celestial bodies, described the simplest method for them: go downhill.