# CFD Simulation of Polydispersed Bubbly Two-Phase Flow Around an Obstacle

Authors: E. Krepper; P. Ruyer; M. Beyer; D. Lucas; H.-M. Prasser; N. Seiler
Journal: Science and Technology of Nuclear Installations
 (2009)

---

## Abstract

This paper concerns the model of a polydispersed bubble population in the frame of an ensemble averaged two-phase flow formulation.
The ability of the moment density approach to represent bubble population size distribution within a multi-dimensional CFD code based on the two-fluid model is studied.
Two different methods describing the polydispersion are presented: (i) a moment density method, developed at IRSN,
to model the bubble size distribution function and (ii) a population balance method considering several different velocity fields of the gaseous phase.
The first method is implemented in the Neptune_CFD code, whereas the second method is implemented in the CFD code ANSYS/CFX.
Both methods consider coalescence and breakup phenomena and momentum interphase transfers related to drag and lift forces. Air-water bubbly flows in a vertical pipe with obstacle of the TOPFLOW experiments series performed at FZD are then used as simulations test cases. The numerical results, obtained with Neptune_CFD and with ANSYS/CFX, allow attesting the validity of the approaches. Perspectives concerning the improvement of the models, their validation, as well as the extension of their applicability range are discussed.


---

## Body

## 1. Introduction

Many flow regimes in Nuclear Reactor Safety Research are characterized
by multiphase flows, with one phase being a liquid and the other phase
consisting of gas or vapor of the liquid phase. The flow regimes found in vertical pipes are dependent on the void fraction of the gaseous phase, which
varies, as void fraction increases, from bubbly flow to slug flow, churn
turbulent flow, annular flow, and finally to droplet flow at highest void
fractions. In the regimes of bubbly and slug flows, a spectrum of different bubble sizes is observed. While dispersed bubbly flows with low gas volume fraction are mostly monodispersed, an increase of the gas volume fraction leads to a broader bubble size distribution due to breakup and coalescence of bubbles. The exchange area for mass, momentum, or heat between continuous and dispersed phases thus cannot be simply modeled
based on the knowledge of just the void fraction and a mean diameter. Moreover,
the forces acting on the bubbles may depend on their individual size  which is the case not only for drag but also for nondrag forces. Among the forces leading to lateral migration of the bubbles, that is, acting
in normal direction with respect to the main drag force, bubble lift force was
found to change the sign as the bubble size varies. Consequently, in the
context of pipe flows, this leads to a radial separation between small and
large bubbles and to further coalescence of large bubbles migrating toward the
pipe center into even larger Taylor bubbles or slugs.An adequate modeling approach must consider all these phenomena. The
paper presents two different approaches both being based on the Eulerian
modeling framework. On one hand, a generalized inhomogeneous multiple size group
(MUSIG) model was applied, for which the dispersed gaseous phase is divided
into N inhomogeneous velocity groups (phases), each of these groups being
subdivided intoMj bubble
size classes. On the other hand, a moment density method is used to model the
evolution of the bubble size distribution function along the flow, where main
statistics of the distribution being described with a reduced set of transport
equations. For both models, bubble breakup and coalescence processes are taken
into account by appropriate models.


## 2. Model for the Local Size Distribution of a Bubble Population

### 2.1. Closure Models for Momentum Exchange and for Bubble Coalescence and Breakup

While modeling a two-phase flow using the Euler/Eulerian approach, the
momentum exchange between the phases has to be considered. Apart from the drag
acting in flow direction, the so-called nondrag forces acting mainly
perpendicularly to the flow direction must be considered. Namely, the lift
force, the turbulence dispersion force, the virtual mass force, and the wall
force play an important role.The turbulent dispersion force acts on smoothing the gas volume fraction
distribution and can be evaluated either from a single expression related to
drag turbulent contribution, [1] or deduced from a model for dispersed phase agitation contribution to the momentum balance, for example, the Tchen model [2]. To avoid the maximum gas fraction at the wall, Tomiyama et al. [3, 4] propose a wall force.The lift forceFL→ considers
the interaction of the bubble with the shear field of the liquid. For a single
bubble, it reads
(1)FL→=−CLρlπd36(wg→−wl→)×rot(wl→),
where ρl is the liquid
density, d the bubble diameter, and  wg→ and  wl→ are the bubble and liquid velocities,
respectively. The classical lift
force, which has a positive coefficient CL, acts in the direction of decreasing liquid velocity. In case of cocurrent upwards pipe flow, this is the direction towards the pipe wall. Numerical [5] and experimental [3] investigations showed that the direction of the lift force changes its sign if a substantial deformation of the bubble occurs. Tomiyama [4] investigated single bubble motion and derived the following correlation for the coefficient of the lift force from these experiments:
(2)CL={min[0.288tanh(0.121Re),f(Eod)],Eod<4f(Eod),for4<Eod<10−0.27Eod>10withf(Eod)=0.00105Eod3−0.0159Eod2−0.0204Eod+0.474.
This coefficient depends on the modified Eötvös number Eod, which is formulated for the maximum horizontal bubble size in
flow direction [6] and on the Reynolds number Re based on bubble size. The sign of the lift
coefficient, and, consequently, the direction of the lift force, depends, thus,
on the bubble size diameter. This behavior, originally found for single bubbles
of air in glycerol, was also established by different experiments for gas-water
polydispersed flows (e.g., [7]) and for different fluids too.For several flow configurations, this bubble size dependency of the lift
force direction can lead to the separation between small and large bubbles.
This effect has been shown to be a key phenomenon for the development of the
flow regime.The bubble coalescence model takes into account the random collision
processes between two bubbles. The applied model is based on the work of Prince
and Blanch [8]. The bubble breakup model considers the collision between a
liquid turbulent eddy of a certain characteristic size and a bubble (see [9]). Coalescence and breakup mechanistic models are common to both methods for taking into account the polydispersion in size of the bubble
population. We introduce corresponding “breakup”  and “coalescence” numerical
coefficients FB and FC used to scale the original
correlations.


### 2.2. Population Balance Approach

#### 2.2.1. The MUSIG Model by Lo

In principle, the Eulerian two-fluid approach can be extended to
simulate a continuous liquid phase and several gaseous dispersed phases solving
the complete set of balance equations for each phase. The investigations,
however, showed that for an adequate description of the gas volume fraction
profile including a population balance model decades of bubble size classes
would be necessary. In a CFD code, such a procedure is limited by the increased
computational effort to obtain converged flow solutions. To solve this problem,
the multiple size group model first implemented by the code developers in CFX-4
solves only one common momentum equation for all bubble size classes
(homogeneous MUSIG model, see [10], Figure 1). Mathematically, the multiple size group model (MUSIG) is based on the population balance method and the two-fluid modeling approach. The dispersed phase is divided into M size fractions. The population balance equation is applied to describe the mass conservation of the size fractions taking into account  the interfraction mass transfer resulting
from bubble coalescence and breakup. This model approach allows a sufficient
number of size fraction groups required for the coalescence and breakup
calculation to be used and has found a number of successful applications to
large-scale industrial multiphase flow problems.Figure 1
Scheme of the standard MUSIG model: all size fractions representing different bubble sizes move with the same velocity field.Nevertheless, the assumption also restricts its applicability to
homogeneous dispersed flows, where the slip velocities of particles are almost
independent of particle size and the particle relaxation time is sufficiently
small with respect to inertial time scales. Thus, the asymptotic slip velocity
can be considered to be attained almost instantaneously. The homogeneous MUSIG
model described above fails to predict the correct phase distribution when
heterogeneous particle motion becomes important. One example is the bubbly flow
in vertical pipes, where the nondrag forces play an essential role on the
bubble motion. The lift force was found to change its sign, when applied for
large deformed bubbles, which are dominated by the asymmetrical wake ([3], see Section 2.1). The lift force in this case has a direction opposite to the shear-induced lift force on a small bubble. For this reason, large bubbles tend to move to the pipe core region resulting in a core void maximum, whereas a near-wall void peak is measured for small bubbles. The
radial separation of small and large bubbles cannot be predicted by the
homogeneous MUSIG model. This has been shown to be a key mechanism for the
establishment of a certain flow regime.


#### 2.2.2. New Strategy: the Inhomogeneous MUSIG Model

A combination of the consideration of different dispersed phases and the
algebraic multiple size group model was proposed to combine both the adequate
number of bubble size classes for the simulation of coalescence and breakup and
a limited number of dispersed gaseous phases to limit the computational effort
[11]. The inhomogeneous MUSIG model was developed in
cooperation with ANSYS CFX and is implemented in CFX since the version CFX-10
([12–14] see Figure 2).Figure 2
Improvement of the MUSIG
approach: the size fractionsMj are
assigned to the velocity field Vj.In the inhomogeneous MUSIG model, the gaseous dispersed phase is divided
into a numberN of the so-called velocity groups (or phases), where each of the
velocity groups is characterized by its own velocity field. Further, the
overall bubble size distribution is represented by dividing the bubble diameter
range within each of the velocity groups j in a number Mj, j=1,…,N,  bubble
subsize fractions. The population balance model, considering bubble
coalescence or bubble breakup is applied to the subsize groups. Hence the mass
exchange between the subsize groups can exceed the size ranges assigned to the
velocity groups resulting in mass transfer terms between the different phases
or velocity groups.The lower and upper boundaries of bubble diameter intervals for the
bubble size fractions can be controlled by either an equal bubble diameter
distribution, an equal bubble mass distribution, or can be based on user
definition of the bubble diameter ranges for each distinct bubble diameter
fraction. The subdivision should be based on the physics of bubble motion for
bubbles of different size, for example, different behavior of differently sized
bubbles with respect to lift force or turbulent dispersion. Extensive model
validation calculations have shown that in most cases,N=2 or 3, velocity
groups are sufficient in order to capture the main phenomena in bubbly or slug
flows [15, 16].


### 2.3. The Moment Density Method Approach

The moment density method proposes an alternative way to model
polydispersion. This method allows to model the time and space evolution of a
realistic distribution with the help of a very reduced number of transport
equations, for example, Kamp et al. [17], Hill [18]. This is an interesting property, especially from a numerical point of view, with regard to alternative methods like population balance methods. In the present study, we restrict our attention to adiabatic bubbly flows of an incompressible gas inside a continuous liquid phase, and only consider a dispersion in bubbles size. We here present the basic formalism of the moment density method.


#### 2.3.1. A Distribution Function Parameterized by Its Statistical Moments

The moment density method requires the use of an approximate
representation of the bubble population thanks to a presumed shape continuous
function for the bubble size distribution functionf, the one being thus totally determined by a finite number of
parameters. The moment densities are intensive physical quantities
characterizing mean properties of the population, like the density of bubble
area or volume, the so-called volumetric interfacial area ai, and the volumetric fraction α of the dispersed phase. Bubbles being assumed as spherical and of diameter d, it yields
(3)ai=∫πd2f∂d,α=∫πd36f∂d.In this work, we introducef as being simply quadratic in the bubble size d,
(4)f={3n(x→,t)d4d13(x→,t)(2d1(x→,t)−d)ifd≤2d1,0elsewhere,
where n is the local density
number of bubbles (m−3) and d1 is the mean diameter of the distribution. It is easy to show that n, d1, ai, and α are related through:
(5)2d1ai=9α,243πnα2=10ai3.
The knowledge of the two moment densities ai and α is thus sufficient to reconstruct the whole
distribution function f. The choice
for these two parameters is dictated by their importance in the description of
the two-phase flow, especially concerning the transfers between the phases.


#### 2.3.2. Transport Equations for the Size Distribution

The transport of the distribution functionf obeys a Liouville-Boltzmann equation from which one can derive
transport equations for the main moment densities as well as for their
transport velocities [19]. The main idea of the method consists in
solving a finite number of moment densities transport equations in the
framework of the two-fluid model. This solving allows to determine the
distribution parameters and thus the whole population evolution. Thus we
consider the solving of the following transport equations:
(6)∂ai∂t+∇(ai〈w→〉ai)=∫2πd(d¯dd¯t)f∂d+∫πd2F˙c,b∂d,∂α∂t+∇⋅(α〈w→〉α)=∫πd22(d¯dd¯t)f∂d,
where F˙c,b is the
source term of the Liouville-Boltzmann equation. The transport velocities 〈w→〉ai and 〈w→〉α of ai and α,
respectively, are assumed to be equal to the mean transport velocity of the
bubbles w→g (no correlation between dispersions in
velocity and in size is taken into account at the current stage of development
of the model). 〈w→〉α is defined by
(7)α〈w→〉α=∫πd36w→f∂d,
where w→ is  the
velocity of a bubble of size d. It
satisfies the transport equation:
(8)∂α〈w→〉α∂t+∇⋅(α〈w→〉α〈w→〉α)=−∇⋅(α〈w→′w→′〉α)+∫πd36(d¯w→d¯t)f∂d.


#### 2.3.3. Closure Relations

The closure of the system of (3)–(8) concerns
models (i) for coalescence and breakup events through F˙c,b, (ii) for the bubble size evolution along its
trajectory d¯d/d¯t, (iii) for the hydrodynamic forces acting on individual
bubbles, through dw→/d¯t, as well as (iv) for the kinetic stress tensor ∇⋅(α〈w→′w→′〉α). The closure issue is based on the mechanistic
model of all these phenomena at the scale of a single bubble. Source terms of
the equations are then evaluated from the integral of individual contributions
over the whole bubble population, whose size distribution function is known.
The simplicity of the analytical expression (4) for f allows to easily derive expressions for the integrals from the
classical models described in Section 2.1. In this study, since the dispersed phase is
considered as incompressible and incondensable, the bubble size evolution along
its trajectory is zero. Detailed expressions for the source terms can be found
in [20]
.It is worth pointing out that solving the system of (3)–(8), coupled with
the system of balance equations for the continuous phase, is sufficient to
describe the evolution of a broad spectrum of bubble sizes. This simplicity of
the formulation is advantageous especially in terms of computational cost.
Further extension of the model, taking into account the dispersion in velocity,
can be considered using a similar formalism. For example, size-related velocity
distribution can be introduced to model size-dependent bubble migration due to
lift force, like it is done for the inhomogeneous MUSIG model. Moreover, this
method has already been successful to consider more general particles velocity
dispersion, for example, the work of Fox [21]
concerning crossing trajectories.


## 2.1. Closure Models for Momentum Exchange and for Bubble Coalescence and Breakup

While modeling a two-phase flow using the Euler/Eulerian approach, the
momentum exchange between the phases has to be considered. Apart from the drag
acting in flow direction, the so-called nondrag forces acting mainly
perpendicularly to the flow direction must be considered. Namely, the lift
force, the turbulence dispersion force, the virtual mass force, and the wall
force play an important role.The turbulent dispersion force acts on smoothing the gas volume fraction
distribution and can be evaluated either from a single expression related to
drag turbulent contribution, [1] or deduced from a model for dispersed phase agitation contribution to the momentum balance, for example, the Tchen model [2]. To avoid the maximum gas fraction at the wall, Tomiyama et al. [3, 4] propose a wall force.The lift forceFL→ considers
the interaction of the bubble with the shear field of the liquid. For a single
bubble, it reads
(1)FL→=−CLρlπd36(wg→−wl→)×rot(wl→),
where ρl is the liquid
density, d the bubble diameter, and  wg→ and  wl→ are the bubble and liquid velocities,
respectively. The classical lift
force, which has a positive coefficient CL, acts in the direction of decreasing liquid velocity. In case of cocurrent upwards pipe flow, this is the direction towards the pipe wall. Numerical [5] and experimental [3] investigations showed that the direction of the lift force changes its sign if a substantial deformation of the bubble occurs. Tomiyama [4] investigated single bubble motion and derived the following correlation for the coefficient of the lift force from these experiments:
(2)CL={min[0.288tanh(0.121Re),f(Eod)],Eod<4f(Eod),for4<Eod<10−0.27Eod>10withf(Eod)=0.00105Eod3−0.0159Eod2−0.0204Eod+0.474.
This coefficient depends on the modified Eötvös number Eod, which is formulated for the maximum horizontal bubble size in
flow direction [6] and on the Reynolds number Re based on bubble size. The sign of the lift
coefficient, and, consequently, the direction of the lift force, depends, thus,
on the bubble size diameter. This behavior, originally found for single bubbles
of air in glycerol, was also established by different experiments for gas-water
polydispersed flows (e.g., [7]) and for different fluids too.For several flow configurations, this bubble size dependency of the lift
force direction can lead to the separation between small and large bubbles.
This effect has been shown to be a key phenomenon for the development of the
flow regime.The bubble coalescence model takes into account the random collision
processes between two bubbles. The applied model is based on the work of Prince
and Blanch [8]. The bubble breakup model considers the collision between a
liquid turbulent eddy of a certain characteristic size and a bubble (see [9]). Coalescence and breakup mechanistic models are common to both methods for taking into account the polydispersion in size of the bubble
population. We introduce corresponding “breakup”  and “coalescence” numerical
coefficients FB and FC used to scale the original
correlations.


## 2.2. Population Balance Approach

### 2.2.1. The MUSIG Model by Lo

In principle, the Eulerian two-fluid approach can be extended to
simulate a continuous liquid phase and several gaseous dispersed phases solving
the complete set of balance equations for each phase. The investigations,
however, showed that for an adequate description of the gas volume fraction
profile including a population balance model decades of bubble size classes
would be necessary. In a CFD code, such a procedure is limited by the increased
computational effort to obtain converged flow solutions. To solve this problem,
the multiple size group model first implemented by the code developers in CFX-4
solves only one common momentum equation for all bubble size classes
(homogeneous MUSIG model, see [10], Figure 1). Mathematically, the multiple size group model (MUSIG) is based on the population balance method and the two-fluid modeling approach. The dispersed phase is divided into M size fractions. The population balance equation is applied to describe the mass conservation of the size fractions taking into account  the interfraction mass transfer resulting
from bubble coalescence and breakup. This model approach allows a sufficient
number of size fraction groups required for the coalescence and breakup
calculation to be used and has found a number of successful applications to
large-scale industrial multiphase flow problems.Figure 1
Scheme of the standard MUSIG model: all size fractions representing different bubble sizes move with the same velocity field.Nevertheless, the assumption also restricts its applicability to
homogeneous dispersed flows, where the slip velocities of particles are almost
independent of particle size and the particle relaxation time is sufficiently
small with respect to inertial time scales. Thus, the asymptotic slip velocity
can be considered to be attained almost instantaneously. The homogeneous MUSIG
model described above fails to predict the correct phase distribution when
heterogeneous particle motion becomes important. One example is the bubbly flow
in vertical pipes, where the nondrag forces play an essential role on the
bubble motion. The lift force was found to change its sign, when applied for
large deformed bubbles, which are dominated by the asymmetrical wake ([3], see Section 2.1). The lift force in this case has a direction opposite to the shear-induced lift force on a small bubble. For this reason, large bubbles tend to move to the pipe core region resulting in a core void maximum, whereas a near-wall void peak is measured for small bubbles. The
radial separation of small and large bubbles cannot be predicted by the
homogeneous MUSIG model. This has been shown to be a key mechanism for the
establishment of a certain flow regime.


### 2.2.2. New Strategy: the Inhomogeneous MUSIG Model

A combination of the consideration of different dispersed phases and the
algebraic multiple size group model was proposed to combine both the adequate
number of bubble size classes for the simulation of coalescence and breakup and
a limited number of dispersed gaseous phases to limit the computational effort
[11]. The inhomogeneous MUSIG model was developed in
cooperation with ANSYS CFX and is implemented in CFX since the version CFX-10
([12–14] see Figure 2).Figure 2
Improvement of the MUSIG
approach: the size fractionsMj are
assigned to the velocity field Vj.In the inhomogeneous MUSIG model, the gaseous dispersed phase is divided
into a numberN of the so-called velocity groups (or phases), where each of the
velocity groups is characterized by its own velocity field. Further, the
overall bubble size distribution is represented by dividing the bubble diameter
range within each of the velocity groups j in a number Mj, j=1,…,N,  bubble
subsize fractions. The population balance model, considering bubble
coalescence or bubble breakup is applied to the subsize groups. Hence the mass
exchange between the subsize groups can exceed the size ranges assigned to the
velocity groups resulting in mass transfer terms between the different phases
or velocity groups.The lower and upper boundaries of bubble diameter intervals for the
bubble size fractions can be controlled by either an equal bubble diameter
distribution, an equal bubble mass distribution, or can be based on user
definition of the bubble diameter ranges for each distinct bubble diameter
fraction. The subdivision should be based on the physics of bubble motion for
bubbles of different size, for example, different behavior of differently sized
bubbles with respect to lift force or turbulent dispersion. Extensive model
validation calculations have shown that in most cases,N=2 or 3, velocity
groups are sufficient in order to capture the main phenomena in bubbly or slug
flows [15, 16].


## 2.2.1. The MUSIG Model by Lo

In principle, the Eulerian two-fluid approach can be extended to
simulate a continuous liquid phase and several gaseous dispersed phases solving
the complete set of balance equations for each phase. The investigations,
however, showed that for an adequate description of the gas volume fraction
profile including a population balance model decades of bubble size classes
would be necessary. In a CFD code, such a procedure is limited by the increased
computational effort to obtain converged flow solutions. To solve this problem,
the multiple size group model first implemented by the code developers in CFX-4
solves only one common momentum equation for all bubble size classes
(homogeneous MUSIG model, see [10], Figure 1). Mathematically, the multiple size group model (MUSIG) is based on the population balance method and the two-fluid modeling approach. The dispersed phase is divided into M size fractions. The population balance equation is applied to describe the mass conservation of the size fractions taking into account  the interfraction mass transfer resulting
from bubble coalescence and breakup. This model approach allows a sufficient
number of size fraction groups required for the coalescence and breakup
calculation to be used and has found a number of successful applications to
large-scale industrial multiphase flow problems.Figure 1
Scheme of the standard MUSIG model: all size fractions representing different bubble sizes move with the same velocity field.Nevertheless, the assumption also restricts its applicability to
homogeneous dispersed flows, where the slip velocities of particles are almost
independent of particle size and the particle relaxation time is sufficiently
small with respect to inertial time scales. Thus, the asymptotic slip velocity
can be considered to be attained almost instantaneously. The homogeneous MUSIG
model described above fails to predict the correct phase distribution when
heterogeneous particle motion becomes important. One example is the bubbly flow
in vertical pipes, where the nondrag forces play an essential role on the
bubble motion. The lift force was found to change its sign, when applied for
large deformed bubbles, which are dominated by the asymmetrical wake ([3], see Section 2.1). The lift force in this case has a direction opposite to the shear-induced lift force on a small bubble. For this reason, large bubbles tend to move to the pipe core region resulting in a core void maximum, whereas a near-wall void peak is measured for small bubbles. The
radial separation of small and large bubbles cannot be predicted by the
homogeneous MUSIG model. This has been shown to be a key mechanism for the
establishment of a certain flow regime.


## 2.2.2. New Strategy: the Inhomogeneous MUSIG Model

A combination of the consideration of different dispersed phases and the
algebraic multiple size group model was proposed to combine both the adequate
number of bubble size classes for the simulation of coalescence and breakup and
a limited number of dispersed gaseous phases to limit the computational effort
[11]. The inhomogeneous MUSIG model was developed in
cooperation with ANSYS CFX and is implemented in CFX since the version CFX-10
([12–14] see Figure 2).Figure 2
Improvement of the MUSIG
approach: the size fractionsMj are
assigned to the velocity field Vj.In the inhomogeneous MUSIG model, the gaseous dispersed phase is divided
into a numberN of the so-called velocity groups (or phases), where each of the
velocity groups is characterized by its own velocity field. Further, the
overall bubble size distribution is represented by dividing the bubble diameter
range within each of the velocity groups j in a number Mj, j=1,…,N,  bubble
subsize fractions. The population balance model, considering bubble
coalescence or bubble breakup is applied to the subsize groups. Hence the mass
exchange between the subsize groups can exceed the size ranges assigned to the
velocity groups resulting in mass transfer terms between the different phases
or velocity groups.The lower and upper boundaries of bubble diameter intervals for the
bubble size fractions can be controlled by either an equal bubble diameter
distribution, an equal bubble mass distribution, or can be based on user
definition of the bubble diameter ranges for each distinct bubble diameter
fraction. The subdivision should be based on the physics of bubble motion for
bubbles of different size, for example, different behavior of differently sized
bubbles with respect to lift force or turbulent dispersion. Extensive model
validation calculations have shown that in most cases,N=2 or 3, velocity
groups are sufficient in order to capture the main phenomena in bubbly or slug
flows [15, 16].


## 2.3. The Moment Density Method Approach

The moment density method proposes an alternative way to model
polydispersion. This method allows to model the time and space evolution of a
realistic distribution with the help of a very reduced number of transport
equations, for example, Kamp et al. [17], Hill [18]. This is an interesting property, especially from a numerical point of view, with regard to alternative methods like population balance methods. In the present study, we restrict our attention to adiabatic bubbly flows of an incompressible gas inside a continuous liquid phase, and only consider a dispersion in bubbles size. We here present the basic formalism of the moment density method.


### 2.3.1. A Distribution Function Parameterized by Its Statistical Moments

The moment density method requires the use of an approximate
representation of the bubble population thanks to a presumed shape continuous
function for the bubble size distribution functionf, the one being thus totally determined by a finite number of
parameters. The moment densities are intensive physical quantities
characterizing mean properties of the population, like the density of bubble
area or volume, the so-called volumetric interfacial area ai, and the volumetric fraction α of the dispersed phase. Bubbles being assumed as spherical and of diameter d, it yields
(3)ai=∫πd2f∂d,α=∫πd36f∂d.In this work, we introducef as being simply quadratic in the bubble size d,
(4)f={3n(x→,t)d4d13(x→,t)(2d1(x→,t)−d)ifd≤2d1,0elsewhere,
where n is the local density
number of bubbles (m−3) and d1 is the mean diameter of the distribution. It is easy to show that n, d1, ai, and α are related through:
(5)2d1ai=9α,243πnα2=10ai3.
The knowledge of the two moment densities ai and α is thus sufficient to reconstruct the whole
distribution function f. The choice
for these two parameters is dictated by their importance in the description of
the two-phase flow, especially concerning the transfers between the phases.


### 2.3.2. Transport Equations for the Size Distribution

The transport of the distribution functionf obeys a Liouville-Boltzmann equation from which one can derive
transport equations for the main moment densities as well as for their
transport velocities [19]. The main idea of the method consists in
solving a finite number of moment densities transport equations in the
framework of the two-fluid model. This solving allows to determine the
distribution parameters and thus the whole population evolution. Thus we
consider the solving of the following transport equations:
(6)∂ai∂t+∇(ai〈w→〉ai)=∫2πd(d¯dd¯t)f∂d+∫πd2F˙c,b∂d,∂α∂t+∇⋅(α〈w→〉α)=∫πd22(d¯dd¯t)f∂d,
where F˙c,b is the
source term of the Liouville-Boltzmann equation. The transport velocities 〈w→〉ai and 〈w→〉α of ai and α,
respectively, are assumed to be equal to the mean transport velocity of the
bubbles w→g (no correlation between dispersions in
velocity and in size is taken into account at the current stage of development
of the model). 〈w→〉α is defined by
(7)α〈w→〉α=∫πd36w→f∂d,
where w→ is  the
velocity of a bubble of size d. It
satisfies the transport equation:
(8)∂α〈w→〉α∂t+∇⋅(α〈w→〉α〈w→〉α)=−∇⋅(α〈w→′w→′〉α)+∫πd36(d¯w→d¯t)f∂d.


### 2.3.3. Closure Relations

The closure of the system of (3)–(8) concerns
models (i) for coalescence and breakup events through F˙c,b, (ii) for the bubble size evolution along its
trajectory d¯d/d¯t, (iii) for the hydrodynamic forces acting on individual
bubbles, through dw→/d¯t, as well as (iv) for the kinetic stress tensor ∇⋅(α〈w→′w→′〉α). The closure issue is based on the mechanistic
model of all these phenomena at the scale of a single bubble. Source terms of
the equations are then evaluated from the integral of individual contributions
over the whole bubble population, whose size distribution function is known.
The simplicity of the analytical expression (4) for f allows to easily derive expressions for the integrals from the
classical models described in Section 2.1. In this study, since the dispersed phase is
considered as incompressible and incondensable, the bubble size evolution along
its trajectory is zero. Detailed expressions for the source terms can be found
in [20]
.It is worth pointing out that solving the system of (3)–(8), coupled with
the system of balance equations for the continuous phase, is sufficient to
describe the evolution of a broad spectrum of bubble sizes. This simplicity of
the formulation is advantageous especially in terms of computational cost.
Further extension of the model, taking into account the dispersion in velocity,
can be considered using a similar formalism. For example, size-related velocity
distribution can be introduced to model size-dependent bubble migration due to
lift force, like it is done for the inhomogeneous MUSIG model. Moreover, this
method has already been successful to consider more general particles velocity
dispersion, for example, the work of Fox [21]
concerning crossing trajectories.


## 2.3.1. A Distribution Function Parameterized by Its Statistical Moments

The moment density method requires the use of an approximate
representation of the bubble population thanks to a presumed shape continuous
function for the bubble size distribution functionf, the one being thus totally determined by a finite number of
parameters. The moment densities are intensive physical quantities
characterizing mean properties of the population, like the density of bubble
area or volume, the so-called volumetric interfacial area ai, and the volumetric fraction α of the dispersed phase. Bubbles being assumed as spherical and of diameter d, it yields
(3)ai=∫πd2f∂d,α=∫πd36f∂d.In this work, we introducef as being simply quadratic in the bubble size d,
(4)f={3n(x→,t)d4d13(x→,t)(2d1(x→,t)−d)ifd≤2d1,0elsewhere,
where n is the local density
number of bubbles (m−3) and d1 is the mean diameter of the distribution. It is easy to show that n, d1, ai, and α are related through:
(5)2d1ai=9α,243πnα2=10ai3.
The knowledge of the two moment densities ai and α is thus sufficient to reconstruct the whole
distribution function f. The choice
for these two parameters is dictated by their importance in the description of
the two-phase flow, especially concerning the transfers between the phases.


## 2.3.2. Transport Equations for the Size Distribution

The transport of the distribution functionf obeys a Liouville-Boltzmann equation from which one can derive
transport equations for the main moment densities as well as for their
transport velocities [19]. The main idea of the method consists in
solving a finite number of moment densities transport equations in the
framework of the two-fluid model. This solving allows to determine the
distribution parameters and thus the whole population evolution. Thus we
consider the solving of the following transport equations:
(6)∂ai∂t+∇(ai〈w→〉ai)=∫2πd(d¯dd¯t)f∂d+∫πd2F˙c,b∂d,∂α∂t+∇⋅(α〈w→〉α)=∫πd22(d¯dd¯t)f∂d,
where F˙c,b is the
source term of the Liouville-Boltzmann equation. The transport velocities 〈w→〉ai and 〈w→〉α of ai and α,
respectively, are assumed to be equal to the mean transport velocity of the
bubbles w→g (no correlation between dispersions in
velocity and in size is taken into account at the current stage of development
of the model). 〈w→〉α is defined by
(7)α〈w→〉α=∫πd36w→f∂d,
where w→ is  the
velocity of a bubble of size d. It
satisfies the transport equation:
(8)∂α〈w→〉α∂t+∇⋅(α〈w→〉α〈w→〉α)=−∇⋅(α〈w→′w→′〉α)+∫πd36(d¯w→d¯t)f∂d.


## 2.3.3. Closure Relations

The closure of the system of (3)–(8) concerns
models (i) for coalescence and breakup events through F˙c,b, (ii) for the bubble size evolution along its
trajectory d¯d/d¯t, (iii) for the hydrodynamic forces acting on individual
bubbles, through dw→/d¯t, as well as (iv) for the kinetic stress tensor ∇⋅(α〈w→′w→′〉α). The closure issue is based on the mechanistic
model of all these phenomena at the scale of a single bubble. Source terms of
the equations are then evaluated from the integral of individual contributions
over the whole bubble population, whose size distribution function is known.
The simplicity of the analytical expression (4) for f allows to easily derive expressions for the integrals from the
classical models described in Section 2.1. In this study, since the dispersed phase is
considered as incompressible and incondensable, the bubble size evolution along
its trajectory is zero. Detailed expressions for the source terms can be found
in [20]
.It is worth pointing out that solving the system of (3)–(8), coupled with
the system of balance equations for the continuous phase, is sufficient to
describe the evolution of a broad spectrum of bubble sizes. This simplicity of
the formulation is advantageous especially in terms of computational cost.
Further extension of the model, taking into account the dispersion in velocity,
can be considered using a similar formalism. For example, size-related velocity
distribution can be introduced to model size-dependent bubble migration due to
lift force, like it is done for the inhomogeneous MUSIG model. Moreover, this
method has already been successful to consider more general particles velocity
dispersion, for example, the work of Fox [21]
concerning crossing trajectories.


## 3. The Experiment–Bubbly Flow Around an Obstacle

In the presented experiment performed at the TOPFLOW facility of
Research center Dresden-Rossendorf (FZD), the large test section with a nominal
diameter of DN200 was used to study the flow field around an asymmetric
obstacle (see Figure3). This is an ideal test case for CFD code validation, since the obstacle creates a pronounced three-dimensional two-phase flow field. Curved stream lines, which form significant angles with the gravity vector, a recirculation zone in the wake and a flow separation at the edge of the obstacle are common in industrial components and installations.Figure 3
Sketch of the movable
obstacle with driving mechanism—a half-moon shaped horizontal plate mounted
on top of a toothed rod.The wire-mesh technology was applied to measure the gas volume fraction
and the gas velocity in different distances up- and downstream the obstacle
[22]. The sensor provides detailed data on the instantaneous flow structure with a high resolution in space and time. In particular, they allow visualizing the structure of the gas-liquid interface [23].The tests were performed both for air/water and for steam/water. In the
current paper, only adiabatic air/water tests were considered. The parameters
are summarized in Table1.Table 1
Water and gas superficial
velocitiesJL  and JG of the tests investigated in this paper.
JG[m/s]—0.03680.0898JL[m/s]1.611—Run n° 0971.017Run n° 074Run n° 096


## 4. Numerical Settings

### 4.1. CFX and Population Balance Method

Pretest calculations using ANSYS/CFX and applying a monodispersed bubble
size approach were performed for the conditions of test run 074 (JL=1.017 m/s, JG=0.0368 m/s) (see [24–26]). In the calculation, a fluid domain was
modeled 1.5 m upstream and downstream the obstacle. Half of the tube including
a symmetry boundary condition set at the xz-plane of the geometry was
simulated. In the present paper, the inhomogeneous MUSIG model approach was
applied to air/water obstacle experiments run 096 (JL=1.017 m/s, JG=0.0898 m/s) and run 097 (JL=1.611 m/s, JG=0.0898 m/s). In the presented calculations for run 096 and run 097, 25, and 20,
respectively, subsize gas fractions representing equidistant bubble sizes up to
25 mm and 20 mm, respectively, were simulated, assigned to 2 dispersed gaseous phases. The first 6 size groups were assigned to the first gaseous phase (or
velocity group) and the remaining size groups were assigned to the second
gaseous phase. The bubble size distribution measured at the largest upstream
position was set as an inlet boundary condition for the calculation.


### 4.2. Neptune_CFD and Moment Density Method

Neptune_CFD code [27], which is based on the two-fluid model formulation, has been used to perform numerical simulations using the moment density method. The test run 074 has been simulated using a calculation domain corresponding to a 1 m long half pipe centered on the obstacle of around 150 000 cells that takes benefit of the symmetry of the experimental setup. Neptune_CFD calculations of test 074 are based on a k-ε model for the liquid flow, the Tchen correlation being used to
model the fluctuations related to bubbles turbulent motion. Bubble
hydrodynamics model include both laminar and turbulent evaluations of drag,
lift, and virtual mass forces, thanks to a drift model for the relative
velocity. Inlet boundary conditions correspond to flat profiles for α,ai,  and the velocities,
allowing to recover the experimental surface averaged values for superficial
velocities and mean bubble diameter. The surface averaged volumetric fraction
and mean diameter are used to reconstruct the inlet bubble size distribution
function. The numerical method is based on a pressure-based method. Mass,
momentum, and energy equations are coupled by an iterative procedure within a
time step. Spatial discretization is based on finite volume framework on
unstructured meshes, all the variables being computed at the center of cells.


## 4.1. CFX and Population Balance Method

Pretest calculations using ANSYS/CFX and applying a monodispersed bubble
size approach were performed for the conditions of test run 074 (JL=1.017 m/s, JG=0.0368 m/s) (see [24–26]). In the calculation, a fluid domain was
modeled 1.5 m upstream and downstream the obstacle. Half of the tube including
a symmetry boundary condition set at the xz-plane of the geometry was
simulated. In the present paper, the inhomogeneous MUSIG model approach was
applied to air/water obstacle experiments run 096 (JL=1.017 m/s, JG=0.0898 m/s) and run 097 (JL=1.611 m/s, JG=0.0898 m/s). In the presented calculations for run 096 and run 097, 25, and 20,
respectively, subsize gas fractions representing equidistant bubble sizes up to
25 mm and 20 mm, respectively, were simulated, assigned to 2 dispersed gaseous phases. The first 6 size groups were assigned to the first gaseous phase (or
velocity group) and the remaining size groups were assigned to the second
gaseous phase. The bubble size distribution measured at the largest upstream
position was set as an inlet boundary condition for the calculation.


## 4.2. Neptune_CFD and Moment Density Method

Neptune_CFD code [27], which is based on the two-fluid model formulation, has been used to perform numerical simulations using the moment density method. The test run 074 has been simulated using a calculation domain corresponding to a 1 m long half pipe centered on the obstacle of around 150 000 cells that takes benefit of the symmetry of the experimental setup. Neptune_CFD calculations of test 074 are based on a k-ε model for the liquid flow, the Tchen correlation being used to
model the fluctuations related to bubbles turbulent motion. Bubble
hydrodynamics model include both laminar and turbulent evaluations of drag,
lift, and virtual mass forces, thanks to a drift model for the relative
velocity. Inlet boundary conditions correspond to flat profiles for α,ai,  and the velocities,
allowing to recover the experimental surface averaged values for superficial
velocities and mean bubble diameter. The surface averaged volumetric fraction
and mean diameter are used to reconstruct the inlet bubble size distribution
function. The numerical method is based on a pressure-based method. Mass,
momentum, and energy equations are coupled by an iterative procedure within a
time step. Spatial discretization is based on finite volume framework on
unstructured meshes, all the variables being computed at the center of cells.


## 5. Comparison of Measured and Calculated Results

### 5.1. The Main Observed Phenomena

Both the steady-state ANSYS CFX calculations applying the inhomogeneous
MUSIG model, and the Neptune_CFD calculations applying the moment density
method, could reproduce all qualitative details of the flow structure of the
two-phase flow field around the diaphragm. The structure of the flow for the
here considered test cases 074, 096, and 097 are essentially similar. These
different tests have been selected for the purpose of investigating certain
phenomena which are more or less pronounced.The numerical results have been compared to three-dimensional wire-mesh
sensor data in Figure4 (ANSYS CFX for the run 096) and Figure 5 (Neptune_CFD for the run 074). The water velocity and the total gaseous void fraction are represented. All qualitative details of the structure of the two-phase flow field around the obstacle could be reproduced.Comparison of time
averaged values  calculated by CFX
(left) and measured (right) up- and downstream of the obstacle in the
air-water test run 096,JL=1.017 m/s, JG=0.0898 m/s, (FB=FC=0.05).
(a)(b)Figure 5
Neptune_CFD numerical
(left-hand side of each pair) versus TOPFLOW experimental (right-hand side of
each pair) results for run 074. The two left-hand side pictures represent the
dispersed phase volume fraction (scale from 0 to6⋅10−2)  across the symmetry plane of the pipe and
the two right-hand side pictures the liquid velocity norm (scale from 0 to 2 ms−1).Shortly, behind the obstacle a strong vortex of the liquid combined with
the accumulation of gas is observed. The measured and calculated shape and
extension of the recirculation area agree very well. Upstream the obstacle, a
stagnation point with lower gas content is seen in experiment and calculation.
Details, like the velocity and void fraction maxima above the gap between the
circular edge of the obstacle and the inner wall of the pipe, are also found in
a good agreement between experiments and calculations. In the unobstructed
cross sectional part of the tube a strong jet is established. Main
discrepancies between experimental and Neptune_CFD results concern the
volumetric fraction upward (below) the obstacle and can be associated to the
flat profiles used as inlet bottom boundary conditions.The structure of the flow is studied in more detail in the following
sections.


### 5.2. Phenomena in the Wake of the Obstacle

#### 5.2.1. Size Distribution of the Bubble Population

More detailed understanding of the flow situation can be gained,
considering the bubble size distribution. According to the applied bubble
breakup model of Luo and Svendsen [9], bubble breakup can be expected in
regions showing high-turbulent eddy dissipation. Figure 6 presents maximum
values of the numerically evaluated turbulent eddy dissipation at the edges of
the obstacle. At the same time, the applied bubble coalescence model of Prince
and Blanch [8] indicates strong importance of coalescence in regions of
bubble accumulation, that is, in the wake behind the obstacle. Both bubble
coalescence (see gas accumulation shown in Figures 4 and 5) and bubble breakup
(see distribution of turbulence dissipation Figure 6), which might partially compensate each other, are expected shortly behind the obstacle.Figure 6
Turbulence eddy dissipation
(run 096) (CFX).Figure7 presents measured cross-sectional averaged bubble size
distributions upstream (z=−0.52 m), shortly behind (z=0.08 m) and
downstream the obstacle (z=0.52 m) for the run 096. The measurements show in
the bubble accumulation zone at z=0.08 m the cross-sectional average shows a
shift toward larger bubbles. The experimentally measured bubble size distribution of run 074 is essentially similar to the run 096. It can be seen on the volumetric fraction maps, Figure 8. Above the obstacle, where turbulent intensity is strong, breakup is not experimentally attested, contrarily to the Luo and Svedsen [9] model prediction, but coalescence occurs, leading to the formation of big bubbles (over 7 mm) that concentrate thereafter in the central part of the pipe.Figure 7
Measured bubble size distribution for run 096.TOPFLOW experimental
volumetric fractions of air bubbles according to different size classes,
extracted from [24].
(a)(b)(c)(d)Both ANSYS CFX and Neptune_CFD calculated bubble size distributions, however, show a shift of the mean bubble diameter toward smaller bubbles
shortly behind the obstacle when both coalescence and breakup are taken into
account. In the calculations, the bubble breakup is overestimated. The
corresponding results are provided by Figure9 for the run 096 and on Figure 10 for run 097. This disagreement was found not solvable by simple changes of breakup or coalescence coefficients, which were set here to FB=FC=0.05. Similar deviations would arise at other
locations of the flow domain. Neptune_CFD results for the run 074 with standard
coefficients FB=FC=1 are shown in Figure 11
(the details can be found in [20]).Figure 9
Bubble size distributions for run 096 (JL=1.017 m/s, JG=0.0898 m/s) (CFX) (FB=FC=0.05).Figure 10
Bubble size distributions
for run 097 (JL=1.611 m/s, JG=0.0898 m/s) (CFX)
(FB=FC=0.05).Figure 11
Neptune_CFD result for
mean Sauter bubble diameter (m), both breakup and coalescence being taken
into account (FB=FC=1) run 074.This suggests to perform computations for which breakup is neglected, but
that still consider the Prince and Blanch [8] coalescence model with standard coefficient. The corresponding Neptune_CFD results show that both bubble size repartition and order of magnitude are consistent with experimental results: big bubbles of more than 7 mm are created and accumulate in the recirculation zone just above the obstacle (see Figure 14). Big bubbles stay in the LHS of the pipe, the unobstructed cross-sectional part of the pipe is free of big bubbles, that is consistent with the experimental observation.As partial conclusions, (i) the Luo and Svendsen breakup model tends to
overestimate the breakup of TOPFLOW experimental tests, whereas (ii) the use of
the Prince and Blanch model for coalescence within the framework of the moment
density method allows to predict satisfying evolution of bubble size across the
flow.


#### 5.2.2.
Bubbles Streamlines

More detailed effects of lateral motion of small and large bubbles can
be revealed by studying bubble streamlines and by analyzing lift forces acting
on bubbles of different size. On one hand, the liquid velocity flow carries the
small bubbles into the region behind the obstacle (see Figures12 and 14, right-hand side for the bubble streamlines). Lateral deviation due to lift
force is illustrated on Figure 13 for the lift force arrows with the population
balance method and Figure 14, right-hand side for the lift coefficient with the moment density method. On the other hand, the air accumulation in the wake region leads to bubble coalescence and the generation of large bubbles as revealed by the analysis of experimental results. This phenomenon is
underestimated in the calculations taking breakup into account.Figure 12
Streamlines for small
(left) and large (right) bubbles (run 096) (CFX).Figure 13
Bubble lift force vectors
for the different gas velocity groups (run 096) (CFX).Neptune_CFD numerical
results for local mean bubble diameter and corresponding streamlines colored
by the equivalent mean lift coefficient (run 074)(FB=0;FC=1).
(a)(b)Caused by the lift force, large bubbles are redirected into the
downstream jet (see Figure13) once they can be formed in the wake by
coalescence. The streamline representation (see Figure 12) clearly shows
this phenomenon for large bubbles already present in the upstream flow.In the actual version of the moment density method used in Neptune_CFD,
dynamics of the bubble population is estimated using a single transport
velocity. The averaged lift contribution takes into account both the bubble
size distribution and the Tomiyama correlation. When a majority of bubbles are
locally above the critical Eötvös number, the lift coefficient changes its
sign. In this case, the direction of the lift force is changed, as it can be
seen on Figure14. As a consequence, the bubble streamlines above the obstacle
deviate to the center of the pipe. This explains the low value of the dispersed
phase volumetric fraction near the left-hand side part of the pipe wall above
the obstacle for Neptune_CFD numerical results (see Figure 5). This result is fully consistent with experimental results.As a partial conclusion concerning the bubble streamlines calculations,
both methods showed their ability to consider the effect of bubble size on the
lateral deviation of bubble streamlines due to lift force. This provides a more
precise understanding and a more accurate prediction of bubbles repartition
across the flow.


### 5.3.
Phenomena in the Jet

In the cross-sectional area beside the obstacle, a strong jet is
established creating strong shear flow. The resulting phenomena are more
pronounced with increasing water velocity, like in  run 097, where the liquid velocity was increased toJL=1.611 m/s. Figure 16 represents measured and ANSYS CFX calculated cross-sectional gas fraction distributions for this run. In the most downstream cross section of the measurements an almost gas bubble free region in the centre of the jet is found. Bubbles are collected at the edges of the jet in regions of largest water velocity gradient.The streamline representation of the ANSYS CFX calculations, however,
(Figure12 for run 096, which is fairly similar to run 097) indicates large bubbles being directed into the jet caused by the lift force. This discrepancy between experiment and ANSYS CFX calculations can possibly be explained by the strong water velocity gradient near the jet. This strong shear flow induces
bubble breakup which is not yet considered in the model of Luo and Svendsen
[9]. In the tests, the big bubbles could migrate toward the jet, but be fragmented at the relatively sharp boarder of this jet. Only a small fraction of the small bubbles created by this breakup process can enter the jet by action of the turbulent dispersion force.The gas distribution resolved by bubble size classes (see Figure15) shows clearly the deviations between measurements and calculations. In the experiment, the gas accumulation behind the obstacle leads to a strong coalescence and the creation of large bubbles. In the CFX-calculations, this
effect is underestimated.Calculated by CFX (left)
and measured (right) gas distributions up- and downstream of the obstacle
resolved to bubble size classes (run 096JL=1.017 m/s, JG=0.0898 m/s, FB=0,FC=0.05).
(a)(b)Figure 16
Calculated by CFX (left)
and measured (right) gas cross-fractional distributions downstream the
obstacle (run 097JL=1.611 m/s, JG=0.0898 m/s, FB=FC=0.05).
Calculations (obstacle shown), distances at
z=0.08 m, 0.16 m, 0.25 m, 0.37 m, and 0.52 m. Measurements (obstacle in the upper left area). distances at  z=0.01 m,
0.015 m, 0.02 m, 0.04 m, 0.08 m, 0.16 m, 0.25 m, and 0.52 m.In Neptune_CFD calculations based on the moment density method the
bubble breakup was neglected. As attested by Figure14, left-hand side, in this
case the bubbles present in the jet are small bubbles. Big bubbles created are
concentrated in the center of the pipe where the largest value of volumetric
fraction can be found (Figure 5). The structure of the gas repartition in the
jet for both calculation and experiments is represented on Figure 17. In the
upper part, the structure of the flow of the numerical calculation is
consistent with experimental results. The most important deviation is observed
in the region just above the obstacle (z=0.08m). A gain in accuracy of these
results could be achieved by considering size-dependent velocity dispersion
within the moment density method, in particular for bubble size-dependent
lateral migration in the regions where there exists a strong mixing of
different bubble sizes (mainly just above the obstacle as far as run 074 is
concerned, see Figure 8).Figure 17
Calculated by
Neptune_CFD (left) and measured (right) gas volumetric fraction run 074 for different elevations (obstacle, symbolized in gray, is atz=0 m and covers the “upper” part of the half-pipe
cross-section representation).Nevertheless, the main trends of the bubble dynamics downward the
obstacle are in agreement with experimental results (see Figure16 for CFX and Figure 17 for Neptune_CFD).


## 5.1. The Main Observed Phenomena

Both the steady-state ANSYS CFX calculations applying the inhomogeneous
MUSIG model, and the Neptune_CFD calculations applying the moment density
method, could reproduce all qualitative details of the flow structure of the
two-phase flow field around the diaphragm. The structure of the flow for the
here considered test cases 074, 096, and 097 are essentially similar. These
different tests have been selected for the purpose of investigating certain
phenomena which are more or less pronounced.The numerical results have been compared to three-dimensional wire-mesh
sensor data in Figure4 (ANSYS CFX for the run 096) and Figure 5 (Neptune_CFD for the run 074). The water velocity and the total gaseous void fraction are represented. All qualitative details of the structure of the two-phase flow field around the obstacle could be reproduced.Comparison of time
averaged values  calculated by CFX
(left) and measured (right) up- and downstream of the obstacle in the
air-water test run 096,JL=1.017 m/s, JG=0.0898 m/s, (FB=FC=0.05).
(a)(b)Figure 5
Neptune_CFD numerical
(left-hand side of each pair) versus TOPFLOW experimental (right-hand side of
each pair) results for run 074. The two left-hand side pictures represent the
dispersed phase volume fraction (scale from 0 to6⋅10−2)  across the symmetry plane of the pipe and
the two right-hand side pictures the liquid velocity norm (scale from 0 to 2 ms−1).Shortly, behind the obstacle a strong vortex of the liquid combined with
the accumulation of gas is observed. The measured and calculated shape and
extension of the recirculation area agree very well. Upstream the obstacle, a
stagnation point with lower gas content is seen in experiment and calculation.
Details, like the velocity and void fraction maxima above the gap between the
circular edge of the obstacle and the inner wall of the pipe, are also found in
a good agreement between experiments and calculations. In the unobstructed
cross sectional part of the tube a strong jet is established. Main
discrepancies between experimental and Neptune_CFD results concern the
volumetric fraction upward (below) the obstacle and can be associated to the
flat profiles used as inlet bottom boundary conditions.The structure of the flow is studied in more detail in the following
sections.


## 5.2. Phenomena in the Wake of the Obstacle

### 5.2.1. Size Distribution of the Bubble Population

More detailed understanding of the flow situation can be gained,
considering the bubble size distribution. According to the applied bubble
breakup model of Luo and Svendsen [9], bubble breakup can be expected in
regions showing high-turbulent eddy dissipation. Figure 6 presents maximum
values of the numerically evaluated turbulent eddy dissipation at the edges of
the obstacle. At the same time, the applied bubble coalescence model of Prince
and Blanch [8] indicates strong importance of coalescence in regions of
bubble accumulation, that is, in the wake behind the obstacle. Both bubble
coalescence (see gas accumulation shown in Figures 4 and 5) and bubble breakup
(see distribution of turbulence dissipation Figure 6), which might partially compensate each other, are expected shortly behind the obstacle.Figure 6
Turbulence eddy dissipation
(run 096) (CFX).Figure7 presents measured cross-sectional averaged bubble size
distributions upstream (z=−0.52 m), shortly behind (z=0.08 m) and
downstream the obstacle (z=0.52 m) for the run 096. The measurements show in
the bubble accumulation zone at z=0.08 m the cross-sectional average shows a
shift toward larger bubbles. The experimentally measured bubble size distribution of run 074 is essentially similar to the run 096. It can be seen on the volumetric fraction maps, Figure 8. Above the obstacle, where turbulent intensity is strong, breakup is not experimentally attested, contrarily to the Luo and Svedsen [9] model prediction, but coalescence occurs, leading to the formation of big bubbles (over 7 mm) that concentrate thereafter in the central part of the pipe.Figure 7
Measured bubble size distribution for run 096.TOPFLOW experimental
volumetric fractions of air bubbles according to different size classes,
extracted from [24].
(a)(b)(c)(d)Both ANSYS CFX and Neptune_CFD calculated bubble size distributions, however, show a shift of the mean bubble diameter toward smaller bubbles
shortly behind the obstacle when both coalescence and breakup are taken into
account. In the calculations, the bubble breakup is overestimated. The
corresponding results are provided by Figure9 for the run 096 and on Figure 10 for run 097. This disagreement was found not solvable by simple changes of breakup or coalescence coefficients, which were set here to FB=FC=0.05. Similar deviations would arise at other
locations of the flow domain. Neptune_CFD results for the run 074 with standard
coefficients FB=FC=1 are shown in Figure 11
(the details can be found in [20]).Figure 9
Bubble size distributions for run 096 (JL=1.017 m/s, JG=0.0898 m/s) (CFX) (FB=FC=0.05).Figure 10
Bubble size distributions
for run 097 (JL=1.611 m/s, JG=0.0898 m/s) (CFX)
(FB=FC=0.05).Figure 11
Neptune_CFD result for
mean Sauter bubble diameter (m), both breakup and coalescence being taken
into account (FB=FC=1) run 074.This suggests to perform computations for which breakup is neglected, but
that still consider the Prince and Blanch [8] coalescence model with standard coefficient. The corresponding Neptune_CFD results show that both bubble size repartition and order of magnitude are consistent with experimental results: big bubbles of more than 7 mm are created and accumulate in the recirculation zone just above the obstacle (see Figure 14). Big bubbles stay in the LHS of the pipe, the unobstructed cross-sectional part of the pipe is free of big bubbles, that is consistent with the experimental observation.As partial conclusions, (i) the Luo and Svendsen breakup model tends to
overestimate the breakup of TOPFLOW experimental tests, whereas (ii) the use of
the Prince and Blanch model for coalescence within the framework of the moment
density method allows to predict satisfying evolution of bubble size across the
flow.


### 5.2.2.
Bubbles Streamlines

More detailed effects of lateral motion of small and large bubbles can
be revealed by studying bubble streamlines and by analyzing lift forces acting
on bubbles of different size. On one hand, the liquid velocity flow carries the
small bubbles into the region behind the obstacle (see Figures12 and 14, right-hand side for the bubble streamlines). Lateral deviation due to lift
force is illustrated on Figure 13 for the lift force arrows with the population
balance method and Figure 14, right-hand side for the lift coefficient with the moment density method. On the other hand, the air accumulation in the wake region leads to bubble coalescence and the generation of large bubbles as revealed by the analysis of experimental results. This phenomenon is
underestimated in the calculations taking breakup into account.Figure 12
Streamlines for small
(left) and large (right) bubbles (run 096) (CFX).Figure 13
Bubble lift force vectors
for the different gas velocity groups (run 096) (CFX).Neptune_CFD numerical
results for local mean bubble diameter and corresponding streamlines colored
by the equivalent mean lift coefficient (run 074)(FB=0;FC=1).
(a)(b)Caused by the lift force, large bubbles are redirected into the
downstream jet (see Figure13) once they can be formed in the wake by
coalescence. The streamline representation (see Figure 12) clearly shows
this phenomenon for large bubbles already present in the upstream flow.In the actual version of the moment density method used in Neptune_CFD,
dynamics of the bubble population is estimated using a single transport
velocity. The averaged lift contribution takes into account both the bubble
size distribution and the Tomiyama correlation. When a majority of bubbles are
locally above the critical Eötvös number, the lift coefficient changes its
sign. In this case, the direction of the lift force is changed, as it can be
seen on Figure14. As a consequence, the bubble streamlines above the obstacle
deviate to the center of the pipe. This explains the low value of the dispersed
phase volumetric fraction near the left-hand side part of the pipe wall above
the obstacle for Neptune_CFD numerical results (see Figure 5). This result is fully consistent with experimental results.As a partial conclusion concerning the bubble streamlines calculations,
both methods showed their ability to consider the effect of bubble size on the
lateral deviation of bubble streamlines due to lift force. This provides a more
precise understanding and a more accurate prediction of bubbles repartition
across the flow.


## 5.2.1. Size Distribution of the Bubble Population

More detailed understanding of the flow situation can be gained,
considering the bubble size distribution. According to the applied bubble
breakup model of Luo and Svendsen [9], bubble breakup can be expected in
regions showing high-turbulent eddy dissipation. Figure 6 presents maximum
values of the numerically evaluated turbulent eddy dissipation at the edges of
the obstacle. At the same time, the applied bubble coalescence model of Prince
and Blanch [8] indicates strong importance of coalescence in regions of
bubble accumulation, that is, in the wake behind the obstacle. Both bubble
coalescence (see gas accumulation shown in Figures 4 and 5) and bubble breakup
(see distribution of turbulence dissipation Figure 6), which might partially compensate each other, are expected shortly behind the obstacle.Figure 6
Turbulence eddy dissipation
(run 096) (CFX).Figure7 presents measured cross-sectional averaged bubble size
distributions upstream (z=−0.52 m), shortly behind (z=0.08 m) and
downstream the obstacle (z=0.52 m) for the run 096. The measurements show in
the bubble accumulation zone at z=0.08 m the cross-sectional average shows a
shift toward larger bubbles. The experimentally measured bubble size distribution of run 074 is essentially similar to the run 096. It can be seen on the volumetric fraction maps, Figure 8. Above the obstacle, where turbulent intensity is strong, breakup is not experimentally attested, contrarily to the Luo and Svedsen [9] model prediction, but coalescence occurs, leading to the formation of big bubbles (over 7 mm) that concentrate thereafter in the central part of the pipe.Figure 7
Measured bubble size distribution for run 096.TOPFLOW experimental
volumetric fractions of air bubbles according to different size classes,
extracted from [24].
(a)(b)(c)(d)Both ANSYS CFX and Neptune_CFD calculated bubble size distributions, however, show a shift of the mean bubble diameter toward smaller bubbles
shortly behind the obstacle when both coalescence and breakup are taken into
account. In the calculations, the bubble breakup is overestimated. The
corresponding results are provided by Figure9 for the run 096 and on Figure 10 for run 097. This disagreement was found not solvable by simple changes of breakup or coalescence coefficients, which were set here to FB=FC=0.05. Similar deviations would arise at other
locations of the flow domain. Neptune_CFD results for the run 074 with standard
coefficients FB=FC=1 are shown in Figure 11
(the details can be found in [20]).Figure 9
Bubble size distributions for run 096 (JL=1.017 m/s, JG=0.0898 m/s) (CFX) (FB=FC=0.05).Figure 10
Bubble size distributions
for run 097 (JL=1.611 m/s, JG=0.0898 m/s) (CFX)
(FB=FC=0.05).Figure 11
Neptune_CFD result for
mean Sauter bubble diameter (m), both breakup and coalescence being taken
into account (FB=FC=1) run 074.This suggests to perform computations for which breakup is neglected, but
that still consider the Prince and Blanch [8] coalescence model with standard coefficient. The corresponding Neptune_CFD results show that both bubble size repartition and order of magnitude are consistent with experimental results: big bubbles of more than 7 mm are created and accumulate in the recirculation zone just above the obstacle (see Figure 14). Big bubbles stay in the LHS of the pipe, the unobstructed cross-sectional part of the pipe is free of big bubbles, that is consistent with the experimental observation.As partial conclusions, (i) the Luo and Svendsen breakup model tends to
overestimate the breakup of TOPFLOW experimental tests, whereas (ii) the use of
the Prince and Blanch model for coalescence within the framework of the moment
density method allows to predict satisfying evolution of bubble size across the
flow.


## 5.2.2.
Bubbles Streamlines

More detailed effects of lateral motion of small and large bubbles can
be revealed by studying bubble streamlines and by analyzing lift forces acting
on bubbles of different size. On one hand, the liquid velocity flow carries the
small bubbles into the region behind the obstacle (see Figures12 and 14, right-hand side for the bubble streamlines). Lateral deviation due to lift
force is illustrated on Figure 13 for the lift force arrows with the population
balance method and Figure 14, right-hand side for the lift coefficient with the moment density method. On the other hand, the air accumulation in the wake region leads to bubble coalescence and the generation of large bubbles as revealed by the analysis of experimental results. This phenomenon is
underestimated in the calculations taking breakup into account.Figure 12
Streamlines for small
(left) and large (right) bubbles (run 096) (CFX).Figure 13
Bubble lift force vectors
for the different gas velocity groups (run 096) (CFX).Neptune_CFD numerical
results for local mean bubble diameter and corresponding streamlines colored
by the equivalent mean lift coefficient (run 074)(FB=0;FC=1).
(a)(b)Caused by the lift force, large bubbles are redirected into the
downstream jet (see Figure13) once they can be formed in the wake by
coalescence. The streamline representation (see Figure 12) clearly shows
this phenomenon for large bubbles already present in the upstream flow.In the actual version of the moment density method used in Neptune_CFD,
dynamics of the bubble population is estimated using a single transport
velocity. The averaged lift contribution takes into account both the bubble
size distribution and the Tomiyama correlation. When a majority of bubbles are
locally above the critical Eötvös number, the lift coefficient changes its
sign. In this case, the direction of the lift force is changed, as it can be
seen on Figure14. As a consequence, the bubble streamlines above the obstacle
deviate to the center of the pipe. This explains the low value of the dispersed
phase volumetric fraction near the left-hand side part of the pipe wall above
the obstacle for Neptune_CFD numerical results (see Figure 5). This result is fully consistent with experimental results.As a partial conclusion concerning the bubble streamlines calculations,
both methods showed their ability to consider the effect of bubble size on the
lateral deviation of bubble streamlines due to lift force. This provides a more
precise understanding and a more accurate prediction of bubbles repartition
across the flow.


## 5.3.
Phenomena in the Jet

In the cross-sectional area beside the obstacle, a strong jet is
established creating strong shear flow. The resulting phenomena are more
pronounced with increasing water velocity, like in  run 097, where the liquid velocity was increased toJL=1.611 m/s. Figure 16 represents measured and ANSYS CFX calculated cross-sectional gas fraction distributions for this run. In the most downstream cross section of the measurements an almost gas bubble free region in the centre of the jet is found. Bubbles are collected at the edges of the jet in regions of largest water velocity gradient.The streamline representation of the ANSYS CFX calculations, however,
(Figure12 for run 096, which is fairly similar to run 097) indicates large bubbles being directed into the jet caused by the lift force. This discrepancy between experiment and ANSYS CFX calculations can possibly be explained by the strong water velocity gradient near the jet. This strong shear flow induces
bubble breakup which is not yet considered in the model of Luo and Svendsen
[9]. In the tests, the big bubbles could migrate toward the jet, but be fragmented at the relatively sharp boarder of this jet. Only a small fraction of the small bubbles created by this breakup process can enter the jet by action of the turbulent dispersion force.The gas distribution resolved by bubble size classes (see Figure15) shows clearly the deviations between measurements and calculations. In the experiment, the gas accumulation behind the obstacle leads to a strong coalescence and the creation of large bubbles. In the CFX-calculations, this
effect is underestimated.Calculated by CFX (left)
and measured (right) gas distributions up- and downstream of the obstacle
resolved to bubble size classes (run 096JL=1.017 m/s, JG=0.0898 m/s, FB=0,FC=0.05).
(a)(b)Figure 16
Calculated by CFX (left)
and measured (right) gas cross-fractional distributions downstream the
obstacle (run 097JL=1.611 m/s, JG=0.0898 m/s, FB=FC=0.05).
Calculations (obstacle shown), distances at
z=0.08 m, 0.16 m, 0.25 m, 0.37 m, and 0.52 m. Measurements (obstacle in the upper left area). distances at  z=0.01 m,
0.015 m, 0.02 m, 0.04 m, 0.08 m, 0.16 m, 0.25 m, and 0.52 m.In Neptune_CFD calculations based on the moment density method the
bubble breakup was neglected. As attested by Figure14, left-hand side, in this
case the bubbles present in the jet are small bubbles. Big bubbles created are
concentrated in the center of the pipe where the largest value of volumetric
fraction can be found (Figure 5). The structure of the gas repartition in the
jet for both calculation and experiments is represented on Figure 17. In the
upper part, the structure of the flow of the numerical calculation is
consistent with experimental results. The most important deviation is observed
in the region just above the obstacle (z=0.08m). A gain in accuracy of these
results could be achieved by considering size-dependent velocity dispersion
within the moment density method, in particular for bubble size-dependent
lateral migration in the regions where there exists a strong mixing of
different bubble sizes (mainly just above the obstacle as far as run 074 is
concerned, see Figure 8).Figure 17
Calculated by
Neptune_CFD (left) and measured (right) gas volumetric fraction run 074 for different elevations (obstacle, symbolized in gray, is atz=0 m and covers the “upper” part of the half-pipe
cross-section representation).Nevertheless, the main trends of the bubble dynamics downward the
obstacle are in agreement with experimental results (see Figure16 for CFX and Figure 17 for Neptune_CFD).


## 6. Summary and Perspectives

In this study, we focused on the model of the bubble size distribution
in the numerical simulation of bubbly flows. More deep understanding of the
flow structure is possible when considering a more accurate characterization of
the polydispersion. For upward two-phase flow in vertical pipes the core peak
in the cross-sectional gas fraction distribution could be reproduced very well
both by the moment density and by the MUSIG population balance methods. For complex
flows, the general three-dimensional structure of the flow could be well
reproduced in the simulations.These test cases of pipe flow with internal obstacle demonstrate the
complicated relationship and interference between size-dependent bubble migration,
bubble coalescence, and breakup effects for real flows. With an appropriate
given distribution function, the numerical effort of the moment density method
is lower compared to the multiple bubble size group method (MUSIG). On the
other hand, applying the MUSIG method the simulation of a flow situation allows
to deal with more general shapes for the distribution function. Both methods
enable to consider the effect of polydispersion in size on the bubble
population dynamics, in particular on the evaluation of interphase momentum
transfers associated with lift and drag. The inhomogeneous population balance
model, using several velocity fields for the bubbly phase, is able to deal with
size separation of a locally polydispersed in size population, whereas the
moment density method accounts for local diversity in bubble hydrodynamics
thanks to a single-averaged contribution of interphase transfers. These two
promising refinements of the two-fluid model for bubbly flows have shown their
ability to recover consistent description of the lift force in the upper part
of the flow, where complex flow structures are observed.While the closure models on bubble forces, which are responsible for the
simulation of bubble migration, allow to explain the bubbles hydrodynamic
behavior observed experimentally, clear deviations occur for bubble coalescence
and breakup. Both methods based on similar coalescence and breakup models lead
to the same conclusion: in the simulations of the TOPFLOW experiments, the Luo
and Svendsen model leads to an overestimation of the breakup that appears as
negligible in the experiments. On the other hand, the coalescence model of
Prince and Blanch seems able to recover the correct bubble size, if used
solely, as attested by corresponding Neptune_CFD calculations. For both
methods, the presently applied models describing bubble breakup and coalescence
could be proven as weak points in numerous CFD analyses. These bubble breakup
and coalescence models depend to a large extent on the turbulence properties of
the two-phase flow, which were not measured and could not be validated in the
pipe flow test cases. Therefore, further investigations are necessary to
determine whether the currently used multiphase flow turbulence models deliver
appropriate and verifiable quantities that can be used for the description of
bubble dynamics processes.Extensions of both the moment density method and the MUSIG method to nucleate
boiling regime numerical simulation are in progress. This includes the
phenomena of compressibility, phase-change, and wall nucleation. To model the
bubble size-dependent lateral migration phenomenon, the moment density method
should also include a model for the bubble velocity distribution. This can be
done using a similar formalism to the present model for bubble size
distribution function.



---

*Source: 320738-2008-10-14.xml*