Finite Element Method Explained: Meshing, Solving & Results

Every modern aircraft, car crash simulation, and structural engineering design depends on FEM. It converts partial differential equations describing continuous physical domains into systems of algebraic equations solvable by computers — by dividing the domain into small elements and approximating the unknown field within each one.

1. From PDE to Algebraic System

Most engineering problems are governed by PDEs. For linear elasticity:

Governing PDE (equilibrium, no body forces): \nabla\cdotσ = 0 (divergence of stress tensor = 0) σ = C : ε (constitutive law: stress = elastic tensor \times strain) ε = ½(\nablau + (\nablau)ᵀ) (strain-displacement relation) Where: u = displacement field (what we solve for) σ = stress tensor (6 independent components in 3D) ε = strain tensor C = 4th-order elasticity tensor (reduced to 2 constants E, ν for isotropic materials) Analytical solutions exist only for simple geometries. FEM solves the weak (variational) form: \int_Ω δε : σ dV = \int_Γ δu \cdot t dA + \int_Ω δu \cdot b dV δu = virtual displacement (test function) t = surface tractions (forces per area) b = body forces (gravity, etc.)

2. Mesh Generation & Quality

The domain is divided into non-overlapping elements. Element types in 2D: tria3/tria6 (triangles with 3 or 6 nodes), quad4/quad8 (quadrilaterals). In 3D: tet4/tet10 (tetrahedra), hex8/hex20 (bricks), wedge6/wedge15 (prisms).

Key quality metrics:

Aspect ratio: Ratio of longest to shortest edge. High aspect ratio (very elongated elements) degrades accuracy. Target < 5 for most elements, < 10 only in boundary layers.
Jacobian ratio: Ratio of maximum to minimum Jacobian det within element. Must be > 0 everywhere (no inverted elements). Collapsed below 0.6 causes poor accuracy.
Angular quality: For triangles/tets, very small or very large angles reduce accuracy. Ideal triangle: equilateral. Delaunay triangulation maximises the minimum angle.
Element size transition: Gradual size transitions (growth factor < 1.5) between fine regions (e.g., stress concentrations) and coarse regions.

Boundary layer meshing (CFD): In fluid dynamics, the velocity changes from zero (wall) to freestream over a thin boundary layer. This requires prismatic layers of very thin, high-aspect-ratio elements near walls — typically 20+ layers with y+ ≈ 1 for direct wall resolution in turbulence modelling. Generating these automatically from complex CAD geometry is among the most challenging mesh generation tasks.

3. Shape Functions

Within each element, the unknown field u is approximated as a weighted sum of shape functions N_i(x), where the weights are the nodal values u_i:

u(x) \approx Σᵢ Nᵢ(x) \cdot uᵢ (within one element) Shape functions for a linear triangle (CST — constant strain): L₁, L₂, L₃ = area coordinates (barycentric) N₁ = L₁, N₂ = L₂, N₃ = L₃ Properties: Nᵢ(xⱼ) = δᵢⱼ (1 at node i, 0 at all other nodes) Σᵢ Nᵢ = 1 (partition of unity) Serendipity quad (Q8) shape functions at a corner node: N₁(ξ,η) = ¼(1+ξ₀)(1+η₀)(ξ₀+η₀-1) where ξ₀ = ξ\cdotξ₁, η₀ = η\cdotη₁ (natural coordinates -1 to +1) Higher-order shape functions (p-refinement): Increasing polynomial degree improves accuracy without smaller elements. p-FEM can achieve exponential convergence for smooth problems.

4. Element Stiffness & Assembly

For each element, the element stiffness matrix K_e is computed by substituting shape functions into the weak form:

Element stiffness matrix: K_e = \int_Ωe Bᵀ C B dV B = strain-displacement matrix = \partialN (derivatives of shape functions) C = elasticity matrix (E, ν constants) Numerical integration (Gauss quadrature): \int_Ωe f dV \approx Σₖ wₖ \cdot f(xₖ) \cdot det(J(xₖ)) xₖ = Gauss points, wₖ = weights, J = Jacobian of mapping Full FEM system: K \cdot u = F K = global stiffness matrix (assembled from Kₑ) u = nodal displacement vector (unknowns) F = force vector (applied loads) For a model with N nodes and 3 DOF/node: Matrix size: 3N \times 3N Typical mesh: 1,000-10,000,000 nodes K is sparse (only connected elements are non-zero) Sparse solvers (PARDISO, MUMPS) required for large problems

5. Boundary Conditions & Solving

Boundary conditions specify known quantities at the domain boundary. For structural analysis:

Dirichlet BC (essential): Prescribed displacements (e.g., fixed support: u = 0). Applied by modifying the stiffness matrix rows/columns.
Neumann BC (natural): Prescribed forces/tractions. Appear naturally in the force vector F.
Contact conditions: Inequality constraints (surfaces cannot interpenetrate). Require iterative solution methods (penalty, Lagrange multiplier, or mortar contact).

For linear problems, K·u = F is solved once using direct (LU decomposition, Cholesky) or iterative (conjugate gradient, GMRES) methods. For nonlinear problems (large deformation, plasticity, contact), Newton-Raphson iteration solves repeatedly until residual < tolerance.

6. Convergence & Error Estimation

FEM produces approximate solutions. Accuracy improves with:

h-refinement: Reducing element size h. For linear elements, displacement error ∝ h², stress error ∝ h.
p-refinement: Increasing polynomial order p. Exponential convergence for smooth problems.
hp-refinement: Combining both. Most efficient for problems with localised singularities (crack tips, reentrant corners).

A-posteriori error estimate (Zienkiewicz-Zhu): η_e = || σ_h* - σ_h ||_Ωe / || σ_h ||_Ω σ_h* = smoothed stress (recovered using SPR or PPR) σ_h = raw FEM stress (discontinuous across elements) η = global error estimator (sum over elements) Target: η < 5-10% for engineering accuracy Use error estimator to drive adaptive mesh refinement: Refine elements where local error > threshold Coarsen elements where local error << threshold

7. Multiphysics Applications

Crash simulation (nonlinear transient FEM): A car crash lasts ~100 ms. Explicit time integration with time steps of ~1 μs. 10–50 million elements. Accounts for large deformation, plasticity, fracture, contact, and airbag inflation. Tools: Abaqus/Explicit, LS-DYNA, PAM-CRASH.
Structural fatigue: FEM provides stress fields for S-N curve and fracture mechanics analysis. Rainflow cycle counting + Miner's rule estimates component life under variable amplitude loading.
Thermal-structural coupling: FEM solves heat transfer (conduction, convection) to get temperature field, then uses thermal strains in structural analysis. Critical for turbine blade design (Δ T = 800°C → thermal stress dominates mechanical load).
Electromagnetic FEM: Solves Maxwell's equations in frequency domain. Used for antenna design, motor performance, transformer core loss, and MRI coil design. Tools: COMSOL, Ansys HFSS/Maxwell.
Biomedical: Patient-specific bone fracture risk from CT scans, stent mechanical design, knee implant contact stresses, cochlear implant fluid-structure interaction.
Geomechanics: Reservoir compaction, fault stability under pore pressure changes, tunnel lining design, dam safety. Coupled hydro-mechanical FEM simulates fluid flow in porous deformable media.