8.2 Discretization methods

Discretization methods overview

Discretization methods convert the continuous partial differential equations (PDEs) governing fluid flow into systems of algebraic equations that a computer can actually solve. Without them, the Navier-Stokes equations and other governing equations of aerodynamics would remain analytically intractable for all but the simplest cases.

These methods work by dividing the computational domain into smaller, discrete elements or volumes, then solving the governing equations at those discrete points. The choice of method depends on the complexity of the geometry, the desired accuracy, and the computational resources available.

Four main approaches dominate CFD:

• Finite difference method (FDM) approximates derivatives directly on a structured grid
• Finite volume method (FVM) enforces conservation laws over small control volumes
• Finite element method (FEM) represents the solution as a weighted sum of basis functions on unstructured elements
• Spectral methods expand the solution in global basis functions like Fourier or Chebyshev polynomials

Each offers different trade-offs between accuracy, computational cost, and geometric flexibility.

Finite difference method

The finite difference method approximates derivatives in the governing PDEs using Taylor series expansions, replacing continuous derivatives with algebraic differences evaluated at discrete grid points.

• Straightforward to implement on structured (rectangular) grids
• Computationally efficient for simple domain shapes
• Accuracy degrades on complex geometries because structured grids can't conform well to curved boundaries
• Struggles with discontinuities (e.g., shocks) unless special treatment is added

FDM is often the first method students encounter because the connection between the Taylor series and the numerical approximation is very direct.

Finite volume method

The finite volume method divides the domain into a set of non-overlapping control volumes and integrates the governing equations over each one. This guarantees that mass, momentum, and energy are conserved at the discrete level, not just approximately.

• Naturally enforces conservation, which is critical for capturing shocks and other flow features correctly
• Handles complex geometries well using unstructured meshes
• The most widely used discretization approach in production CFD codes (ANSYS Fluent, OpenFOAM, etc.)
• Accuracy depends heavily on the flux interpolation scheme chosen

Finite element method

FEM discretizes the domain into elements (triangles, quadrilaterals in 2D; tetrahedra, hexahedra in 3D) and approximates the solution as a weighted sum of basis functions defined on each element.

• Excels at handling complex and irregular geometries
• Supports adaptive mesh refinement, where you concentrate resolution only where it's needed
• Computationally more expensive per degree of freedom than FDM or FVM
• Dominant in structural analysis and increasingly used in aeroacoustics and high-order CFD

Spectral methods

Spectral methods represent the solution as a linear combination of global basis functions, typically Fourier series (for periodic problems) or Chebyshev polynomials (for non-periodic problems).

• Achieve exponential convergence for smooth solutions, meaning accuracy improves dramatically as you add more basis functions
• Require far fewer grid points than FDM or FVM for the same accuracy on smooth problems
• Poorly suited to complex geometries or flows with discontinuities (shocks cause Gibbs-type oscillations)
• Most commonly used in direct numerical simulation (DNS) of turbulence in simple geometries

Finite difference formulations

Forward, backward, and central differences

These are the building blocks of finite difference approximations. Each estimates the derivative $f'(x)$ using nearby function values separated by a grid spacing $h$:

• Forward difference: $f'(x) \approx \frac{f(x+h) - f(x)}{h}$ — first-order accurate, $O(h)$
• Backward difference: $f'(x) \approx \frac{f(x) - f(x-h)}{h}$ — first-order accurate, $O(h)$
• Central difference: $f'(x) \approx \frac{f(x+h) - f(x-h)}{2h}$ — second-order accurate, $O(h^2)$

Central differences are more accurate because the leading error terms from the Taylor expansion cancel. Higher-order stencils (using more neighboring points) can be derived the same way to achieve $O(h^4)$ or better accuracy.

Explicit vs. implicit schemes

When advancing the solution in time, you have two fundamental choices:

Explicit schemes calculate the solution at the next time step using only known values from the current time step.

• Simple to code and cheap per time step
• Conditionally stable: the time step must satisfy a stability limit (see CFL condition below), which can force very small steps on fine grids

Implicit schemes solve a system of equations that couples unknown values at the next time step together.

• Unconditionally stable for many formulations, allowing much larger time steps
• Each time step is more expensive because you must solve a linear (or nonlinear) system
• Often more efficient overall for steady-state problems or stiff equations where explicit time steps would be prohibitively small

Accuracy and stability considerations

• Spatial accuracy depends on the order of the finite difference stencil and the grid resolution
• Temporal accuracy depends on the order of the time integration scheme (e.g., forward Euler is first-order, Runge-Kutta 4 is fourth-order) and the time step size
• Stability is governed by the interplay between scheme type, grid spacing, and time step

The Courant-Friedrichs-Lewy (CFL) condition is the key stability criterion for explicit schemes:

$CFL = \frac{u \, \Delta t}{\Delta x} \leq CFL_{max}$

where $u$ is the characteristic velocity, $\Delta t$ is the time step, and $\Delta x$ is the grid spacing. For many explicit schemes, $CFL_{max} = 1$. Violating this condition causes the numerical solution to blow up.

Finite volume formulations

Control volume approach

The finite volume method works by:

1. Dividing the domain into non-overlapping control volumes (cells)
2. Integrating the governing equations in their conservation (integral) form over each control volume
3. Evaluating fluxes across each cell face so that what flows out of one cell flows into its neighbor

This flux-balance structure is what guarantees discrete conservation. It's the main reason FVM dominates aerodynamic CFD.

Conservation laws

The three fundamental conservation equations in differential form:

Mass (continuity): $\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \vec{u}) = 0$

Momentum: $\frac{\partial (\rho \vec{u})}{\partial t} + \nabla \cdot (\rho \vec{u} \otimes \vec{u}) = -\nabla p + \nabla \cdot \overline{\overline{\tau}}$

Energy: $\frac{\partial (\rho E)}{\partial t} + \nabla \cdot (\rho H \vec{u}) = \nabla \cdot (k \nabla T) + \nabla \cdot (\overline{\overline{\tau}} \cdot \vec{u})$

In FVM, these are integrated over each control volume and converted (via the divergence theorem) into surface integrals of fluxes across cell faces.

Flux evaluation and interpolation

Since solution values are stored at cell centers but fluxes are needed at cell faces, interpolation is required:

• Upwind schemes use values from the upstream direction of the flow. First-order upwind is very stable but diffusive; second-order upwind improves accuracy while retaining the directional bias.
• Central schemes use a symmetric stencil around the face. They're more accurate on smooth flows but can produce oscillations near shocks.
• Gradient reconstruction methods (Green-Gauss, least-squares) compute gradients at cell centers, which are then used to interpolate face values for higher-order accuracy.

The choice of flux scheme has a huge impact on solution quality. First-order upwind smears shocks and boundary layers; second-order or higher schemes resolve them more sharply but may need limiters to prevent spurious oscillations.

Boundary conditions treatment

Boundary conditions are enforced by modifying the fluxes at domain boundaries. Common types in aerodynamics:

• No-slip wall: $\vec{u} = 0$ at the surface (fluid velocity matches the wall)
• Inflow: specified velocity, pressure, or total conditions
• Outflow: zero-gradient (extrapolated) or specified static pressure
• Symmetry: zero normal velocity and zero normal gradient for scalar quantities

Implementation typically uses ghost cells (fictitious cells outside the boundary whose values are set to enforce the desired condition) or direct extrapolation from interior cells.

Finite element formulations

Weak form of governing equations

FEM doesn't solve the PDEs directly. Instead, it converts them into a weak form:

1. Multiply the governing PDE by a test function $w(x)$
2. Integrate the product over the domain
3. Apply integration by parts to shift derivatives from the solution onto the test function, reducing the smoothness requirements on the solution
4. Incorporate boundary conditions into the resulting boundary integral terms

The weak form is equivalent to the original PDE for smooth solutions but allows a broader class of approximate solutions.

Element types and shape functions

The domain is divided into elements, and within each element the solution is approximated using shape functions (also called basis functions):

• Element types: triangles and quadrilaterals (2D), tetrahedra and hexahedra (3D)
• Shape function order: linear, quadratic, or higher-order polynomials
• Common families: Lagrange (nodal) basis functions, hierarchical basis functions

Shape functions are constructed so that the solution is continuous across element boundaries. Higher-order shape functions improve accuracy but increase the number of unknowns per element.

Galerkin method

The Galerkin method is the standard FEM approach where the test functions are chosen to be the same as the shape functions. This choice:

• Produces a symmetric system of equations for self-adjoint problems (e.g., diffusion)
• Provides optimal convergence rates in the appropriate mathematical norm
• Can produce oscillations for convection-dominated flows, which is why stabilized variants (SUPG, Petrov-Galerkin) are often used in aerodynamics

Assembly and solution procedures

Once element-level equations are formed, they must be combined and solved:

1. Assembly: Element matrices and vectors are assembled into a global system $[K]\{u\} = \{F\}$
2. Storage: Sparse matrix formats (CSR, CSC) store only the nonzero entries, since most entries in the global matrix are zero
3. Solution: The linear system is solved using either:
   • Direct methods (LU decomposition, Cholesky factorization): exact to machine precision but scale poorly with problem size (memory and time grow rapidly)
   • Iterative methods (Conjugate Gradient, GMRES): scale much better for large problems but require good preconditioners to converge efficiently

Spectral method formulations

Fourier and Chebyshev polynomials

The two most common spectral basis function families:

Fourier series are used for periodic domains: $f(x) = \sum_{k=-\infty}^{\infty} \hat{f}_k e^{ikx}$

Fourier transforms can be computed efficiently using the Fast Fourier Transform (FFT), which operates in $O(N \log N)$ time. This makes Fourier-based spectral methods extremely fast for periodic problems like homogeneous turbulence.

Chebyshev polynomials are used for non-periodic domains: $T_n(x) = \cos(n \arccos(x))$

Chebyshev grids cluster points near the domain boundaries, which is actually beneficial for resolving boundary layers in aerodynamic flows.

Collocation and Galerkin approaches

Two strategies for enforcing the governing equations:

• Collocation: Require the equations to be satisfied exactly at a set of collocation points. Simpler to implement and computationally cheaper, but can be less stable for nonlinear problems and may need filtering to suppress aliasing errors.
• Galerkin (spectral): Project the equations onto the basis functions by integrating against them. More stable and accurate for nonlinear problems, but computationally more expensive because of the required numerical integration (quadrature).

Advantages and limitations

• Exponential convergence for smooth solutions: adding a few more modes can drop the error by orders of magnitude
• Far fewer degrees of freedom needed compared to FDM or FVM for equivalent accuracy on smooth problems
• Poorly suited to complex geometries (the basis functions are defined on simple domains)
• Discontinuities (shocks) cause Gibbs phenomena, producing oscillations that pollute the entire solution
• Most effective for DNS of turbulence in canonical geometries (channels, periodic boxes)

Discretization error analysis

Truncation and round-off errors

Two distinct sources of error arise in any numerical simulation:

Truncation error is the difference between the exact PDE and its discretized approximation. It comes from dropping higher-order terms in Taylor series or truncating basis function expansions. Truncation error decreases as you refine the grid or increase the polynomial order.

Round-off error is the difference between the exact solution of the discrete equations and what the computer actually computes, due to finite-precision floating-point arithmetic. Round-off error accumulates over many operations and can become significant in long simulations or ill-conditioned systems. Using double precision (64-bit) instead of single precision (32-bit) or employing compensated summation algorithms helps mitigate this.

There's a practical tension: making the grid extremely fine reduces truncation error but increases the total number of operations, which can amplify round-off error.

Convergence and order of accuracy

Convergence means the numerical solution approaches the exact solution as the grid is refined ($h \to 0$) or the polynomial order is increased ($p \to \infty$).

Order of accuracy describes how fast the error decreases:

• For finite difference/volume methods: $\text{error} = O(h^n)$, where $n$ is the scheme order. A second-order scheme ($n=2$) cuts the error by a factor of 4 when you halve the grid spacing.
• For spectral methods: $\text{error} = O(e^{-\alpha N})$, where $N$ is the number of modes. This exponential convergence is dramatically faster for smooth solutions.

Higher-order schemes converge faster but can be less robust and more expensive per grid point.

Grid refinement studies

A grid refinement study systematically verifies that your simulation is converging at the expected rate:

1. Run the simulation on at least three successively refined grids (e.g., coarse, medium, fine)
2. Compare a key output quantity (lift coefficient, drag, pressure at a point) across the grids
3. Use Richardson extrapolation to estimate the grid-independent (exact) solution and the observed order of accuracy
4. Compute the Grid Convergence Index (GCI) to quantify the uncertainty band around your finest-grid result

This process is essential for establishing credibility in any CFD result. Reporting a single-grid answer without a refinement study gives no indication of how much discretization error remains.

Discretization scheme selection

Problem-specific considerations

No single method is best for all problems. The right choice depends on:

• Geometry complexity: FVM and FEM handle complex shapes (aircraft, engine inlets) with unstructured meshes. FDM and spectral methods work best on simple or rectangular domains.
• Solution smoothness: Spectral methods are optimal for smooth flows (attached boundary layers, low-speed aerodynamics). FVM with shock-capturing schemes is needed for transonic and supersonic flows with discontinuities.
• Boundary conditions: Spectral methods are most efficient for periodic conditions. FVM and FEM accommodate the full range of wall, inflow, outflow, and far-field conditions used in external aerodynamics.

Computational cost vs. accuracy trade-offs

• Higher-order schemes give better accuracy per degree of freedom but cost more per grid point
• For a given accuracy target, spectral methods need the fewest unknowns (if the solution is smooth), followed by high-order FEM, then FVM, then low-order FDM
• Adaptive mesh refinement (AMR) concentrates grid resolution only where it's needed (near shocks, boundary layers, vortex cores), reducing total cell count while maintaining accuracy in critical regions

Hybrid and adaptive methods

Modern CFD increasingly combines methods to exploit their respective strengths:

Hybrid methods:

• Spectral element method combines the geometric flexibility of finite elements with the high accuracy of spectral polynomial bases within each element
• Discontinuous Galerkin (DG) method uses a finite element framework with discontinuous basis functions at element interfaces, enabling sharp shock capturing and straightforward parallel implementation with adaptive refinement

Adaptive methods dynamically adjust resolution based on error estimates:

• h-adaptivity: refines or coarsens the mesh while keeping the polynomial order fixed
• p-adaptivity: increases or decreases the polynomial order within elements while keeping the mesh fixed
• hp-adaptivity: combines both, using h-refinement near discontinuities (where high-order polynomials would oscillate) and p-refinement in smooth regions (where adding polynomial order is more efficient than adding cells)