Difference between revisions of "BubbleFoam"

Revision as of 04:28, 12 February 2010

Valid versions:

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 κ − ε 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

4.3.2 Implementation of the phase continuity equation

4.3.3 Implementation of the pressure equation

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.

@@ Line 161: / Line 161: @@
 #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
+#Correct the velocities with the new pressure ﬁeld, and update the phase fractions
-the phase fractions
 #Solve the transport equations for the turbulence quantities