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

*without Boussinesq approximation.*

**buoyantPimpleFoam**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:

21 22 23 24 25 26 27 28 29 30 31 32 |
solve ( UEqn == fvc::reconstruct ( ( - ghf*fvc::snGrad(rhok) - fvc::snGrad(p_rgh) )*mesh.magSf() ) ); |

The Eq. \eqref{eq:final} is calculated using the * fvc::reconstruct* function from the face flux fields and the variable

*is defined in createFields.H as \(\rho/\rho_0\):*

**rhok**
55 56 57 58 59 60 61 62 63 64 65 |
// Kinematic density for buoyancy force volScalarField rhok ( IOobject ( "rhok", runTime.timeName(), mesh ), 1.0 - beta*(T - TRef) ); |

Reference |

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

Please wait till the next update!

Dear Fumiya Nozaki,

Firstly, let me congratulate you on your blog. I am a foreign researcher living in Japan and working on OpenFOAM and I am one of the followers of your blog for a long time. I guess, this is the clearest and the most explanatory article on buoyantBoussinesqSimpleFoam. It was very easy for me to follow up line by line with original OpenFOAM source code. Thank you so much!

I just wanted to ask your opinion on the first lines of TEqn.H file. For instance, ” alphat = turbulence->nut()/Prt;” or “volScalarField alphaEff(“alphaEff”, turbulence->nu()/Pr + alphat);” parts of the source code. How can we implement these parts into math? I guess it requests turbulence file from somewhere in the source code and apply the required equation but I just want to sure, how these lines implements those properties and using which file. I am looking forward to hearing from you.

Dear Fumiya Nozaki,

This is really an informative blog. I read your posts on foamcalc. I wonder is there any utlity to find time averaging of fields after simulations. Say I ran my simulation from 0-2000s. Now after everthing is done I would like to do time averaging of velocity filed from 500-2000s. Is there any utility (not in run time using function objects) to create an average field (Umean , UPrime2Mean etc) in foamcalc.

Could you also provide info on more post processing tools available that can be used in conjunction with openfoam.

Someting like energy spectrum analysis, Proper orthogonal decomposition, etc…