Lesson 4 of 13 25 min

Finite Volume Method

The Finite Volume Method (FVM) is the backbone of commercial CFD. Unlike finite differences that approximate derivatives, FVM starts from the integral form of conservation laws. This ensures that mass, momentum, and energy are conserved exactly at the discrete level — a property that makes FVM robust for flows with discontinuities, shocks, and complex geometries.

Why Not Just Use Finite Differences?

Finite differences have limitations for CFD:

IssueFD ProblemFVM Solution
ConservationNot inherently conservativeConservative by construction
Complex geometryRequires structured gridsWorks with unstructured meshes
DiscontinuitiesMay produce oscillationsNatural handling with upwinding
Physical interpretationMathematical abstractionClear flux balance

FVM's philosophy: What goes in must come out — enforced cell by cell.

The Control Volume Concept

Integral Conservation Laws

Consider the general conservation equation:

$$\frac{\partial \phi}{\partial t} + \nabla \cdot (\mathbf{u}\phi) = \nabla \cdot (\Gamma \nabla \phi) + S$$

Where:

  • $\phi$ = conserved quantity (mass, momentum, energy)
  • $\mathbf{u}$ = velocity field
  • $\Gamma$ = diffusion coefficient
  • $S$ = source term

Integrating over a control volume $V$:

$$\int_V \frac{\partial \phi}{\partial t} dV + \int_V \nabla \cdot (\mathbf{u}\phi) dV = \int_V \nabla \cdot (\Gamma \nabla \phi) dV + \int_V S\, dV$$

Applying the divergence theorem:

$$\boxed{\frac{d}{dt}\int_V \phi\, dV + \oint_S (\mathbf{u}\phi) \cdot \mathbf{n}\, dA = \oint_S (\Gamma \nabla \phi) \cdot \mathbf{n}\, dA + \int_V S\, dV}$$

This is the integral conservation law: the rate of change of $\phi$ in the volume equals the net flux through the boundary plus sources.

Visualize a finite volume cell with face fluxes. Click faces to see flux calculations.

Discretizing the Domain

Cell-Centered Storage

In FVM, we divide the domain into non-overlapping control volumes (cells). Values are stored at cell centers:

  • $\phi_P$ = value at cell center P
  • $\phi_E, \phi_W, \phi_N, \phi_S$ = values at neighboring cells

The integral equation becomes:

$$\frac{(\phi V)_P^{n+1} - (\phi V)_P^n}{\Delta t} + \sum_f F_f = \sum_f D_f + S_P V_P$$

Where:

  • $F_f$ = convective flux through face $f$
  • $D_f$ = diffusive flux through face $f$
  • $S_P V_P$ = source term contribution

The Flux Balance

See how fluxes balance across cell faces. Toggle terms to visualize convection, diffusion, and sources.

For a 2D cell with faces e, w, n, s:

$$\frac{d(\phi V)_P}{dt} + (F_e - F_w) + (F_n - F_s) = (D_e - D_w) + (D_n - D_s) + S_P V_P$$

Key insight: The flux leaving one cell automatically enters the neighboring cell. This guarantees conservation at the discrete level.

Computing Face Fluxes

Convective Flux

The convective flux through face $e$ (east face):

$$F_e = (\rho u)_e \phi_e A_e$$

Where:

  • $(\rho u)_e$ = mass flux through face (from velocity field)
  • $\phi_e$ = face value of $\phi$ (needs interpolation!)
  • $A_e$ = face area
The challenge: We store $\phi$ at cell centers, but need it at faces. This leads to convection schemes (upwind, central, QUICK) — covered in Lesson 7.

Diffusive Flux

The diffusive flux through face $e$:

$$D_e = \Gamma_e \frac{\phi_E - \phi_P}{\delta_{PE}} A_e$$

Where:

  • $\Gamma_e$ = diffusivity at face (often averaged)
  • $\delta_{PE}$ = distance between P and E cell centers
  • $A_e$ = face area

Diffusive flux is straightforward because it naturally uses cell-center values.

Mass Flux and Velocity Interpolation

The mass flux $\dot{m}_e = (\rho u)_e A_e$ requires the velocity at the face. Options:

  • Linear interpolation: $u_e = \frac{1}{2}(u_P + u_E)$
  • Rhie-Chow interpolation: Prevents checkerboard oscillations (see Lesson 8)

The FVM Discretized Equation

After computing all face fluxes, we get an algebraic equation for each cell:

$$a_P \phi_P = a_E \phi_E + a_W \phi_W + a_N \phi_N + a_S \phi_S + S_P$$

Or in matrix form:

$$\mathbf{A}\boldsymbol{\phi} = \mathbf{b}$$

Where:

  • $a_P, a_E, a_W, ...$ = coefficients from flux discretization
  • The matrix $\mathbf{A}$ is typically sparse (few non-zeros per row)

Coefficient Interpretation

CoefficientPhysical Meaning
$a_E, a_W, ...$Influence of neighbors through flux
$a_P$Self-contribution (usually $\sum a_{nb}$ + source)
$S_P$Source term contribution

Properties of Good Discretization

  • Conservation: $\sum$ coefficients = 0 (for pure convection)
  • Boundedness: Coefficients same sign → no unphysical oscillations
  • Transportiveness: Upwinding for convection-dominated flows

FVM for Different Cell Types

Structured (Hexahedral) Meshes

Regular grid with clear i,j,k indexing:

  • Easy to implement
  • Efficient memory access
  • Limited to simple geometries

Unstructured (Tetrahedral/Polyhedral) Meshes

General connectivity:

  • Handles complex geometry
  • More complex data structures
  • Same FVM principles apply

Face-Based Data Structures

Modern FVM codes store:

  • List of cells with their volumes
  • List of faces with areas, normals, owner/neighbor cells
  • Face-cell connectivity

This allows the same flux routines for any cell type.

Conservation Guarantee

Why FVM is Conservative

Consider two adjacent cells P and E sharing face $e$:

$$\text{Cell P:} \quad ... + F_{Pe} = ... $$

$$\text{Cell E:} \quad ... - F_{eE} = ... $$

Since $F_{Pe} = F_{eE}$ (same flux, opposite sign), summing all cells:

$$\sum_{\text{cells}} \frac{d(\phi V)}{dt} = \sum_{\text{boundary faces}} F_f$$

Result: Internal fluxes cancel! Only boundary fluxes affect total conservation.

What Can Go Wrong

FVM guarantees discrete conservation, but:

  • Truncation error still exists
  • Non-conservative source terms can break balance
  • Boundary condition implementation matters

Time Discretization

Explicit Methods

$$\phi_P^{n+1} = \phi_P^n + \Delta t \cdot RHS^n$$

  • Simple to implement
  • CFL-limited time step
  • Good for transient accuracy

Implicit Methods

$$\phi_P^{n+1} = \phi_P^n + \Delta t \cdot RHS^{n+1}$$

  • Requires solving linear system
  • Unconditionally stable (for some schemes)
  • Larger time steps possible

Crank-Nicolson (Second-Order)

$$\phi^{n+1} = \phi^n + \frac{\Delta t}{2}(RHS^n + RHS^{n+1})$$

Combines accuracy with stability for smooth transients.

FVM vs FEM vs FDM

AspectFDMFVMFEM
BasisPoint valuesCell averagesBasis functions
ConservationNot inherentBuilt-inWeak form
GeometryStructured onlyAny meshAny mesh
Dominant useResearchCFDStructures
ComplexitySimpleModerateComplex

FVM dominates CFD because:

  • Conservation is paramount for fluids
  • Complex geometries are common in engineering
  • Well-suited for hyperbolic/parabolic PDEs

Practical Implementation

Steps to Discretize a PDE with FVM

  • Write conservation form: $\frac{\partial \phi}{\partial t} + \nabla \cdot \mathbf{F} = S$
  • Integrate over control volume
  • Apply divergence theorem: Surface integrals
  • Approximate face fluxes: Interpolation + gradients
  • Assemble linear system: Sparse matrix + RHS
  • Solve iteratively: Jacobi, Gauss-Seidel, or GMRES

Code Structure (Pseudocode)

FOR each cell P:
    Initialize aP = 0, source = 0

    FOR each face f of cell P:
        Compute face flux F_f
        Add flux contribution to aP and neighbor coefficients
        Add to source if boundary face

    Store equation: aP * phi_P = sum(a_nb * phi_nb) + source

Solve linear system A * phi = b

Key Takeaways

  • FVM integrates conservation laws over control volumes, ensuring discrete conservation
  • Face fluxes are the core concept — what leaves one cell enters the neighbor
  • Cell-centered storage requires interpolation to get face values
  • Convective flux needs schemes (upwind, central) to handle face interpolation
  • Diffusive flux naturally uses cell-center differences
  • Sparse linear system results from discretization, solved iteratively
  • Conservation guarantee: Internal fluxes cancel exactly, only boundaries affect total

What's Next

The mesh is where FVM discretization happens. The next lesson covers Mesh Generation & Quality — how to create good meshes, what quality metrics to check, and why poor meshes lead to poor solutions.