Conservation Laws and the Reynolds Transport Theorem

Conservation Laws and the Reynolds Transport Theorem :: Fundamental Conservation Principles in CFD

This note covers the mathematical foundation of all computational fluid dynamics. These foundations are the conservation laws. Every CFD simulation uses a discrete approximation of these conservation statements. The simulations range from incompressible pipe flow to compressible hypersonic re-entry. The Reynolds Transport Theorem (RTT) bridges the classical continuum mechanics formulation and the control volume formulation. The formulation is system-based in classical mechanics. The formulation is domain-based in CFD code implementation.

**Table of Contents**

• Reynolds Transport Theorem :: The bridge between system and control volume formulations
• Conservation of Mass :: The continuity equation
• Conservation of Momentum :: Newton's second law in continuum form
• Conservation of Energy :: First law of thermodynamics
• Conservation of Species :: Mass fraction transport
• General Conservation Form :: The universal transport equation
• Well-Posedness :: When are the conservation laws complete?
• Global Conservation Error :: What happens when a scheme fails to conserve?
• See Also
• References

Reynolds Transport Theorem :: The Bridge Between System and Control Volume Formulations

The Reynolds Transport Theorem relates the time derivative of an extensive property for a system to the time derivative over a control volume. A system is a closed, material volume. The system moves and deforms with the flow. A control volume is a fixed region in space. Mass flows through a control volume.

Let B be an extensive property. Examples are mass M, momentum Mv, and energy E. Define the intensive property as beta equals partial B partial m. Beta represents property per unit mass. For a system S(t):

$$$$ \frac{dB_{sys}}{dt} = \frac{d}{dt} \int_{V_{sys}(t)} \rho \beta \, dV $$$$

The Reynolds Transport Theorem states:

$$$$ \frac{dB_{sys}}{dt} = \frac{\partial}{\partial t} \int_{CV} \rho \beta \, dV + \int_{CS} \rho \beta (\vec{v} \cdot \hat{n}) \, dA $$$$

CV is the control volume. CS is the control surface. ^n is the outward unit normal. The first term on the right-hand side is the local rate of change within the CV. This is also called the temporal rate of change. The second term is the net convective flux of B across the control surface. A positive value means B leaves the CV.

RTT is the single most important theorem in CFD. Every numerical scheme implements a discrete version of this balance. The schemes include finite volume, finite element, and discontinuous Galerkin. The finite volume method builds explicitly on RTT. Each cell is a control volume. The conservation statement for that cell reads:

rate of accumulation equals net inflow plus source.

Or:

$$$$ \frac{\partial}{\partial t} \int_{V_i} \rho \beta \, dV = - \int_{\partial V_i} \rho \beta (\vec{v} \cdot \hat{n}) \, dA + \int_{V_i} S_\beta \, dV $$$$

The minus sign reflects that positive ^n points outward. A positive flux integral denotes outflow, which is a loss. The source term S_beta represents volumetric production or destruction.

The FVM integrates the differential conservation equation over each cell. The integral form is the conservation statement. FVM is conservative by construction. The flux out of one cell is the flux into its neighbor. Finite difference schemes do not discretize the integral form. Those schemes can fail to conserve. This is one reason FVM dominates CFD.

**Why the distinction matters:** The control volume is fixed in practice. This is the Eulerian description. RTT maps the material derivative to Eulerian partial derivatives exactly. The material derivative follows a fluid particle. Without RTT, you cannot write any conservation equation in a form suitable for a fixed-domain numerical scheme. The material derivative relationship is Db/Dt equals partial b partial t plus v dot nabla b for intensive property b. This is the pointwise form. RTT generalizes this form to the integral statement. The two forms are equivalent under smoothness assumptions. The continuity equation must hold.

**Practical implication for simulation:** You set up a control volume in any CFD code. OpenFOAM, Fluent, and Star-CCM+ all evaluate RTT for each cell. The cell volume integral is the property at the cell centroid times the cell volume. The surface integral is a sum over faces. The sum is rho_f times beta_f times v_f dot n_f times A_f. The accuracy of the conservation balance depends on two factors. The face fluxes must match the physics. The mesh quality matters. Skewness and non-orthogonality of face normals degrade accuracy.

Conservation of Mass :: The Continuity Equation

Mass is conserved in any fluid flow process. Standard Navier-Stokes systems have no mass sources or sinks. Apply RTT with beta equal to one. Mass per unit mass is unity:

$$$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \vec{v}) = 0 $$$$

In open notation, this is:

$$$$ \frac{D\rho}{Dt} + \rho \nabla \cdot \vec{v} = 0 $$$$

The second form separates the material derivative of density from the velocity divergence. For incompressible flow, the material derivative of density is zero. The equation collapses to a significant constraint:

$$$$ \nabla \cdot \vec{v} = 0 $$$$

The divergence of the velocity field must vanish everywhere. This condition is stronger than it seems. The flow is isochoric. Volume elements of a fluid do not expand or contract. A discretization that conserves mass locally prevents spurious pressure oscillations. Local conservation means cell-by-cell conservation. Local conservation ensures that the momentum balance does not generate or destroy mass through boundary inconsistencies.

The continuity equation is first-order in time and space. Numerical treatment is usually trivial. No solver is required for this equation. The role in determining well-posedness of the full system is enormous. Compressible flow couples the density equation to the energy equation via an equation of state. The ideal gas law is one example. The system is not closed until the energy equation is introduced.

**Practical note:** The divergence-free constraint creates the pressure-velocity coupling problem. See Pressure-Velocity Coupling. Pressure does not appear in the continuity equation itself. This constraint causes staggered grids, Rhie-Chow interpolation, and the SIMPLE algorithm. These methods enforce the divergence-free constraint on a numerical level. Pressure is not a thermodynamic variable in the incompressible limit.

Conservation of Momentum :: Newton's Second Law in Continuum Form

Momentum is mass times velocity. The momentum equation states that the rate of change of momentum in a control volume equals the sum of forces acting on that volume:

$$$$ \frac{\partial (\rho \vec{v})}{\partial t} + \nabla \cdot (\rho \vec{v} \otimes \vec{v}) = \nabla \cdot \mathbf{\sigma} + \rho \vec{g} $$$$

Sigma is the stress tensor. For a Newtonian fluid:

$$$$ \mathbf{\sigma} = -p \mathbf{I} + \boldsymbol{\tau} $$$$

The deviatoric stress tau is:

$$$$ \boldsymbol{\tau} = \mu \left[ \nabla \vec{v} + (\nabla \vec{v})^T \right] + \lambda (\nabla \cdot \vec{v}) \mathbf{I} $$$$

For incompressible flow, lambda equals negative two-thirds mu by Stokes hypothesis. The volumetric term vanishes due to incompressibility:

$$$$ \frac{\partial \vec{v}}{\partial t} + \vec{v} \cdot \nabla \vec{v} = -\frac{1}{\rho} \nabla p + \nu \nabla^2 \vec{v} + \vec{g} $$$$

Nu equals mu over rho is the kinematic viscosity. This is the Navier-Stokes equation in its most commonly taught form. The convective term v dot nabla v is non-linear. The non-linearity creates the complexity of turbulent flow. The pressure gradient is a constraint force. The pressure gradient enforces divergence-free velocity. The viscous term provides the diffusion mechanism.

**Conservation implications:** The conservative form uses nabla dot rho v tensor v instead of rho v dot nabla v. The momentum equation ensures that the flux of momentum across an interface is equal and opposite for two neighboring cells. This property is critical in shock-capturing schemes. If the convective flux is not in conservative form, shock speed predictions via Rankine-Hugoniot conditions will be incorrect. The momentum jump across a discontinuity depends on pressure and viscous stress balance. A non-conservative discretisation predicts wrong shock positions.

Compressible flow uses the same conservative structure for the convective term. This structure links the momentum equation to the continuity equation. Density-based solvers solve continuity, momentum, and energy as a coupled system.

**ANSYS comparison:** Fluent supports conservative and non-conservative forms of the convective term. The default is conservative. CFX uses an auto-differencing formulation. CFX preserves conservation by construction. The CFX method uses finite elements but has an equivalent integral form. Star-CCM+ uses a polyhedral finite volume method. Star-CCM+ uses conservative flux reconstruction. OpenFOAM always uses the conservative form.

Conservation of Energy :: First Law of Thermodynamics

The energy equation expresses conservation of total energy. Total energy includes internal energy, kinetic energy, and potential energy:

$$$$ \frac{\partial (\rho E)}{\partial t} + \nabla \cdot (\rho E \vec{v}) = \nabla \cdot (k \nabla T) + \nabla \cdot (\mathbf{\sigma} \cdot \vec{v}) + \rho \vec{g} \cdot \vec{v} + S_E $$$$

E equals e plus one-half v squared plus Phi. The terms represent internal energy plus kinetic energy plus potential energy per unit mass. The terms are:

• Heat conduction: nabla dot (k nabla T). Fourier's law applies. Thermal conductivity is k.
• Work by surface stresses: nabla dot (sigma dot v). This term includes viscous dissipation and pressure work.
• Body force work: rho g dot v.
• Volumetric source: S_E. Chemical reaction heat release and radiation source are examples.

Constant properties and negligible viscous dissipation simplify the energy equation for incompressible flow. The equation becomes a scalar transport equation for temperature:

$$$$ \frac{\partial T}{\partial t} + \vec{v} \cdot \nabla T = \alpha \nabla^2 T + \frac{\dot{q}}{\rho c_p} $$$$

Alpha equals k over (rho c_p) is the thermal diffusivity. Buoyancy matters in the Boussinesq approximation. Temperature couples back to momentum through the buoyancy term. The buoyancy term is rho(T) g, which approximates to rho_0 times [1 minus beta_T times (T minus T_0)] g.

Compressible flow requires the full energy equation. Compressible solvers solve rho E as the primary variable. Temperature and internal energy are reconstructed from E via the equation of state. Density-based solvers use this approach. Examples are Fluent pressure-based solver in compressible mode. OpenFOAM rhoPimpleFoam and sonicFoam are examples too.

A common pitfall in compressible simulations is solving the internal energy equation instead of the total energy equation. Energy is created or destroyed because the kinetic energy component of E is not advected. This error shows as a slow drift in total enthalpy. The drift is especially visible in steady-state runs. Initial transients do not fully decay.

ANSYS Fluent uses total enthalpy as the default energy variable for compressible flows. The variable is h_0 equals h plus one-half v squared. This approach equals solving for total energy. OpenFOAM solves for total energy directly in rhoPimpleFoam.

Conservation of Species :: Mass Fraction Transport

Reactive flows or multi-component mixing transport mass fraction Y_k for species k:

$$$$ \frac{\partial (\rho Y_k)}{\partial t} + \nabla \cdot (\rho \vec{v} Y_k) = - \nabla \cdot \vec{j}_k + \dot{\omega}_k $$$$

j_k is the species diffusion flux. Fick's law gives j_k. The law is j_k equals negative rho D_k nabla Y_k. D_k is the binary diffusion coefficient. omega_k is the net production rate by chemical reactions. The rate is zero for non-reacting flows.

Conservation of total mass requires the sum of all Y_k to equal one. If any species equation is solved, the remaining mass fraction is determined by difference. Reaction source terms ω_k must be stoichiometrically consistent. The net mass production must sum to zero. Multiphase reactive flows can have mass transfer between phases. A non-normalized chemistry mechanism breaks local mass conservation. Spurious mass accumulates. Negative species fractions appear.

**LES-DNS note:** Large-eddy simulation with combustion filters the species transport equation. The filtered equation includes a subgrid-scale flux term. The term is rho-bar times u-tilde_k Y-tilde_k minus rho-bar times u-tilde_k times Y-tilde_k. This term is unclosed. The term is analogous to the subgrid stress tensor. This unclosed term is one of the most challenging problems in turbulent combustion modeling. Flamelet models, transported PDF methods, and linear-Est models are closure approaches.

ANSYS Fluent offers species transport with detailed chemistry. CHEMKIN-format mechanisms are supported. The Eddy-Dissipation Model, Eddy-Dissipation Concept, and Partially-Stirred Reactor (PaSR) models are available. OpenFOAM supporting solvers are reactingFoam, reactingTwoPhaseEulerFoam, and flameSolver. These solvers offer similar capabilities.

General Conservation Form :: The Universal Transport Equation

Any conserved quantity phi can be put in the general form:

$$$$ \frac{\partial (\rho \phi)}{\partial t} + \nabla \cdot (\rho \vec{v} \phi) = \nabla \cdot (\Gamma \nabla \phi) + S_\phi $$$$

Gamma is the diffusion coefficient. S_phi is the source term. Table below shows the mapping for each conservation equation:

All CFD discretisation and all solvers handle this exact equation. Different phi, Gamma, and S_phi values give different physics. The FVM integrates the equation per cell. The convective flux is evaluated at faces. The diffusive flux is evaluated at faces. Only the meaning of phi, Gamma, and S_phi changes.

**Key insight for scheme selection:** The convective term is non-linear when phi depends on a quantity that is itself solved. Momentum depends on velocity and velocity appears in the flux. This non-linearity requires iterative solution methods. Newton-Raphson, SIMPLE, and Picard are examples. Diffusion is linear for constant Gamma. Diffusion becomes non-linear when Gamma depends on phi. Temperature-dependent thermal conductivity is an example.

The general form is what every fvSchemes dictionary in OpenFOAM operates on. The div entry in fvSchemes specifies discretisation for the convective term nabla dot rho v phi. The laplacian entry specifies discretisation for nabla dot (Gamma nabla phi). The same discretisation scheme applies to any phi. The same linearUpwind or TVD scheme works for momentum, temperature, turbulent kinetic energy, and species fractions without modification.

Well-Posedness and Initial-Boundary Conditions

A system of PDEs is well-posed if three conditions hold (Hadamard, 1923):

1. A solution exists 2. The solution is unique 3. The solution depends continuously on the data. Small changes in input produce small changes in output.

The Navier-Stokes equations require these conditions for well-posedness:

1. One initial condition (IC) for each time derivative in the system. The NS equations are first-order in time. You need v(x, 0), p(x, 0), rho(x, 0) for compressible flow, and all passive scalars at t = 0. 2. The correct number of boundary conditions (BCs) at each boundary: 3. Elliptic equations (steady-state or each implicit time step of a transient) require two BCs per boundary face. This equals the order of the PDE. Each second-order equation needs one BC at each boundary face. 4. Hyperbolic equations (inviscid compressible flow at each time step) require a number of BCs equal to the number of characteristic incoming waves at each boundary. 5. The type of BC matters. Dirichlet specifies the value. Neumann specifies the derivative. Robin mixes both. Dirichlet BCs are more common at inlets for velocity, pressure, and temperature. Neumann BCs appear at outlets. zeroGradient handles outflow.

Well-posedness determines solver stability and physical correctness. An ill-posed problem produces divergent, oscillatory, or physically meaningless results. The BC choice at an outlet strongly determines well-posedness. Specifying pressure at the outlet for incompressible flow (Dirichlet on pressure) is almost always correct. Specifying zero-gradient for pressure together with zero-gradient for velocity causes pressure drift and non-unique solutions.

Incompressible flow requires velocity (Dirichlet) or traction (Neumann) at all boundaries. Pressure is defined up to an additive constant. You need at least one pressure Dirichlet BC to fix the absolute pressure level. Without it, the pressure Poisson system is singular. This causes linear solver failures in OpenFOAM. If no pressure reference is set, the pressure matrix has a null space. A fixedValue at one patch and zeroGradient at all others causes this failure.

No numerical scheme can save an ill-posed simulation. Setting BCs correctly is the most important part of any CFD setup. No amount of mesh refinement or scheme order improvement fixes a poorly posed problem.

Global Conservation Error :: What Happens When a Scheme Fails to Conserve

Conservation is essential in CFD, not optional. A non-conservative scheme converges to the wrong solution as the mesh refines. The global conservation error for any conserved quantity is:

$$$$ E_{cons} = \left| \sum_{CV} \left[ \frac{\partial}{\partial t} \int_{V_i} \rho \beta \, dV + \int_{\partial V_i} \rho \beta (\vec{v} \cdot \hat{n}) \, dA - \int_{V_i} S_\beta \, dV \right] \right| $$$$

A globally conservative scheme has E_cons equal zero. Roundoff error is acceptable. An approximately conservative scheme has E_cons decrease at O(Delta x^n) as Delta x goes to zero. The convergence rate is degraded. The asymptotic solution can be wrong.

Incorrect shock speed in compressible flow is the most visible manifestation of non-conservation. A shock moves with speed sigma. The Rankine-Hugoniot jump conditions relate the states upstream and downstream:

$$$$ [\rho(v - \sigma) \cdot \hat{n}] = 0 $$$$ [\rho(v - \sigma)(v \cdot \hat{n}) + p \hat{n}] = 0 $$$$ [(\rho E + p)(v - \sigma) \cdot \hat{n}] = 0 $$$$

These jump conditions require that the numerical scheme is in conservative form. A non-conservative scheme predicts the wrong sigma. The shock appears in the wrong place. The post-shock state is wrong.

Incompressible flow shows local conservation error as mass imbalance at outflow. Fluid is generated or lost. The imbalance creates unphysical pressure gradients and drift. The effect shows as a slow increase or decrease in total mass within the domain. The domain BCs specify zero net mass flow. The error accumulates because each continuity balance is off by a small amount. The error propagates through the pressure-velocity coupling.

**Practical detection:** Use the postProcess function object in OpenFOAM. The command is postProcess -func "fieldValues (massFlow (inlet outlet))". This command computes net mass flow at inlet and outlet patches. A net flux imbalance indicates a conservation error. The most common causes are:

1. Poor mesh quality. Skewness and non-orthogonality degrade flux reconstruction. 2. Inappropriate convection scheme. Linear schemes are not bounded. Oscillations break conservation at large amplitudes. 3. Inconsistent time stepping. Large delta t causes CFL violations and flux overshoots. 4. Inconsistent BCs. Inlet flow rate does not match outlet treatment.

**The practical takeaway:** Choose a conservative discretisation. Finite volume is conservative by construction. Use it. Finite element can be conservative. Use mixed methods or mass-lumped discretisation. Spectral element is conservative if the flux uses a Riemann solver at faces. Do not use a first-order upwind scheme and expect total energy conservation. The effective diffusivity of FOU is one-half u delta x per cell. This value acts as artificial mass, momentum, and energy generation or destruction.

See Also :: Related Notes

• Governing Equations -- derives the conservation system and closure relations
• Incompressible Flow -- the divergence-free constraint and its implications
• Compressible Flow -- hyperbolic structure, Rankine-Hugoniot conditions, characteristic analysis
• Numerical Schemes -- conservative versus non-conservative discretisation, scheme families
• FVM Discretization -- the integral form and cell-by-cell balance
• Pressure-Velocity Coupling -- enforcing incompressibility
• Turbulence Models -- conservation equations for k and omega
• OpenFOAM Case Setup -- fvSchemes configuration for conservation
• OpenFOAM Numerical Schemes -- div and laplacian entries in fvSchemes

References

1. Pope, S.B. (2000). Turbulent Flows. Cambridge University Press. Chapter 3 covers conservation equations. Chapter 5 covers RANS equations. 2. Versteeg, H.K. and Malalasekera, W. (2007). An Introduction to Computational Fluid Dynamics: The Finite Volume Method (2nd ed.). Pearson. Chapter 3 covers the conservation form of governing equations. 3. Hirsch, C. (1989). Numerical Computation of Internal and External Flows (Vol. 1-2). Wiley. Vol. 1, Chapter 7 covers conservation laws. 4. Anderson, J.D. (1995). Computational Fluid Dynamics (McGraw-Hill). Chapter 3 covers Euler equations. Chapter 8 covers conservation forms. 5. Barth, T.J. and Frederickson, P. (1990). "High-order entropy-stable finite-difference schemes." Journal of Computational Physics, 89(1), 1-34. This paper covers conservation and entropy conditions. 6. Toro, E.F. (2009). Riemann Solvers and Numerical Methods for Fluid Dynamics (3rd ed.). Springer. Part I covers basic theory. The chapter covers Rankine-Hugoniot conditions. 7. Hadamard, J. (1923). Über die Grenzwerte von Integralien, welche von mehreren veränderlichen abhängen. Acta Mathematica, 47, 51-78. This paper gives the original well-posedness definition.