How the Finite Element Method Works

Jump to Examples

In Structural Statics I calculated beams, trusses and columns by hand, and at the end a continuous beam with the stiffness method in Python. Real structures are rarely made only of beams: slabs, walls, openings, a staircase, a strange corner where three things meet. For all of that, almost every structural program today uses the finite element method – FEM. The more I learn about it, the more I think the most important skill is not clicking through the program but knowing what it is really doing.

The colourful pictures make it look easy. They are also the reason why FEM is sometimes called "the method that turns wrong input into beautiful results". This note explains what happens inside, in plain words and with one small program, and lists the most common mistakes.

A Bit of History

The idea of cutting a complicated body into simple pieces is old. In 1941 Alexander Hrennikoff modelled a plate as a lattice of small bars, and in 1943 the mathematician Richard Courant solved a torsion problem by splitting the area into triangles with a simple function on each one – the mathematical core of FEM, although nobody had computers to use it yet.

The practical push came from aircraft. Swept and delta wings of the 1950s were too complicated for the old beam formulas. At Boeing, M. J. Turner, Ray Clough, H. C. Martin and L. J. Topp published in 1956 how to assemble a wing from small triangular panels, each with its own stiffness matrix. At the same time John Argyris in London (later Stuttgart) developed the matrix methods from energy principles. In 1960 Clough, by then a professor at Berkeley, gave the method its name: the finite element method.

The rest is the story of computers. Olgierd Zienkiewicz wrote the first textbook in 1967, NASA's general program NASTRAN appeared at the end of the 1960s, and Ed Wilson's SAP from Berkeley around 1970 – the ancestor of many structural programs still in use. With personal computers in the 1980s FEM moved to every design office, and today a laptop solves a model with a million unknowns in minutes.

The Idea: Many Small Springs

Start with the simplest structural element there is: a spring. Pull it with a force F and it stretches by u, with F = k · u, where k is the stiffness. A steel bar under tension is exactly such a spring, with k = E · A / L.

Now connect several springs. Each connection point is a node, and each node can move. For every node you can write one equilibrium equation: the forces from all springs connected to it must balance the external load. Put all the equations together and you get one matrix equation:

K · u = F

K ... stiffness matrix of the whole structure (how strongly it resists each movement)
u ... the unknown displacements of all nodes
F ... the loads at the nodes

That is the whole method. A real structure is cut into many small pieces – the elements – and each element behaves like a slightly more complicated spring, described by its own small stiffness matrix. The computer adds all of them into the big matrix K (the assembly), removes the movements blocked by the supports, solves the system of equations, and then goes back into each element to calculate strains and stresses from the displacements of its nodes.

The trick that makes it work: inside an element we assume a simple shape of the displacement – linear, or a parabola – described by the displacements of its nodes (the shape functions). The real displacement in a slab or a wall is not linear, but if the elements are small enough, many small straight pieces follow any curve well. That is also why FEM is always an approximation, and why the mesh matters.

Types of Elements

The rule of thumb is to use the simplest element that still describes the behaviour. A building frame is modelled with beams and shells, not with millions of solid elements: the result is the same, it is faster, and – most importantly – you get internal forces (moments, shear forces) that you can check and design with. Solid elements give stresses, which are hard to check by hand.

Elements also come in linear and quadratic versions. Linear elements have nodes only at the corners and assume straight-line displacements; quadratic ones have extra nodes in the middle of the edges and can bend. In bending, linear elements are too stiff – the effect has a name, shear locking – and Example 2 shows how much it matters.

The Mesh and Convergence

Because each element can only follow a simple shape, a coarse mesh gives a structure that is too stiff: smaller deflections and usually smaller stresses than in reality. As the mesh gets finer, the results move towards the true solution – the model converges. Two practical consequences:

• Make a mesh study. Calculate the model with a mesh, then with a mesh about twice as fine. If the results you care about change by only a few percent, the mesh is fine enough. If they change a lot, keep refining.
• Displacements converge faster than stresses. Stresses come from the derivatives of displacements, and a derivative of an approximation is always rougher than the approximation itself. A mesh that gives good deflections may still give poor stresses.

And a warning about the beautiful colours: programs average the stresses at nodes and draw smooth colour maps, which hide how rough the real results are. Look at the unaveraged results once in a while – big jumps between neighbouring elements mean the mesh is too coarse there.

Singularities: When Stress Goes to Infinity

In some places the theory of elasticity gives an infinite stress: under a point load, at a point support, at a sharp inside corner of an opening, at the edge of a fully fixed support. A real structure is never like that – the load is spread over a bearing plate, the corner has a radius, concrete cracks a little, steel yields. But the model does exactly what you told it. The finer the mesh, the higher the peak; it never converges.

So a red spot at a corner or under a support is often not a real problem but a singularity. Recognise it by refining the mesh: if the peak keeps growing while the results a little further away stay the same, it is a singularity. Design from the results away from it, or from integrated values – the total force in a section, the moment per metre of a slab – which converge even when the peak does not. Example 2 shows this very clearly.

The Classic Mistakes

• Supports. The most common source of wrong results. A fully fixed support where in reality there is a pin makes the structure much stiffer; a missing support makes it a mechanism (the program complains about a "singular matrix" or gives huge displacements). Think about what the real connection can do.
• Units. Millimetres and metres, N and kN, MPa and kN/m² in the same model. One mixed unit and the result is off by a factor of 1000. FreeCAD, for example, works in mm, N and MPa.
• Loads. Self-weight forgotten, or counted twice (once by the program, once by hand). Loads in the wrong direction. Point loads where a distributed load would be correct (see singularities).
• Element choice. Linear tetrahedra in bending are far too stiff. Use quadratic elements, or enough elements through the thickness.
• Trusting the colours. A model that runs without errors is not a correct model. The program checks the mathematics, not the engineering.
• Forgetting what is not in the model. A linear elastic model knows nothing about cracking of concrete, yielding of steel, buckling (unless you ask for it), or the creep of a slab over the years.

Examples

Example 1: A Stiffness Matrix by Hand

A steel rod hangs from the ceiling. The upper part (2 m) has a cross-section of 4 cm², the lower part (2 m) 2 cm². At the bottom it carries 20 kN. Three nodes (1 at the ceiling, 2 in the middle, 3 at the bottom), two elements, one unknown per node: the vertical displacement.

k1 = E·A1/L1 = 210 000 000 kN/m² · 0.0004 m² / 2 m = 42 000 kN/m
k2 = E·A2/L2 = 210 000 000 kN/m² · 0.0002 m² / 2 m = 21 000 kN/m

Each element adds its stiffness  [ k  -k ]  to the rows and columns of its two nodes:
                                 [ -k  k ]

        node 1     node 2       node 3
K = [   k1        -k1           0     ]   [ 42 000   -42 000        0 ]
    [  -k1       k1 + k2       -k2    ] = [-42 000    63 000  -21 000 ]
    [   0         -k2           k2    ]   [      0   -21 000   21 000 ]

Node 1 is fixed (u1 = 0), so remove row and column 1:

[ 63 000  -21 000 ] [u2]   [  0 ]
[-21 000   21 000 ] [u3] = [ 20 ]

u2 = 0.476 mm,  u3 = 1.429 mm

Check by hand: each part stretches by P/k – 20/42 000 = 0.476 mm and 20/21 000 = 0.952 mm – and together 1.429 mm. Exactly the same. The reaction at the ceiling comes from the removed first row: −42 000 · 0.000476 = −20 kN, so the support pulls up with 20 kN. A real program does exactly this, only with millions of rows.

Example 2: A Cantilever from Small Elements, in Python

Now a real 2D problem: a cantilever 2 m long, 20 cm high and 10 cm thick, made of a concrete-like material (E = 30 000 MPa, ν = 0.2), fixed on the left, with 10 kN at the free end. Beam theory gives a deflection of P·L³ / (3·E·I) = 13.33 mm (plus about 0.1 mm from shear), a bending stress of 30 MPa at the support and 15 MPa at mid-span.

The program below models it with four-node plane elements (Q4) – a little rectangle that knows only two things: its corners can move horizontally and vertically, and inside it everything changes linearly. It is the same six steps as in Structural Statics: mesh, element matrices, assembly, loads, supports, solve.

import numpy as np

# A cantilever made of 4-node plane elements (Q4, plane stress)
L, h, t = 2.0, 0.2, 0.1          # length, height, thickness [m]
E, nu = 30e6, 0.2                # concrete-like material, E in kN/m2 (= 30 000 MPa)
P = 10.0                         # load at the free end [kN], downwards
D = E / (1 - nu**2) * np.array([[1, nu, 0], [nu, 1, 0], [0, 0, (1 - nu) / 2]])

def B_matrix(xy, xi, eta):
    """Strain-displacement matrix of a Q4 element at the point (xi, eta)."""
    dN = 0.25 * np.array([[-(1 - eta), 1 - eta, 1 + eta, -(1 + eta)],
                          [-(1 - xi), -(1 + xi), 1 + xi, 1 - xi]])
    J = dN @ xy                                   # Jacobian: element -> real size
    dNx = np.linalg.solve(J, dN)
    B = np.zeros((3, 8))
    B[0, 0::2] = dNx[0]; B[1, 1::2] = dNx[1]
    B[2, 0::2] = dNx[1]; B[2, 1::2] = dNx[0]
    return B, np.linalg.det(J)

def solve(nx, ny):
    # 1. mesh: nodes on a regular grid, node number = i*(ny+1) + j
    node = lambda i, j: i * (ny + 1) + j
    xy_all = np.array([[x, y] for x in np.linspace(0, L, nx + 1)
                              for y in np.linspace(-h / 2, h / 2, ny + 1)])
    n_dof = 2 * len(xy_all)
    K, F = np.zeros((n_dof, n_dof)), np.zeros(n_dof)
    elements = {}
    # 2. element stiffness matrices, assembled into the global matrix
    for i in range(nx):
        for j in range(ny):
            en = [node(i, j), node(i + 1, j), node(i + 1, j + 1), node(i, j + 1)]
            dofs = np.ravel([[2 * n, 2 * n + 1] for n in en])
            ke = np.zeros((8, 8))
            for xi in (-0.57735, 0.57735):        # 2 x 2 Gauss integration
                for eta in (-0.57735, 0.57735):
                    B, detJ = B_matrix(xy_all[en], xi, eta)
                    ke += B.T @ D @ B * detJ * t
            K[np.ix_(dofs, dofs)] += ke
            elements[(i, j)] = (xy_all[en], dofs)
    # 3. load: P shared between the nodes at the free end
    for j in range(ny + 1):
        F[2 * node(nx, j) + 1] -= P / ny * (0.5 if j in (0, ny) else 1.0)
    # 4. supports: all nodes at x = 0 fixed in both directions
    fixed = [d for j in range(ny + 1) for d in (2 * node(0, j), 2 * node(0, j) + 1)]
    free = np.setdiff1d(np.arange(n_dof), fixed)
    # 5. solve K u = F
    u = np.zeros(n_dof)
    u[free] = np.linalg.solve(K[np.ix_(free, free)], F[free])
    # 6. results: deflection, and stress sigma_x on the top edge
    tip = -u[2 * node(nx, ny // 2) + 1]
    def sigma_x(i, xi):                           # top edge of element column i
        xy, dofs = elements[(i, ny - 1)]
        B, _ = B_matrix(xy, xi, 1.0)
        return (D @ B @ u[dofs])[0] / 1000        # kN/m2 -> MPa
    return tip, sigma_x(nx // 2, -1.0), sigma_x(0, -1.0), len(free)

I = t * h**3 / 12
print(f"beam theory: deflection {P * L**3 / (3 * E * I) * 1000:.2f} mm (+ about 0.1 mm from shear)")
print(f"beam theory: stress at mid-span {P * L / 2 / (t * h**2 / 6) / 1000:.1f} MPa, at the support {P * L / (t * h**2 / 6) / 1000:.1f} MPa")
for nx, ny in [(4, 1), (10, 2), (20, 2), (40, 4), (80, 8), (160, 16)]:
    tip, s_mid, s_corner, n = solve(nx, ny)
    print(f"mesh {nx:3d} x {ny:<2d} ({n:5d} unknowns): deflection {tip * 1000:5.2f} mm, "
          f"stress mid-span {s_mid:5.1f} MPa, at the fixed corner {s_corner:5.1f} MPa")

Output:

beam theory: deflection 13.33 mm (+ about 0.1 mm from shear)
beam theory: stress at mid-span 15.0 MPa, at the support 30.0 MPa
mesh   4 x 1  (   16 unknowns): deflection  3.68 mm, stress mid-span   3.2 MPa, at the fixed corner   7.5 MPa
mesh  10 x 2  (   60 unknowns): deflection  9.38 mm, stress mid-span   9.6 MPa, at the fixed corner  20.6 MPa
mesh  20 x 2  (  120 unknowns): deflection 12.02 mm, stress mid-span  13.0 MPa, at the fixed corner  27.0 MPa
mesh  40 x 4  (  400 unknowns): deflection 13.03 mm, stress mid-span  14.4 MPa, at the fixed corner  30.2 MPa
mesh  80 x 8  ( 1440 unknowns): deflection 13.32 mm, stress mid-span  14.8 MPa, at the fixed corner  32.2 MPa
mesh 160 x 16 ( 5440 unknowns): deflection 13.39 mm, stress mid-span  14.9 MPa, at the fixed corner  34.5 MPa

Three lessons in one table:

• A coarse mesh is dangerously stiff. With 4 elements in a row the cantilever bends only 3.7 mm instead of 13.4 – almost four times too little. Linear elements cannot bend, so they resist with fake shear stiffness (shear locking). A program would happily show this result in nice colours.
• The deflection converges to 13.39 mm – beam theory with the shear part gives 13.43 mm – and the stress at mid-span converges to 15 MPa. Away from the support, the model and the hand calculation agree.
• The stress at the fixed corner does not converge. It passes the 30 MPa from beam theory and keeps growing with every refinement: 30.2, 32.2, 34.5 MPa… That is a singularity – the edge is fixed completely, including the sideways movement that the Poisson effect wants. In a real structure that peak would not exist in that form; for design, use the bending moment at the support (10 kN × 2 m = 20 kNm), not the peak stress.

The program uses a full (dense) matrix to stay short, so it gets slow beyond a few thousand unknowns. Real programs store only the non-zero numbers (sparse matrices) and solve millions of unknowns.

Example 3: The Same Cantilever in FreeCAD

The same check with a real program. The names of the buttons change a little between FreeCAD versions; these are from FreeCAD 1.0:

1. Geometry. In the Part workbench make a box (Cube): length 2000 mm, width 100 mm, height 200 mm. FreeCAD works in millimetres, newtons and MPa.
2. Analysis. Switch to the FEM workbench and create an Analysis container – everything else goes into it.
3. Material. Add a solid material and set Young's modulus 30 000 MPa and Poisson's ratio 0.2.
4. Support. Add a Fixed boundary condition and select the end face at x = 0.
5. Load. Add a Force load on the other end face: 10 000 N, direction downwards.
6. Mesh. Create a mesh with Gmsh. Choose second-order elements (quadratic tetrahedra) and a maximum element size of about 20 mm. With first-order tetrahedra you will see shear locking again.
7. Solve. Add the CalculiX solver, write the input file and run it.
8. Check. The vertical displacement at the free end should be close to the 13.4 mm from Example 2. The stress near the fixed edge will show a red peak above 30 MPa – now you know why.

Then repeat with a coarser and a finer mesh, and with first-order elements, and watch how the results move. Fifteen minutes of that teaches more about FEM than a semester of theory.

My Checklist for Every Model

• Sum of the reactions = sum of the loads. The first thing I look at. If they differ, a load is missing, doubled or in the wrong direction.
• The deformed shape makes sense. Exaggerate it and look: does the structure move the way you would expect? Do the supports hold where they should?
• A hand calculation for the order of magnitude. A simple beam or a strip of slab. If FEM and the hand result differ by a factor of two, one of them is wrong – and it is not always the hand calculation.
• A mesh study for the results you design with.
• Peaks at corners, point loads and supports are checked as possible singularities before anybody panics.
• Units checked once more.

A Few More Thoughts

The finite element method is one of the great tools of engineering: the same idea – many small springs and one big system of equations – calculates wings, bridges, dams, car crashes, the flow of heat through a wall and the vibration of a guitar body. Free programs like FreeCAD, CalculiX, OpenSees or Code_Aster can do things that cost a fortune twenty years ago.

But the method has not changed one thing: the computer solves the model you give it, not the structure you meant. Good FEM users first draw the result by hand – roughly how it will bend, where the moments will be – and only then look at the colours. If the picture surprises them, they look for the mistake before they believe the result. That habit is worth more than any software licence.