TVD Limiters -- Sweby diagram, bounded high-order convection

The boundedness problem

Central-differencing (linear interpolation) is second-order accurate for convective flux. It is unbounded. Near steep gradients or discontinuities, it produces non-physical undershoots and overshoots. For a scalar phi that is between 0 and phi_{max}, the value should represent a physical quantity. Examples: mass fraction, temperature, void fraction, turbulent kinetic energy. An undershoot below zero is non-physical. An overshoot above phi_{max} is non-physical. Negative k breaks the k-epsilon model. Negative density breaks all equations that depend on it.

First-order upwind is monotonic. It never creates new extrema. If phi_P is between 0 and 1 for all cells P, then phi_f is between 0 and 1 for all faces f. The downside is first-order accuracy. The truncation error is O(Delta x). This is systematic numerical diffusion. It smears gradients regardless of mesh fineness. Second-order accuracy is required for fidelity. The TVD (Total Variation Diminishing) framework resolves this dichotomy.

Total variation diminishing

Harten formalized the TVD criterion in 1983. It formalizes what bounded and high-order actually means. The total variation of a discrete field phi in 1D is:

$$TV(phi) = sum_i |phi_{i+1} - phi_i|$$

A scheme is TVD if TV(phi^{n+1}) is less than or equal to TV(phi^n) for all n. The scheme does not create new local maxima or minima. It does not amplify existing extrema. This is equivalent to a discrete maximum principle. The concept extends naturally to higher dimensions.

The practical consequence: TVD schemes are adaptive. Where the solution is smooth (low variation), they operate at second-order accuracy. Where they detect a steep gradient or discontinuity (high variation), they automatically degrade to first-order locally. This prevents oscillations.

The Sweby diagram

Consider a 1D advection problem on a uniform grid. Upwind-biased second-order interpolation from upwind cell U, current cell C, and downstream cell D:

$$phi_f = phi_C + 1/2 psi(r) (phi_D - phi_{C'})$$

r equals (phi_C - phi_U) over (phi_D - phi_C) is the ratio of successive gradients. psi(r) is the limiter function. For central differencing (unbounded second-order), psi(r) equals 1 for all r. For first-order upwind, psi(r) equals 0. The two values are extremes. A TVD limiter blends between them.

The Sweby diagram (Sweby, 1986) defines the region in the (r, psi) plane where a scheme is TVD:

1. psi(r) greater than or equal to 0. The limiter never flips the sign of the corrected gradient. 2. psi(r) less than or equal to 1. Never amplifies beyond linear interpolation. 3. psi(r) less than or equal to 2r. TVD constraint on expansion waves. 4. psi(r) greater than or equal to r minus 1 on r greater than 1.

These four conditions define a bounded region. The TVD region is a diamond-shaped area with vertices at (0, 0), (1, 1), (2, 0), and the line psi equals 1 for 0 less than or equal to r less than or equal to 1. Central differencing with psi equals 1 lies on the top edge. It is bounded for monotonic profiles but unbounded at discontinuities. Upwind with psi equals 0 lies on the bottom edge. It is monotonic but diffusive. The goal is a limiter that hugs the line psi equals 1 where possible (second-order accuracy). It drops inside the TVD region near discontinuities.

Limiter functions

minMod

The name minMod is descriptive. It is the most diffusive TVD limiter:

$$psi(r) = max(0, min(1, r))$$

This produces a V-shaped curve from (0, 0) rising linearly to (1, 1), then horizontal at psi equals 1. The scheme stays completely within the TVD region. It is robust. It degrades to upwind quickly even for modestly varying gradients (r less than 1). It is the safe choice for highly discontinuous flows. It is too diffusive for wall-resolved RANS.

OpenFOAM: Gauss minMod. Fluent: minMod MUSCL.

van Leer

vanLeer is smooth and a good all-rounder. It is the classic TVD limiter:

$$psi(r) = (r + |r|) / (1 + r) = 0 for r less than or equal to 0, and 2r / (1 + r) for r greater than 0$$

This is continuous and differentiable everywhere in the TVD region (except at r equals 0 where it has a cusp). It hugs the upper bound psi equals 1 more tightly than minMod for r greater than or equal to 1 (compression waves, shocks). It still dampens at r less than 1 (expansion waves). It is the most widely used TVD limiter in industrial CFD.

OpenFOAM: Gauss vanLeer. Fluent: MUSCL vanLeer. Star-CCM+: MUSCL vanLeer (default).

SuperBee

SuperBee is the most compressive TVD limiter. It preserves sharp profiles and discontinuities better than van Leer. It can over-sharpen slightly:

$$psi(r) = max(0, min(1, 2r), min(2, r))$$

The Sweby plot shows a triangular region rising steeply from (0, 0) to (0.5, 1), then horizontal at psi equals 2 for 0.5 <= r <= 1. SuperBee compresses gradients instead of smearing them. It can make discontinuous profiles sharper than they should be. It can produce slight overshoots at shocks. It creates thin boundary layers. This is excellent for shock-capturing and interface tracking. It can cause non-physical overshoots in smooth regions near steep gradients.

OpenFOAM: Gauss SuperBee. Fluent: MUSCL SuperBee.

Venkatakrishnan

venkatakrishnan is a modern differentiable limiter. It avoids the cusp at r equals 1 found in van Leer. This makes it better suited for Newton-like iterative solvers:

$$psi(r) = (r + sign(r, epsilon)) / (1 + sign(r, epsilon))$$

sign(r, epsilon) equals (r^2 + epsilon r) over (r^2 + epsilon). epsilon is a small parameter typically between 10^{-6} and 10^{-3}. As epsilon goes to 0, this converges to van Leer. For finite epsilon, it is smooth everywhere. The smoothness benefits the inner Newton iteration loops of coupled solvers. It does not appear in base OpenFOAM. It is available in newer versions and in the foam-extend branch.

Comparing limiters: what they do in practice

The following qualitative comparison comes from standard test cases from the Wolf Dynamics FVM training materials and standard benchmarks (Sod shock tube, advection of a sine wave, Taylor-Green vortex).

The recommendation from Wolf Dynamics training is:

• Default for RANS: vanLeer or venkatakrishnan. Good accuracy with low oscillation risk. Works well on polyhedral or hex-dominant meshes.
• Compressible and shock flows: SuperBee for shock resolution. Pair it with a compression scheme (Riemann solvers for density-based solvers).
• Getting a diverging run to stabilize: switch to minMod or upwind. Achieve convergence. Then switch back to vanLeer for the production run. This is the start robust, finish accurate discipline.

CFL constraints and limiters

A TVD limiter is necessary but not sufficient for boundedness. The time step must also satisfy the Courant-Friedrichs-Lewy (CFL) condition:

$$Co = |U| delta_t / delta_x \in Co_{max}$$

For explicit time stepping, Co_{max} is less than or equal to 1. The limit can be tighter depending on the scheme. For implicit time stepping (as used in OpenFOAM PIMPLE and SIMPLE solvers), Co_{max} can exceed 1. Boundedness is not guaranteed at high Co regardless of the limiter. The practical limit for boundedness with TVD schemes is Co about 2 to 5 for PIMPLE with correctors. The limit depends on mesh and flow regime. For pure PISO, maintain Co less than 1.

See the OpenFOAM Numerics note for implementation details on maxCo in controlDict. See the Pressure-velocity coupling note for the relationship between corrector loops and permissible CFL numbers.

Barth-Jespersen limiters

The Barth-Jespersen (1987) approach limits solution values instead of flux reconstruction. It ensures that phi_P at any cell P never exceeds the maximum or falls below the minimum of phi among all neighboring cells:

$$phi_P^{new} = min(phi_{P,max}, max(phi_{P,min}, phi_P^{high-order}))$$

phi_{P,min} and phi_{P,max} are taken from the stencil (neighbors plus current). This is a cell-wise limiter. It operates after the high-order solution is computed. It clamps extrema.

Barth-Jespersen is less aggressive than flux limiters. It preserves more of the high-order solution. It is not TVD in the strict Sweby sense. Fluent uses it as its default bounded second-order scheme. OpenFOAM provides it via the barthJespersenLimited gradient limiter. Star-CCM+ uses a similar cell-wise limiting strategy.

Implementation in OpenFOAM

In OpenFOAM, TVD convection schemes are specified in fvSchemes:

divSchemes
{
    default         none;
    div(phi,U)      Gauss vanLeer;
    div(phi,k)      Gauss limitedLinear 1;
    div(phi,omega)  Gauss SuperBee;
}

The limitedLinear scheme requires a second parameter. The parameter is the limiter coefficient psi in range 0 to 1. This coefficient controls the degree of limitation. limitedLinear 1 is full linear with no limitation. limitedLinear 0 is full upwind. limitedLinear 0.5 is half-linear and half-upwind. This is a practical compromise.

vanLeer, minMod, and SuperBee take no parameters. The limiter function is fully specified. venkatakrishnan takes an optional epsilon parameter.

default none is a deliberate discipline from Wolf Dynamics training. It forces the user to explicitly select a scheme for every convective term. If a term is missing, the solver errors out instead of silently using an inappropriate default. OpenFOAM has no default discretization. Nothing is more dangerous in CFD than an unexamined default.

Limiters and turbulence equations

Turbulence transport equations (k, epsilon, omega, nu_t) are particularly sensitive to the convection scheme. Negative k or omega is non-physical and breaks turbulence models immediately.

• k less than or equal to 0 means sqrt(k) equals 0 or is imaginary. Turbulent viscosity breaks.
• omega less than or equal to 0 means epsilon equals beta^* k omega less than or equal to 0. Dissipation becomes negative. This causes exponential blow-up.

For turbulence equations, prefer schemes that are guaranteed bounded. vanLeer, minMod, or limitedLinear with coefficient less than or equal to 0.5 are good choices. Avoid pure linear (central differencing) on turbulence equations unless the mesh is excellent (orthogonality less than 50 degrees, skewness less than 1.0) and the flow is mildly varying. This is not a preference. The solver will crash if omega goes negative. This happens regardless of whether you are solving RANS or LES.

See the OpenFOAM Turbulence Modelling note for turbulence model configuration. See the OpenFOAM Numerics note for specific fvSchemes recommendations for RANS versus LES.

Cross-references

• FVM Overview -- convective flux, collocated grids
• Numerical Dispersion -- truncation error from upwind
• Compression schemes -- TVD limiters in compressible flow
• OpenFOAM Numerics -- fvSchemes implementation details
• OpenFOAM Turbulence -- why turbulence equations need bounded schemes
• Sweby, P. K. (1986). Towards the Limit Convergence of Limiter Functions. SIAM Journal on Numerical Analysis, 23(4), 756-765.
• Harten, A. (1983). High Resolution Schemes for Hyperbolic Conservation Laws. Journal of Computational Physics, 49(3), 357-393.
• Barth, T. J., and Jespersen, D. C. (1987). The Design and Application of Optimized Monotone Reconstruction Schemes. AIAA Paper 83-1924.