Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Numerical study of bubbly flows using direct quadrature method of moments
S.C.P. Cheung1, G.H. Yeoh2, J.Y. Tu1, E. Krepper3 and D. Lucas3
School of Aerospace, Mechanical and Manufacturing Engineering,
RMIT University, VIC 3083, Australia
chipok.cheung@rmit.edu.au, jiyuan.tu@rmit.edu.au
Australian Nuclear Science and Technology Organization (ANSTO),
NSW 2234, Australia
ghy@ansto.gov.au
Institute of Safety Research, Forschungszentrum Rossendrof e.V.,
P.O. Box 510 119, 01314 Dresden, Germany
e.krepper@fzd.de, d.lucas@fzd.de
Keywords: Population Balance, Method of Moment, Bubbly Flow
Abstract
To model the macroscopic bubble interactions (e.g. coalescence and breakage), the twofluid model in conjunction with the
population balance equation (PBE) approach has been considered as a practical and accurate formulation of handling bubbly
flow systems. Recently, the MUltiple SIze Group (MUSIG) model appears to be one of the most direct solution methods which
solves the PBE with discrete class approach and fuses seamlessly with the computational fluid dynamics (CFD) framework.
Nonetheless, for complex bubbly flow structures with wide range of bubble size, large number of classes must be used to
attain sufficient resolution for the bubble size distribution (BSD). This poses severe limitations on both computational time and
resources. This paper focuses on introducing an alternative direct quadrature method of moments (DQMOM) (Marchisio and
Fox, 2005) where the BSD is tracked through its moments by integrating out the internal coordinate. The main advantage of
DQMOM is that the number of scalars to be solved is very small (i.e. usually 46). To assess the performance of DQMOM in
measureup with the MUSIG model, predictions of both models are validated against two experimental data by Hibiki et al.
(2001) and Lucas et al. (2005). In general, the model predictions compared very well against the measured data. Associated
numerical issues and drawbacks for the DQMOM model are also discussed.
Introduction
Twophase bubbly flow structures are featured in various
industrial applications such as chemical, petroleum, mining,
food and pharmaceutical. For safety analysis or design
optimization, reliable predictions of the void fraction
distribution and other twophase flow parameters are of
paramount importance. With such mounting industrial
interests, population balance modelling has recently become
the centre of focus attempting to synthesize the behaviour of
the bubble mechanisms and it dynamical evolution with
respect to different flow regimes.
In the present stateoftheart, twofluid model can be
considered as one of the most practical and accurate
macroscopic formulations to model the
thermalhydrodynamics of twophase flow systems. Within
the field equations, which are expressed by the conservation
of mass, momentum and energy for each phase, interfacial
Nomenclature
aif interfacial area concentration
a coalescence rate
b
BC, Bb
DC, Db
breakage rate
birth rate due to coalescence and breakup
death rate due to coalescence and breakup
f bubble size distribution function
F, total interfacial force
g gravitational acceleration
j superficial velocity
M mass scale of gas phase (bubble)
N number density of gas phase (bubble)
P pressure
p number of fragments/daughter bubbles
t physical time
u velocity
We Weber number
Greek Sybmol
a void fraction
6 Dirac's delta function
e turbulence kinetic energy dissipation
Ye effective viscosity
5 internal space vector of the PBE
p density
a surface tension
qi weighted abscissas
Subscripts
g gas
i Index of abscissas or gas/liquid phase
I Liquid
min Minimum operator
max Maximum operator
Paper No
transfer terms appear in each of the equations. These terms
require essential closure relations and should be modelled
accurately. Interfacial transfer terms in the twofluid model
are strongly related to the local transfer mechanisms such as
the degree of turbulence near the interfaces and the
interfacial area concentration. Theoretically speaking, the
interfacial area concentration (,, ) is a geometrical parameter
of the local interfacial structure which describes the
available area for the interfacial mass, momentum and
energy transport. All the above interfacial transport
mechanisms between phases are proportional the local
interfacial area concentration. However, the closure
relations for the interfacial transfer terms remain far from
resolution and they still represent the weakest link in the
twofluid model.
Since the interfacial area concentration represents the key
parameter that links the interaction of the phases, much
attention have been concentrated towards better
understanding the coalescence and breakage effects due to
interactions among bubbles and between bubbles and
turbulent eddies for gasliquid bubbly flows. The primary
objective is to better describe the temporal and spatial
evolution of the twophase geometrical structure. Population
balance approaches (Cheung et al., 2008, Wang et al., 2005,
Chen et al., 2004) and volumetric interfacial area transport
equation (Hibiki and Ishii, 2000, Yao and Morel, 2004, Sun
et al., 2004, Wu et al., 1998) have been proposed to predict
the interfacial area concentration.
Benefitting from the early introduction to commercial
package (Lo, 1996), the population balance approach based
on the MUSIG model has been frequently employed to
predict the nonuniform bubble size distribution in a
gasliquid mixture by solving a range of bubble classes.
Although encouraging results have been reported (Chen et
al., 2004, Cheung et al., 2007a), in case of wide range of
bubble sizes in a complex twophase flow system were
being considered, a substantial number of equations might
be required to adequately track the range of bubble sizes.
For flows where large bubbles could exist, especially in
large diameter pipe, computational resource for solving such
large number of transport equations could be extremely
excessive. This model drawback is fundamentally caused by
the fact that it adopts class method to discretize the bubble
size distribution where the pivot size or abscissa of each
class is fixed. In practical calculations where the number of
bubble classes is limited, bubble size distribution cannot be
adequately represented.
In this paper, an alternative approach to predict gasliquid
bubbly flows is presented by the consideration of Method of
Moments (MOM). Here, the bubble size distribution is
tracked through its moments by integrating out the internal
coordinates. The main advantage of MOM is its numerical
economy that condenses the problem substantially by only
tracking the evolution of a small number of moments (i.e.
usually 46). As aforementioned, this becomes rather critical
in modeling complex flow problems when the bubble
dynamics is strongly coupled with already timeconsuming
calculations of turbulence multiphase flows. However, due
to the difficulties related with expressing transport equations
in terms of the moments themselves, the Direct Quadrature
Method of Moments (DQMOM) is applied instead which
essentially involves the direct solution of the transport
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
equations for weights and abscissas of the quadrature
approximation. As will become clearer later, each node of
the quadrature approximation can be treated as a distinct gas
phase. DQMOM, similar to MUSIG, thus offers a powerful
approach for describing polydisperse bubbly flows
undergoing coalescence and breakage processes in the
context of Computational Fluid Dynamics simulations.
Model Description
TwoFluid model for Gasliquid flow
The threedimensional twofluid model solves the
ensembleaveraged of mass and momentum transport
equations governing each phase. These equations can be
written as:
a(pii) + V.(p ) 0
at+V(p, )=0 (1)
Bt
Sava, .)
(pau) + V (piaiuu,)
N.J
aiVP + apg
S+V.x[j, (Vui +(VufY)]+F,
where g is the gravity acceleration vector and P is the
pressure. From the above equation, it is noted that closure
law is required to determine the momentum transfer of the
total interfacial force. This force strongly governs the
distribution of the liquid and gas phases within the flow
volume. On the L.H.S of equation (2), Fi represents the total
interfacial force which is composed of the drag force, lift
force, wall lubrication force and the turbulent dispersion
force respectively. Numerical details on handling these
interfacial forces can be found in Cheung et al. (2007b) and
references therein. For handling the turbulence effects, the
Shear Stress Transport (SST) model is adopted for the liquid
phase (Menter, 1994), while the Sato's bubbleinduced
turbulent viscosity model (Sato et al., 1981) was employed
for the gas phase.
Direct Quadrature Method of Moments (DQMOM)
For the PBE, an integrodifferential form describing the local
Bubble Size Distribution (BSD) can be written as:
af( t)
+ V (( t) f(, t)= S (, t) (3)
at
where f ( t) is the function of bubble size distribution
dependent on the internal space vector whose
components could be characteristics dimensions, surface
area, volume and so on. t is the external variables
representing the physical time in external coordinate
respectively. u (, t) is velocity vector in external space.
The R.H.S of equation (3) is the net source or sink term
of the PBE which denotes the birth and death rates of
bubbles due to coalescence and breakage processes defined
by:
S(<, t) = f a( ', ')f t)f (', t)d(
2 J0 s s 0 '
f (,t)Ofa( ', ') ( t)d '
+ y(5')b(')p(1 f ')f(',t)ds
b()f(, t) (4)
Paper No
Here, the first and second terms denote birth and death rate
of bubble of space vector due to coalescence processes;
the third and fourth terms account for the birth and death
rate caused by the breakage processes
respectively. a(, ') is the coalescence rate between
bubbles of size E and ,'. Conversely, b()) is the breakage
rate of bubbles of size E y({') is the number of
fragments/daughter bubbles generated from the breakage of
a bubble of size ' and p((/ I') represents the probability
density function for a bubble of size E to be generated by
breakage of a bubble of size E'.
Like all method of moment approaches, the basic idea of
DQMOM founded upon a transforming the problem into
lowerorder moments of the size distribution where the
integral of the BSD function is approximated by a finite set
of Dirac's delta functions (McGraw, 1997). Taking the
bubble mass, M, as the internal coordinate, the BSD can
then be expressed as:
f(M,t) Nk(t)(M Mk(t)) (5)
kl1
where Nk presents the number density or weight of the ith
class and consists of all bubbles per unit volume with a
pivot size or abscissa Mk. Obviously, the quadrature method
is brought down to solving 2N unknowns, Nk and Mk. A
number of approaches in the specific evaluation of the
quadrature abscissas and weights have been proposed. With
the aim to solve multidimensional problems, Marchisio and
Fox (2005) extended the method by developing the
DQMOM where the quadrature abscissas and weights are
formulated as transport equations. The main idea of the
method is to keep track of the primitive variables appearing
in the quadrature approximation, instead of moments of the
BSD. As a result, the evaluation of the abscissas and
weights are obtained using matrix operations. More details
of the method can be found in above reference.
In the present study, in order to be consistent with the
variables used in the twofluid model, the weights and
abscissas can be related to the size fraction of the dispersed
phase (fk) and a variable defined as /k = fk /Mk. As a
preliminary study, bubbles are assumed to travel with the
gas velocity, the size fraction of fk is related to the
weights and abscissas by:
Pgofk N,Mi = (k (6)
Using the above expression, the transport equations become:
al tgafk) = bk (7)
gagk c ) + V.,gagg k )=ak (8)
The moment transform of the coalescence and breakup of
the term Sk can then be expressed as:
Sk = B D + B D (9)
where the terms B and D represent the birth and death rates
of the coalescence and breakup of bubbles which is
equivalent to S(,t)in equation (4). On the basis of the
approximation given in equation (6), the birth and death
rates can be written as:
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Bc = ZN iNj( M +M Ya(MiM ) (10)
2 i j
Dc = Mi a(M,M iNj (11)
i j
B = Z M: b(M ,M ,)I (12)
i j
Dk = ZM b(Mi, Mj) (13)
i j
Based on equation (1013), through a series of matrix
operations, the source terms ai and bi can then be
determined and the weights N, and M, can be also evaluated
according to its definition in terms of f, and V, As a
preliminary generic study, the birth and death rates are
evaluated according to the widelyadopted coalescence
kernel by Prince and Blanch (1990) and the breakup
mechanism of Luo and Svendsen (1996).
Experimental Details
The twophase flow experiment conducted by Hibiki et al.
(2001) has been performed at the ThermalHydraulics and
Reactor Safety Laboratory in Purdue University. The test
section was a round tube made of acrylic with an inner
diameter (D) of 50.8 mm and a length (L) of 3061 mm. The
temperature of the apparatus was kept at a constant
temperature (200C) within the deviation of +0.20C by a heat
exchanger installed in a water reservoir. Local flow
measurements using the double sensor and hotfilm
anemometer probes were performed at three axial (height)
locations of z/D = 6.0, 30.3 and 53.5 and 15 radial locations
of r/R = 0 to 0.95. The schematic diagram of the
experimental arrangement is shown in Figure la. A range of
superficial liquid velocities jl and superficial gas velocities
j, have been performed, which covered mostly the bubbly
flow region, including finely dispersed bubbly flow and
50.8mm
3061mm
:53.5
3500mm
z/D= 6.0
51.2mm
 I
z/D=
zID:
IeE I IhF
Bubble Injected from Bubble Injected from
Spurger designated capillaries
Figure 1: Schematic diagram of experimental
arrangement of (a) Hibikit et al. (2001) and (b)
MTLOOP of Lucas et al. (2005)
Paper No
bubblytoslug transition flow regions. Area averaged
superficial gas velocity was obtained from local void
fraction and gas velocity measured by the double sensor
probe, whereas area averaged superficial liquid velocity
was obtained from local void fraction measured by the
double sensor probe and local liquid velocity measured by
the hotfilm anemometry. More details regarding the
experimental setup can be found in Hibiki et al. (2001). In
this paper, numerical predictions have been compared
against local measurements at two flow conditions: Case 1
with =0.491 m/s and =0.0556 m/s; Case 2
=0.986 m/s and =0.113 m/s. The inlet void fractions
are 5% and 10% respectively. The bubble size from the air
injection is 2.5 mm for both cases.
Figure lb depicts the schematic diagram of the MTLOOP
experiment designed by Lucas et al. (2005). Similar to the
experiment by Hibiki et al. (2001), gasliquid bubbly flow
was measured in a vertical medium size cylindrical pipe
with a height of 3500 mm and an inner diameter of 51.2 mm.
Water with a constant temperature of 30C was circulated
from the bottom to the top of the pipe. Air was injected via
an injection device which consists of 19 capillaries with an
inner diameter of 0.8 mm. These capillaries were equally
distributed over the cross section area of the pipe. In order
to measure the local instantaneous gas fraction void fraction
as well as bubble size distribution, an electrode wiremesh
sensor is placed above a certain distance from the injection
device, which can be varied from 30 mm to 3030 mm (i.e.
corresponding to dimensionless axial location Z/D = 0.660).
Gas and liquid flow rates were closely controlled via the
compressed air inlet and circulated pump with maximum
superficial gas and liquid velocities, jg = 14 m/s and jf = 4
m/s, respectively. In this paper, numerical predictions have
been compared against local measurements at two flow
conditions: Case 1 with =1.017 m/s and =0.140 m/s;
Case 2 =1.017 m/s and =0.219 m/s.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
velocities, void fractions for both phases and uniformly
distributed bubble size in accordance with the flow
conditions described above. At the pipe outlet, a relative
average static pressure of zero was specified. A
threedimensional mesh containing hexahedral elements
was generated resulting in a total of 12,000 elements (i.e. 20
radial, 20 circumferential and 40 axial) covering the entire
pipe domain.
For the MTLOOP experiment, a uniform gas volume
fraction was specified at the inlet boundary. On the other
hand, 12 equally spaced point sources of the gas phase were
placed at the circumference of the 600 radial sector to
represent the wall injection method in TOPFLOW. Gas
injection rate at each point source was assumed to be
identical. Details of the boundary conditions are
summarized in Table 1. Based on grid sensitivity test
performed for the MTLOOP and TOPFLOW experiments,
grid independent solutions have revealed that computational
meshes which consisted of 18,223 elements did not
appreciably change even though finer computational meshes
were tested. For all flow conditions, reliable convergence
was achieved within 600 iterations for a satisfied
convergence criterion based on the RMS (Root Mean
Square) residuals of 1.0 x 104 and for a physical time scale
of the fully implicit solution of 0.008s.
Results and discussion
Experiment of Hibiki et al. (2001)
Local radial profiles of the void fraction, interfacial area,
Sauter mean diameter, gas and liquid velocities at two
measuring axial locations of z/D = 6.0 and 53.5 were
predicted through the twofluid model and DQMOM. The
computed results are compared against the measured data of
Hibiki et al. (2001). In order to assess its predictive
capability, additional comparison is also carried out against
Numerical details
For the twofluid model, two sets of equations governing the
conservation of mass and momentum were solved via the
ANSYS Inc, CFX11 computer code. For DQMOM, two
sets of transport equations governing four weights and four
abscissas were chosen to predict the bubble size distribution
of which the evaluation of the source terms ai and bi in
equations (7) and (8) were determined through matrix
operations carried out by an inhouse external subroutine.
Both breakage and coalescence calibration factors, FB and
Fc, were adjusted to 0.15 and 0.05 respectively. Comparing
with our previous study (Cheung et al., 2007b), FB and Fc
were specified to 1.0 and 0.05 in the MUSIG model based
on experimental calibrations. Such discrepancy of
calibration factor between both approaches could be
attributed to the additional flexibility of the DQMOM. As
both weights and abscissas are variables within the
DQMOM, calculation of bubble size distribution could be
very sensitivity to strength of coalescence and breakage
sources. As a result, a different set of calibration factor are
adopted in this study. Radial symmetry was assumed, so that
the numerical simulations were performed on a 600 radial
sector of the pipe with symmetry boundary conditions at
both sides. Inlet conditions were assumed to be
homogeneous in regards to the superficial liquid and gas
05
0 a
Exp
04   MUSIG
 DQMOM
o03 =0.491 m/s
< =6, ILi m/s
02 D=6.0
02
0 
0 02 04 06 08
Radial Location[]
05
Exp
04   MUSIG
 DQMOM
03 =0.986 m/s
=0.113 m/s
z/D=6.0
02 
02
0 02 04 06 08
Radial Location [
05
SExp
04   MUSIG
DQMOM
03 =0.491 m/s
S < =ii,2 m/s
SzD=53.5
02 
01 
0 02 04 06 08
RadialLocation[]
05
Exp
04   MUSIG
DQMOM
S03 =0.986 m/s
S=0.113 m/s
z/D=53.5
02 
01
* .
0 02 04 06 08 1
Radial Location []
Figure 2: Predicted and measured gas volume fraction
profiles at z/D = 6.0 and 53.5 for both flow conditions
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
the predicted MUSIG results obtained through our previous
work in Cheung et al. 12' 1 .i
Figure 2 shows the void fraction distributions of both flow
conditions at the two axial locations for the measured data
and computer results of the DQMOM and MUSIG. In
isothermal gasliquid bubbly flows, Serizawa and Kataoka
(1990) classified the phase distribution patterns into four
basic types of distributions: wall peak, intermediate peak,
core peak and transition. The void fraction peaking near the
pipe wall represented the flow phase distributions caused by
the typical wall peak behaviour. In both flow condition, it
was observed that the wall peaking profile started to develop
at the axial locations of z/D = 6.0 (near the inlet) and
become well established at z/D = 53.5 (near the exit). Model
predictions of both MUSIG and DQMOM captured the
radial void fraction distributions considerably well at the
two locations. Nevertheless, it appeared that DQMOM gave
slightly better predictions especially at the welldeveloped
wall peaking characteristic at z/D = 53.5 in both test cases.
Figure 3 illustrates the Interfacial Area Concentration (IAC)
distributions of both flow conditions at the respective two
axial locations for the measurements and the two model
predictions. The measured data followed the similar profile
as the void fraction distribution as stipulated in Figure 2.
Here again, predictions from both MUSIG and DQMOM
models were in well agreement with measurements. This
further ascertains the predictability of the DQMOM in
comparison with MUSIG model. The Sauter mean bubble
diameter distributions are exemplified in Figure 4. At z/D =
53, good agreement was achieved for DQMOM near the
pipe center while MUSIG underpredicted the bubble sizes
there. For the flow condition of =0.491 m/s and
=0.0556 m/s, DQMOM marginally overpredicted the
bubble sizes but followed similar trend with the
experimental distribution.
: Exp
4   MUSIG
S DQMOM
S35
3 
25
S5 =0.491 m/s
1 < = s m/s
o5 z/D=6.0
0 02 014' 016' 018'
Radial Locaon []
5
45
SExp
4 MUSIG
__ DQMOM
S35
S3
25
=0.986 m/s
=0.113 m/s
z/D=6.0
o0 I I I I ,
0 02 04 06 08
Radial Locaon []
Exp
4   MUSIG
__ DQMOM
t 35 
I  
3
25 0 
2 **
s5 =0.491 m/s
1 < = 6n5 m/s
o5 D=53.5
0 02 04 06 08
Radial Locaon [
4   MUSIG
 DQMOM
35
3
25
I^ .., ., n 
=0.986 m/s
=0.113 m/s
z/D=53.5
0 02 04 06 08
Radial Locaton []
Figure 4: Local predicted and measured sauter mean
bubble diameter distributions at z/D = 6.0 and 53.5 for
both flow conditions
Experiment of Lucas et al. (2005)
Figure 5 shows the predicted volume fraction distributions
obtained from the DQMOM and MUSIG model in
comparison with the MTLOOP measurements at
dimensionless axial position (i.e. Z/D=4.5 and 60) for the
flow conditions of case M107 and M118. In general, it is
600
500 Exp
  MUSIG
 DQMOM
g400 
3300
=0.491 m/s
20 < =6ii6i m/s
200
100
0 02 04 06 08
Radial Locaton []
02 04 06 08
Radial Location [
600 L
500
0 Exp
  MUSIG
 DQMOM
),uu
S =0.491 m/s
200 < 'li m/s
200 zD=53.5
100 
0 02 04 06 08
Radial Locaon []
600
500 xp
S  MUSIG
 DQMOM
400 =0.986 m/s
=0.113 m/s
300 zD=53.5 0
200  
100
0 02 04L 016 08 1
Radial Lcation[
35 Exp
 MUSIG
 DOMOM
30 
25 =1.017 m/s
=0.140 m/s
20 zD=4.5
is 
10 
to
40
25
15
10 =0.219 m/s
s zID=4.5
02 04 06 08
Radial Location []
35 a Exp
 DQMOM
, 30
25
20c 
1 =1.017 m/s
10 =0.219 m/s *
5 z/D=4.5 0
0 02 04 06 08
Radial Location []
S EP.
 MUSIG
SDQOMOM
=1.017 m/s
=0.140 m/s
z/D=60
***... *
Ra o ]
02 04 06 08
Radial Location []
02 04 06 08
Radial Location []
Figure 3: Local predicted and measured interfacial area
concentration profiles at z/D = 6.0 and 53.5 for both flow
conditions
Figure 5: Predicted and measured gas volume fraction
profiles at z/D = 4.5 and 60 for both flow conditions in
MTLOOP experiment
Paper No
Paper No
observed that both population balance models are capable to
capture the transition process from "wall peak" to "core
peak" of the gas volume fraction distribution. This clear
attested that the adopted interfacial force models (especially
for the lift force model) successfully predicted the trend of
the gas phase lateral motions and its sequential demixing
behaviour of bubbles.
The dynamical changes of bubble size distribution dictates
fundamental interfacial area between gas and liquid phase;
which is closely coupled with the interfacial momentum
transfer and the gas volume fraction profile, a close
investigation of both population balance approaches in
predicting the dynamical changes of bubble size is essential
in examining its bubble coalescence and breakage kernels.
Figure 6 shows the predicted crosssectional averaged
bubble size distribution from the DQMOM and MUSIG
model in comparison with the MTLOOP experimental data
at the axial location Z/D=60. Based on our measurement,
close to the gas injection capillaries; where coalescence
processes were started, the initial bubble size distributions
were found to considerably narrow (i.e. ranging from 0 to
9mm). Nonetheless, the bubble size distribution was
significantly broadened after a series of merging procedures.
Especially for the Case M118, the bubble size ranged from 0
to 18mm which is almost twice of its initial size range. In
general, the predicted bubble size distributions by the
MUSIG model are well agreed with the measurements.
On the other hand, the DQMOM model overpredicted the
bubble size distribution. This could be attributed to the
insufficient number of moment adopted in the simulation.
More moments could be tested in future which might better
capture the evolution of bubble size distribution throughout
the system.
8
S6
4
2
1 0 1 2 3 4 5 6 7 8 9101112131415161718
Bubble Diameter [mmun]
O 2 4 6 8 10 12 14 16 18
Bubble Diameter [mm]
Figure 6: Predicted and measured gas volume fraction
profiles at z/D = 6.0 and 53.5 for both flow conditions
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Conclusion
A twofluid model coupled with a population balance model
is presented in this paper to handle isothermal gasliquid
bubbly flows. The DQMOM was implemented in the CFD
code ANSYS Inc., CFX11 to determine the temporal and
spatial geometrical changes of the gas bubbles. Computed
results by the DQMOM twofluid model were assessed
against experimental data by Hibiki et al. (2001) and Lucas
et al. (2005) as well as the computed results from the
MUSIG model based on the class method. Reasonably good
agreements for the void fraction, interfacial area
concentration, bubble Sauter mean diameter and bubble size
distribution have been achieved.
Acknowledgement
The financial support provided by the Australian Research
Council (ARC project ID DP0877734) is gratefully
acknowledged.
References
Chen, P., Sanyal, J. & Dudukovic, M. P. CFD modeling of
bubble columns flows: implementation of population
balance. Chemical Eng. Sci., Vol. 59, 52015207 (2004).
C li'iin:. S. C. P., Yeoh, G. H. and Tu, J. Y. On the numerical
study of isothermal vertical bubbly flow using two
population balance approaches, Chemical Eng. Sci., Vol. 62,
46594674,1 i.I ,
(liLciiin: S. C. P., Yeoh, G. H. & Tu, J. Y. On the modelling
of population balance in isothermal vertical bubbly flows 
Average bubble number density approach. Chemical Eng.
and Proc., Vol. 46, 742756 (2007b).
( lciiil. S. C. P., Yeoh, G. H. & Tu, J. Y. Population balance
modeling of bubbly flows considering the hydrodynamics
and thermomechanical processes. AIChE J., Vol. 54,
16891710 (2008).
Hibiki, T. & Ishii, M. Onegroup interfacial area transport of
bubbly flows in vertical round tubes. Int. J. of Heat and
Mass Transfer, Vol. 43, 27112726 (2000).
Hibiki, T., Ishii, M. & Xiao, Z. Axial interfacial area
transport of vertical bubbly flows. Int. J. of Heat and Mass
Transfer, Vol. 44, 18691888 (2001).
Lo, S. M. Application of Population Balance to CFD
Modelling of Bubbly Flow via the MUSIG model. AEA
Technology, AEAT1096 (1996).
Lucas, D., Krepper, E. & Prasser, H.M. Development of
cocurrent airwater flow in a vertical pipe. Int. J. of
Multiphase Flow, Vol. 31, 13041328 (2005).
Luo, H. & Svendsen, H. F Theoretical model for drop and
bubble breakup in turbulent dispersions. AIChE J., Vol. 42,
12251233 (1996).
Marchisio, D. L. & Fox, R. O. Solution of population
balance equations using the direct quadrature method of
moments. J. of Aerosol Sci., Vol. 36, 4373 (2005).
McGraw, R. Description of aerosol dynamics by the
quadrature method of moments. Aerosol Sci. and Tech., Vol.
27, 255265 (1997).
Menter, F.R. Twoequation eddy viscosity turbulence
models for engineering applications. AIAA J., Vol. 32,
15981605 (1994).
Prince, M. J. & Blanch, H. W. Bubble Coalescence and
Breakup in AirSparged BubbleColumns. AIChE J., Vol.
Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
36, 14851499 (1990).
Sato, Y., Sadatomi, M. & Sekoguchi, K. Momentum and
heat transfer in twophase bubbly flow I. Int. J. of
Multiphase Flow, Vol. 7, 167178 (1981).
Serizawa, A. & Kataoka, I. Turbulence Suppression in
Bubbly 2Phase Flow. Nuclear Eng. and Design, Vol. 122,
116 (1990).
Sun, X. D. Kim, S., Ishii, M. & Beus, S. G. "Modeling of
bubble coalescence and disintegration in confined upward
twophase flow, Nuclear Eng.and Design, Vol. 230, 326,
(2004).
Wang, T. F., Wang, J. F & Jin, Y. Theoretical prediction of
flow regime transition in bubble columns by the population
balance model. Chemical Eng. Sci., Vol. 60, 61996209
(2005).
Wu, Q., Kim, S., Ishii, M. & Beus, S. G. Onegroup
interfacial area transport in vertical bubbly flow. Int. J. of
Heat and Mass Transfer, Vol. 41, 11031112 (1998).
Yao, W. & Morel, C. Volumetric interfacial area prediction
in upward bubbly twophase flow. Int. J. of Heat and Mass
Transfer, Vol. 47, 307328 (2004).
