OpenFOAM Numerical Schemes — fvSchemes Deep Dive

#+CATEGORY: openfoam # :PROPERTIES: # :ID: uuid-openfoam-numerical-schemes # :END:

OpenFOAM Numerical Schemes :: fvSchemes Configuration Reference

This note provides a reference-grade walkthrough of OpenFOAM's discretisation schemes in fvSchemes. Every scheme listed here corresponds to theory documented in Numerical Schemes, TVD Limiters, Compression Schemes, and Numerical Dissipation.

The fvSchemes file is where you translate numerical theory into solver configuration. Fluent achieves the same effect through its Solution Methods panel — but Fluent abstracts away the mathematical detail. OpenFOAM exposes every choice, which is both its greatest strength and its most common source of user errors.

**Table of Contents**

• ddtSchemes :: Temporal discretisation
• gradSchemes :: Cell and face gradients
• divSchemes :: Convective fluxes (the most critical choices)
• laplacianSchemes :: Diffusive fluxes
• snGradSchemes :: Face-normal gradients for non-orthogonal correction
• interpolationSchemes :: Face value interpolation
• TVD Schemes :: Sweby-limited second-order
• Finding Keywords :: Auto-generate fvSchemes
• See Also
• References

Time Discretisation (ddtSchemes) :: From First-Order Euler to Crank-Nicolson

The choice of temporal scheme determines how accurately the transient term =\partial \phi / \partial t/ is approximated. First-order is stable but introduces numerical diffusion in time. Second-order (Crank-Nicolson) is accurate but can oscillate for discontinuous initial conditions.

ddtSchemes
{
    default         Euler;           // 1st-order implicit — default for steady-state
    // for transient:
    default         CrankNicolson 0.9;  // 2nd-order with 10% dissipation
}

The Crank-Nicolson scheme with blending factor =\alpha \in [0, 1]/:

= \frac{\partial \phi}{\partial t} \approx (1-\alpha) \frac{\phi^{n+1} - \phi^n}{\Delta t} + \alpha \frac{\phi^{n+1} - \phi^{n-1}}{2\Delta t}

For =\alpha = 0/ (pure Euler, 1st-order): unconditionally stable, numerical dissipation =O(\Delta t)/. For =\alpha = 1/ (pure Crank-Nicolson, 2nd-order): bounded for smooth solutions, oscillates for discontinuities. The typical value =\alpha = 0.9/ provides 90% second-order accuracy and 10% numerical dissipation (preventing oscillations).

For LES, Crank-Nicolson is standard (2nd-order temporal accuracy minimises energy cascade errors). For transient RANS (e.g., vortex shedding), Euler or CrankNicolson 0.9 are both acceptable — the turbulence model's inherent dissipation dominates. For DNS (direct numerical simulation), CrankNicolson is mandatory — any temporal dissipation contaminates the smallest scales.

Gradient Schemes (gradSchemes) :: Gauss, Least-Squares, and Green-Gauss

The gradient scheme computes =\nabla \phi/ at cell centres from cell-centred values (collocated grid). Two options:

gradSchemes
{
    default         Gauss linear;      // Green-Gauss + linear interpolation
    // or: leastSquares      // cell-based least-squares (unlimited stencil)
    grad(p)         Gauss linear corrected;     // non-orthogonal correction
    grad(U)         Gauss linear corrected;
}

Green-Gauss gradient: =\nabla \phi_P = \frac{1}{V_P} \sum_f \phi_f \hat{n}_f A_f/. The face value =\phi_f = (\phi_P + \phi_N) / 2/ is linear interpolation from the two neighbouring cell centres. The correction term adds: =(\nabla \phi)_f = ( \phi_f - \phi_{H P}) / |\vec{d}_{P \to f}|/ where =\phi_{HP}/ is the projected value along the cell-centre line (the orthogonal component removed).

Least-squares gradient: minimises =\sum_i w_i (\phi_i - \phi_P - \nabla \phi_P \cdot \vec{r}_{iP})^2/ over all neighbours i of cell P. This is more robust on highly skewed meshes because it uses the actual neighbour geometry rather than the orthogonal-correction approximation. However, it is more expensive (larger stencil, more arithmetic).

For quality meshes (orthogonality > 0.8), Gauss linear corrected is standard. For poor-quality meshes (orthogonality < 0.6, high skewness), leastSquares is more robust but slower. The choice between the two does not significantly affect solution accuracy when the mesh is good — the gradient error is sub-dominant to other discretisation errors in that case.

Divergence (divSchemes) :: Convective Fluxes

This is where the discretisation philosophy matters most. The convective term is =\nabla \cdot (\vec{v} \phi)/, evaluated at cell faces as =\sum_f \rho_f \phi_f (\vec{v}_f \cdot \hat{n}_f) A_f/. The key question is: how is the face value =\phi_f reconstructed?

divSchemes
{
    default         none;
    // Momentum — 2nd-order upwind-biased
    div(phi,U)      bounded Gauss linearUpwind Grad U;

    // Passive scalars — 1st-order upwind (more stable)
    div(phi,k)      bounded Gauss upwind;

    // VOF — MULES with interface compression
    div(phi,alpha)  MULES interIsoAdvector;

    // Viscous term
    div((nuEff*dev2(T(grad(U)))))    Gauss linear;

    // Custom TVD
    div(phi,psi)    bounded Gauss TVD linear;
}

The standard recipe for RANS:

Momentum: linearUpwind (2nd-order) — this is critical. 1st-order upwind for momentum artificially diffuses the momentum boundary layer, leading to incorrect separation points and incorrect drag coefficients.

Turbulent scalars (k, omega, epsilon): upwind (1st-order) — turbulent scalars are very sensitive to numerical dissipation. 1st-order is more stable but under-predicts turbulent kinetic energy and over-predicts dissipation (in high-gradient regions).

Passive scalars (temperature, species): linear or TVD — temperature is passively transported by the flow, so numerical dissipation can smear temperature gradients. Use TVD for sharp gradients (e.g., flame fronts). Linear central differences for smooth scalars.

Volume fraction (VOF): MULES — strictly bounded, multidimensionally limited. Without MULES, the VOF field becomes unbounded.

The linearUpwind scheme uses the gradient (from the gradSchemes file) to extrapolate from cell-centre to face value:

= \phi_f = \phi_P + (\vec{r}_f - \vec{r}_P) \cdot \nabla \phi_P /

In the upwind direction (direction of flow). This is a second-order upwind-biased scheme: it uses information from the upwind direction (like upwind) but with a second-order correction (like central differencing).

For TVD schemes, the flux is corrected by a Sweby limiter function =\phi(r)/ to ensure boundedness:

= \phi_f = \phi_{\text{minmod}} + \phi(r) (\phi_{\text{central}} - \phi_{\text{minmod}}) /

where =r/ is the ratio of successive gradients. The Sweby limiting ensures =0 \leq \phi(r) \leq 2/ in the upwind direction. The TVD limiters are described in TVD Limiters, with the common ones: minmod, vanLeer, superbee, venkatakrishnan (differentiable, for Newton iteration).

For compressible flow (rhoPimpleFoam, sonicFoam), the convective fluxes use limitedLinear or MUSCL — schemes designed for shock-capturing. This maps to Compression Schemes.

Laplacian (laplacianSchemes) :: Diffusive Fluxes

Diffusion terms are =\nabla \cdot (\Gamma \nabla \phi)/, evaluated at faces as =\sum_f \Gamma_f (\nabla \phi \cdot \hat{n}_f) A_f/.

laplacianSchemes
{
    default         none;
    laplacian(nu,U)         Gauss linear corrected;      // viscous stress
    laplacian(nuEff,U)      Gauss linear corrected;      // turbulent + laminar
    laplacian(1,rho,p)      Gauss linear corrected;       // pressure diffusion
    laplacian(Dt,T)         Gauss linear corrected;       // thermal
    laplacian(D,T)          Gauss linear corrected;       // species
}

The linear scheme computes the face gradient from the two neighbouring cell centres. The "corrected" option adds the non-orthogonal correction (described above under gradSchemes). For meshes with orthogonality > 0.9, the correction is small; for meshes with orthogonality ~ 0.5–0.7, the correction is critical.

For highly non-orthogonal meshes (skewness > 0.9), you can increase the number of non-orthogonal correctors in fvSolution (see Pressure-Velocity Coupling for how this affects the SIMPLE algorithm's non-orthogonal correction loop):

PIMPLE
{
    nNonOrthogonalCorrectors  2;     // default is 0; 1–2 for non-orthogonal meshes
}

Each non-orthogonal corrector re-solves the laplacian with updated non-orthogonal corrections. Zero correctors (default) means the face gradient uses only the orthogonal component — equivalent to treating the non-orthogonal part as a source term. One corrector applies the full non-orthogonal correction in one step; two correctors iteratively improve the correction. For most snappyHexMesh-generated polyhedral meshes, 1–2 correctors are needed.

Interpolation Schemes :: Face Value Reconstruction

The interpolation scheme controls how cell-centred values are interpolated to cell faces:

interpolationSchemes
{
    default         linear;
    interpolation(p)        linear;
    interpolation(U)        linear;
    interpolation(rho)      linear;
}

Linear interpolation is standard: =\phi_f = (\phi_P + \phi_N) / 2/. Upwind interpolation: =\phi_f = \phi_U/ (upwind cell value). Upwind is only used when combined with a TVD scheme (for the convective flux) — it defines the zero-order reconstruction that TVD blends with the central linear interpolation.

TVD Schemes :: Sweby-Limited Second-Order

OpenFOAM implements TVD through the =TVD/ scheme wrapper:

divSchemes
{
    div(phi,psi)    bounded Gauss TVD vanLeer;   // bounded, 2nd-order
}

Supported limiters: vanLeer (smooth, 2nd-order), superbee (sharp, but non-differentiable), minmod (diffusive, very safe), venkatakrishnan (differentiable — for Newton-Raphson), SMART (hybrid van Leer/minmod). The Sweby diagram from TVD Limiters applies.

The venkatakrishvan limiter is unique because it introduces a small artificial smoothness parameter =\epsilon/ to make the limiter function differentiable:

= \phi_{vkr}(r) = \frac{(r^2 + \epsilon + r\epsilon) \cdot |r - 1| + (2r - r^2 - 1)\epsilon + 2r\epsilon + r^2\epsilon^2}{(r+1)^2 + \epsilon + 2\epsilon(r+1) + \epsilon^2}

This is critical for implicit solvers using Newton-Raphson: the limiter must be differentiable or the Jacobian is discontinuous and Newton's method will not converge reliably. superbee is non-differentiable at r = 1 and r = 2 — it causes convergence problems in implicit solvers. vanLeer and venkatakrishvan are safe.

Banana Method for fvSchemes Keywords

OpenFOAM includes a utility that auto-generates fvSchemes keyword lists for all fields in your case:

banana 0          # list all keyword entries for fields in 0/
banana 0 -print   # print fvSchemes entries for all fields in 0/
banana            # list all keywords found in 0/ across all files

The output is a complete set of fvSchemes divSchemes and laplacianSchemes entries for every field in your case. For example, if you add a new scalar =T/ to your case, =banana 0 -print/ outputs:

div(SomethingT)    bounded Gauss linearUpwind Grad T;
laplacian(Dt,T)    Gauss linear corrected;
grad(T)            Gauss linear corrected;
interpolation(T)   linear;

This is the fastest way to discover which keywords you need. It works by scanning the source code for all fields and identifying all divergence and laplacian operators used. If banana reports entries for fields you don't recognise, check your solver source code to identify which operators are activated.

See Also :: Related Notes

. Numerical Schemes — upstream theory . TVD Limiters — Sweby diagram, limiter functions . Compression Schemes — shock-capturing schemes . OpenFOAM Case Setup — fvSchemes in context . OpenFOAM Solver Selection — which solver needs which schemes . ANSYS Fluent Overview — Fluent's scheme selection comparison

References

. OpenFOAM User Guide. ESI OpenCFD. (fvSchemes reference.) . OpenFOAM Programming Guide. ESI OpenCFD. (TVD scheme implementation.) . Versteeg, H.K. & Malalasekera, W. (2007). An Introduction to CFD. Chapter 5 (Discretisation schemes). . Hirsch, C. (1989). Numerical Computation of Internal and External Flows. Vol. 1, Chapter 11.