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)
\begin{equation}
c_p(T) = \sum_{i=0}^7 c_i T^i, \tag{1}
\end{equation}
the temperature in the j-th cell \(T_j\) is calculated from the following equation
\begin{equation}
\displaystyle \int_{T_{std}}^{T_j} \left( \sum_{i=0}^7 c_i T^i \right) dT = h_j, \tag{2}
\end{equation}
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.

Summary
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!

Tensor Operations in OpenFOAM


Tensor operations are frequently used in turbulence modeling because the strain rate and vorticity tensors can be expressed in terms of the velocity gradient tensor.

In this blog post, I will pick out some typical operations and give brief explanations of them with some usage examples in OpenFOAM.

Keywords
strain rate tensor, vorticity tensor, Q-criterion, Hodge dual

OpenFOAM Version
OpenFOAM-dev

Gradient of a Vector Field | fvc::grad(u)

The gradient of a velocity vector \(\boldsymbol{u}\) returns a velocity gradient tensor (second rank tensor).
\begin{eqnarray}
\nabla \boldsymbol{u} &\equiv& \partial_{i} u_j \\
&=& \left(
\begin{matrix}
\partial u_1/\partial x_1 & \partial u_2/\partial x_1 & \partial u_3/\partial x_1 \\
\partial u_1/\partial x_2 & \partial u_2/\partial x_2 & \partial u_3/\partial x_2 \\
\partial u_1/\partial x_3 & \partial u_2/\partial x_3 & \partial u_3/\partial x_3
\end{matrix} \right)
\end{eqnarray}

Symmetric Part of a Second Rank Tensor | symm(T)

\begin{eqnarray}
{\rm symm}(\boldsymbol{T}) &\equiv& \frac{1}{2} (\boldsymbol{T} + \boldsymbol{T}^T) \\
&=& \frac{1}{2} \left(
\begin{matrix}
2T_{11} & T_{12} + T_{21} & T_{13} + T_{31} \\
T_{21} + T_{12} & 2T_{22} & T_{23} + T_{32} \\
T_{31} + T_{13} & T_{32} + T_{23} & 2T_{33}
\end{matrix} \right)
\end{eqnarray}

Twice the Symmetric Part of a Second Rank Tensor | twoSymm(T)

\begin{eqnarray}
{\rm twoSymm}(\boldsymbol{T}) &\equiv& \boldsymbol{T} + \boldsymbol{T}^T \\
&=& \left(
\begin{matrix}
2T_{11} & T_{12} + T_{21} & T_{13} + T_{31} \\
T_{21} + T_{12} & 2T_{22} & T_{23} + T_{32} \\
T_{31} + T_{13} & T_{32} + T_{23} & 2T_{33}
\end{matrix} \right)
\end{eqnarray}

Skew-symmetric Part of a Second Rank Tensor | skew(T)

\begin{eqnarray}
{\rm skew}(\boldsymbol{T}) &\equiv& \frac{1}{2} (\boldsymbol{T} – \boldsymbol{T}^T) \\
&=& \frac{1}{2} \left(
\begin{matrix}
0 & T_{12} – T_{21} & T_{13} – T_{31} \\
T_{21} – T_{12} & 0 & T_{23} – T_{32} \\
T_{31} – T_{13} & T_{32} – T_{23} & 0
\end{matrix} \right)
\end{eqnarray}

The symm and skew operations of the velocity gradient tensor field frequently appear in the source code. The following is a typical example.

In the above code, the symmetric and antisymmetric parts of the velocity gradient tensor \(\partial u_j/\partial x_i\) are defined as follows:
\begin{eqnarray}
S_{ij} &=& \frac{1}{2} \left( \frac{\partial u_j}{\partial x_i} + \frac{\partial u_i}{\partial x_j} \right), \\
\Omega_{ij} &=& \frac{1}{2} \left( \frac{\partial u_j}{\partial x_i} – \frac{\partial u_i}{\partial x_j} \right),
\end{eqnarray}
where \(S_{ij}\) is the strain rate tensor and \(-\Omega_{ij}\) is the vorticity (spin) tensor.

Hodge Dual | *T

\begin{equation}
*T = \left( T_{23},\;-T_{13},\;T_{12} \right)
\end{equation}

The vorticity vector \(\boldsymbol{\omega}\) is calculated as the Hodge dual of the skew-symmetric part of the velocity gradient tensor.
\begin{eqnarray}
\boldsymbol{\omega} &=& 2 \times \left( * \Omega_{ij} \right) \\
&=& 2 \times \left( \Omega_{23},\;-\Omega_{13},\;\Omega_{12} \right) \\
&=& \left( \frac{\partial u_3}{\partial x_2}-\frac{\partial u_2}{\partial x_3},\;\frac{\partial u_1}{\partial x_3}-\frac{\partial u_3}{\partial x_1},\;\frac{\partial u_2}{\partial x_1}-\frac{\partial u_1}{\partial x_2} \right)
\end{eqnarray}

The operator * in front of the skew is the Hodge dual operator.

Inner Product of Two Second Rank Tensors | T & S

\begin{equation}
P_{ij} = T_{ik} S_{kj}
\end{equation}

Double Inner Product of Two Second Rank Tensors | T && S

\begin{align}
s = T_{ij}S_{ij} &= T_{11}S_{11} + T_{12}S_{12} + T_{13}S_{13} \\
&+ T_{21}S_{21} + T_{22}S_{22} + T_{23}S_{23} \\
&+ T_{31}S_{31} + T_{32}S_{32} + T_{33}S_{33}
\end{align}

Trace of a Second Rank Tensor | tr(T)

\begin{equation}
{\rm tr}(\boldsymbol{T}) \equiv T_{11} + T_{22} + T_{33} = \boldsymbol{T} {\rm \&}{\rm \&} \boldsymbol{I}
\end{equation}

The Q-criterion can be used to identify vortex cores:
\begin{align}
Q &= \frac{1}{2} \left[ (\nabla \cdot \boldsymbol{u})^2 – {\rm tr}(\nabla \boldsymbol{u}^2)\right] \\
&= \frac{1}{2} \left[ (\nabla \cdot \boldsymbol{u})^2 + \Omega_{ij}\Omega_{ij} – S_{ij}S_{ij} \right].
\end{align}
For incompressible flows, it can be simplified as follows:
\begin{equation}
Q = \frac{1}{2} \left[ \Omega_{ij}\Omega_{ij} – S_{ij}S_{ij} \right].
\end{equation}

References

[1] OpenFOAM Programmer’s Guide
[2] CFD Direct | Tensor Mathematics

Recent changes in the basics of the OpenFOAM solvers

constrainHbyA

Comment by Henry G. Weller:

The boundary conditions of HbyA are now constrained by the new “constrainHbyA” function which applies the velocity boundary values for patches for which the velocity cannot be modified by assignment and pressure extrapolation is not specified via the new
“fixedFluxExtrapolatedPressureFvPatchScalarField”.

The member function assignable() is defined in the abstract base class of boundary conditions (fvPatchField):

It is overridden in the basic boundary condition classes (fixedValue, mixed, directionMixed and sliced):

fvc::flux and dotInterpolate

Comment by Henry G. Weller:

surfaceInterpolationScheme: Added dotInterpolate member-function

dotInterpolate interpolates the field and “dots” the resulting face-values with the vector field provided which removes the need to create a temporary field for the interpolate. This reduces the peak storage of OpenFOAM caused by the divergence of the gradient of vector fields, improves memory management and under some conditions decreases run-time.

This development is based on a patch contributed by Paul Edwards, Intel.

constrainPressure

In the case of compressible flow, the following relation \eqref{eq:1} holds for the velocity on the face:
\begin{equation}
\rho \boldsymbol{U}_f \cdot \boldsymbol{S}_f = \left(\frac{\rho \boldsymbol{H}}{A_p}\right)_f \cdot \boldsymbol{S}_f\;-\;\left(\frac{\rho}{A_p}\right)_f (\nabla p)_f \cdot \boldsymbol{S}_f. \tag{1} \label{eq:1}
\end{equation}
Thus
\begin{equation}
(\nabla p)_f \cdot \boldsymbol{n}_f = \frac{1}{|\boldsymbol{S}_f|\left(\frac{\rho}{A_p}\right)_f} \left(\left(\frac{\rho \boldsymbol{H}}{A_p}\right)_f \cdot \boldsymbol{S}_f\;-\;\rho \boldsymbol{U}_f \cdot \boldsymbol{S}_f \right). \tag{2} \label{eq:2}
\end{equation}
The fixedFluxPressure boundary condition calculates the boundary normal gradient of the pressure field using the Eq. \eqref{eq:2} so that the boundary fluxes satisfy the Eq. \eqref{eq:1}.

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:

\begin{align}
\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}
\end{align}

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:

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

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

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

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)\):

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

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

conduction_pb1
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
conduction_pb1_comparison
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).

t_heatsource

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:

\begin{align}
\frac{\partial \rho}{\partial t} + \nabla \cdot \left( \rho \boldsymbol{u} \right) = 0, \tag{1} \label{eq:mass}
\end{align}
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:

\begin{align}
\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}
\end{align}

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:

\begin{align}
-\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}
\end{align}

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)

\begin{align}
&\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}
\end{align}

  • Internal energy (sensibleInternalEnergy)

\begin{align}
&\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}
\end{align}

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:
\begin{align}
\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}
\end{align}

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

\begin{align}
\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}
\end{align}

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:

References

[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

flangeT

Convection
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.

References

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:

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

\begin{align}
\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}
\end{align}

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:

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

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:

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

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:

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

\begin{align}
\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}
\end{align}

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

\begin{align}
\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}
\end{align}

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

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

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

\begin{align}
\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}
\end{align}

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:

\begin{align}
&-\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}
\end{align}

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

\begin{align}
\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}
\end{align}

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\):

Reference

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

Please wait till the next update!