4. Spatial Discretization Methods

Introduction and context

The governing equations of fluid mechanics are normally expressed as partial differential equations in continuous space. A computational fluid dynamics (CFD) solution, however, is obtained at a finite number of points, cells, or elements. Spatial discretization is the process of converting spatial derivatives and volume or surface integrals in the governing equations into algebraic expressions that can be solved numerically.

This conversion is not merely a coding step. It determines whether a method conserves mass, momentum, and energy; how accurately it represents gradients and shocks; whether it produces excessive numerical diffusion; and how its computational cost grows as the mesh is refined. The discretization method therefore influences both the credibility of a simulation and the engineering decisions based on it.

Three major families dominate practical CFD: finite difference methods (FDM), finite volume methods (FVM), and finite element methods (FEM). Spectral methods form an important fourth family, particularly for smooth flows and fundamental research. Within each family, the choice of interpolation, gradient reconstruction, and convection scheme can be as important as the underlying governing equations.

This chapter develops these methods from a common conservation-law framework and examines their accuracy, consistency, conservation properties, stability-related behavior, and practical suitability.

Learning objectives

After studying this chapter, the reader should be able to:

1. Explain the role of spatial discretization in CFD and distinguish consistency, conservation, order of accuracy, and numerical diffusion.

2. Derive common finite difference approximations and assess their truncation errors.

3. Formulate a finite volume discretization from an integral conservation law.

4. Describe the Galerkin finite element formulation and the significance of weak forms.

5. Explain the basic principles and limitations of spectral methods.

6. Compare interpolation, gradient, and convection schemes for different flow regimes.

7. Select an appropriate spatial discretization strategy for engineering applications.

8. Diagnose common errors associated with mesh quality, boundary treatment, and excessive numerical diffusion.

From continuous equations to discrete algebra

Consider a scalar conservation equation,

∂φ/∂t + ∇·(uφ) = ∇·(Γ∇φ) + Sφ,

where φ is a transported variable, u is velocity, Γ is a diffusion coefficient, and Sφ is a source term. The same structure represents mass, momentum components, energy, or a passive scalar after suitable definitions of φ, Γ, and Sφ.

A spatial discretization replaces the continuous domain Ω by a set of discrete entities. In FDM these are usually nodal points; in FVM they are control volumes; in FEM they are elements connected through interpolation functions. The result is a system of ordinary differential equations in time, or—after temporal discretization—a system of algebraic equations:

Aφ = b.

The matrix A contains coefficients generated by diffusion, convection, source terms, and boundary conditions. Its properties affect convergence and boundedness. For example, a discretization that generates strong positive off-diagonal coefficients may produce oscillatory or unbounded solutions even when the physical problem is well posed.

Several concepts must be separated. A scheme is consistent if its discrete equation approaches the governing differential equation as grid spacing tends to zero. Its order of accuracy describes the rate at which the discretization error decreases with grid spacing. A conservative scheme ensures that flux leaving one control volume enters its neighbor with equal and opposite magnitude. Stability concerns the bounded response of the numerical solution to perturbations. Consistency and stability together are central to convergence, although practical nonlinear CFD also depends on solver behavior, boundary conditions, and mesh quality.

Finite difference methods: derivatives at discrete points

Finite difference methods approximate derivatives directly at selected points. For a uniform one-dimensional grid with spacing Δx, Taylor expansion gives

φ(x+Δx) = φ(x) + Δx φ′(x) + (Δx²/2)φ″(x) + (Δx³/6)φ‴(x) + ···.

Combining expansions produces the familiar formulas.

The forward difference is

φ′i ≈ (φi+1 − φi)/Δx,

with first-order truncation error. The backward difference is

φ′i ≈ (φi − φi−1)/Δx,

also first order. The central difference is

φ′i ≈ (φi+1 − φi−1)/(2Δx),

with second-order truncation error. For the second derivative,

φ″i ≈ (φi+1 − 2φi + φi−1)/Δx²,

which is second order on a uniform grid.

For a one-dimensional steady convection–diffusion equation,

ρu dφ/dx = Γ d²φ/dx²,

central differencing gives a balance involving φi−1, φi, and φi+1. It is accurate for smooth solutions, but when convection dominates diffusion, the resulting solution can oscillate. A first-order upwind difference instead uses the upstream value according to the sign of u. It is generally bounded and robust, but introduces artificial diffusion.

Nonuniform grids require geometrical distances and, often, coordinate transformations. Applying uniform-grid formulas without accounting for unequal spacing is a common source of inconsistency. In multiple dimensions, finite difference methods are straightforward on structured Cartesian or logically rectangular grids, but become more complicated on highly irregular geometries.

FDM is attractive because derivative formulas are simple, matrices are often compact, and high-order approximations are relatively easy to construct. Its principal limitation is geometric flexibility. Complex industrial domains, moving boundaries, and locally refined unstructured meshes are usually handled more naturally by FVM or FEM.

Finite volume methods: conservation through flux balance

The finite volume method begins with the integral form of the conservation law. For a control volume V with boundary surface S,

∂/∂t ∫V φ dV + ∮S (uφ)·n dA = ∮S Γ∇φ·n dA + ∫V Sφ dV.

The domain is divided into non-overlapping control volumes. Applying the divergence theorem converts volume divergence terms into surface fluxes. For a two-dimensional cell P with neighboring faces f, the semi-discrete equation is

(d/dt)(φP VP) + Σf Ff φf = Σf Df(φN − φP) + SP VP,

where Ff = (u·n)f Af is the mass flux through face f and Df represents a diffusive conductance. The key property is local conservation: a face flux is shared by two adjacent cells, so the flux leaving one cell enters the other with the opposite sign.

The finite volume formulation is especially powerful for compressible flow, multiphase flow, and flows containing shocks because conservation is enforced over each cell even on an irregular mesh. Global conservation follows from cancellation of internal face fluxes, leaving only boundary fluxes and sources.

Face interpolation and variable locations

Cell-centered FVM stores primary variables at cell centers, whereas fluxes are required at faces. Face interpolation is therefore fundamental. Linear interpolation on a uniform grid gives

φf = (φP + φN)/2.

On a nonuniform grid, a distance-weighted expression is more appropriate. For convective fluxes, the face value is chosen according to a convection scheme. In first-order upwind, φf equals the value in the upstream cell. Central differencing uses a centered interpolation. Higher-resolution schemes blend upwind robustness with higher-order accuracy.

The choice of face velocity is also important. In collocated arrangements, pressure and velocity are stored at the same locations, and interpolation can create pressure–velocity decoupling. Pressure–velocity coupling procedures such as Rhie–Chow-type momentum interpolation are commonly used to prevent checkerboard pressure fields.

Diffusion and non-orthogonal meshes

For an orthogonal mesh, the diffusive flux through a face is approximated by

(Γ∇φ·n)f Af ≈ Γf Af (φN − φP)/dPN.

On non-orthogonal meshes, the vector connecting neighboring cell centers is not aligned with the face normal. The gradient is decomposed into an orthogonal contribution and a non-orthogonal correction. The correction may be treated explicitly, iteratively, or through a more elaborate reconstruction. Ignoring it can reduce accuracy and, for severely distorted meshes, compromise solution quality.

Finite element methods: weighted residuals and weak formulations

Finite element methods approximate the unknown field within each element using shape functions:

φ(x) ≈ Σj Nj(x) φj.

The governing differential equation is multiplied by a weighting function w and integrated over the domain. For diffusion, integration by parts gives a weak form such as

∫Ω Γ ∇w·∇φ dΩ − ∫∂Ω w Γ∇φ·n dA + ∫Ω w u·∇φ dΩ = ∫Ω w Sφ dΩ.

The derivatives acting on φ are reduced in order, which permits solutions that are piecewise smooth and makes natural flux boundary conditions arise naturally from the boundary integral. With Galerkin FEM, the weighting functions are usually the same as the shape functions. Element matrices are calculated locally and assembled into a global sparse system.

FEM offers exceptional geometric flexibility, systematic order elevation through higher-degree shape functions, and a strong mathematical foundation for error estimation and adaptive refinement. It is widely used for diffusion-dominated flows, fluid–structure interaction, porous media, and multiphysics problems.

Pure Galerkin FEM can, however, generate oscillations in convection-dominated problems because the standard formulation has insufficient directional bias. Stabilized methods such as streamline upwind/Petrov–Galerkin, Galerkin least squares, and discontinuity-capturing formulations introduce controlled stabilization. These methods improve robustness but must be designed carefully to avoid excessive artificial diffusion or loss of accuracy near layers.

Spectral and high-order methods

Spectral methods represent a solution using globally smooth basis functions, for example Fourier series on periodic domains or Chebyshev polynomials on finite intervals. Derivatives are obtained analytically from the basis representation or through differentiation matrices. For sufficiently smooth solutions, the error can decrease very rapidly as the number of modes increases; this is often called spectral or exponential convergence.

The strength of spectral methods is their high accuracy per degree of freedom for smooth, well-resolved flows. Their weaknesses arise when the solution contains shocks, discontinuities, complex geometry, or localized nonsmooth features. Global basis functions can produce oscillations near discontinuities, known as Gibbs phenomena. Spectral-element methods combine high-order local elements with spectral-like basis functions and are more adaptable to complex domains.

Modern high-order discontinuous Galerkin methods also combine element-local polynomial representations with numerical fluxes at element interfaces. They offer conservation, high-order accuracy, and parallel scalability, but require careful flux, limiter, and time-step design.

Assessing a discretization

Consistency and truncation error

To test consistency, the exact smooth solution is substituted into the discrete equation. The residual that remains is the truncation error. For example, central differencing of a first derivative has an error proportional to Δx², whereas first-order upwind has an error proportional to Δx. A formal order does not guarantee equal practical accuracy: a second-order scheme on a poor-quality mesh may be less accurate than a first-order scheme on a well-designed mesh.

Conservation

Conservation is essential when integral quantities matter. FVM provides local conservation by construction when fluxes are shared consistently. FDM may conserve globally or locally only if specifically formulated in flux-difference form. FEM conservation depends on the formulation; standard continuous Galerkin methods do not automatically ensure local conservation, although mixed, discontinuous, or locally conservative formulations can provide it.

Numerical diffusion and dispersion

Numerical diffusion, or artificial viscosity, smooths gradients and broadens shear layers, thermal fronts, and shocks. It is especially associated with first-order upwind convection. Numerical dispersion produces phase errors and may appear as oscillations or wiggles, particularly in central or high-order schemes applied to under-resolved convection-dominated flows.

A modified-equation analysis makes these effects visible. A first-order upwind approximation to advection behaves like the original equation plus an artificial diffusion term whose magnitude is related to |u|Δx. Thus, refining the grid reduces numerical diffusion, while increasing velocity or using a larger cell increases it. Central schemes have lower leading-order diffusion but can possess dispersive error and may require stabilization.

Boundedness and monotonicity

A bounded solution should not create new maxima or minima when the physical equation does not permit them. Monotone schemes are valuable for shocks and steep fronts. High-order schemes can violate boundedness unless limited. Flux limiters and nonlinear reconstruction methods attempt to use high-order behavior in smooth regions and revert toward upwind behavior near strong gradients.

Mesh quality and grid convergence

Skewness, non-orthogonality, aspect ratio, abrupt size changes, and inadequate boundary-layer resolution affect the observed error. A credible study should include mesh refinement and compare integral outputs, local profiles, and conservation residuals. If a result changes substantially under refinement, the original mesh was not sufficient for the reported conclusion.

Interpolation, gradient, and convection schemes

Gradient reconstruction is required for diffusion, pressure correction, turbulence models, and higher-order convection. Common approaches include Green–Gauss reconstruction, which estimates gradients from face values and face normals, and least-squares reconstruction, which fits a local gradient to neighboring values. Least-squares methods are often more robust on irregular meshes, while Green–Gauss methods can be effective on high-quality meshes.

Central differencing is formally second order and has low numerical diffusion, making it suitable for smooth, diffusion-dominated flows and well-resolved turbulence. It can oscillate in convection-dominated or shock-containing flows. First-order upwind is strongly robust and bounded but smears important features. Second-order upwind reconstructs a gradient from the upstream cell and improves accuracy, though it may still require limiting.

Blended schemes combine central and upwind contributions. The blending factor should not be regarded as a substitute for mesh refinement; excessive upwind weighting can conceal inadequate resolution. QUICK and related biased higher-order schemes can improve accuracy on structured meshes, but may require limiters on irregular grids. Total variation diminishing and flux-limited schemes offer a practical compromise by reducing spurious oscillations without imposing first-order diffusion everywhere.

The appropriate scheme depends on the flow. For a low-Reynolds-number diffusion-dominated problem, central or higher-order schemes are usually suitable. For a high-Reynolds-number flow with strong separation, a bounded second-order method is often a reasonable engineering choice. For shocks, conservation, monotonicity, and an appropriate Riemann or numerical flux are more important than formal smooth-flow order.

An illustrative applied case: discretizing flow through a sudden expansion

Background

Consider steady incompressible flow through a two-dimensional channel that expands abruptly. The problem contains a recirculation region downstream of the expansion, strong velocity gradients near separation, and a pressure recovery zone.

Situation and problem

An initial simulation uses a coarse, skewed mesh and first-order upwind differencing for momentum. The pressure drop is plausible, but the recirculation bubble is shorter than expected and the velocity profile is overly smooth. A second simulation uses central differencing without sufficient mesh refinement and produces oscillatory velocities near the reattachment region.

Analysis

The first result is dominated by numerical diffusion. Upwind differencing damps the separated shear layer and transports excessive momentum across it. The second result has less diffusion but insufficient boundedness and resolution. Mesh skewness also degrades gradient reconstruction, especially near the expansion corner.

Approach and solution

The domain is remeshed with local refinement at the expansion, in the separated shear layer, and near the walls. Non-orthogonal correction is included in the diffusive flux. A least-squares gradient reconstruction is selected, and a bounded second-order convection scheme is used after an initial first-order solution establishes a stable flow field. The calculation is repeated on progressively finer meshes, while pressure drop, reattachment length, and mass imbalance are monitored.

Outcome

The refined solutions show a longer and more stable recirculation region, while the pressure drop approaches a mesh-independent value. The final result is accepted only after the change in key integral quantities between the two finest meshes becomes small and the mass imbalance is negligible.

Learning points

The case demonstrates that formal scheme order cannot compensate for inadequate mesh resolution. First-order schemes can be useful for initialization, but their results should not automatically be used for final engineering decisions. A bounded higher-order scheme, suitable gradients, mesh refinement, and conservation checks must be considered together.

Practical and professional considerations

In industrial CFD, discretization choices are technical decisions with managerial consequences. A highly accurate method may increase computational cost and require specialist verification, whereas a robust lower-order method may be appropriate for rapid design screening. The selection should reflect the purpose of the simulation: preliminary design, certification, safety assessment, optimization, or fundamental analysis.

Professional practice requires documenting mesh topology, cell count, quality measures, discretization schemes, convergence criteria, boundary treatment, and grid-sensitivity results. Reproducibility is particularly important when simulation results support safety, environmental compliance, or expensive equipment decisions. Numerical settings should never be adjusted merely to obtain a visually attractive flow field.

Digital CFD platforms increasingly automate mesh adaptation, error indicators, and scheme selection. Machine-learning tools can assist surrogate modeling and adaptive sampling, but they do not remove the need for conservation checks and physical validation. Automated refinement based only on residuals may miss an important quantity of interest; refinement should also consider gradients, vorticity, wall shear, heat flux, and uncertainty in the engineering output.

Challenges and common mistakes

Several errors recur in spatial discretization practice:

• Applying uniform-grid finite differences to strongly nonuniform meshes without deriving the appropriate coefficients.
• Treating a low residual as proof of accuracy without checking conservation and mesh independence.
• Using first-order upwind for final predictions of shear layers, wakes, or heat-transfer coefficients.
• Applying an unbounded central or high-order scheme to a shock or under-resolved front.
• Ignoring non-orthogonal correction in a distorted finite volume mesh.
• Using gradients reconstructed from poor neighboring cells near boundaries or sharp corners.
• Comparing two schemes while simultaneously changing the mesh, solver tolerances, and boundary conditions.
• Reporting only contour plots rather than quantitative errors, integral balances, and uncertainty.

No discretization is universally superior. Accuracy, boundedness, conservation, geometric flexibility, implementation complexity, and computational cost must be evaluated together.

Current developments and best practices

High-order finite volume, discontinuous Galerkin, spectral-element, and flux-reconstruction methods are being developed for large-eddy simulation, aeroacoustics, turbomachinery, and wave propagation. Adaptive mesh refinement is becoming more goal-oriented, focusing computational effort on specified outputs rather than uniformly reducing cell size. Hybrid meshes and high-order boundary representations are also improving the treatment of complex geometry.

A sound workflow is to begin with a conservative and robust scheme, verify boundary conditions and balances, then increase spatial order as the mesh and solution become reliable. Perform at least a basic grid-convergence study, inspect local gradients and integral outputs, and select schemes according to the physics rather than software defaults. For convection-dominated problems, use bounded high-resolution methods; for smooth problems, exploit higher-order accuracy; for shocks, prioritize conservation and monotonicity.

Chapter Summary

Spatial discretization transforms continuous fluid equations into solvable algebraic systems. Finite difference methods approximate derivatives at points and are efficient on structured grids. Finite volume methods integrate conservation laws over control volumes and provide a natural framework for local conservation on complex meshes. Finite element methods use weighted residuals and weak forms, offering geometric flexibility and systematic high-order approximation. Spectral methods achieve exceptional accuracy for smooth solutions but are less suited to discontinuities and complex nonsmooth flows.

Consistency, truncation error, order of accuracy, conservation, boundedness, numerical diffusion, and numerical dispersion are distinct criteria that must be considered together. Interpolation and gradient reconstruction determine how cell or nodal information is converted into fluxes and derivatives. Central schemes reduce diffusion but may oscillate; upwind schemes improve robustness but smear gradients; bounded higher-resolution schemes seek a practical balance. Reliable CFD requires appropriate mesh quality, conservation checks, refinement studies, and scheme selection based on the flow physics and engineering purpose.

Key Terms

Consistency: The property that a discrete equation approaches the governing differential equation as the mesh spacing tends to zero.

Control volume: A finite region over which a conservation law is integrated in the finite volume method.

Convection scheme: A numerical rule used to approximate transported quantities at faces or interfaces.

Diffusion: Physical or numerical smoothing caused by molecular transport or discretization error.

Discontinuous Galerkin method: A finite element-type method using discontinuous polynomial fields and numerical interface fluxes.

Finite difference method: A method that approximates derivatives using values at discrete points.

Finite element method: A method that approximates a field with element shape functions and derives equations through weighted residuals.

Finite volume method: A method based on integral conservation over non-overlapping control volumes.

Flux: The rate at which a conserved quantity crosses a surface.

Galerkin method: A weighted residual method in which weighting and approximation functions are selected from the same space.

Gradient reconstruction: Estimation of spatial derivatives from neighboring discrete values.

Mesh convergence: The approach of a computed result toward a limiting value as the mesh is refined.

Monotonicity: The property of avoiding new numerical maxima or minima in situations where the physical solution is bounded.

Numerical dispersion: Phase or wave-propagation error that can create oscillatory numerical behavior.

Numerical diffusion: Artificial smoothing introduced by a discretization scheme.

Order of accuracy: The power governing how discretization error decreases with mesh spacing.

Spectral method: A method representing the solution with global smooth basis functions such as Fourier modes or polynomials.

Truncation error: The residual produced when an exact solution is substituted into a discrete equation.

Upwind scheme: A convection discretization that uses information from the upstream direction.

Weak form: An integral formulation in which derivatives are transferred partly from the solution to weighting functions.

Review Questions

Short Answer Questions

1. What is meant by consistency in a spatial discretization scheme?

2. Why does the finite volume method provide local conservation?

3. State the formal order of accuracy of central and first-order upwind differencing for a first derivative.

4. What is numerical diffusion, and how does it affect a shear layer?

5. Why can standard Galerkin finite elements oscillate in convection-dominated flow?

Long Answer Questions

1. Derive the central difference approximation for a first derivative using Taylor series and identify its truncation error.

2. Explain the finite volume formulation of a convection–diffusion equation, including face fluxes and conservation.

3. Compare finite difference, finite volume, finite element, and spectral methods with respect to geometry, conservation, accuracy, and application.

4. Discuss the roles of interpolation, gradient reconstruction, and convection schemes in a cell-centered CFD solver.

5. Explain how mesh quality, numerical diffusion, boundedness, and grid convergence influence the credibility of a CFD result.

Analytical / Application Questions

1. A high-Reynolds-number channel simulation produces a wake that is much shorter with first-order upwind than with a bounded second-order scheme. Explain the numerical cause and propose a verification strategy.

2. Design a discretization strategy for compressible flow containing a shock and a smooth downstream boundary layer. Justify your choices of conservation property, convection treatment, gradient reconstruction, and mesh refinement.

3. A finite volume mesh has severe non-orthogonality near a curved wall. Analyze the likely errors if non-orthogonal corrections are omitted and recommend corrective actions.

22. Multiple Choice Questions with Answers

Question 1

Which feature is most directly guaranteed by a standard conservative finite volume flux formulation?

A. Exponential convergence for smooth solutions

B. Exact local cancellation of shared internal face fluxes

C. Absence of numerical diffusion

D. Second-order accuracy on every mesh

Correct Answer: B. Exact local cancellation of shared internal face fluxes

Question 2

What is the formal truncation order of the central difference approximation for a first derivative on a uniform grid?

A. First order

B. Second order

C. Third order

D. Zeroth order

Correct Answer: B. Second order

Question 3

A first-order upwind scheme is often robust in convection-dominated flow primarily because it

A. eliminates all truncation error

B. introduces stabilizing artificial diffusion

C. uses a global spectral basis

D. guarantees exact pressure recovery

Correct Answer: B. introduces stabilizing artificial diffusion

Question 4

Which situation is most favorable for a classical spectral method?

A. A flow containing a strong shock and discontinuous material interface

B. A highly irregular industrial geometry with sharp corners

C. A smooth, well-resolved flow in a simple or periodic domain

D. A flow with unresolved contact discontinuities

Correct Answer: C. A smooth, well-resolved flow in a simple or periodic domain

Question 5

In a finite element weak formulation, integration by parts primarily

A. removes all boundary conditions

B. reduces the derivative order imposed on the trial solution

C. converts every element into a control volume

D. guarantees monotonicity

Correct Answer: B. reduces the derivative order imposed on the trial solution

Question 6

Which numerical behavior is most associated with an unbounded central scheme applied to an under-resolved convection-dominated problem?

A. Spurious oscillations

B. Exact local conservation in every formulation

C. Automatic mesh adaptation

D. Complete suppression of gradients

Correct Answer: A. Spurious oscillations

Question 7

On a severely skewed finite volume mesh, which improvement is most directly relevant to estimating spatial derivatives?

A. Ignoring neighboring cells

B. Using a suitable least-squares or corrected gradient reconstruction

C. Replacing convection with a source term

D. Reducing the physical Reynolds number

Correct Answer: B. Using a suitable least-squares or corrected gradient reconstruction

Question 8

A grid-refinement study is primarily used to determine whether

A. the governing equations contain time derivatives

B. the result is approaching a mesh-independent value

C. the fluid is incompressible

D. the solver uses a direct matrix inversion

Correct Answer: B. the result is approaching a mesh-independent value

Question 9

Which statement best describes numerical diffusion?

A. A physical property that is always independent of mesh size

B. Artificial smoothing that broadens gradients and is often increased by first-order upwinding

C. A phase error that necessarily creates oscillations

D. A boundary condition for heat transfer

Correct Answer: B. Artificial smoothing that broadens gradients and is often increased by first-order upwinding

Question 10

For a shock-containing flow, which combination is generally the most appropriate priority?

A. Global smoothness and unrestricted central differencing

B. Conservation, boundedness, and an appropriate shock-capturing flux

C. Very low polynomial order with no mesh refinement

D. Ignoring discontinuities during gradient reconstruction

Correct Answer: B. Conservation, boundedness, and an appropriate shock-capturing flux