67  Computational methods for engineering models

Meshes, weak forms, and conservation in computation

The governing equations may be right and still produce a useless simulation. A coarse mesh smears a boundary layer. A bad time step destabilises a model. A solver converges to something, but not to the physics you meant. Upper-year engineering lives in that gap between the equation on paper and the model you actually trust.

This chapter is about the mathematical move that closes that gap. A continuum equation is not yet computable. A computer needs finitely many unknowns, finitely many equations, and an update rule or linear system it can actually solve. Discretisation supplies that finite structure.

Every discretisation keeps some truths and sacrifices others. One method preserves conservation neatly. Another handles curved geometry more naturally. Another is easy to implement but unstable for large time steps. Knowing which sacrifices the physics can tolerate is what makes a simulation trustworthy.

67.1 What this chapter helps you do

Symbols to keep handy

These are the bits of notation you'll see a lot. If a line of symbols feels like a fence, read it out loud once, then keep going.

  • _i: phi sub i — a basis or shape function

  • r: r — the stability ratio, D delta-t divided by (delta-x) squared

  • _: integral over omega — an integral over the domain

  • A = : A u equals b — the algebraic system after discretisation

  • u_i^n: u i n — the approximation at grid point i and time step n

Definitions to keep handy

These are the words we keep coming back to. If one feels slippery, come back here and steady it before you push on.

  • discretisation: Turning a continuous model into a finite computation (a grid, a mesh, and a finite set of unknowns).

  • mesh / grid: A geometric partition of the domain that the computer can work with.

  • weak form: A rewritten form of the PDE designed to behave well under approximation (especially in finite elements).

  • stencil: The local pattern of grid points used by a finite difference formula.

  • sparse system: A big linear system with mostly zeros, which needs special storage and solvers.

Here is the main move this chapter is making, in plain terms. You do not need to be fast. You just need to keep the thread.

  • Coming in: Governing equations become useful in upper-year engineering only when they can be computed on real geometries and real data.

  • Leaving with: Discretisation is not bookkeeping. It is the mathematical act that turns a continuum model into a trustworthy simulation.

67.2 From equation to algebraic system

Discretisation is translation

If you have not seen this before, think of it as translating from “continuous language” into “spreadsheet language.”

  • a field becomes a table of values
  • derivatives become differences, flux balances, or basis-function integrals
  • boundary conditions become constraints on a finite set of unknowns

The point is not to throw away calculus. The point is to keep the physics while making the problem computable.

Take the one-dimensional diffusion equation

\frac{\partial u}{\partial t} = D\frac{\partial^2 u}{\partial x^2}

defined on a grid with spatial spacing \Delta x and time step \Delta t. The continuous field u(x,t) becomes a table of approximate values u_i^n, read “u at grid point i and time level n.”

Here the indices are doing the organisational work: i counts where you are on the grid, and n counts when you are in the time-marching loop.

Replace \partial u/\partial t with a forward difference (u_i^{n+1}-u_i^n)/\Delta t and \partial^2 u/\partial x^2 with a central difference (u_{i+1}^n - 2u_i^n + u_{i-1}^n)/(\Delta x)^2:

\frac{u_i^{n+1} - u_i^n}{\Delta t} \approx D\frac{u_{i+1}^n - 2u_i^n + u_{i-1}^n}{(\Delta x)^2}

Rearrange:

u_i^{n+1} = u_i^n + \frac{D\Delta t}{(\Delta x)^2} \left(u_{i+1}^n - 2u_i^n + u_{i-1}^n\right)

The PDE has become a recurrence. Once the spatial grid is fixed, the whole problem is now algebra.

This is the central translation of computational engineering:

  • fields become arrays
  • derivatives become finite differences, flux balances, or basis-function integrals
  • boundary conditions become algebraic constraints
  • the PDE becomes either a time-marching update or a linear/nonlinear system to solve

Use the selector to see how the same diffusion law changes shape as it is translated into three different computational languages.

67.3 Three dominant viewpoints

Upper-year engineering simulation is often organised around three families of discretisation.

Finite difference

Replace derivatives directly with local difference quotients on a grid. Conceptually simple. Especially natural on structured meshes and textbook geometries.

Finite volume

Integrate the balance law over control volumes, then approximate fluxes at faces. Especially attractive when conservation matters strongly, as in fluid flow, transport, and environmental modelling.

Finite element

Represent the solution using basis functions \phi_i over elements and solve a weak or variational form. The weak form replaces the requirement that the PDE holds exactly at every point with the requirement that a weighted integral of the residual is zero. This relaxed condition is easier to satisfy on irregular meshes, which makes finite element methods especially useful on complicated geometries and in structural mechanics.

Weak form, in words

The “strong form” says: the PDE must be true at every point.

The weak form says: if you test the PDE against a family of weight functions and integrate, the error averages to zero.

That sounds abstract, but it is the practical reason finite elements work on curved geometry: the method enforces the governing law in an integral sense that is stable under approximation.

These are not three unrelated subjects. They are three ways of answering the same question: how should a continuum law be turned into finite equations?

67.4 Stability, consistency, and convergence

A discretisation is not good merely because it can be coded.

Three recurring tests matter:

  • consistency: does the discrete equation approach the continuous one as the mesh is refined?
  • stability: do errors remain controlled instead of blowing up?
  • convergence: does the computed solution approach the true solution as the grid is refined?

Students often meet these separately. In practice they are linked by the Lax equivalence theorem. For a linear, well-posed PDE — one where small changes in the input produce small changes in the solution — and a consistent discretisation, stability is both necessary and sufficient for convergence. A consistent method that is unstable may be useless. A stable method that is too diffusive may hide the physics you care about.

For the explicit diffusion update above, the ratio

r = \frac{D\Delta t}{(\Delta x)^2}

controls stability. For the 1D explicit diffusion scheme the stability bound is

r = \frac{D\Delta t}{(\Delta x)^2} \leq \frac{1}{2}

This means mesh spacing and timestep cannot be chosen independently. Refining space without tightening time can destroy stability.

If r exceeds \frac{1}{2}, the numerical solution can oscillate or blow up even when the physical problem is perfectly well behaved. This bound means halving the grid spacing forces the time step to shrink by a factor of four — a practical cost of explicit time-stepping.

Implicit methods remove or relax this restriction at the cost of solving a linear system at each time step. Backward Euler is unconditionally stable (any \Delta t is permitted) but introduces more numerical diffusion than a comparable explicit step. Crank-Nicolson averages the current and future spatial operators, recovering second-order accuracy in time while remaining stable for all step sizes. For fine spatial grids or stiff problems, solving a system at each step is almost always the correct exchange.

NoteDiscretisation always introduces error

Replacing a continuous field with finitely many unknowns introduces error automatically. The right question is never “is there error?” but “what kind of error, how large, and does it destroy the physical meaning of the model?”

Computational mathematics is therefore not separate from modelling. The numerical method is part of the model.

67.5 The core method

A first pass through a computational engineering model usually goes like this:

  1. Write the governing equation and the boundary/initial conditions.
  2. Decide which discrete viewpoint matches the physics and geometry best.
  3. Define the grid, mesh, or basis functions.
  4. Derive the algebraic update or system.
  5. Check the main stability or conditioning constraint.
  6. Interpret the numerical output in physical terms and test mesh sensitivity.

The last step matters as much as the derivation. If halving \Delta x changes the answer substantially, the previous mesh was telling you more about the discretisation than about the system.

67.6 Worked example 1: explicit diffusion update

Take

\frac{\partial u}{\partial t} = D\frac{\partial^2 u}{\partial x^2}

with D = 1, \Delta x = 0.1, and \Delta t = 0.002.

Then

r = \frac{D\Delta t}{(\Delta x)^2} = \frac{1(0.002)}{0.01} = 0.2

For the standard explicit 1D diffusion scheme, the stability requirement is r \leq \tfrac{1}{2}. Since 0.2 < 0.5, this time step is within the stable range and is plausibly acceptable for this grid.

Now change only the time step to \Delta t = 0.008. Then

r = \frac{0.008}{0.01} = 0.8

and the scheme is no longer stable — 0.8 > 0.5. Same PDE. Same material parameter. Different numerical behaviour because the discretisation changed.

This is why simulation judgement is mathematical judgement. The equation on the page did not fail. The computational choice did.

67.7 Worked example 2: a finite-volume mass balance

Suppose a pollutant concentration is tracked over a control volume of width \Delta x. The half-integer subscripts label the cell faces: F_{i-1/2} is the flux at the left face of cell i (between cells i-1 and i), and F_{i+1/2} is the flux at the right face (between cells i and i+1). Then a finite-volume update has the form

\frac{c_i^{n+1} - c_i^n}{\Delta t} = -\frac{F_{i+1/2}^n - F_{i-1/2}^n}{\Delta x} + s_i^n

This is the discrete balance law in its most honest form:

  • future storage minus current storage
  • equals net flux imbalance
  • plus any local source

The attraction of the finite-volume approach is that conservation is built into the discretisation. If what leaves one cell enters the next, the algebra itself remembers that.

That is why finite volume is so common in transport and environmental systems. The method echoes the underlying physics instead of merely approximating the derivatives.

67.8 Worked example 3: why finite elements dominate structural analysis

In a structural problem, geometry is often awkward: holes, supports, joints, changing thickness, mixed materials. A regular grid is not always natural.

Finite elements instead approximate the displacement field using basis or shape functions:

u_h(x) = \sum_i U_i \phi_i(x)

The unknowns are the nodal coefficients U_i. The governing equations become a matrix system

K\mathbf{U} = \mathbf{f}

where K is the stiffness matrix — a different K from the control gain in Chapter 1 (Feedback control and state-space models), but the same letter is standard in both fields.

The important mathematical point is not the assembly details. It is that the continuum problem has been projected onto a finite-dimensional space spanned by the \phi_i. Geometry and material behaviour are therefore represented through the choice of mesh and basis functions, not just through a derivative stencil.

A second important point is that K is almost always sparse: each basis function \phi_i overlaps only with its neighbours, so most entries of K are zero. A mesh with millions of nodes still produces a matrix with only a few non-zeros per row. This sparsity is why large FEM problems are computationally tractable — but solver choice matters.

Direct sparse solvers (such as UMFPACK or CHOLMOD) factorize K in one shot. They work well at moderate problem sizes but produce fill-in: new non-zero entries created during factorization that were zero in the original K. Fill-in grows faster than the problem does, which limits how far direct methods can scale.

Iterative solvers — Krylov methods such as conjugate gradient (for symmetric positive-definite K, meaning K = K^T and \mathbf{x}^T K \mathbf{x} > 0 for all nonzero \mathbf{x}, which holds for many standard structural problems) or GMRES (for non-symmetric cases such as convection-dominated flow) — access K only through matrix-vector products. They can reach millions of unknowns where direct methods cannot.

Both solver families operate on K stored in compressed sparse row (CSR) or an equivalent format. At scale, understanding these choices is inseparable from using finite element analysis.

67.9 Where this goes

This chapter draws on two prior threads: the numerical methods sequence (finite differences, ODE solvers, linear systems) and continuum transport (the conservation laws that finite volume discretises directly). If the stability analysis or the flux-balance notation felt unfamiliar, those are the chapters to revisit.

The next natural continuation is Nonlinear optimisation for design and operations. Once a simulation model can be trusted to represent the system, design questions become computational questions: which shape, parameter, schedule, or control setting produces the best feasible outcome?

This chapter also deepens what you should mean by “model.” In upper-year work, the numerical method is not downstream implementation detail. It is part of the model’s claim to credibility.

TipApplications
  • heat-transfer simulation in machine components
  • finite-element structural analysis
  • finite-volume flow and transport models
  • weather, hydrology, and terrain-based environmental models
  • sparse linear systems in scientific computing
  • mesh-sensitivity and timestep-sensitivity studies

67.10 Exercises

These are project-style exercises. Do not stop at the formula. Explain what the numerical choice means for the model.

67.10.1 Exercise 1

For the explicit diffusion scheme with D = 2, \Delta x = 0.05, and \Delta t = 0.0004, compute

r = \frac{D\Delta t}{(\Delta x)^2}

and decide whether the time step is plausibly stable for the standard explicit scheme.

67.10.2 Exercise 2

A finite-volume cell has inflow flux F_{i-1/2} = 3, outflow flux F_{i+1/2} = 5, source term s_i = 1, and width \Delta x = 2.

Write the sign of the net flux contribution and explain whether storage tends to increase or decrease before the source term is included.

67.10.3 Exercise 3

Write a short comparison note explaining when you would prefer:

  1. finite difference
  2. finite volume
  3. finite element

Use one example application for each method.

67.10.4 Exercise 4

Choose one simulation problem from your own field and prepare a one-page model brief naming:

  1. the governing equation
  2. the boundary and initial conditions
  3. the likely discretisation family
  4. one likely stability or conditioning concern
  5. one mesh- or timestep-sensitivity check you would run before trusting the result