Difference between revisions of "BubbleFoam"

Revision as of 05:18, 12 February 2010

Valid versions:

This document is incomplete. The Author(s) is/are still working on it and the available content is not final and might not be proofread.

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} \rho_{\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 [10]. 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{b}},$

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 \textrm{U},$

while the virtual mass force term is evaluated as

$\mathbf{M}_{\varphi,\textrm{vm}} = \alpha_{\textrm{a}} \alpha_{\textrm{b}} \left( \alpha_{\textrm{b}} C_{\textrm{a,vm}} \rho_{\textrm{b}} + \alpha_{\textrm{a}} C_{\textrm{b,vm}} \rho_{\textrm{a}} \right) \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 [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 brieﬂy summarized in this chapter.

4.1 Numerical methodology

4.1.1 Phase momentum equation

4.1.2 Phase continuity equation

4.1.3 Pressure equation

4.2 Solution procedure

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

Solve the phase continuity equations
Update lift, drag and virtual mass coeﬃcients
Construct the momentum equation matrix
Predict the phase velocity ﬁelds, without considering the pressure gradient at this stage
Solve the pressure equation
Correct the velocities with the new pressure ﬁeld, and update the phase fractions
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
 
    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 ﬂ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)
    {
        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

5.2 Initial conditions

5.3 Solver configuration

5.3.1 The environmentalProperties dictionary

5.3.2 The RASProperties dictionary

5.3.3 The transportProperties dictionary

5.4 Solver controls

5.4.1 The controlDict dictionary

5.4.2 The fvSchemes dictionary

5.4.3 The fvSolution dictionary

6 References

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

--Alberto 05:10, 12 February 2010 (UTC)

@@ Line 24: / Line 24: @@
 *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 <math>\kappa − \varepsilon</math> model and standard wall functions.
+*Turbulence modelling through <math>\kappa - \varepsilon</math> model and standard wall functions.
 === bubbleFoam limitations ===