From OpenFOAMWiki

Valid versions: OF version 16.png OF version 17.png

1 Introduction and applications

The bubbleFoam solver is a two-phase solver based on the Euler-Euler two-fluid methodology [1, 2, 3, 4, 10], suitable to compute dispersed gas-liquid and liquid-liquid flows. In the Euler-Euler two-fluid 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-fluid 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-fluid equations derived in [10, 8] for the simulation of gas-liquid flows. The model undergoes the following assumptions:

  • Phases are incompressible
  • The dispersed phase particle diameter is constant
  • The flow 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 flows with strong density ratio
  • Robust solution algorithm, able to deal with complete flow 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 coefficient is computed as a blend of the drag coefficients 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-fluid equations solved by the bubbleFoam solver are described in this chapter. The general two-fluid 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-fluid 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 flow 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 defined, 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 fluid 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, defined 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 [10]. In particular, indicating with a and b the two phases considered in the two-fluid 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 coefficients 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),


 \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}}


\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 flows [5]. 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 briefly 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 [10] to overcome the instability arising when the phase volume fraction tends to become zero, which is a frequent situation is fully separated flows. 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}}


\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}},


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},


\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}}}

on cell faces, leading to

\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 coefficients
  3. Construct the momentum equation matrix
  4. Predict the phase velocity fields, without considering the pressure gradient at this stage
  5. Solve the pressure equation
  6. Correct the velocities with the new pressure field, 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 clarified 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
    volTensorField Rca = -nuEffa*(fvc::grad(Ua)().T());
    Rca = Rca + (2.0/3.0)*sqr(Ct)*I*k - (2.0/3.0)*I*tr(Rca);
    surfaceScalarField phiRa =
      - fvc::interpolate(nuEffa)
       *mesh.magSf()*fvc::snGrad(alpha)/fvc::interpolate(alpha + 
    // 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::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)
    volTensorField Rcb = -nuEffb*fvc::grad(Ub)().T();
    Rcb = Rcb + (2.0/3.0)*I*k - (2.0/3.0)*I*tr(Rcb);
    surfaceScalarField phiRb =
       - fvc::interpolate(nuEffb)
        *mesh.magSf()*fvc::snGrad(beta)/fvc::interpolate(beta + 
    UbEqn = 
        (scalar(1) + Cvm*rhob*alpha/rhob)*
          + 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)

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 = "
        << (
        << 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::div(phi, alpha, scheme)
           + fvm::div(-fvc::flux(-phir, beta, scheme), alpha, scheme)
        // Under-relaxation is applied to the dispersed phase continuity equation
        // The dispersed phase continuity equation is solved
        // 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::div(phi, beta, scheme)
            + fvm::div(-fvc::flux(phir, scalar(1) - beta, scheme), 
               beta, scheme)
        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 flows 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)
        if (isA<zeroGradientFvPatchScalarField>(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 + 
    //  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
        if (nonOrth == nNonOrthCorr)
            // The surface normal pressure gradient is computed
            surfaceScalarField SfGradp = pEqn.flux()/Dp;
            // The phase fluxes and the total flux are corrected to account for the updated
            // pressure field. 
            phia -= rUaAf*SfGradp/rhoa;
            phib -= rUbAf*SfGradp/rhob;
            phi = alphaf*phia + betaf*phib;
            // Pressure is explicitly relaxed before correcting the phase velocity fields
            //  The pressure surface normal gradient is recomputed.
            SfGradp = pEqn.flux()/Dp;
            // Velocities are not updated using directly the pressure, as done conventionally, 
            // but reconstructing the velocity update from the fluxes, which are the primary
            // variable [8]
            Ua += (fvc::reconstruct(phiDraga - rUaAf*SfGradp/rhoa));
            //Ua += rUaA*(fvc::reconstruct(phiDraga/rUaAf - 
            //      SfGradp/rhoa));
            Ub += (fvc::reconstruct(phiDragb - rUbAf*SfGradp/rhob));
            //Ub += rUbA*(fvc::reconstruct(phiDragb/rUbAf - 
	    //      SfGradp/rhob));
            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.

|-- 0
|   |-- Ua
|   |-- Ub
|   |-- alpha
|   |-- epsilon
|   |-- k
|   `-- p
|-- constant
|   |-- RASProperties
|   |-- g
|   |-- polyMesh
|   |   |-- blockMeshDict
|   |   `-- boundary
|   `-- transportProperties
`-- system
    |-- controlDict
    |-- fvSchemes
    `-- fvSolution
The directory
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

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

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

    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;
    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

6 References

  1. D. A. Drew. Averaged equations for two-phase flows. Studies in Applied Mathematics, L(3):205 – 231, 1971.
  2. D. A. Drew. Continuum modeling of two-phase flows. In R. Meyer, editor, Theory of dispersed multiphase flow, pages 173 – 190. Academic Press, 1983.
  3. H. Enwald, E. Peirano, and A. E. Almstedt. Eulerian two-phase flow theory applied to fluidization. International Journal of Multiphase Flow, 22:21 – 66, 1996.
  4. D. P. Hill. The computer simulation of dispersed two-phase flow. PhD thesis, Imperial College of Science, Technology and Medicine, London, U.K., 1998.
  5. J. P. Oliveira and R. I. Issa. Numerical aspects of an algorithm for Eulerian simulation of two-phase flows. International Journal for Numerical Methods in Fluids, 43:1177 – 1198, 2003.
  6. OpenCFD. OpenFOAM - The Open Source CFD Toolbox - Programmer’s Guide. OpenCFD Ltd., United Kingdom, 1.6 edition, 2009.
  7. OpenCFD. OpenFOAM - The Open Source CFD Toolbox - User’s Guide. OpenCFD Ltd., United Kingdom, 1.6 edition, 2009.
  8. H. Rusche. Computational fluid dynamics of dispersed two-phase flows at high phase fractions. PhD thesis, Imperial College of Science, Technology and Medicine, London, 2002.
  9. L. Schiller and A. Naumann. Uber die grundlegenden Berechnungen bei der Schwerkraftaufbereitung. Zeitschrift des Vereins deutscher Ingenieure, 77(12):318, 1933.
  10. H. G. Weller. Derivation, modelling and solution of the conditionally averaged two-phase flow equations. Technical report, OpenCFD Ltd., United Kingdom, 23 February 2005.

--Alberto 03:29, 17 February 2010 (UTC)