Solvers & Convergence

After discretizing the governing equations, we have a massive system of linear equations. For a 3D mesh with 10 million cells solving for 5 variables, that's 50 million unknowns. Direct solution methods like Gaussian elimination have $O(n^3)$ complexity — utterly impractical. Instead, CFD relies on iterative solvers that progressively refine an initial guess.

The Linear System

Matrix Form

The discretized equations form:

$$\mathbf{A}\mathbf{x} = \mathbf{b}$$

Where:

• $\mathbf{A}$ = coefficient matrix (sparse, mostly zeros)
• $\mathbf{x}$ = solution vector (velocities, pressure, etc.)
• $\mathbf{b}$ = right-hand side (source terms, boundary conditions)

Sparsity Pattern

For a 5-point stencil in 2D (or 7-point in 3D), each equation involves only a few neighbors:

$$a_P \phi_P = a_E \phi_E + a_W \phi_W + a_N \phi_N + a_S \phi_S + b_P$$

The matrix has at most 7 non-zeros per row (in 3D). Storing only non-zeros is essential.

Basic Iterative Methods

Jacobi Method

The simplest iterative approach: solve for each unknown using values from the previous iteration.

$$\phi_P^{n+1} = \frac{1}{a_P}\left( a_E \phi_E^n + a_W \phi_W^n + a_N \phi_N^n + a_S \phi_S^n + b_P \right)$$

• Naturally parallelizable (all updates independent)
• Slow convergence
• Requires diagonal dominance: $|a_P| > \sum |a_{nb}|$

Gauss-Seidel Method

Use the latest available values immediately:

$$\phi_P^{n+1} = \frac{1}{a_P}\left( a_E \phi_E^n + a_W \phi_W^{n+1} + a_N \phi_N^n + a_S \phi_S^{n+1} + b_P \right)$$

• Faster convergence than Jacobi (about 2x)
• Sequential updates reduce parallelism
• Sweep direction matters (influences which neighbors are "new")

Successive Over-Relaxation (SOR)

Accelerate Gauss-Seidel with an over-relaxation factor $\omega$:

$$\phi_P^{n+1} = (1-\omega)\phi_P^n + \omega \phi_P^{GS}$$

Where $\phi_P^{GS}$ is the Gauss-Seidel update.

• $\omega = 1$: Standard Gauss-Seidel
• $1 < \omega < 2$: Over-relaxation (faster)
• $\omega \approx 1.5 - 1.9$: Typical for CFD
• $\omega \geq 2$: Unstable

Residual Monitoring

What Is a Residual?

The residual measures how well the current solution satisfies the equations:

$$r_P = b_P - a_P \phi_P + \sum a_{nb} \phi_{nb}$$

For a perfect solution, $r_P = 0$ everywhere.

Scaled Residuals

Raw residuals depend on variable magnitudes. Scaled residuals are normalized:

$$R_\phi = \frac{\sum_{cells} |r_P|}{\sum_{cells} |a_P \phi_P|}$$

Alternatively, normalize by inlet flux or other reference values.

Typical Convergence Criteria

Monitoring Beyond Residuals

Residuals alone don't guarantee accuracy. Also monitor:

• Mass balance — Net mass flow should be zero (closed domain) or match inlet (open)
• Key quantities — Drag coefficient, pressure drop, exit temperature
• Point monitors — Values at specific locations should stabilize

Convergence Patterns

Healthy Convergence

• Residuals drop monotonically
• 3-4 orders of magnitude reduction
• Key quantities stabilize

Warning Signs

Multigrid Methods

The Problem with Smooth Errors

Basic iterative methods eliminate high-frequency errors (cell-to-cell oscillations) quickly but struggle with low-frequency errors (smooth, gradual variations). On fine meshes, smooth errors look like constants locally.

The Multigrid Idea

Low-frequency errors on a fine grid become high-frequency on a coarser grid. The multigrid strategy:

• Smooth on fine grid (removes high-frequency errors)
• Restrict residual to coarser grid
• Solve (approximately) on coarse grid
• Prolongate correction back to fine grid
• Smooth again on fine grid

V-Cycle Algorithm

multigrid(level):
    if level == coarsest:
        solve_directly()
    else:
        smooth(level)                    // Pre-smoothing
        residual = compute_residual()
        coarse_residual = restrict(residual)
        coarse_correction = multigrid(level+1)
        correction = prolongate(coarse_correction)
        apply_correction(correction)
        smooth(level)                    // Post-smoothing

Why Multigrid Is Optimal

For a grid with $N$ cells:

• Basic iterative methods: $O(N^2)$ operations to converge
• Multigrid: $O(N)$ operations — optimal scaling

This is why multigrid is essential for large-scale CFD.

Algebraic Multigrid (AMG)

• Works on unstructured meshes
• Automatic coarsening
• Used in most commercial solvers

Pressure Equation Challenges

The pressure correction equation (Poisson-type) is often the hardest to solve:

• Elliptic nature — Information propagates everywhere
• No diagonal dominance — Near-singular for some BCs
• Slow convergence — Low-frequency modes dominate

Solver Selection

Equation Type Matters

Krylov Subspace Methods

Modern CFD often uses Krylov methods:

• Conjugate Gradient (CG) — Symmetric positive definite only
• BiCGSTAB — General non-symmetric systems
• GMRES — Robust, needs memory for basis vectors

Combined with preconditioners (ILU, AMG), these are state-of-the-art.

Under-Relaxation Revisited

Why It Helps Convergence

The SIMPLE algorithm over-estimates corrections. Under-relaxation stabilizes:

$$\phi^{new} = \phi^{old} + \alpha (\phi^{calc} - \phi^{old})$$

Effect on Convergence

Adaptive Under-Relaxation

Some solvers adjust $\alpha$ based on residual behavior:

• Increase if converging smoothly
• Decrease if oscillating

Practical Convergence Strategy

Starting a Difficult Case

• Initialize well — Use potential flow or steady boundary values
• Start conservative — Low under-relaxation, first-order schemes
• Ramp up gradually — Increase $\alpha$, switch to higher-order
• Monitor everything — Not just residuals

Convergence Checklist

• [ ] Residuals dropped 3+ orders of magnitude
• [ ] Mass imbalance < $10^{-6}$ of inlet flux
• [ ] Key quantities stabilized (< 0.1% change)
• [ ] Physical plausibility verified
• [ ] Results independent of further iterations

Steady vs. Transient

Steady-State Iterations

• Goal: Find time-independent solution
• Use SIMPLE, under-relaxation
• Many outer iterations (hundreds to thousands)

Transient Time-Stepping

• Goal: Resolve time evolution
• Use PISO or dual time-stepping
• Few iterations per time step (3-10)
• Residuals must converge each step

Common Convergence Issues

Issue: Residuals Stuck at $10^{-3}$

• Mesh quality issues (high skewness)
• Inconsistent boundary conditions
• Physics not reaching steady state

• Improve mesh quality
• Check boundary condition values
• Consider transient solution

Issue: One Variable Won't Converge

• Check outlet BC (needs pressure reference)
• Verify mass balance at boundaries
• Increase pressure under-relaxation iterations

Issue: Periodic Oscillations

Key Takeaways

• Iterative solvers are essential for large CFD systems — direct methods don't scale
• Gauss-Seidel converges faster than Jacobi but is less parallelizable
• Multigrid achieves optimal $O(N)$ convergence by eliminating smooth errors on coarse grids
• Residuals measure equation imbalance — target 3-4 orders reduction
• Monitor key quantities — residuals alone don't guarantee accuracy
• Under-relaxation stabilizes coupled solvers like SIMPLE
• AMG enables multigrid on unstructured meshes automatically

What's Next

With the numerical machinery in place, we turn to one of CFD's greatest challenges: Turbulence Modeling. The next lesson covers Reynolds averaging, RANS equations, and popular turbulence models like k-ε and k-ω SST.