Numerical Dispersion and Diffusion -- truncation error, artificial viscosity

The two sins of convection discretisation

When you approximate the convective flux at a face using a finite-order scheme, you introduce a discretisation error. This error is the difference between the exact infinite-order flux and your truncated approximation. The error has two components:

1. Numerical diffusion (also called artificial viscosity or truncation error). This is a smearing of gradients. It is always present in upwind or lower-order schemes. 2. Numerical dispersion. This is an oscillatory error. It produces spurious waves (overshoots and undershoots) near steep gradients.

The name dispersion comes from wave mechanics. A true solution to the advection equation u_t plus a u_x equals zero preserves the shape of an initial profile. A discretised solution with dispersion errors produces spurious oscillations. These oscillations disperse away from the true solution. They look like ripples behind a boat.

Taylor series error expansion

Consider first-order upwind advection in 1D on a uniform grid with constant Delta x. Flow moves from left to right.

$$phi_f approx phi_P$$

The exact face value by Taylor series is:

$$phi_f = phi_P + (Delta x / 2) partial phi partial x\big|_P + (Delta x^2 / 8) partial^2 phi partial x^2\big|_P + O(Delta x^3)$$

The upwind scheme truncates after the zeroth order. The leading truncation error is:

$$epsilon_{trunc} = (Delta x / 2) partial phi partial x + (Delta x^2 / 8) partial^2 phi partial x^2 + O(Delta x^3)$$

For the convective term nabla dot (U phi), substituting the upwind approximation into the exact expression reveals that the leading truncation term equals a diffusion term. The diffusion has magnitude:

$$Gamma_{num} = |U| Delta x / 2$$

This is artificial viscosity. It is not physical. It depends only on grid spacing and velocity. It does not depend on fluid properties.

For a second-order (central-difference or linear) scheme, the leading error is O(Delta x^2). The error has the form of a dispersion term (third derivative):

$$epsilon_{trunc}^{linear} = (Delta x^2 / 6) partial^3 phi partial x^3 + O(Delta x^4)$$

This odd-derivative term produces oscillations instead of smoothing. Linear (upwind minus upwind) interpolation is non-diffusive by construction. Its error is purely dispersive. This is why it produces wiggles near gradients instead of smearing them.

The lesson: first-order schemes smear. Second-order schemes wiggle. The TVD framework seeks a combination of both. It provides second-order accuracy on smooth profiles for low dispersion error. It uses first-order behavior at discontinuities to prevent oscillations.

Richardson extrapolation and mesh sensitivity

Richardson extrapolation is the standard technique for estimating discretisation error from a discrete solution. Given a solution on three uniformly refined grids (coarse H, medium H/r, fine H/r^2) with refinement ratio r (typically r equals 2):

$$p = (1 / ln r) ln | (f_H - f_m) / (f_m - f_f) |$$

The value p is the observed order of accuracy. The values f_H, f_m, f_f are quantities of interest (drag coefficient, pressure coefficient at a probe) on the three grids. The extrapolated infinite-refinement value is:

$$f_{ex} = f_f + (f_f - f_m) / (r^p - 1)$$

The discretisation error estimate on the fine grid is:

$$epsilon_f = f_{ex} - f_f$$

The discretisation uncertainty (95% confidence) is:

$$U_f = 1.25 epsilon_f$$

The factor 1.25 is the grid-convergence index factor for three grids.

Practical mesh sensitivity guideline

For RANS simulations using second-order schemes:

The AIAA Committee on Standards recommends computing drag and lift on at least three grids with r equal 1.3 to 1.5. Report p and the Grid Convergence Index. If p is less than 1.5, the solution is not in the asymptotic range. The error estimate is unreliable. The mesh is still too coarse.

Numerical diffusion in practice

Numerical diffusion is the single most common source of error in engineering CFD. It manifests as:

• Overprediction of mixing: numerical smearing of species or temperature interfaces makes mixing appear faster than it physically is.
• Underprediction of shear layer growth: artificial viscosity thickens shear layers before the physical instability has grown.
• Overprediction of boundary layer thickness: artificial diffusion adds momentum to the near-wall region. The boundary layer appears mushier than it is.
• Suppression of small-scale structures: LES and DNS are particularly sensitive. The subgrid model expects the physical dissipation rate. It does not expect physical plus numerical dissipation.

Numerical diffusion magnitude is proportional to |U| Delta x over 2 for first-order schemes. It is proportional to O(Delta x^2) for second-order schemes. This means:

1. First-order schemes produce numerically-dominant diffusion on meshes that would otherwise be fine enough for second-order accuracy. 2. Near walls where Delta x is small, numerical diffusion is less significant. Near walls, the y+ discipline and log-law wall function matter more. 3. In the far-field where Delta x is large, numerical diffusion from coarse meshes can swamp physical diffusion entirely.

Mesh density far from regions of interest still matters. Numerical diffusion from a coarse far-field can propagate into the region of interest. This effect occurs in recirculating flows and internal flows.

Numerical dispersion in practice

Numerical dispersion is the oscillatory error of central-differencing schemes. It is less common than numerical diffusion. It is more dramatic when it occurs. It manifests as:

• Pressure oscillations near shocks or sharp gradients. This is the classic wiggly pressure that can crash a solver with negative pressures.
• Overshoots in scalar fields. Mass fractions exceed 100%. Temperature drops below absolute zero. Turbulent kinetic energy becomes negative.
• Spurious vorticity in shear layers where the grid resolves the kinks in the advected profile.

Numerical dispersion is proportional to Delta x^2 for central differencing. It is proportional to the limiter function value for TVD schemes.

Strategies to control discretisation error

Recommended default strategy for production CFD

1. Mesh: verify grid convergence on at least three grids using the GCI method 2. Scheme: second-order linear or TVD (vanLeer) for all convection terms. Never use first-order for the production run. 3. Non-orthogonal correction: use at least one non-orthogonal corrector for meshes with maximum non-orthogonality greater than 30 degrees. 4. Time stepping: verify time-step independence separately from mesh independence. Use delta t to delta t/2 refinement. 5. Report: always report mesh resolution (total cells, y+ range, max non-orthogonality), scheme order, and grid-convergence order of accuracy p.

See the FVM Overview note for the discretisation scheme taxonomy. See the OpenFOAM Numerics note for three ready-to-use fvSchemes sets (accurate, balanced, robust).

Cross-references

• FVM Overview -- discretisation error from convection and diffusion
• TVD Limiters -- how TVD reduces dispersion error
• Compression Schemes -- shock smearing as numerical diffusion
• Mesh Quality -- how mesh quality affects truncation error
• AIAA Committee on Standards (2009). Verification and Validation Terminology, Definitions, and Guidelines. AIAA G-077-1998 (R2009).
• Roache, P. J. (1998). Verification and Validation in Computational Science and Engineering. Hermosa Publishers.