Home / Courses / Computational Fluid Dynamics / Solvers

Computational Fluid Dynamics · Chapter 6 of 10 · Intermediate

Solving the Discretized Equations

Discretization leaves a huge sparse system of algebraic equations. Solving it directly is too costly, so CFD iterates, and watches the residual fall to know when to stop.

01

Readiness check

This chapter solves the assembled system. Tick only what you can do closed-notes.

  • Recognise a tridiagonal matrix.
  • Do forward elimination and back substitution.
  • Rearrange an equation to solve for one unknown.
  • Recall the cell equation aPφP = Σanbφnb + b.
  • Compute a residual as b − Ax.
0 or 1 weak itemsContinue with this chapter.
2 weak itemsRevisit the assembled system in Chapter 3.
3 or more weak itemsReview linear solvers in Math: Numerical Methods.
02

The core idea

Each cell equation links a node to a few neighbours, so the global matrix is large but sparse. A direct solver suits the tridiagonal one-dimensional case; everything bigger is solved by iteration, stopped when the residual is small.

aPφP = Σ anbφnb + bTDMA: forward sweep, back substitutionresidual R = b − Aφ → 0

Because a control volume only talks to its immediate neighbours, the coefficient matrix has nonzeros on just a few diagonals. In one dimension it is tridiagonal, and the Thomas algorithm (TDMA) solves it exactly in a single forward elimination and back substitution, far cheaper than general elimination. In two and three dimensions the matrix is sparse but not tridiagonal, and direct solution is too expensive, so point-iterative methods like Gauss-Seidel sweep through the cells, updating each from its neighbours, repeating until the residual, the amount each equation is unsatisfied, drops below tolerance. Watching that residual fall is how every CFD run knows it has converged.

The skill works when: you use TDMA for tridiagonal systems and iterate with a residual stopping test otherwise.
The skill breaks down when: a run is stopped before the residual has fallen, or a direct solve is attempted on a huge 3D matrix.
The concept. Each cell couples only to its neighbours, so the matrix is sparse. In one dimension it is tridiagonal (shaded band), which the Thomas algorithm solves directly; larger systems are iterated.
03

The skills, taught in order

Five skills cover the sparse system, the direct-versus-iterative choice, the TDMA, point iteration, and convergence.

6.1 The sparse linear system

Assembling every cell equation gives Aφ = b, where A is large but mostly zero because each row has only a handful of neighbour entries. Exploiting this sparsity is what makes CFD tractable; a dense solve would be impossibly expensive.

6.2 Direct versus iterative

Direct methods solve the system exactly in a fixed number of operations but scale badly for large matrices. Iterative methods improve a guess step by step and stop at a tolerance, which suits the millions of cells in a real mesh. CFD is almost always iterative.

6.3 The TDMA

For a tridiagonal system the Thomas algorithm (TDMA) is a direct solver in disguise: a forward sweep eliminates the sub-diagonal, then a back substitution recovers every unknown. It is exact and cheap, and it is applied line by line inside two- and three-dimensional solvers.

MethodTypeBest for
TDMA (Thomas)directtridiagonal, one-dimensional lines
Gauss-Seideliterativelarge sparse systems
Jacobiiterativesimple, parallel
Multigriditerativevery large, fast convergence

6.4 Point-iterative methods

Gauss-Seidel sweeps the cells, setting each φP = (Σanbφnb + b)/aP using the latest neighbour values. It converges faster than Jacobi because it reuses updates immediately. Multigrid accelerates convergence by solving on a hierarchy of coarse and fine grids.

6.5 Residuals and convergence

The residual R = b − Aφ measures how far each equation is from being satisfied. A run is converged when the normalised residual falls below a tolerance, often three to five orders of magnitude. Residual plots are the standard convergence diagnostic in every code.

Engineering connection: these solvers run inside the SIMPLE loop, solving the momentum, pressure-correction, and scalar equations on every outer iteration.

04

Worked example 1: solving a tridiagonal system with the TDMA

A discretized one-dimensional problem gives −φi−1 + 2φi − φi+1 = 2 for four interior nodes, with boundary values φ0 = φ5 = 0. Solve the tridiagonal system with the TDMA.

Figure 1. The discrete Poisson problem with a uniform source. The TDMA returns the exact nodal values in one forward and one backward sweep.
  1. ProblemSolve the four-node tridiagonal system in Figure 1.
  2. Given / find−φi−1 + 2φi − φi+1 = 2, φ0 = φ5 = 0. Find φ1 to φ4.
  3. AssumptionsTridiagonal system; the TDMA is exact for it.
  4. ModelWrite the four equations, then apply the forward elimination and back substitution of the Thomas algorithm.
  5. Equations2φ1 − φ2 = 2 (φ0 = 0)−φi−1 + 2φi − φi+1 = 2 (interior)
  6. SolveThe forward sweep gives modified coefficients, and back substitution yields φ = [4, 6, 6, 4].
  7. CheckSubstituting back: 2(4) − 6 = 2, −4 + 2(6) − 6 = 2, all four equations hold. The values match the analytic parabola φ = −x² + 5x at x = 1, 2, 3, 4.
  8. ConclusionThe TDMA solved the system exactly and cheaply. The same line solver is swept through a multidimensional mesh row by row.
Result. φ = [4, 6, 6, 4], satisfying every equation exactly.
05

Worked example 2: Gauss-Seidel iteration

Solve the same system by Gauss-Seidel iteration, starting from φ = [0, 0, 0, 0]. Carry out three sweeps and show the values approaching the exact solution.

Figure 2. Gauss-Seidel improves the estimate each sweep, marching the nodal values toward the exact solution. The residual falls correspondingly, the universal CFD convergence signal.
  1. ProblemIterate the system in Figure 2 with Gauss-Seidel for three sweeps.
  2. Given / findSame equations, start φ = [0, 0, 0, 0]. Find the values after three sweeps.
  3. AssumptionsSweep in order 1 to 4, using the latest available neighbour values.
  4. ModelUpdate each node as φi = (2 + φi−1 + φi+1)/2, with boundary zeros.
  5. Equationsφi = (b + φi−1 + φi+1)/aP
  6. SolveSweep 1: [1, 1.5, 1.75, 1.875]. Sweep 2: [1.75, 2.75, 3.31, 2.66]. Sweep 3: [2.38, 3.84, 4.25, 3.13], all climbing toward the exact [4, 6, 6, 4].
  7. CheckEach sweep moves the values closer to the exact solution and shrinks the residual; many more sweeps would reach [4, 6, 6, 4] to tolerance. Gauss-Seidel converges because the system is diagonally dominant.
  8. ConclusionIteration trades exactness per step for cheapness, reaching the solution gradually. The residual, not a fixed step count, decides when to stop.
Result. After three sweeps φ ≈ [2.38, 3.84, 4.25, 3.13], converging toward [4, 6, 6, 4].
06

Misconceptions and diagnostics

MistakeSymptomDiagnostic questionCorrection
Direct solve on a huge matrixRuns out of memory or time"How many cells are there?"Use iterative solvers for large 2D and 3D systems.
Stopping too earlyResult still changing"Has the residual dropped enough?"Iterate until the normalised residual is below tolerance.
Jacobi when Gauss-Seidel is fineNeedlessly slow convergence"Can I reuse updated values immediately?"Gauss-Seidel reuses updates and converges faster.
Ignoring diagonal dominanceIteration diverges"Is aP at least the sum of neighbours?"Ensure the coefficients keep the matrix diagonally dominant.
07

Practice ladder

Level 1 · Direct skill

One Gauss-Seidel update is φ2 = (2 + φ1 + φ3)/2. With φ1 = 1 and φ3 = 0, find φ2.

Show answer

φ2 = (2 + 1 + 0)/2 = 1.5, the first-sweep value from the worked example.

Level 2 · Mixed concept

After the worked example's three sweeps, by roughly how much is φ3 = 4.25 still short of the exact 6?

Show answer

It is short by 1.75, about 29 percent. Several more sweeps are needed to drive the residual below a typical tolerance.

Level 3 · Independent problem

For 2φ1 − φ2 = 1 and −φ1 + 2φ2 = 1, do two Gauss-Seidel sweeps from zero and compare to the exact solution.

Show answer

Exact: φ1 = φ2 = 1. Sweep 1: φ1 = 0.5, φ2 = 0.75. Sweep 2: φ1 = 0.875, φ2 = 0.938. Converging to 1.

Level 4 · Transfer to real engineering

Explain how you would read a residual plot from a real CFD run to decide whether the solution has converged.

What good work looks like

Residuals dropping several orders of magnitude and leveling off, plus a check that an integrated quantity (drag, mass balance) has stopped changing.

08

Working with AI, and proving it yourself

Use AI as an examiner, not a solver

"Check that my Gauss-Seidel updates use the latest neighbour values."
"Give me three matrices; I will say whether they are diagonally dominant."
"Solve this system." Doing the TDMA or sweeps yourself is the skill.
"Has it converged?" Reading the residual is the point.

Portfolio task

Solve a small discretized system both ways: the TDMA for an exact answer and Gauss-Seidel for a few sweeps, and plot the residual against iteration.

Must include: a TDMA solution, several Gauss-Seidel sweeps, and a residual that falls toward zero.
09

Retrieval and spaced review

Closed notes. Answer out loud, then reveal.

1. Why is the CFD matrix sparse?

Each cell equation links only to its neighbours, so most entries are zero.

2. What is the TDMA good for?

Tridiagonal systems, solved exactly by a forward sweep and back substitution.

3. How does Gauss-Seidel update a node?

φP = (Σanbφnb + b)/aP, using the latest neighbour values.

4. What is the residual?

R = b − Aφ, how far each equation is from being satisfied.

5. When is a run converged?

When the normalised residual falls below tolerance, often by several orders of magnitude.

TodayFinish this quiz and Levels 1 and 2 of the ladder.
+1 dayRe-derive the TDMA result and two sweeps from a blank page.
+3 daysIterate two new systems.
+7 daysCarry the solvers into unsteady flow, Chapter 7.
+30 daysReuse the residual as the stopping test in any CFD run.
10

Textbook mapping

ItemMapping
Primary sourceVersteeg and Malalasekera, An Introduction to Computational Fluid Dynamics, Chapter 7 (Solution of Discretised Equations)
Cross-referencePatankar, Ch. 4 · Math: Numerical Methods
Core topics6.1 Sparse system · 6.2 Direct vs iterative · 6.3 TDMA · 6.4 Point iteration · 6.5 Residuals and convergence
Engineering connectionThese solvers run inside the SIMPLE loop on every iteration.
Read nextChapter 7: Unsteady Flows and Time Stepping.