BubbleFoam

Valid versions:

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),$

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

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

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 + 
       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();
 
 
    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 + 
	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 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 + 
	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
            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
            p.relax();
 
            //  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));
            Ua.correctBoundaryConditions();
 
            Ub += (fvc::reconstruct(phiDragb - rUbAf*SfGradp/rhob));
            //Ub += rUbA*(fvc::reconstruct(phiDragb/rUbAf - 
	    //      SfGradp/rhob));
            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

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)