BubbleFoam

1 Introduction and applications

The bubbleFoam solver is a two-phase solver based on the Euler-Euler two-ﬂuid methodology [1, 2, 3, 4, 10], suitable to compute dispersed gas-liquid and liquid-liquid ﬂows. In the Euler-Euler two-ﬂuid approach, the phases are treated as interpenetrating continua, which are capable of exchanging properties, like momentum, energy and mass one with the other. Typical applications of the two-ﬂuid approach, as implemented in bubbleFoam, are:

• bubble columns
• stirred tank reactors
• static mixers

2 bubbleFoam capabilities and limitations

2.1 bubbleFoam capabilities

The bubbleFoam solver implements the two-ﬂuid equations derived in [10, 8] for the simulation of gas-liquid ﬂows. The model undergoes the following assumptions:

• Phases are incompressible
• The dispersed phase particle diameter is constant
• The ﬂow is isothermal
• Only momentum exchange is accounted for in the momentum transport equations

The main features of the solver are the following:

• Capability to solve for dispersed two-phase ﬂows with strong density ratio
• Robust solution algorithm, able to deal with complete ﬂow separation
• Turbulence modelling through $\kappa - \varepsilon$ model and standard wall functions.

2.2 bubbleFoam limitations

The bubbleFoam solver currently has the following limitations:

• Only one dispersed phase and a continuous phase can be described. It is not possible to account for multiple dispersed phases (i.e. represent a dispersed phase diameter distribution)
• The diameter of the particles1 constituting the dispersed phase is assumed to be constant. Aggragation, breakage and coalescence phenomena are not accounted for
• The drag coeﬃcient is computed as a blend of the drag coeﬃcients evaluated for each phase on the basis of the phase fractions, and no alternative drag models are available
• The interaction between the phases happens only through the momentum exchange term in the corresponding momentum equations:
• It is not possible to model the heat transfer between the phases
• It is not possible to model the mass transfer between the phases
• No chemical reaction model is available.

3 bubbleFoam theory

3.1 Overview

The two-ﬂuid equations solved by the bubbleFoam solver are described in this chapter. The general two-ﬂuid continuity and momentum equations are introduced, the closure expressions for the phase stress tensor and the momentum transfer term are presented, and the transport equations of the turbulence model are summarized.

3.2 Governing equations

In the two-ﬂuid approach, a continuity and a momentum equation are solved for each phase present in the system. These equations can be derived by conditionally averaging the single phase ﬂow equations. The reader interested in the theoretical background is invited to refer to, for example, [1, 2, 3, 10]. The continuity equations for each phase $\varphi$ has the form $\frac{\partial}{\partial t} \left( \alpha_{\varphi} \rho_{\varphi} \right) +\nabla \cdot \left( \alpha_{\varphi} \rho_{\varphi} \mathbf{U}_{\varphi} \right) = 0$

where $\alpha_\varphi$ represents the phase fraction of phase $\varphi$, $\rho_\varphi$ is the density of the material constituting the same phase, and $\mathbf{U}_\varphi$ is the phase velocity. The phase momentum equation is $\frac{\partial}{\partial t} \left( \alpha_{\varphi} \rho_{\varphi} \mathbf{U}_\varphi \right) + \nabla \cdot \left( \alpha_{\varphi} \rho_{\varphi} \mathbf{U}_\varphi \mathbf{U}_\varphi \right) + \nabla \cdot \alpha_{\varphi} \boldsymbol{\tau}_{\varphi} + \nabla \cdot \left( \alpha_{\varphi} \mathbf{R}_{\varphi} \right) = - \alpha_{\varphi} \nabla p + \alpha_{\varphi} \rho_{\varphi} \mathbf{g} + \mathbf{M}_{\varphi},$

in which $\tau_\varphi$ is the phase laminar stress tensor, assumed to be Newtonian, $\mathbf{R}_{\varphi}$ is the phase Reynolds stress tensor, $p$ is the pressure, $\mathbf{g}$ is the gravitational acceleration vector, and $\mathbf{M}_{\varphi}$ is the momentum exchange term. The laminar stress tensor is deﬁned, for each phase, as $\boldsymbol{\tau}_{\varphi} = - \rho_{\varphi} \nu_{\varphi} \left[\nabla \mathbf{U}_{\varphi} + \nabla^{\textrm{T}} \mathbf{U}_{\varphi} \right] + \frac{2}{3}\rho_{\varphi}\nu_{\varphi} \left( \nabla \cdot \mathbf{U}_{\varphi} \right) \mathbf{I}$

where $\nu_\varphi$ is the molecular kinematic viscosity of the ﬂuid constituting phase $\varphi$, and $\mathbf{I}$ is the identity matrix. The phase Reynolds stress tensor is given by $\mathbf{R}_{\varphi} = - \rho_{\varphi} \nu_{\varphi,\textrm{t}} \left[ \nabla \mathbf{U}_{\varphi} +\nabla^{\textrm{T}} \mathbf{U}_{\varphi} \right] + \frac{2}{3} \rho_{\varphi} \nu_{\varphi,\textrm{t}} \left( \nabla \cdot \mathbf{U}_{\varphi} \right) \mathbf{I} + \frac{2}{3} \rho_{\varphi} \kappa_{\varphi} \mathbf{I}$

where $\kappa_\varphi$ is the phase turbulent kinetic energy (Note that in the bubbleFoam implementation, it assumed that $\kappa_a = \kappa_b$ . In other words, the turbulent kinetic energy is identical for both the phases present in the system), and $\nu_{\varphi, t}$ is the phase turbulent kinematic viscosity, deﬁned as $\nu_{\varphi,\textrm{t}} = C_{\mu} \frac{\kappa_{\varphi}^{2}}{\varepsilon_{\varphi}},$

in which $C_\nu$ is a constant, and $\varepsilon_\varphi$ is the phase turbulent dissipation rate. The phase effective viscosity is calculated as the sum of the phase molecular viscosity and of the phase turbulent viscosity, as $\nu_{\varphi,\textrm{eff}} = \nu_{\varphi} + \nu_{\varphi,\textrm{t}}.$

The momentum exchange term can be decomposed in a drag contribution, a lift force contribution and a virtual mass contribution $\mathbf{M}_{\varphi} = \mathbf{ \mathbf{M}_{\varphi,\textrm{drag}}} + \mathbf{\mathbf{M}_{\varphi,\textrm{lift}}} + \mathbf{\mathbf{M}_{\varphi,\textrm{vm}}}$

where the terms are modelled according to the mixture model . In particular, indicating with a and b the two phases considered in the two-ﬂuid model, where a represents the dispersed phase, the drag term is described by equation $\mathbf{M}_{\textrm{a},\textrm{drag}} = \frac{3}{4} \alpha_{\textrm{a}} \alpha_{\textrm{b}} \left( \alpha_{\textrm{b}} \frac{C_{\textrm{D,a}}\rho_{\textrm{b}}}{\textrm{d}_{\textrm{a}}} + \alpha_{\textrm{a}} \frac{C_{\textrm{D}\textrm{b}}\rho_{\textrm{a}}}{\textrm{d}_{\textrm{b}}}\right)\left|\mathbf{U}_{\textrm{r}}\right|\mathbf{U}_{\textrm{r}},$

in which $\textrm{d}_{\textrm{a}}$ and $\textrm{d}_{\textrm{b}}$ are the phase particle diameters, $\mathbf{U}_\textrm{r} = \mathbf{U}_\textrm{a} - \mathbf{U}_\textrm{b}$ is the relative velocity vector, and $C_{\textrm{D,a}}$ and $C_{\textrm{D,b}}$ are the drag coeﬃcients computed with respect to each phase, according to $C_{\textrm{D},\varphi} = \frac{24} {\operatorname{Re}_{\varphi}} \left( 1 + 0.15 \operatorname{Re}_{\varphi}^{0.687} \right),$

where $\operatorname{Re}_{\varphi}=|\mathbf{U}_{\varphi}|\textrm{d}_{\varphi}/\nu_{\varphi}$.

The lift for term is modelled as (Notice that $\mathbf{U}$ should actually be $\mathbf{U}_\textrm{b}$. We report $\mathbf{U}$ for consistency with the code implementation.} $\mathbf{M}_{\varphi,\textrm{lift}} = \alpha_{\textrm{a}} \alpha_{\textrm{b}} \left( \alpha_{\textrm{b}} C_{\textrm{a,lift}} \rho_{\textrm{b}} + \alpha_{\textrm{a}} C_{\textrm{b,lift}} \rho_{\textrm{a}} \right) \mathbf{U}_{\textrm{r}} \times \nabla \times \mathbf{U},$

while the virtual mass force term is evaluated as $\mathbf{M}_{\varphi,\textrm{vm}} = \alpha_{\textrm{a}} \alpha_{\textrm{b}} C_{\textrm{vm}} \rho_{\textrm{b}} \left( \left.\frac{\textrm{d}\mathbf{U}_{\textrm{b}}}{\textrm{d} t} \right|_{\textrm{b}} - \left. \frac{\textrm{d} \mathbf{U}_{\textrm{a}}}{\textrm{d} t} \right|_{\textrm{a}} \right),$

where $\left.\frac{\textrm{d} \mathbf{U}_{\textrm{a}}}{\textrm{d} t} \right|_{\textrm{a}} = \frac{\partial \mathbf{U}_{\textrm{a}}}{\partial t} + \mathbf{U}_{\textrm{a}} \cdot \nabla \mathbf{U}_{\textrm{a}}$

and $\left.\frac{\mathbf{d} \mathbf{U}_{\textrm{b}}}{\textrm{d} t} \right|_{\textrm{b}} = \frac{\partial \mathbf{U}_{\textrm{b}}}{\partial t} + \mathbf{U}_{\textrm{b}} \cdot \nabla \mathbf{U}_{\textrm{b}}$.

3.2.1 Turbulence model

The bubbleFoam solver uses a two-equation $\kappa-\varepsilon$ turbulence model for the continuous phase, and accounts for the influence of the turbulence on the dispersed phase by scaling the dispersed phase turbulent viscosity. The equation for the turbulent kinetic energy of the continuous phase ( $\textrm{b}$) reads $\frac{\partial}{\partial t} \left( \alpha_{\textrm{b}} \kappa_{\textrm{b}} \right) + \nabla \cdot \left( \alpha_{\textrm{b}} \mathbf{U}_{\textrm{b}} \kappa_{\textrm{b}} \right) - \nabla^2 \left( \sigma_{\kappa} \nu_{\textrm{b},\textrm{eff}} \kappa_{\textrm{b}} \right) = \alpha_{\textrm{b}} G - \alpha_{\textrm{b}} \varepsilon_{\textrm{b}}$

where $G$ is the producton of the turbulent kinetic energy, given by $G = 2 \nu_{\textrm{b},\textrm{t}} \left[ \nabla \mathbf{U}_{\textrm{b}} \cdot \operatorname{dev} \left( \nabla \mathbf{U}_{\textrm{b}} + \nabla^{\textrm{T}} \mathbf{U}_{\textrm{b}} \right) \right],$ and $\sigma_{\kappa}$ is the turbulent Schmidt number. The turbulent dissipation rate is determined by solving the transport equation $\frac{\partial}{\partial t} \left( \alpha_{\textrm{b}} \varepsilon_{\textrm{b}} \right) + \nabla \cdot \left( \alpha_{\textrm{b}} \mathbf{U} \varepsilon_{\textrm{b}} \right) - \nabla^2 \left( \sigma_{\varepsilon} \nu_{\textrm{b},\textrm{eff}} \varepsilon_{\textrm{b}}\right) = C_1 \alpha_{\textrm{b}} G \frac{\varepsilon_{\textrm{b}}}{\kappa_{\textrm{b}}} - C_2 \alpha_{\textrm{b}} \frac{\varepsilon_{\textrm{b}}^2}{\kappa_{\textrm{b}}}$

where $C_{1}$ and $C_{2}$ are constant of the turbulence model. The turbulence viscosity of the continuous phase is calculated from its definition $\nu_{\textrm{b},\textrm{t}} = C_{\mu} \frac{\kappa^2}{\varepsilon},$

while the turbulent viscosity of the dispersed phase is evaluated as $\nu_{\textrm{a},\textrm{t}} = C_{\textrm{t}}^2 \nu_{\textrm{b},\textrm{t}},$

being the coefficient $C_{\mu}$ and the turbulence response coefficient $C_{\textrm{t}}$ constants of the model. Standard wall functions are adopted to treat the zone next to the wall.

4 bubbleFoam implementation

The numerical solution of the two-phase equations relies on a segregated algorithm based on the PISO procedure extended to two-phase ﬂows . The momentum equations are manipulated to stabilize the system of equation at the limits of the range of volume fractions, to avoid singularities, as suggested in [5, 10]. The details of the numerical methodology are brieﬂy summarized in this chapter.

4.1 Numerical methodology

4.1.1 Phase momentum equation

The numerical method used in the bubbleFoam solver relies on the phase-intensive formulation of the phase momentum equations proposed in  to overcome the instability arising when the phase volume fraction tends to become zero, which is a frequent situation is fully separated ﬂows. According to this approach, assuming the phase density $\rho_\varphi$ to be constant, the momentum equation is rewritten in non-conservative form to extract the volume fraction from the transport terms, leading to $\frac{\partial \mathbf{U}_{\varphi}}{\partial t} + \mathbf{U}_{\varphi} \cdot \nabla \mathbf{U}_{\varphi} + \nabla \cdot \left( \frac{\boldsymbol{\tau}_{\varphi}}{\rho_{\varphi}} + \mathbf{R}_{\varphi} \right) + \frac{\nabla \alpha_{\varphi}}{\alpha_{\varphi}} \cdot \left( \frac{\boldsymbol{\tau}_{\varphi}}{\rho_{\varphi}} + \mathbf{R}_{\varphi} \right) = -\frac{\nabla p}{\rho_{\varphi}} + \mathbf{g} + \frac{\mathbf{M}_{\varphi}}{\alpha_{\varphi} \rho_{\varphi}}.$

By renaming the total stress tensor as $\mathbf{R}_{\varphi,\textrm{eff}} = \frac{\mathbf{\tau}_{\varphi}}{\rho_{\varphi}} + \mathbf{R}_{\varphi},$

and decomposing it in a diffusive component and a correction term $\mathbf{R}_{\varphi,\textrm{eff}} = \mathbf{R}_{\varphi,\textrm{eff}}^{\textrm{D}} + \mathbf{R}_{\varphi,\textrm{eff}}^{\textrm{C}}$

with $\mathbf{R}_{\varphi,\textrm{eff}}^{\textrm{D}} = -\nu_{\varphi,\textrm{eff}} \nabla \mathbf{U}_{\varphi},$ $\mathbf{R}_{\varphi,\textrm{eff}}^{\textrm{C}} = \mathbf{R}_{\varphi,\textrm{eff}} + \nu_{\varphi,\textrm{eff}} \nabla \mathbf{U}_{\varphi},$

we obtain $\frac{\partial \mathbf{U}_{\varphi}}{\partial t} + \mathbf{U}_{\varphi} \cdot \nabla \mathbf{U}_{\varphi} - \nabla \cdot \left( \nu_{\varphi,\textrm{eff}} \nabla \mathbf{U}_{\varphi} \right) + \nabla \cdot \mathbf{R}_{\varphi,\textrm{eff}}^{\textrm{C}} - \nu_{\varphi,\textrm{eff}} \frac{\nabla \alpha_{\varphi}}{\alpha_{\varphi}} \cdot \nabla \mathbf{U}_{\varphi} + \frac{\nabla \alpha_{\varphi}}{\alpha_{\varphi}} \cdot \mathbf{R}_{\varphi,\textrm{eff}}^{\textrm{C}} = - \frac{\nabla p}{\rho_{\varphi}} + \mathbf{g} + \frac{\mathbf{M}_{\varphi}}{\alpha_{\varphi} \rho_{\varphi}}.$

Finally, introducing the total phase velocity $\mathbf{U}_{\varphi}^{\textrm{T}} = \mathbf{U}_{\varphi} - \nu_{\varphi,\textrm{eff}} \frac{\nabla \alpha_{\varphi}}{\alpha_{\varphi}},$

introducing $A_{\textrm{d}} = \frac{3}{4} \left( \alpha_{\textrm{b}} \frac{C_{\textrm{D,a}}\rho_{\textrm{b}}}{\textrm{d}_{\textrm{a}}} + \alpha_{\textrm{a}} \frac{C_{\textrm{D}\textrm{b}}\rho_{\textrm{a}}}{\textrm{d}_{\textrm{b}}}\right)\left|\mathbf{U}_{\textrm{r}}\right|,$ $A_{\textrm{l}} = \left( \alpha_{\textrm{b}} C_{\textrm{a,lift}} \rho_{\textrm{b}} + \alpha_{\textrm{a}} C_{\textrm{b,lift}} \rho_{\textrm{a}} \right) \mathbf{U}_{\textrm{r}} \times \nabla \times \mathbf{U},$ $A_{\textrm{vm}} = C_{\textrm{vm}} \rho_{\textrm{b}}\,$

and considering, as an example, phase a, the momentum equation becomes $\frac{\partial \mathbf{U}_{\textrm{a}}}{\partial t} + \mathbf{U}_{\textrm{a}}^{\textrm{T}} \cdot \nabla \mathbf{U}_{\textrm{a}} - \nabla \cdot \left( \nu_{\textrm{a},\textrm{eff}} \nabla \mathbf{U}_{\varphi} \right) + \nabla \cdot \mathbf{R}_{\varphi,\textrm{eff}}^{\textrm{C}} + \frac{\nabla \alpha_{\textrm{a}}}{\alpha_{\textrm{a}}} \cdot \mathbf{R}_{\textrm{a},\textrm{eff}}^{\textrm{C}} = - \frac{\nabla p}{\rho_{\textrm{a}}} + \mathbf{g} + \frac{\alpha_{\textrm{b}}}{\rho_{\textrm{a}}} \left[ A_{\textrm{d}} \left( \mathbf{U}_{\textrm{b}} - \mathbf{U}_{\textrm{a}} \right) + A_{\textrm{l}} + A_{\textrm{vm}} \left( \left.\frac{\textrm{d}\mathbf{U}_{\textrm{b}}}{\textrm{d} t} \right|_{\textrm{b}} - \left. \frac{\textrm{d} \mathbf{U}_{\textrm{a}}}{\textrm{d} t} \right|_{\textrm{a}} \right) \right ]$

The momentum equation can be furtherly rewritten by introducing the total convective momentum fluxes $\phi_{\textrm{a}}^{\textrm{T}} = \phi_{\textrm{a}} - \nu_{\textrm{a},\textrm{eff},\textrm{f}} \frac{S_{\textrm{f}} \nabla_{\textrm{f}}^{\perp}\alpha_{\textrm{a}}}{\alpha_{\textrm{a},\textrm{f}} + \delta},$

and $\phi_{\textrm{b}}^{\textrm{T}} = \phi_{\textrm{b}} - \nu_{\textrm{b},\textrm{b},\textrm{f}} \frac{S_{\textrm{f}} \nabla_{\textrm{f}}^{\perp}\alpha_{\textrm{b}}}{\alpha_{\textrm{b},\textrm{f}} + \delta},$

where $\nabla_{\textrm{f}}^{\perp}\alpha_{\varphi}$ is the surface normal gradient to the surface $S_{\textrm{f}}$ of the phase fraction, the subscript $\textrm{f}$ indicates values computed at cell faces, and $\delta$ is a small number. Once this substitution is performed, the momentum equation can be discretized. In particular

• the unsteady, the convective and diffusive terms are treated fully implicitly
• the stress term correction is handles explicitly as well as the lift terms
• the drag and virtual mass terms are handled semi-implicitly, treating implicitly the part containing the velocity of the phase the equation is written for, and moving it to the pressure equation, and explicitly the other part
• the pressure gradient is not included in the momentum equation directly, but its effect is accounted for only when the phase velocities are corrected
• the same procedure adopted for the pressure gradient is used for the gravity, which does not appear in the momentum equation in the code, but its effect is accounted for in the pressure equation

The discretized equations, can be represented as $\mathbf{U}_{\textrm{a}} = \frac{\mathbf{H}_{\textrm{a}}}{A_{\textrm{a}}} - \frac{\nabla p}{\rho_{\textrm{a}} A_{\textrm{a}}},$ $\mathbf{U}_{\textrm{b}} = \frac{\mathbf{H}_{\textrm{b}}}{A_{\textrm{b}}} - \frac{\nabla p}{\rho_{\textrm{b}} A_{\textrm{b}}},$

where $A_{\varphi}$ is the diagonal part of the matrix originating from the discretization of the phase momentum equation and $\mathbf{H}_{\textrm{a}}$ is the remaining part. These equations will be used to correct the velocity after the pressure field is updated.

4.1.2 Phase continuity equation

The phase continuity equations have to be solved ensuring that the phase fraction of each phase is kept between zero and one. To obtain this result, the dispersed phase continuity equation is rewritten as a function of the mean and relative velocity. Rewriting the phase velocity in terms of the relative velocity $\mathbf{U}_{\textrm{r}} = \mathbf{U}_{\textrm{a}} - \mathbf{U}_{\textrm{b}},$

and the mean velocity weighted on the phase fractions $\mathbf{U} = \alpha_{\textrm{a}} \mathbf{U}_{\textrm{a}} + \alpha_{\textrm{b}} \mathbf{U}_{\textrm{b}},$

we find $\mathbf{U}_{\textrm{a}} = \mathbf{U} + \alpha_{\textrm{b}} \mathbf{U}_{\textrm{r}}.$

Substituting into the phase continuity equation, a new expression is obtained $\frac{\partial \alpha_{\textrm{a}}}{\partial t} + \nabla \cdot \left( \alpha_{\textrm{a}} \mathbf{U} \right) + \nabla \cdot \left( \alpha_{\textrm{a}} (1-\alpha_{\textrm{a}}) \mathbf{U}_{\textrm{r}}\right) = 0,$

which, if iteratively (Note that the equation is not linear) solved in a fully implicit manner, provides a bounded solution for the phase fraction field.

The volume fraction of the continuous phase can be computed as $\alpha_{\textrm{b}} = 1 - \alpha_{\textrm{a}}$. However in some cases, this approach could not be satisfactory from a numerical point of view, and the solution of a phase continuity equation also for the continuous phase is recommended, in order to increase the convergence rate. In such a case, after solving a continuity equation for phase b, the boundness of the dispersed phase volume fraction is ensured by evaluating the new value of the dispersed phase fraction as $\alpha_{\textrm{a},\textrm{bounded}} = \frac{1}{2} \left[ 1 - \left( 1 - \alpha_{\textrm{a}} \right)^2 + \left( 1 - \alpha_{\textrm{b}} \right)^2 \right],$

and recomputing the continuous phase fraction as $\alpha_{\textrm{b},\textrm{bounded}} = 1 - \alpha_{\textrm{a}, \textrm{bounded}}$.

4.1.3 Pressure equation

The pressure equation is obtained imposing that the divergence of the mixture flux is zero $\nabla \cdot \phi = \nabla \cdot \left( \alpha_{\textrm{a},\textrm{f}} \phi_{\textrm{a}} + \alpha_{\textrm{b},\textrm{f}} \phi_{\textrm{b}} \right) = 0,$.

The phase fluxes are obtained by interpolating equations $\mathbf{U}_{\textrm{a}} = \frac{\mathbf{H}_{\textrm{a}}}{\textbf{A}_{\textrm{a}}} - \frac{\nabla p}{\rho_{\textrm{a}} A_{\textrm{a}}}$ $\mathbf{U}_{\textrm{b}} = \frac{\mathbf{H}_{\textrm{b}}}{\textbf{A}_{\textrm{b}}} - \frac{\nabla p}{\rho_{\textrm{b}} A_{\textrm{b}}}$ $\phi_{\textrm{a}} = \phi_{\textrm{a}}^* - \left. \frac{1}{\rho_{\textrm{a}} A_{\textrm{a}}}\right|_{\textrm{f}} S_{\textrm{f}} \nabla_{\textrm{f}}^{\perp} p + \phi_{\textrm{dg},\textrm{a}}$ $\phi_{\textrm{b}} = \phi_{\textrm{b}}^* - \left. \frac{1}{\rho_{\textrm{b}} A_{\textrm{b}}}\right|_{\textrm{f}} S_{\textrm{f}} \nabla_{\textrm{f}}^{\perp} p + \phi_{\textrm{dg},\textrm{b}}$ $\phi_{\textrm{a}}^* = \left. \frac{\mathbf{H}_{\textrm{a}}}{\mathbf{A}_{\textrm{a}}}\right|_{\textrm{f}} \cdot \mathbf{S}_{\textrm{f}},$ $\phi_{\textrm{b}}^* = \left. \frac{\mathbf{H}_{\textrm{b}}}{\mathbf{A}_{\textrm{b}}}\right|_{\textrm{f}} \cdot \mathbf{S}_{\textrm{f}},$

The effect of the gravity and of the explicit part of the drag are included in the fluxes as $\phi_{\textrm{dg},\textrm{a}} = \left. \frac{\alpha_{\textrm{b}} C_{\textrm{D}}}{\rho_{\textrm{a}} A_{\textrm{a}}} \right|_{\textrm{f}} \phi_{\textrm{b}}^* + \frac{\mathbf{g} \cdot \mathbf{S}_{\textrm{f}}}{A_{\textrm{a}}},$ $\phi_{\textrm{dg},\textrm{b}} = \left. \frac{\alpha_{\textrm{a}} C_{\textrm{D}}}{\rho_{\textrm{b}} A_{\textrm{b}}} \right|_{\textrm{f}} \phi_{\textrm{a}}^* + \frac{\mathbf{g} \cdot \mathbf{S}_{\textrm{f}}}{A_{\textrm{b}}}$.

4.2 Solution procedure

The solution procedure adopted in the bubbleFoam solver is summed up as follows:

1. Solve the phase continuity equations
2. Update lift, drag and virtual mass coeﬃcients
3. Construct the momentum equation matrix
4. Predict the phase velocity ﬁelds, without considering the pressure gradient at this stage
5. Solve the pressure equation
6. Correct the velocities with the new pressure ﬁeld, and update the phase fractions
7. Solve the transport equations for the turbulence quantities

4.3 Code representation

4.3.1 Implementation of the phase momentum equation

The implementation of the phase momentum equations, in its phase-intensive form is better clariﬁed in the following commented code snippet.

// The fvVectorMatrix for the two phase momentum equations are declared,
// setting the correct dimensional units

fvVectorMatrix UaEqn(Ua, Ua.dimensions()*dimVol/dimTime);
fvVectorMatrix UbEqn(Ub, Ub.dimensions()*dimVol/dimTime);

{
// The dispersed phase stress tensor is computed

Rca = Rca + (2.0/3.0)*sqr(Ct)*I*k - (2.0/3.0)*I*tr(Rca);

surfaceScalarField phiRa =
- fvc::interpolate(nuEffa)
scalar(0.001));

// The momentum predictor is set up bringing the virtual mass term
// on the LHS for convenience. The drag term is managed semi-implicitly, and the
// explicit part is treated at cell faces, moving it to the pressure equation. The implicit
// part of the drag is treated as a source term. The buoyancy term is moved to the
// pressure equation and managed at faces too
UaEqn =
(
(scalar(1) + Cvm*rhob*beta/rhoa)*
(
fvm::ddt(Ua)
+ fvm::div(phia, Ua, "div(phia,Ua)")
- fvm::Sp(fvc::div(phia), Ua)
)

- fvm::laplacian(nuEffa, Ua)
+ fvc::div(Rca)

+ fvm::div(phiRa, Ua, "div(phia,Ua)")
- fvm::Sp(fvc::div(phiRa), Ua)
+ (fvc::grad(alpha)/(fvc::average(alpha) + scalar(0.001)) & Rca)
==
//  g                   // Buoyancy term transfered to p-equation
- fvm::Sp(beta/rhoa*dragCoef, Ua)
//+ beta/rhoa*dragCoef*Ub// Explicit drag transfered to p-equation
- beta/rhoa*(liftCoeff - Cvm*rhob*DDtUb)
);

UaEqn.relax();

Rcb = Rcb + (2.0/3.0)*I*k - (2.0/3.0)*I*tr(Rcb);

surfaceScalarField phiRb =
- fvc::interpolate(nuEffb)
scalar(0.001));

UbEqn =
(
(scalar(1) + Cvm*rhob*alpha/rhob)*
(
fvm::ddt(Ub)
+ fvm::div(phib, Ub, "div(phib,Ub)")
- fvm::Sp(fvc::div(phib), Ub)
)

- fvm::laplacian(nuEffb, Ub)
+ fvc::div(Rcb)

+ fvm::div(phiRb, Ub, "div(phib,Ub)")
- fvm::Sp(fvc::div(phiRb), Ub)

+ (fvc::grad(beta)/(fvc::average(beta) + scalar(0.001)) & Rcb)
==
//  g                    // Buoyancy term transfered to p-equation
- fvm::Sp(alpha/rhob*dragCoef, Ub)
//+ alpha/rhob*dragCoef*Ua // Explicit drag transfered to p-eqn.
+ alpha/rhob*(liftCoeff + Cvm*rhob*DDtUa)
);

UbEqn.relax();
}

4.3.2 Implementation of the phase continuity equation

{
// A word is defined to store the discretisation scheme label used to compute the
// divergence of the dispersed phase volume fraction. The text between " " is searched
// for in the fvSchemes dictionary
word scheme("div(phi,alpha)");

// The relative flux is calculated
surfaceScalarField phir = phia - phib;

// The Courant number is computed on the base of the relative velocity
Info<< "Max Ur Courant Number = "
<< (
max
(
mesh.surfaceInterpolation::deltaCoeffs()*mag(phir)
/mesh.magSf()
)*runTime.deltaT()
).value()
<< endl;

// The phase fraction correction loop is started
for (int acorr=0; acorr<nAlphaCorr; acorr++)
{
// The dipersed phase continuity equation object is created using
// implicit discretization for all its terms
fvScalarMatrix alphaEqn
(
fvm::ddt(alpha)
+ fvm::div(phi, alpha, scheme)
+ fvm::div(-fvc::flux(-phir, beta, scheme), alpha, scheme)
);

// Under-relaxation is applied to the dispersed phase continuity equation
alphaEqn.relax();

// The dispersed phase continuity equation is solved
alphaEqn.solve();

// The continuous phase continuity equation is defined as done for the dispersed
// phase. Notice that by default this equation is commented out because only the
// equation for the dispersed phase is solved.

/*
fvScalarMatrix betaEqn
(
fvm::ddt(beta)
+ fvm::div(phi, beta, scheme)
+ fvm::div(-fvc::flux(phir, scalar(1) - beta, scheme),
beta, scheme)
);
betaEqn.relax();
betaEqn.solve();

alpha = 0.5*(scalar(1) + sqr(scalar(1) - beta)
- sqr(scalar(1) - alpha));
*/

// The continuous phase volume fraction is computed as the complement to one of
// the dispersed phase fraction.
beta = scalar(1) - alpha;
}

Info<< "Dispersed phase volume fraction = "
<< alpha.weightedAverage(mesh.V()).value()
<< "  Min(alpha) = " << min(alpha).value()
<< "  Max(alpha) = " << max(alpha).value()
<< endl;
}

rho = alpha*rhoa + beta*rhob;

4.3.3 Implementation of the pressure equation

The details of the implementation of the pressure equation are commented in the code snippet below, where the implementation of the pseudo-staggered algorithm developed by Weller [10, 8] to obtain a stable solution in completely separated ﬂows is adopted.

{
// The phase fraction fields are interpolated at cell faces
surfaceScalarField alphaf = fvc::interpolate(alpha);
surfaceScalarField betaf = scalar(1) - alphaf;

//  The central coeffcients are computed
volScalarField rUaA = 1.0/UaEqn.A();
volScalarField rUbA = 1.0/UbEqn.A();

// The values of the central coefficients are interpolated at cell faces
surfaceScalarField rUaAf = fvc::interpolate(rUaA);
surfaceScalarField rUbAf = fvc::interpolate(rUbA);

//  The phase velocities are predicted
Ua = rUaA*UaEqn.H();
Ub = rUbA*UbEqn.H();

// Accounting for the explict drag term and buoyancy effect. Both
// these terms are managed at cell faces and their effect is included in the phase flux.
// Note that the whole drag term is interpolated and then multiplied by the flux,
// instead of reusing the previously computed value of the central coefficient at cell
// faces (rUaAf).

surfaceScalarField phiDraga =
fvc::interpolate(beta/rhoa*dragCoef*rUaA)*phib +
rUaAf*(g & mesh.Sf());
surfaceScalarField phiDragb =
fvc::interpolate(alpha/rhob*dragCoef*rUbA)*phia +
rUbAf*(g & mesh.Sf());

forAll(p.boundaryField(), patchi)
{
{
phiDraga.boundaryField()[patchi] = 0.0;
phiDragb.boundaryField()[patchi] = 0.0;
}
}

// The phase fluxes and the total flux are computed
phia = (fvc::interpolate(Ua) & mesh.Sf()) +
fvc::ddtPhiCorr(rUaA, Ua, phia) + phiDraga;
phib = (fvc::interpolate(Ub) & mesh.Sf())
+ fvc::ddtPhiCorr(rUbA, Ub, phib) + phiDragb;

phi = alphaf*phia + betaf*phib;

//  The mixture central coefficient is computed
surfaceScalarField Dp("(rho*(1|A(U)))", alphaf*rUaAf/rhoa +
betaf*rUbAf/rhob);

//  The non-ortogonal correction loop is started
for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
//  The pressure equation is set up
fvScalarMatrix pEqn
(
fvm::laplacian(Dp, p) == fvc::div(phi)
);

// References for the pressure are set
pEqn.setReference(pRefCell, pRefValue);

//  The pressure equation is solved
pEqn.solve();

if (nonOrth == nNonOrthCorr)
{
// The surface normal pressure gradient is computed

// The phase fluxes and the total flux are corrected to account for the updated
// pressure field.

phi = alphaf*phia + betaf*phib;

// Pressure is explicitly relaxed before correcting the phase velocity fields
p.relax();

//  The pressure surface normal gradient is recomputed.

// Velocities are not updated using directly the pressure, as done conventionally,
// but reconstructing the velocity update from the fluxes, which are the primary
// variable 

//Ua += rUaA*(fvc::reconstruct(phiDraga/rUaAf -
Ua.correctBoundaryConditions();

//Ub += rUbA*(fvc::reconstruct(phiDragb/rUbAf -
Ub.correctBoundaryConditions();

U = alpha*Ua + beta*Ub;
}
}
}

#include "continuityErrs.H"

5 bubbleFoam setup and solution strategy

5.1 Case structure

The structure of a typical bubbleFoam case is represented in the following directory tree, which lists the directories and the files contained in the bubbleColumn tutorial.

bubbleColumn
|-- 0
|   |-- Ua
|   |-- Ub
|   |-- alpha
|   |-- epsilon
|   |-- k
|   -- p
|-- constant
|   |-- RASProperties
|   |-- g
|   |-- polyMesh
|   |   |-- blockMeshDict
|   |   -- boundary
|   -- transportProperties
-- system
|-- controlDict
|-- fvSchemes
`-- fvSolution
The directory
0
contains six files named as the corresponding variables in the solver:
• Ua
dispersed phase velocity
• Ub
continuous phase velocity
• alpha
dispersed phase fraction
• epsilon
turbulent dissipation rate
• k
turbulent kinetic energy
• p
pressure.

Each of these files contain the initial and boundary conditions for the corresponding property field.

The
constant
directory contains the physical model settings and the mesh files, in particular
• RASProperties
contains the turbulence model settings
• g
contains the gravity vector definition
• polyMesh/blockMeshDict
is the mesh dictionary for blockMesh
• polyMesh/boundary
is the file containing the boundary specifications to construct the mesh
• transportProperties
contains the phase properties.
Finally, the
system
directory contains the solver settings:
• controlDict
contains information on the time stepping, the desired simulation time, and allows the setup of the flags for data storage
• fvSchemes
contains the definitions of the numerical schemes used during the simulation
• fvSolution
contains the linear solver settings and tolerances, and the definition of the PISO parameters.

Each dictionary specific to the solver is described in detail in the following sections. The reader is invited to refer to the official documentation for further information on how to build computational meshes with blockMesh, and on the numerical schemes for the different terms of the transport equations.

5.2 Initial conditions

Initial conditions are specified following the general rules adopted in OpenFOAM. The reader is invited to refer to the User's Guide, and the setFields utility documentation.

5.3 Solver configuration

The configuration of the physical model of a case to be run with the bubbleFoam solver is set up using the following dictionary files:

• environmentalProperties
• RASProperties
• transportProperties

5.3.1 The environmentalProperties dictionary

The environmentalProperties dictionary only contains the specification of the gravitational acceleration vector:

dimensions      [0 1 -2 0 0 0 0];
value           ( 0 -9.81 0 );

The user has to specify the dimensional units, the magnitude and the orientation of the gravity vectory, with respect to the reference frame in which the mesh of the system under consideration was built.

5.3.2 The RASProperties dictionary

The RASProperties dictionary contains all the information regarding the turbulence model. The first option the user has to provide is which turbulence model has to be used. This is selected by specifying the option

RASModel            laminar;

The possible choices are

• laminar
No turbulence model is used during the simulation,
• kEpsilon}
The $\kappa-\varepsilon$ mixture model is used, and wall functions are adopted to describe the zone near the wall.

The next option in the dictionary is a switch

turbulence          off;

which determines if the turbulence model is used or not during the simulation. If the turbulence model has to be used, the value of the switch has to be set to on, while setting it to off is equivalent to specify the laminar model.

The option

printCoeffs         off;

determines if the coefficients of the turbulence model are printed when the simulation is started (on), or not (off).

5.3.3 The transportProperties dictionary

The transportProperties coefficient contains the physical properties of the fluid that constitute the two phases and the coefficients used in the momentum transfer term in the momentum equation, as shown in the example below.

rhoa            rhoa [1 -3 0 0 0 0 0] 1;

rhob            rhob [1 -3 0 0 0 0 0] 1000;

nua             nua [0 2 -1 0 0 0 0] 1.6e-05;

nub             nub [0 2 -1 0 0 0 0] 1e-06;

da              da [0 1 0 0 0 0 0] 0.003;

db              db [0 1 0 0 0 0 0] 0.0001;

Cvm             Cvm [0 0 0 0 0 0 0] 0.5;

Cl              Cl [0 0 0 0 0 0 0] 0;

Ct              Ct [0 0 0 0 0 0 0] 1;

The keywords of the transportProperties dictionary are

• rhoa
material density of phase a
• rhob
material density of phase b
• nua
molecular kinematic viscosity of phase a
• nub
molecular kinematic viscosity of phase b
• da
phase diameter of phase a
• db
phase diameter of phase b
• Cvm
virtual mass coefficient
• Cl
lift force coefficient
• Ct
turbulence response time coefficient

5.4 Solver controls

The bubbleFoam solver is controlled by the standard dictionaries controlDict, fvSchemes and fvSolution. The users is invited to refer to the User's Guide for the syntax and the options generally available in these dictionaries.

However, the fvSolution dictionary requires some clarification, since it contains options specific to the solver. An example of the fvSolution dictionary is reported in the code snipped below.

solvers
{
p PCG
{
preconditioner   DIC;
tolerance        1e-10;
relTol           0;
};
Ua PBiCG
{
preconditioner   DILU;
tolerance        1e-05;
relTol           0;
};
Ub PBiCG
{
preconditioner   DILU;
tolerance        1e-05;
relTol           0;
};
alpha PBiCG
{
preconditioner   DILU;
tolerance        1e-10;
relTol           0;
};
beta PBiCG
{
preconditioner   DILU;
tolerance        1e-10;
relTol           0;
};
k PBiCG
{
preconditioner   DILU;
tolerance        1e-05;
relTol           0;
};
epsilon PBiCG
{
preconditioner   DILU;
tolerance        1e-05;
relTol           0;
};
}

PISO
{
nCorrectors     2;
nNonOrthogonalCorrectors 0;
nAlphaCorr      2;
correctAlpha    no;
pRefCell        0;
pRefValue       0;
}

The first part of the dictionary containes the usual settings for the linear solvers of the different equations. The PISO parameters, however, need to be clarified, being some of them specific to the solver:

• nCorrectors
specifies the number of PISO loops (recommended values: 2 or 3)
• nNonOrthogonalCorrectors
specified the number of corrector steps to account for non-orthogonality of the mesh
• nAlphaCorr
specifies the number of corrections to execute on the dispersed phase fraction per each PISO corrector step
• correctAlpha
specifies if the dispersed phase correction has to be corrected or not. Possible values are yes and no
• pRefCell
specifies the cell taken as reference for the pressure. The position of the cell defines the reference point for the pressure in the computational domain
• pRefValue
specified the value of the pressure at the reference cell