Finite Element Method (FEM) Explained Simply

Engineers design bridges, jet engines, and biomedical implants on computers before building them — all thanks to FEM. It turns impossible-to-solve differential equations into large but solvable systems of linear equations by chopping a domain into small pieces and approximating the answer on each piece.

1. The Big Idea

Most physics problems boil down to partial differential equations (PDEs): stress in a beam, heat through a wall, airflow around a wing. These equations rarely have closed-form solutions for complex geometries. FEM solves them numerically by:

Discretize

Split the domain into small elements (triangles, tetrahedra)

Approximate

Assume the solution varies polynomially within each element

Assemble

Combine all elements into a global system of equations

Solve

Use linear algebra to find unknown values at every node

More elements means better accuracy but more computation. A modern automotive crash simulation uses 5–10 million elements and takes hours on a cluster of GPUs.

2. Strong Form → Weak Form

Consider a 1D bar under axial load. The strong (differential) form of equilibrium is:

-d/dx [EA \cdot du/dx] = f(x) for x \in [0, L] E = Young's modulus, A = cross-section area, u = displacement, f = distributed load

The strong form requires u to be twice differentiable everywhere — too strict for numerical methods. We multiply by a test function v(x) and integrate by parts to get the weak form:

\int₀ᴸ EA \cdot (du/dx) \cdot (dv/dx) dx = \int₀ᴸ f\cdotv dx + [EA\cdot(du/dx)\cdotv]₀ᴸ This reduces the smoothness requirement: u and v only need first derivatives (C⁰ continuity), not second.

The weak form is the mathematical foundation of FEM. It says: find u such that the equation holds for all admissible test functions v. In FEM, both u and v are approximated by the same polynomial shape functions — the Galerkin method.

3. Meshing the Domain

The domain (a beam, a turbine blade, a skull) is divided into non-overlapping elements:

1D: Line segments (bars, beams)
2D: Triangles or quadrilaterals (plates, shells)
3D: Tetrahedra, hexahedra, prisms, or pyramids

Mesh quality matters enormously. Elements with extreme aspect ratios or very obtuse angles degrade accuracy. Adaptive meshing refines the mesh where gradients are large (e.g., around a crack tip or hole) and coarsens it where the field varies slowly.

Mesh convergence: As elements shrink (h → 0), the FEM solution converges to the exact PDE solution. The convergence rate depends on element type: linear elements converge O(h²) in displacement, quadratic elements converge O(h³). Always run a mesh convergence study before trusting results.

4. Shape Functions & Interpolation

Within each element, the unknown field u is approximated as a weighted sum of shape functions N_i:

u(x) = Σᵢ Nᵢ(x) \cdot uᵢ where uᵢ are the nodal values (unknowns) and Nᵢ(x) are polynomial shape functions: - Linear (2-node): N₁ = (1-ξ)/2, N₂ = (1+ξ)/2 (ξ \in [-1,1]) - Quadratic (3-node): adds a mid-node - Isoparametric: same functions map geometry and field

Shape functions have the partition of unity property: Σ Nᵢ = 1 everywhere. Each Nᵢ equals 1 at its own node and 0 at all other nodes. This ensures continuity across element boundaries (C⁰ for Lagrange elements).

5. The Stiffness Matrix

Substituting the shape function approximation into the weak form produces, for each element, a local equation:

kᵉ \cdot uᵉ = fᵉ Element stiffness matrix: kᵉᵢⱼ = \int EA \cdot (dNᵢ/dx) \cdot (dNⱼ/dx) dx Element force vector: fᵉᵢ = \int f \cdot Nᵢ dx + boundary terms

For a 2-node 1D linear element of length Lᵉ:

kᵉ = (EA/Lᵉ) \cdot [ 1 -1 ] [-1 1 ] This is the "spring stiffness" matrix — the basic building block.

For 2D and 3D elements, the integrals are evaluated numerically using Gaussian quadrature — summing the integrand at specific sample points with specific weights. A quadrilateral element typically uses 2×2 Gauss points; a hexahedron uses 2×2×2 = 8 Gauss points.

6. Assembly & Solving

The global stiffness matrix K is assembled by mapping each element's local matrix into the global DOF numbering using a connectivity table. The result is a large, sparse, symmetric positive-definite system:

K \cdot U = F K = N\timesN global stiffness matrix (N = total DOFs) U = vector of all unknown nodal displacements F = vector of applied forces and boundary conditions For a million unknowns: K has ~10⁶ rows but only ~50 non-zeros per row \to sparse matrix storage, iterative solvers (CG, GMRES)

After applying boundary conditions (e.g., u=0 at supports), the system is solved. Post-processing computes derived quantities: stress σ = E·ε, strain ε = du/dx, von Mises equivalent stress for failure prediction.

7. FEM in the Real World

Structural analysis: Every car body, aircraft fuselage, and bridge is FEM-analysed for static loads, fatigue, and crash impact. ANSYS, Abaqus, NASTRAN are industry standard.
Heat transfer: FEM solves the heat equation to predict thermal gradients in electronics, furnaces, and re-entry heat shields.
Fluid dynamics (CFD): FEM variants (SUPG, DG) solve Navier-Stokes for complex geometries. OpenFOAM (finite volume) and COMSOL (FEM) are widely used.
Electromagnetics: FEM solves Maxwell's equations for antenna design, motor optimisation, and MRI coil design.
Biomedical: Patient-specific bone stress analysis from CT scans (mesh generated directly from voxel data), heart valve simulation, dental implant design.

FEM in the browser: JavaScript libraries like FEA.js and Three.js can visualise FEM meshes and results. WebAssembly makes it feasible to run small FEM solvers entirely client-side — our own stress analysis simulation uses this approach.