Temperature calculation from energy variables in OpenFOAM

The energy conservation equation is expressed in terms of internal energy \(e\) or enthalpy \(h\) in OpenFOAM and the temperature field is calculated from these solution variables. More precisely, we can specify the energy variable in the energy entry in the thermophysicalProperties file and the available options are the followings:

  • specific sensible enthalpy \(h {\rm [J/kg]}\) (sensibleEnthalpy)
  • specific sensible internal energy \(e {\rm [J/kg]}\) (sensibleInternalEnergy)

If we select the hePsiThermo thermophysical model, the temperature field \(T\) is calculated from the solved energy variable in the following function, where the compressibility \(\psi\), dynamic viscosity \(\mu\) and thermal diffusivity \(\alpha\) are also calculated.

Source Code

The calculation procedure of the temperature field depends on the selected energy variable, so the following THE function accordingly switches the called method.

If we choose the sensible enthalpy as the energy variable, the following THs function is called to calculate the temperature from the sensible enthalpy.

The calculation of the temperature is done iteratively using the Newton-Raphson method.

If the specific heat capacity at constant pressure \(c_p\) is expressed in the form of temperature polynomial function (hPolynomial)
c_p(T) = \sum_{i=0}^7 c_i T^i, \tag{1}
the temperature in the j-th cell \(T_j\) is calculated from the following equation
\displaystyle \int_{T_{std}}^{T_j} \left( \sum_{i=0}^7 c_i T^i \right) dT = h_j, \tag{2}
where \(h_j\) is the sensible enthalpy value in the j-th cell. In general the equation will be nonlinear, the iterative solution technique is implemented.

Compressible Flow Solvers

The above function calculate() is called by rhoPimpleFoam, rhoSimpleFoam and sonicFoam etc. from the line “thermo.correct()” after solving the energy conservation equation EEqn:

After updating the compressibility field \(\psi\), the pressure Poisson equation pEqn is constructed and solved.

The density field \(\rho\) is calculated with updated pressure and compressibility fields.

Energy variable Function used to calculate temperature
sensible enthalpy THs
absolute enthalpy THa
sensible internal energy TEs
absolute internal energy TEa

Study with great curiosity and let’s see the physical phenomena from a HIGHER point of view!

Introduction to laplacianFoam and simple validation calculation

In this blog post, I will try to give a description of the governing equation of the laplacianFoam in OpenFOAM that solves a simple Laplace equation, e.g. for thermal diffusion in a solid.

Governing Equation

The heat conduction equation is given by the following equation:

\frac{\partial T}{\partial t} = \frac{1}{\rho c_p} \nabla \cdot \left(k \nabla T\right) + \frac{q}{\rho c_p}, \tag{1} \label{eq:conductionEqn}

where \(T\;{\rm [ K]}\) is the absolute temperature field, \(\rho\;{\rm [kg/m^3]}\) is the density field, \(q\;{\rm [W/m^3]}\) is the rate of energy generation per unit volume, \(k\;{\rm [W/(m\cdot K)]}\) is the thermal conductivity and \(c_p\;{\rm [J/(kg\cdot K)]}\) is the specific heat at constant pressure.

If the heat capacity \(\rho c_p\) is spatially uniform, the Eq. \eqref{eq:conductionEqn} can be transformed into the following form irrespective of whether the thermal conductivity \(k\) is spatially uniform or not:

\frac{\partial T}{\partial t} = \nabla \cdot \left(\alpha \nabla T\right) + \frac{q}{\rho c_p}, \tag{2} \label{eq:conductionEqn2}

where \(\alpha = k/\rho c_p\;{\rm [m^2/s]}\) is the thermal diffusivity.

The laplacianFoam (in OpenFOAM-4.x and earlier versions) doesn’t consider the heat generation and the implemented equation is

\frac{\partial T}{\partial t} = \nabla \cdot \left(\alpha \nabla T\right), \tag{3} \label{eq:laplacianFoam}

but the solver in the latest development version supports the fvOptions so that we can solve \eqref{eq:conductionEqn2} and specify a volumetric heat source.

We can find the Eq. \eqref{eq:conductionEqn2} solved in laplacianFoam.C.

The variable DT represents the thermal diffusivity \(\alpha\) and it is specified in the constant/transportProperties file.

Comparison with Analytical Solution

We consider the steady state problem of source-free heat conduction in a concentric cylinder whose inner and outer walls are maintained at constant temperature of 400 K and 300 K respectively. This problem can be analytically solved and the following relation holds for the temperature distribution in the radial direction \(T(r)\):

\frac{T_1 – T(r)}{T_1 – T_2} = \frac{\ln{(r/r_1)}}{\ln{(r_2/r_1)}}, \tag{4} \label{eq:cylinderT}

where the subscripts 1 and 2 indicate the values at the inner and outer walls respectively shown in Figure 1.

Fig. 1 Problem setting

The temperature distribution obtained using laplacianFoam is shown in Figure 2 and the comparison between the numerical and analytical solutions is shown in Figure 3. The agreement with the analytical solution is good.

Fig. 2
Fig. 2 Temperature distribution
Fig. 3 Comparison of numerical and analytical solutions
Introduction of Source Term

In the next post, I’ll deal with a simple example case to introduce how to specify a heat generation source term in laplacianFoam using the fvOptions functionality (scalarSemiImplicitSource).


buoyantPimpleFoam and buoyantSimpleFoam in OpenFOAM

In this blog post, I will try to give a description of the governing equations of the following solvers in OpenFOAM for buoyant, turbulent flow of compressible fluids:

  • buoyantPimpleFoam (Transient solver)
  • buoyantSimpleFoam (Steady-state solver)

The solvers in OpenFOAM are named properly and comprehensibly so that the users can guess the meanings from their names 🙂

Mass Conservation

The mass conservation (continuity) equation of buoyantPimpleFoam is given by the following equation:

\frac{\partial \rho}{\partial t} + \nabla \cdot \left( \rho \boldsymbol{u} \right) = 0, \tag{1} \label{eq:mass}
where \(\boldsymbol{u}\) is the velocity field and \(\rho\) is the density field. For the steady-state solver, the time derivative term is omitted.

Momentum Conservation

The momentum conservation equation of buoyantPimpleFoam is given by the following equation:

\frac{\partial \left(\rho \boldsymbol{u} \right)}{\partial t} + \nabla \cdot \left(\rho\boldsymbol{u}\boldsymbol{u}\right) = &-\nabla p + \rho \boldsymbol{g} + \nabla \cdot \left(2 \mu_{eff} D \left( \boldsymbol{u} \right) \right) \\
&-\nabla \left( \frac{2}{3}\mu_{eff}\left(\nabla \cdot \boldsymbol{u} \right) \right), \tag{2} \label{eq:momentum}

where \(p\) is the static pressure field and \(\boldsymbol{g}\) is the gravitational acceleration. The effective viscosity \(\mu_{eff}\) is the sum of the molecular and turbulent viscosity and the rate of strain (deformation) tensor \(D\left(\boldsymbol{u}\right)\) is defined as \(D\left(\boldsymbol{u}\right) = \frac{1}{2}\left( \nabla \boldsymbol{u} + \left(\nabla \boldsymbol{u}\right)^{T} \right)\). As is the case in the continuity equation, the time derivative term is omitted in buoyantSimpleFoam.

In terms of the implementation in OpenFOAM, the pressure gradient and gravity force terms are rearranged in the following form:

-\nabla p + \rho \boldsymbol{g}
&= -\nabla \left( p_{rgh} + \rho \boldsymbol{g} \cdot \boldsymbol{r} \right) + \rho \boldsymbol{g} \tag{3a}\\
&= -\nabla p_{rgh} – \left( \boldsymbol{g} \cdot \boldsymbol{r} \right) \nabla \rho \:- \rho \boldsymbol{g} + \rho \boldsymbol{g} \tag{3b} \\
&= -\nabla p_{rgh} – \left( \boldsymbol{g} \cdot \boldsymbol{r} \right) \nabla \rho, \tag{3c} \label{eq:prgh}

where \(p_{rgh} = p \;- \rho \boldsymbol{g} \cdot \boldsymbol{r}\) and \(\boldsymbol{r}\) is the position vector.

In the above code, MRF stands for Multiple Reference Frame, which is one method for solving the problems including the rotating parts with the static mesh. This option is activated if constant/MRFProperties file exists. The line 11(fvOptions(rho, U)) is to add the user-specific source terms using fvOptions functionality.

Energy Conservation

We can choose either internal energy \(e\) or enthalpy \(h\) as the energy solution variable [3]. This selection is made according to energy keyword in thermophysicalProperties file.

For each variable, the energy conservation equations of buoyantPimpleFoam are given by the following equations:

  • Enthalpy (sensibleEnthalpy)

&\frac{\partial \left(\rho h \right)}{\partial t} + \nabla \cdot \left( \rho \boldsymbol{u} h \right) + \frac{\partial \left(\rho K \right)}{\partial t} + \nabla \cdot \left( \rho \boldsymbol{u} K \right) -\frac{\partial p}{\partial t} \\
&= \nabla \cdot \left( \alpha_{eff} \nabla h \right) + \rho \boldsymbol{u} \cdot \boldsymbol{g} \tag{4} \label{eq:energyH}

  • Internal energy (sensibleInternalEnergy)

&\frac{\partial \left(\rho e \right)}{\partial t} + \nabla \cdot \left( \rho \boldsymbol{u} e \right) + \frac{\partial \left(\rho K \right)}{\partial t} + \nabla \cdot \left( \rho \boldsymbol{u} K \right) + \nabla \cdot \left( p \boldsymbol{u} \right) \\
&= \nabla \cdot \left( \alpha_{eff} \nabla e \right) + \rho \boldsymbol{u} \cdot \boldsymbol{g} \tag{5} \label{eq:energyE}

where \(K \equiv |\boldsymbol{u}|^2/2\) is kinetic energy per unit mass and the enthalpy per unit mass \(h\) is the sum of the internal energy per unit mass \(e\) and the kinematic pressure \(h \equiv e + p/\rho\). We can get the following relations from this definition:
\frac{\partial \left(\rho e \right)}{\partial t} = \frac{\partial \left(\rho h \right)}{\partial t} \;-\; \frac{\partial p}{\partial t}, \tag{6a} \\
\nabla \cdot \left( \rho \boldsymbol{u} e \right) = \nabla \cdot \left( \rho \boldsymbol{u} h \right) \;-\; \nabla \cdot \left( p \boldsymbol{u} \right). \tag{6b}

The effective thermal diffusivity \(\alpha_{eff}/\rho\) is the sum of laminar and turbulent thermal diffusivities:

\alpha_{eff} = \frac{\rho \nu_t}{Pr_t} + \frac{\mu}{Pr} = \frac{\rho \nu_t}{Pr_t} + \frac{k}{c_p}, \tag{7} \label{eq:alphaEff}

where \(k\) is the thermal conductivity, \(c_p\) is the specific heat at constant pressure, \(\mu\) is the dynamic viscosity, \(\nu_t\) is the turbulent (kinematic) viscosity, \(Pr\) is the Prandtl number and \(Pr_t\) is the turbulent Prandtl number. The expression of the laminar thermal diffusivity changes depending on the selected thermodynamics package. As is the case in the momentum equation, the time derivative terms are omitted in buoyantSimpleFoam.

The line 21 represents the contributions from the radiative heat transfer. The effective thermal diffusivity \(\alpha_{eff}\) \eqref{eq:alphaEff} is calculated in heThermo.C:

When solving for enthalpy \(h\), the pressure-work term \(dp/dt\) can be excluded by setting the dpdt option to no in thermophysicalProperties file [3]. The volScalarField \(dp/dt\) is defined and uniformly initialized to zero in createFields.H and updated in pEqn.H if this option is activated.

When the mesh is static, the fvc::absolute function does nothing as shown below:


[1] Energy Equation in OpenFOAM | CFD Direct
[2] OpenFOAM User Guide: 7.1 Thermophysical models | CFD Direct
[3] OpenFOAM 2.2.0: Thermophysical Modelling

Conduction, convection and radiation in OpenFOAM (Under construction)

This page is under construction. I will update by adding a description of each solver and model.



Conduction + Convection (Conjugate Heat Transfer)
  • chtMultiRegionSimpleFoam (Steady)
  • chtMultiRegionFoam (Transient)
+ Radiation

All the above solvers but laplacianFoam are able to deal with the radiative heat transfer. There are the following five (virtually four) models available in OpenFOAM. Their source code is located in src/thermophysicalModels/radiation/radiationModels and we can see the brief descriptions in the header file of each radiation class.

  • P1

Keywords: optical thickness [2]

  • fvDOM (Finite Volume Discrete Ordinates Method)

  • viewFactor

  • opaqueSolid

  • none

The settings of the radiation models are described in constant/radiationProperties file.


Solvers for heat transfer problems in OpenFOAM – buoyantBoussinesqPimpleFoam

There are several solvers for solving the heat transfer problems in OpenFOAM. In this post, I will try to give a description of how to derive the governing equations of the following two solvers:

  • buoyantBoussinesqSimpleFoam
  • buoyantBoussinesqPimpleFoam

In the presence of gravity body force, the mass and momentum conservation equations are:

\frac{\partial \rho}{\partial t} + \nabla \cdot \left(\rho \boldsymbol{u}\right) = 0, \tag{1} \label{eq:massEqn}

\frac{\partial \left(\rho \boldsymbol{u} \right)}{\partial t} + \nabla \cdot \left(\rho\boldsymbol{u}\boldsymbol{u}\right) = &-\nabla p + \rho \boldsymbol{g} + \nabla \cdot \left(2 \mu_{eff} D \left( \boldsymbol{u} \right) \right) \\
&-\nabla \left( \frac{2}{3}\mu_{eff}\left(\nabla \cdot \boldsymbol{u} \right) \right), \tag{2} \label{eq:momEqn}

where \(\boldsymbol{u}\) is the velocity field, \(p\) is the pressure field, \(\rho\) is the density field, and \(\boldsymbol{g}\) is the gravitational acceleration. The effective viscosity \(\mu_{eff}\) is the sum of the molecular and turbulent viscosity and the rate of strain tensor \(D\left(\boldsymbol{u}\right)\) is defined as \(D\left(\boldsymbol{u}\right) = \frac{1}{2}\left( \nabla \boldsymbol{u} + \left(\nabla \boldsymbol{u}\right)^{T} \right)\).

If both density and gravitational acceleration are constant, the gravitational force can be expressed in terms of gradient:

\rho \boldsymbol{g} = \nabla \left( \rho \boldsymbol{g} \cdot \boldsymbol{r} \right), \tag{3} \label{eq:Fgravity}

where \(\boldsymbol{r}\) is the position vector, so that the pressure gradient and gravity force can be lumped together as shown in the following equation:

\nabla p \;-\; \rho \boldsymbol{g} = \nabla \left( p \;-\; \rho \boldsymbol{g} \cdot \boldsymbol{r} \right). \tag{4} \label{eq:conversion}

Boussinesq approximation

Now we consider the case when the density is not constant. The Boussinesq approximation is valid when the variation of the density induced by the temperature change is small, as Ferziger and Peric (2001, pp.14-15) states that:

In flows accompanied by heat transfer, the fluid properties are normally functions of temperature. The variations may be small and yet be the cause of the fluid motion. If the density variation is not large, one may treat the density as constant in the unsteady and convection terms, and treat it as variable only in the gravitational term. This is called the Boussinesq approximation.

Hereafter, we denote the reference density by \(\rho_0\) at the reference temperature \(T_0\). If we replace \(\rho\) by \(\rho_0\) in the Eq.\eqref{eq:massEqn}\eqref{eq:momEqn} except for the gravitational term, then we get:

\nabla \cdot \boldsymbol{u} = 0, \tag{5} \label{eq:massEqn2}

\frac{\partial \left(\rho_0 \boldsymbol{u} \right)}{\partial t} + \nabla \cdot \left(\rho_0\boldsymbol{u}\boldsymbol{u}\right) = -\nabla p + \rho \boldsymbol{g} + \nabla \cdot \left(2 \mu_{eff} D \left( \boldsymbol{u} \right) \right). \tag{6} \label{eq:momEqn2}

Then divide both sides of Eq.\eqref{eq:momEqn2} by \(\rho_0\) and we get:

\frac{\partial \boldsymbol{u}}{\partial t} + \nabla \cdot \left(\boldsymbol{u}\boldsymbol{u}\right) = -\frac{1}{\rho_0} \left( \nabla p \;-\; \rho\boldsymbol{g} \right) + \nabla \cdot \left(2 \nu_{eff} D \left( \boldsymbol{u} \right) \right). \tag{7} \label{eq:momEqn3}

Here, the density \(\rho\) in the gravitational term is expressed by the linear function of the temperature \(T\):

\rho \approx \rho_0\left[ 1 \;-\; \beta\left( T \;-\; T_0 \right) \right]. \tag{8} \label{eq:rhoEqn}

This relation is derived from the definition of the volumetric thermal expansion coefficient \(\beta\):

\beta \equiv -\frac{1}{\rho}\frac{\partial \rho}{\partial T}
\approx -\frac{1}{\rho_0}\frac{\rho \;-\; \rho_0}{T \;-\; T_0}. \tag{9} \label{eq:beta}

We have to know when these approximations are valid. Ferziger and Peric (2001, p.15) states that:

This approximation introduces errors of the order of 1% if the temperature differences are below e.g. 2℃ for water and 15℃ for air. The error may be more substantial when temperature differences are larger; the solution may even be qualitatively wrong.

If temperature differences are larger, we should use buoyantSimpleFoam or buoyantPimpleFoam without Boussinesq approximation.

In terms of the implementation, the pressure gradient and gravity force terms are rearranged in the following form:

&-\nabla \left(\frac{p}{\rho_0}\right) + \left(\frac{\rho}{\rho_0}\right)\boldsymbol{g} \\
=&-\nabla \left(\frac{p \;- \rho \boldsymbol{g} \cdot \boldsymbol{r}}{\rho_0} + \frac{\rho \boldsymbol{g} \cdot \boldsymbol{r}}{\rho_0} \right) + \left(\frac{\rho}{\rho_0}\right)\boldsymbol{g} \tag{10a}\\
=&-\nabla p_{rgh} – \left( \boldsymbol{g} \cdot \boldsymbol{r} \right) \nabla \left(\frac{\rho}{\rho_0}\right), \tag{10b} \label{eq:final}

where we define a new variable \(p_{rgh} = \left(p-\rho\boldsymbol{g}\cdot\boldsymbol{r}\right)/\rho_0\) and finally we get:

\frac{\partial \boldsymbol{u}}{\partial t} + \nabla \cdot \left(\boldsymbol{u}\boldsymbol{u}\right) =\; &- \nabla p_{rgh} – \left( \boldsymbol{g} \cdot \boldsymbol{r} \right) \nabla \left(\frac{\rho}{\rho_0}\right) \\ &+ \nabla \cdot \left(2 \nu_{eff} D \left( \boldsymbol{u} \right) \right). \tag{11} \label{eq:momEqn4}

This momentum equation \eqref{eq:momEqn4} is implemented in UEqn.H as shown in the following code:

The Eq. \eqref{eq:final} is calculated using the fvc::reconstruct function from the face flux fields and the variable rhok is defined in createFields.H as \(\rho/\rho_0\):


[1] Joel H. Ferziger and Milovan Peric (2001): Computational Methods for Fluid Dynamics

Please wait till the next update!