OpenFOAM Numerics and Discretization

Overview

Numerics in OpenFOAM are configured entirely through system/fvSchemes (discretization schemes) and system/fvSolution (linear solver settings, under-relaxation, coupling parameters). There is no GUI to toggle discretization — every scheme choice is a line in a text file, and every line is a decision that will either produce physical results or produce results that look physical. The distinction often comes down to whether your divergence scheme satisfies the boundedness condition on a face-by-face basis.

See Numerical Schemes for the general FVM discretization theory underlying fvSchemes, Numerical Dissipation for the relationship between scheme choice and artificial viscosity, and Case Setup for how discretization interacts with boundary conditions and mesh quality.

The fvSchemes Configuration File

The fvSchemes file specifies discretization schemes for edifferential operator appearing in the governing equations. It is organized into categories matching the mathematical operators:

Each category contains one or more scheme lines in the format:

divSchemes
{
    default         Gauss linear;
    div(phi,U)      Gauss linearUpwind grad(U);
}

The general form is Gauss <faceScheme> [optionalCorrection]: the integration method (Gauss, upwind), the face reconstruction scheme, and optionally a gradient-based correction.

Time Discretization (ddtSchemes)

The CrankNicolson scheme uses $\delta \in [0, 1]$ to blend between first-order upwind ($\delta=0$) and second-order centered ($\delta=1$). The default crankNicolson 0.9 in most cases is empirically justified: $\delta=1$ can produce oscillatory solutions for advective-dominated problems, while $\delta=0.9$ adds just enough numerical diffusion to maintain stability without sacrificing the second-order accuracy on smooth solutions.

Convection Discretization (divSchemes)

The convection scheme is the most consequential choice in fvSchemes. It determines how the convective term $\nabla \cdot (\phi \mathbf{U})$ is discretized on each face of the finite volume.

The TVD (Total Variation Diminishing) schemes use a Sweby limiter function. For a cell-centered scheme with face value $\phi_f$, the TVD formulation uses the ratio of successive differences:

$$r = \frac{\phi'_{upstream}}{\phi'_{downstream}} = \frac{\phi_C - \phi_U}{\phi_D - \phi_C}$$

The Sweby limiter function $\psi(r)$ constrains the slope limiter so that the scheme satisfies the TVD condition:

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

This is the classic Sweby TVD region ($0 \le \psi \le 2$ for $0 \le r \le 2$, and $\psi = 0$ otherwise) — a sufficient condition for monotonicity preservation. vanLeer gives the smoothest interpolation ($\psi(r) = \max(r, 1)$); Minmod is the most diffusive ($\psi(r) = \max(0, r)$); SuperBee is the least diffusive ($\psi(r) = \max(0, \min(2r, 1), \min(r, 2))$).

Gradient Schemes (gradSchemes)

Gradient computation is needed for second-order convection schemes, diffusion terms, and wall function formulations. The choice affects both accuracy and stability.

For orthogonal meshes, Gauss linear gives the exact gradient (the discrete divergence theorem is exact for linear fields on orthogonal grids). For non-orthogonal meshes, the non-orthogonal correction must be applied, which is why Gauss linear corrected is the default in most cases: it splits the gradient into orthogonal and non-orthogonal parts:

$$\int_{\partial V} \phi \mathbf{n} \, dA = \sum_{f \in \partial V} \phi_f \mathbf{S}_f = \sum_f \left( \phi_C + (\nabla \phi)_C \cdot \mathbf{d}_{CF} + (\nabla \phi)_N \right) \mathbf{S}_f$$

The non-orthogonal correction is iterated over in the linear solver (controlled by nNonOrthogonalCorrectors in fvSolution).

Laplacian Schemes (laplacianSchemes)

The Laplacian operator $\nabla \cdot (\Gamma \nabla \phi)$ appears in diffusion terms (viscous stress, turbulent diffusion, heat conduction).

The corrected variant handles non-orthogonal meshes by iterating on the non-orthogonal correction term. On a highly non-orthogonal mesh, the =uncorrected~ variant reduces accuracy to first order. The choice between corrected and limited is: limited applies a limiter to the non-orthogonal correction itself (preventing blow-up), while corrected iterates to convergence on the correction.

Interpolation Schemes (interpolationSchemes)

Face values are obtained by interpolating from cell centers:

CFL Number Theory

The Courant-Friedrichs-Lewy (CFL) condition governs the maximum stable time step for explicit (and, more weakly, implicit) time integration. The local cell Courant number is:

$$Co = \frac{|\mathbf{U}| \Delta t}{\Delta x}$$

and the global (maximum) Courant number is:

$$Co_{max} = \max_i \frac{|\mathbf{U}_i| \Delta t}{\Delta x_i}$$

For explicit schemes, $Co \le 1$ is necessary for stability. For implicit schemes (SIMPLE, PISO, PIMPLE), the CFL can be significantly larger (up to $Co \approx 100$ for steady-state RANS, $Co \approx 1-5$ for transient LES/DNS).

The =controlDict~ specifies the maximum Courant number for the solver to use:

// controlDict - dynamic time stepping
adjustTimeStep  yes;
maxCo           0.5;
maxDeltaT       1e-3;

This allows the solver to automatically reduce the time step when the local Courant number exceeds maxCo. For transient RANS (PIMPLE), a common setting is =maxCo 0.5~. For LES, more conservative (=maxCo 0.3~ or even lower). For steady-state RANS (SIMPLE), =maxCo~ is irrelevant because the time derivative is discarded.

Recommended Scheme Sets

Three typical scheme configurations for fvSchemes. Use the appropriate set based on your simulation type.

Accurate Set (for final production runs)

ddtSchemes
{
    default         CrankNicolson 0.9;
}
divSchemes
{
    default         none;
    div(phi,U)      Gauss linearUpwind grad(U);
    div(phi,k)      Gauss limitedLinear 1;
    div(phi,epsilon) Gauss limitedLinear 1;
    div(phi,omega)  Gauss limitedLinear 1;
    div(phi,nuTilda) Gauss limitedLinear 1;
    div(phi,alpha)  Gauss interfaceCompression 1;
    div(phir,alpha) Gauss interfaceCompression 1;
    div(phiTheta)   Gauss limitedLinear 1;
    div(phi,Theta)  Gauss limitedLinear 1;
}
laplacianSchemes
{
    default         none;
    laplacian(nu,U) Gauss linear corrected;
    laplacian(nuTilda,U) Gauss linear corrected;
    laplacian(1*alpha,alphacf) Gauss linear corrected;
    laplacian(Dp,D) Gauss linear corrected;
    laplacian(Dk,k) Gauss linear corrected;
    laplacian(Depsilon,epsilon) Gauss linear corrected;
    laplacian(Domega,omega) Gauss linear corrected;
    laplacian(DTheta,Theta) Gauss linear corrected;
    laplacian(Dphi,phi) Gauss linear corrected;
    laplacian(lam,T) Gauss linear corrected;
}
gradSchemes
{
    default         Gauss linear corrected;
}
snGradSchemes
{
    default         corrected;
}
interpolationSchemes
{
    default         linear;
}
snGradRemapSchemes
{
    default         leastSquares;
}

Balanced Set (for convergence to stable solution)

ddtSchemes
{
    default         Euler;
}
divSchemes
{
    default         none;
    div(phi,U)      Gauss linearUpwind grad(U);
    div(phi,phi)    Gauss linear;
    div(phi,k)      Gauss linear;
    div(phi,epsilon) Gauss linear;
    div(phi,omega)  Gauss linear;
    div(phi,nuTilda) Gauss linear;
    div(phi,alpha)  Gauss upwind;
}
laplacianSchemes
{
    default         none;
    laplacian(nu,U) Gauss linear corrected;
    laplacian(Dk,k) Gauss linear corrected;
}
gradSchemes
{
    default         Gauss linear corrected;
}
snGradSchemes
{
    default         corrected;
}
interpolationSchemes
{
    default         linear;
}
snGradRemapSchemes
{
    default         leastSquares;
}

Robust Set (for challenging cases)

ddtSchemes
{
    default         Euler;
}
divSchemes
{
    default         none;
    div(phi,U)      Gauss upwind;
    div(phi,phi)    Gauss upwind;
    div(phi,k)      Gauss upwind;
    div(phi,epsilon) Gauss upwind;
    div(phi,omega)  Gauss upwind;
    div(phi,nuTilda) Gauss upwind;
    div(phi,alpha)  Gauss upwind;
}
laplacianSchemes
{
    default         none;
    laplacian(nu,U) Gauss linear corrected;
}
gradSchemes
{
    default         Gauss linear corrected;
}
snGradSchemes
{
    default         corrected;
}
interpolationSchemes
{
    default         linear;
}
snGradRemapSchemes
{
    default         linear bounded Gauss 0.33;
}

ANSYS Comparison: Discretization Schemes

Fluent's discretization schemes map to OpenFOAM fvSchemes keywords as follows:

Fluent's PISO-based "Coupled" pressure-based solver is conceptually closer to OpenFOAM's PIMPLE (PISO + SIMPLE under-relaxation) than to Fluent's Segregated PISO. Fluent's density-based solver maps to OpenFOAM's =rhoSimpleFoam~ / =rhoPimpleFoam~ / =rhoCentralFoam~ family.

fvOptions and Source Terms vs ANSYS

OpenFOAM's =fvOptions~ framework provides a unified interface for adding source terms to any equation. This is the OpenFOAM analog of ANSYS Fluent's Source Terms panel (momentum sources, viscous dissipation, species source terms, porous media resistance, etc.); both serve the same purpose but differ in implementation philosophy.

The =fvOptions~ format in system/fvOptions:

// system/fvOptions
darcyForchheimer
{
    type            darcyForchheimer;
    selectionMode   cellZone;
    cellZone        porousZone;
    d               [0 0 -1 0 0 0 0] (1e7 1e7 1e7);
    f               [0 1 -2 0 0 0 0] (1e4 1e4 1e4);
    active          yes;
}

This adds a momentum sink term $S_i = -D_{ij} U_j - \frac{1}{2} F_{ij} |U_j| U_j$ to the momentum equation, exactly analogous to Fluent's =viscous-resistance~ and =inertial-resistance~ parameters in the porous media zone definition.

The coded option allows inline C++ for bespoke source terms:

codedFvOptions
{
    type            codedFvOption;
    ...
    codeOptions
    {
        "-I../lnInclude";
    }
    libs            ["libmyFoamOptions.so"];
    name            mySource;
    codeCorrect     "#{codeContent;}";
    codeAddSupply   "#{codeContent;}";
}

This is OpenFOAM's approach to what ANSYS handles via UDFs (DEFINE_SOURCE, DEFINE_MZMS etc.). The key difference: OpenFOAM compiles the code into the solver at runtime (via the C++ preprocessor macros in the =code~ block), while Fluent UDFs are compiled as shared libraries loaded at runtime. Both achieve the same goal — custom source terms — but the OpenFOAM approach integrates more deeply with the solver's data structures.

Cross-Links

References

• Weller, H. G. et al. (1998). "A tensorial approach to computational hydrodynamics." Computers in Physics, 12(6), 620-631.
• Hirsch, C. (2007). Numerical Computation of Internal and External Flows (2nd ed.). Wiley.
• Versteeg, H. & Malalasekera, W. (2007). An Introduction to Computational Fluid Dynamics: The Finite Volume Method (2nd ed.). Pearson.
• OpenFOAM User Guide: fvSchemes and fvSolution. openfoam.com/documentation/guide.
• OpenFOAM Programming Guide: fvOptions framework. openfoam.com/documentation/guide.
• ANSYS Fluent Theory Guide, Release 2024 R2. ANSYS Inc.
• OpenFOAM Foundation v14 Modular Solver Documentation. foam-extend.org.