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