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!