Finite Volume Method -- general transport equation, Gauss theorem, discretisation schemes

General transport equation

All incompressible and compressible solvers derive from the general transport equation, describing conservation of an arbitrary scalar quantity phi per unit volume.

$$\underbrace{\frac{\partial}{\partial t}(\rho \phi)}_{\text{transient}}

• \underbrace{\nabla \cdot (\rho \mathbf{U} \phi)}_{\text{convective}}

= \underbrace{\nabla \cdot (\Gamma \nabla \phi)}_{\text{diffusive}}

• \underbrace{S_\phi}_{\text{source}}$$

Substitute phi with the physical quantity of interest to recover the governing equations.

All of CFD reduces to this single equation with different substitutions. The second-order spatial diffusive term means all discretisation approximations should be at least second-order accurate. First-order-upwind schemes aid robustness but systematically bias results with numerical diffusion. See Numerical Dispersion.

Gauss divergence theorem

The finite volume method starts by integrating the general transport equation over a control volume V bounded by a closed surface A:

$$\int_V \frac{\partial}{\partial t}(\rho \phi) \, dV

• \int_V \nabla \cdot (\rho \mathbf{U} \phi) \, dV

= \int_V \nabla \cdot (\Gamma \nabla \phi) \, dV

• \int_V S_\phi \, dV$$

The critical step applies Gauss's divergence theorem to convert volume integrals of divergences into surface integrals. After this step, the fundamental FVM statement becomes: rate of change within the cell plus convective flux across faces equals diffusive flux across faces plus source. Every finite-volume solver converts this balance. The platforms differ in how they compute face values phi_f and normal gradients nabla phi dot S, which is the entire craft of CFD discretisation.

Collocated vs staggered variable arrangement

The general transport equation treats phi as a continuous field. The discretisation must decide where scalar values are stored relative to the velocity components:

Collocated arrangement (OpenFOAM, ANSYS Fluent, Star-CCM+)

All variables are stored at the same cell-center location, offering data-structure simplicity and easy handling of complex geometry. The disadvantage is pressure-velocity checkerboard decoupling, which requires careful handling through Rhie-Chow interpolation and under-relaxation.

OpenFOAM resolves checkerboard mode through Rhie-Chow interpolation, which introduces a fourth-order diffusion term. This term suppresses oscillations without significant numerical diffusion, combined with under-relaxation and the SIMPLE/PISO/PIMPLE algorithm.

Staggered arrangement (older CFD codes, some incompressible solvers)

Velocity components reside at cell faces and pressure at cell centres, coupling pressure and velocity inherently without checkerboard mode. This arrangement complicates data structures and multi-domain coupling, so modern industrial codes rarely use it. The complexity outweighs the decoupling benefit.

OpenFOAM, Fluent, and Star-CCM+ all use collocated grids. The Rhie-Chow correction is effective. Data-structure simplicity is essential for complex-mesh applications. CFX takes a different approach with a node-centred arrangement, distinctive among commercial codes.

Convective flux interpolation

The convective term requires phi_f, the scalar value at each face, interpolated from neighbouring cell centres. This interpolation is the single most important approximation in FVM, balancing accuracy against stability.

First-order upwind is bounded, stable, and numerically diffusive. Central differencing is accurate but unbounded. Higher-order bounded schemes like TVD and QUICK offer both stability and accuracy, making scheme choice critical.

Linear (central differencing)

Linear interpolation is second-order accurate by Taylor series but produces non-physical oscillations near steep gradients, including undershoots below the minimum and overshoots above the maximum. The bounded convection scheme requires a monotonicity principle, covered in the TVD Limiters note.

First-order upwind

First-order upwind is fully bounded and stable, never creating new extrema. It is first-order accurate with leading truncation error as numerical diffusion proportional to grid spacing and absolute velocity. This false diffusion smears gradients and is the most ubiquitous source of error in CFD practice.

QUICK (Quadratic Upstream Interpolation for Convective Kinematics)

QUICK is a quadratic interpolation using upstream and downstream upstream nodes for second-order accuracy on monotonic profiles. Fluent used it widely in the 1990s and 2000s. TVD-limited schemes have superseded it in the modern canon.

TVD / Sweby-scheme family

TVD schemes combine second-order accuracy with first-order upwind stability. Limiter functions control behavior: second-order on smooth profiles, degraded to first-order locally at steep gradients or discontinuities to prevent oscillations. OpenFOAM implements vanLeer, minMod, SuperBee, and Venkatakrishnan limiters. Fluent also implements MUSCL and AUSM-family schemes for compressible flow. See Compression Schemes.

Diffusive flux and non-orthogonal correction

The diffusive term requires surface-normal gradient nabla phi dot S_f at each face. Perfectly orthogonal meshes allow direct computation from cell-centre values. Real meshes are not orthogonal, so all platforms use standard decomposition splitting the gradient into orthogonal and non-orthogonal correction terms.

Multiple non-orthogonal correctors iteratively reduce non-orthogonal error, controlled by nNonOrthogonalCorrectors in fvSolution or PIMPLE settings. Two to three correctors suffice for most meshes. Maximum non-orthogonality must stay below 70 degrees. Beyond that, it is a mesher problem, not a solver problem.

Gradient computation

Gradient computation is required by linearUpwind convection, non-orthogonal correction, and the diffusion term. Gauss theorem provides a straightforward formula but faces a bootstrapping problem: it needs face values, which linear interpolation needs. Upwind is the main alternative for gradient resolution.

Least-squares gradient reconstruction is more robust than Gauss linear on skewed or irregular meshes. It fits a local polynomial to cell centres without relying on face-area vectors. Star-CCM+ and Fluent default to least squares. OpenFOAM uses Gauss linear with cellLimited or cellMDLimited options.

Scheme comparison across platforms

All major platforms implement the same core discretisation choices with different default names. Fluent and Star-CCM+ default to second-order for most convection terms, while OpenFOAM requires explicit selection in fvSchemes. The Wolf Dynamics training enforces explicit defaults. Fluent users may not realize they run second-order until they need first-order for a diverging transient case. OpenFOAM users confront the accuracy-stability trade-off explicitly, a pedagogical advantage.

Temporal discretisation

The transient term uses time-stepping schemes. Euler implicit is unconditionally stable and the safe default. Backward is second-order and unconditionally stable. Crank-Nicolson is second-order but conditionally stable. Explicit Euler is rarely used in CFD.

Crank-Nicolson with blending factor theta in range 0 to 1 gives pure Euler at theta 0 and pure Crank-Nicolson at theta 1. The sweet spot for most transient RANS is theta 0.7 to 0.9 for nearly second-order accuracy with Euler damping preventing oscillation. Pure Crank-Nicolson can wiggle near steep gradients or discontinuities, a common pitfall in compressible flow and multiphase interface tracking.

Summary: the discretisation trade-space

Discretisation design reduces to three axes: accuracy (second-order versus first-order, noting that boundedness is stability with a safety net, not accuracy), stability (bounded TVD schemes versus unbounded central difference, since formal accuracy does not matter if the solution is invalid), and robustness (handling mesh distortion, steep gradients, and multi-dimensional flow, where least-squares gradients beat Gauss linear on bad meshes and vanLeer is more robust than SuperBee).

The professional discipline is to start robustly and finish accurately: use first-order on a hard case to get convergence, then restart with second-order schemes for the production answer. Decades of production CFD runs validate this strategy.

Cross-references

• TVD Limiters -- Sweby diagram, limiter functions, bounded high-order
• Numerical Dispersion -- truncation error origins
• Boundary Conditions -- how BCs interact with scheme selection
• OpenFOAM Numerics -- fvSchemes implementation, scheme sets
• OpenFOAM Turbulence -- turbulence transport equations and their schemes