1 Introduction

This article provides information on the equation describing conservation of energy relevant to fluid dynamics and computational fluid dynamics (CFD). It first assembles an equation for combined mechanical and thermal energy, i.e. total energy, in terms of material derivatives. It then presents an equation for thermal, or internal, energy. The total energy equation is then provided in terms of local (partial) derivatives, both in terms of internal energy and enthalpy. The implementation of the energy equation in solvers in OpenFOAM is then described.

Some of the information in this article is also presented in the book Notes on CFD: General Principles.

2 Total Energy

The law of conservation of energy states that the total energy of an isolated system remains constant, i.e. it is conserved over time and energy is not created or destroyed but is transformed from one form to another. Here we consider only mechanical and thermodynamic energy, the contributions of which are described in the following sections, using usual notation of tensor algebra and calculus, including D∕Dt  \relax \special {t4ht= representing the material derivative.

2.1 Mechanical Power

The rate of change of mechanical, or kinetic, energy is:

 DK--    D-(|U-|2∕2)    DU--
ρ Dt  ≡ ρ   Dt     ≡ ρ Dt  ∙U,
\relax \special {t4ht=
where U  \relax \special {t4ht= is velocity, specific kinetic energy (kinetic energy per unit mass)         2
K  ≡ |U| ∕2  \relax \special {t4ht= and ρ  \relax \special {t4ht= is mass density. The power flux, or rate of change of strain energy, is
∇ ∙(σ∙U ) ≡ (∇ ∙σ)∙U  + σ∙∇U,
\relax \special {t4ht=
where σ   \relax \special {t4ht= is the mechanical stress tensor. The stress tensor may be decomposed σ \relax \special {t4ht= into a scalar thermodynamic pressure p  \relax \special {t4ht= and viscous stress tensor τ   \relax \special {t4ht= by σ = τ − pI  \relax \special {t4ht=, where I  \relax \special {t4ht= is the identity tensor. The power source, or rate of change of potential energy, is
ρg ∙U,
\relax \special {t4ht=
where g  \relax \special {t4ht= is a body acceleration, e.g. gravity.

2.2 Thermodynamic Power

The rate of change of thermal, or internal, energy is

ρDt ,
\relax \special {t4ht=
where e  \relax \special {t4ht= is specific internal energy (internal energy per unit mass). The heat flux is
− ∇ ∙q,
\relax \special {t4ht=
where q  \relax \special {t4ht= is the heat flux vector, defined positive inwards. The heat source is
\relax \special {t4ht=
where r  \relax \special {t4ht= is any specific heat source.

2.3 Conservation of Energy

The rate of change of total energy for a particle of material must equal the input of mechanical and thermodynamic power from fluxes and sources acting on the particle. In the limit where particle size is infinitesimally small

 De     DK
ρ--- + ρ---- =  − ∇ ∙q + ρr+ ∇ ∙(σ ∙U )+ ρg ∙U
◟◝D◜t◞   ◟-D◝t◜ ◞  ◟----◝◜----◞  ◟-------◝◜-------◞
thermo   mech      thermo             mech  \relax \special {t4ht=

3 Internal Energy

An equation for internal energy is produced by simplifying the mechanical contributions which, expressed as

  DU--∙       ∙   ∙      ∙         ∙
ρ Dt   U − (∇  σ)  U − σ ∙∇U  − ρg  U,
\relax \special {t4ht=
combines with the momentum equation
ρ DU--= ∇ ∙σ + ρg,
\relax \special {t4ht=
to reduce the mechanical contributions to   ∙
σ ∙∇U  \relax \special {t4ht=. Equation 7 can then be expressed as
ρ--- = − ∇ ∙q+  ρr + σ ∙∙∇U
\relax \special {t4ht=
The σ∙∙∇U  \relax \special {t4ht= term represents the contribution of mechanical power to internal energy, and thus, random motion of particles. The expression (∇ ∙σ)∙U  \relax \special {t4ht= must then represent a power due bulk motion of particles.

4 Total Energy/Enthalpy, local derivatives

We can express our equations in terms of the local derivative (or partial derivative, spatial derivative, …) ∂∕∂t  \relax \special {t4ht=, where D∕Dt  ≡ ∂∕∂t+  U ∙∇ \relax \special {t4ht=. Applying conservation of mass, the following relationship holds for any tensor Q  \relax \special {t4ht=:

 DQ     ∂ρQ
ρ-Dt- ≡ -∂t--+ ∇ ∙(ρUQ  )  \relax \special {t4ht=

Combining equations 7 and 11, and decomposing the stress tensor σ = τ − pI  \relax \special {t4ht=, gives:

∂ρe             ∂ρK
----+∇ ∙(ρUe )+ ----+ ∇ ∙(ρUK  )+ ∇ ∙(Up ) = − ∇ ∙q+ ∇ ∙(τ∙U )+ ρr+ρg ∙U
 ∂t              ∂t\relax \special {t4ht=

Enthalpy is the sum of internal energy and kinematic pressure, i.e. h ≡ e + p∕ρ  \relax \special {t4ht=. Combining this with equation 12 gives:

∂ρh-+ ∇ ∙(ρUh )+ ∂ρK-+ ∇ ∙(ρUK  )− ∂p-= − ∇ ∙q+ ∇ ∙(τ∙U )+ ρr+ ρg∙ U
 ∂t               ∂t               ∂t \relax \special {t4ht=

Total energy can be defined as E ≡ e + K  \relax \special {t4ht=. Combining this with equation 12 gives:

∂t  + ∇ ∙(ρUE  )+ ∇ ∙(Up ) = − ∇ ∙q + ∇ ∙(τ∙ U) + ρr + ρg ∙U  \relax \special {t4ht=

5 Energy Equation in OpenFOAM Solvers

The solution of the energy equation is included in several solvers in OpenFOAM for compressible flow, combustion, heat transfer, multiphase flow and particle tracking. The source code can be found for these solvers within files in sub-directories of the $FOAM_SOLVERS directory of OpenFOAM (including the compressible, combustion, heatTransfer, multiphase and lagrangian sub-directories).

The energy equation is generally implemented in the form of total energy expressed in equations 12 and 13, without the mechanical sources ∇ ∙(τ∙U )  \relax \special {t4ht= and ρg ∙U  \relax \special {t4ht=. A heat flux q = − αeff∇e  \relax \special {t4ht= is assumed, where the effective thermal diffusivity αeff   \relax \special {t4ht= is the sum of laminar and turbulent thermal diffusivities. The implementation of each energy equation contains thermal source terms ρr  \relax \special {t4ht= relevant to the particular solver.

For example, the sonicFoam solver contains the following implementation of the energy equation from equation 12.

fvScalarMatrix EEqn
    fvm::ddt(rho, e) + fvm::div(phi, e)
  + fvc::ddt(rho, K) + fvc::div(phi, K)
  + fvc::div(fvc::absolute(phi/fvc::interpolate(rho), U), p, "div(phiv,p)")
  - fvm::laplacian(turbulence->alphaEff(), e)
    fvOptions(rho, e)

sonicFoam solves equations sequentially, so solves the momentum equation for U  \relax \special {t4ht= before updating the specific kinetic energy field K = |U |2∕2  \relax \special {t4ht= for the energy equation above. More commonly, the energy equation is implemented in terms of both internal energy e  \relax \special {t4ht= and enthalpy h  \relax \special {t4ht=, as both equations 12 and 13, allowing the user to choose the solution variable, e  \relax \special {t4ht= or h  \relax \special {t4ht=, at run time. For example, the rhoPimpleFoam solver has the following implementation:

volScalarField& he = thermo.he();

fvScalarMatrix EEqn
    fvm::ddt(rho, he) + fvm::div(phi, he)
  + fvc::ddt(rho, K) + fvc::div(phi, K)
  + ( == "e"
      ? fvc::div
            fvc::absolute(phi/fvc::interpolate(rho), U),
      : -dpdt
  - fvm::laplacian(turbulence->alphaEff(), he)
    fvOptions(rho, he)

Here, “he” represents either h  \relax \special {t4ht= or e  \relax \special {t4ht=. The 5th term switches between ∇ ∙(Up )  \relax \special {t4ht= and − ∂p∕∂t  \relax \special {t4ht= depending on the solution variable chosen by the user.

The rhoCentralFoam solver includes an implementation of an energy equation best represented by equation 14 that includes the mechanical source ∇ ∙(τ∙ U)  \relax \special {t4ht=.

6 Total vs Internal Energy

The choice of energy equation has a significant on some solutions particularly across shocks. In the well known 1D shockTube tutorial example (Sod’s problem), the initial discontinuity causes a shock to propagate into the low pressure region and an expansion wave to propagate upstream. The figure below shows the temperature after 0.007 s, with simulation results compared with the analytical solution. Using the version of sonicFoam prior to OpenFOAM v2.2.0 that solves a thermal energy equation, the temperature difference across the shock is badly predicted. Using sonicFoam from v2.2.0 onwards that solves a total energy equation, conservation of total energy ensures the temperature difference is predicted accurately.

epsfig\relax \special {t4ht=

Figure 1: Shock tube problem, solution at 0.007s

Energy Equation in OpenFOAM