Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Analysis of Developing Gas/Liquid TwoPhase Flows
Elena A. Tselishcheva*, Michael Z. Podowski*, Steven P. Antal*
Donna Post Guillent, Matthias Beyer*, Dirk Lucas*
(*) Center for Multiphase Research, MANE Department, Rensselaer Polytechnic Institute, Troy, NY 12180 USA
(t) Idaho National Laboratory Idaho Falls, ID 83415, USA
() Forschungszentrum DresdenRossendorf e.V Institute of Safety Research, Dresden, Germany
tselie@ rpi.edu, podowm@rpi.edu, antals@rpi.edu, Donna.Guillen@inl.gov, M.Beyer fzd.de, d.lucas fzd.de
Keywords: developing gas/liquid twophase flow, Computational Multiphase Fluid Dynamics, bubbly and churnturbulent
flow regimes, NPHASECMFD code
Abstract
The objective of this work has been to study the mechanisms governing flow and phase distributions in developing gas/liquid
two phase flows in general, and the evolution of different size bubbles in an adiabatic vertical pipe in particular. Flow
regimes from bubbly to churturbulent have been accounted for. The main emphasis of the work has been on the modeling
of various interfacial forces between the dispersed bubbles and the continuous liquid, as well as of bubble/bubble interactions
(coalescence and breakup).
The proposed modeling concept uses a complete set of transport equations for each field, such as the continuous liquid and
dispersed bubble fields. The overall model has been implemented in a stateoftheart computational multiphase fluid
dynamics code, NPHASECMFD. This threedimensional fourfield model, including the continuous liquid field and three
dispersed gas fields representing bubbles of different sizes, has been carefully tested for numerical convergence and accuracy,
and then validated against the TOPFLOW experimental results.
The NPHASECMFD simulations were aimed at demonstrating the capability of the proposed modeling concepts to predict
the evolution of bubble concentration from channel inlet to nearequilibrium (fullydeveloped) conditions downstream. Along
with several interfacial closure laws, the effect of elevation on air density has also been included in the model.
Introduction
Computational Multiphase Fluid Dynamics (CMFD) has
become a very useful tool for the analysis and design
optimization of a large class of multicomponent flow
systems. However, the accuracy of numerical predictions
for gas/liquid twophase flows using CMFD methods
strongly depends on the proper formulation of models
governing the underlying local physical phenomena.
The purpose of this project has been to develop, test and
validate a multifield model of adiabatic gas/liquid flows at
intermediate gas concentrations (e.g., bubbly and
churnturbulent flow regimes), in which multiplesize
bubbles are divided into a specified number of groups,
each representing a prescribed range of sizes. In particular,
the ability of the proposed modeling concept has been
investigated to predict the evolution of bubble
concentration from channel inlet to nearequilibrium
(fullydeveloped) conditions downstream along an
adiabatic vertical pipe.
The simulations were performed using a stateoftheart
computational multiphase fluid dynamics code,
NPHASECMFD (Antal et al., 2000). A complete
fourfield model, including the continuous liquid field and
three dispersed gas fields representing bubbles of different
sizes, was first carefully tested for numerical convergence
and accuracy, and then used to reproduce the TOPFLOW
experimental data (Prasser et al, 2007).
Nomenclature
ij,k
L
P
v
x,y,z
m
M
g
Numerical indexes
Superficial velocity (ms1)
Length (m)
Pressure (Pa)
Velocity vector (ms1)
Coordinates
Volumetric mass transfer between fields
representing the same phase
Interfacial force per unit volume (Nm3)
Acceleration of gravitation (ms2)
Paper No
A "' Interfacial area density (m'1)
urel Relative velocity of the dispersed field (m s')
Re Reynolds Number
D Diameter (m)
Db Diameter of bubble (mm)
CD Drag force coefficient
CL Lift force coefficient
C, Wall force coefficient
CD Turbulent dispersion coefficient
k Turbulent kinetic energy (m2S2)
K, Measure of the probability of coalescence or
breakup
Pc Shear production rate
P,, Relative probability of coalescence
ReB Bubble relative Reynolds number
f, Frequency of collisions
y, Distance from the wall (m)
Greek letters
a Volumetric gas fraction
p Dynamic viscosity (kg ms 1')
r Shear stress (kg m's2)
p Density (kg m3)
v Kinematic viscosity (m2s1)
e Turbulent dissipation rate (m2s3)
qki Measure of the effect of departure of large
bubbles' shape from spherical
a Surface tension
Subscripts
l,v Indexes for liquid and gas respectively
t Turbulent
b Bubble
Model formulation for churnturbulent flows
The multifield modeling concept of multiphase/
multicomponent flows is based on ensemble averaging the
governing equations for each component fluid. Such
modeling, with appropriate closure laws, is capable of
capturing flow regimes from bubbly through
churturbulent and slug, to annular flow. A typical form
of the resultant conservation equations for adiabatic flows
is shown below.
Mass conservation
0(akPk) + ( ) M (1)
+V.(apkvk= (1)
Momentum conservation
' +V(k PkVkV,)=k
(2)
tkVp +V.(k tk__k)+kPkg+ M, +Emv,kVm,
J
where mk= m, .
m
In Eqs.(1) and (2), ak is the volumetric fraction of fieldk,
mk is the volumetric mass transfer term into fieldk from
other fields representing the same phase, M' is the
interfacial momentum transfer per unit time interfaciall
force) between fields k and j, = + R is the total
y is th oa
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
shear stress term, the subscript 'i' refers to interfacial
variables, and the remaining notation is conventional.
The interfacial interactions between the individual bubble
fields are specified by mechanistic models for both drag
and non drag forces. The total interfacial force on phasek
can be expressed as a superposition of several component
forces (Drew, 1992):
M' = MD +1"I' +MTD +M +" (3)
where Mf is the drag force, 1I' is the virtual mass
force, MD is the turbulent dispersion force, M' is the
lift force and I" is the wall force.
Numerous closure laws can be found in literature to model
the various interfacial forces, usually based on theoretical
or experimental studies on single bubble. In the axial
direction, the drag force and virtual mass force (Drew,
1992) play major roles. The drag force model is given by
M= D VI,( VI,)A" (4)
8
Several models have been developed for the flow regime
dependent drag coefficient, CD. In the present study, the
expression proposed by Wallis (1976) has been used
S24 1+0.1Re
C .= Re, L
for Reb 1000
[0.44 for Re> 1000
where Re is the bubble Reynolds number relative to the
surrounding it liquid
Reb l (6)
In the radial direction, the nondrag interfacial forces are
dominant and control the gas profile. One of them is the
turbulent dispersion force. For small bubbles, this force
can be evaluated from the model proposed by Podowski
(2009). The radial component of this force is
MD7 CTP,ciV[(1 a)v vc]
(7)
=C pVa v'c 'c+CTDp,a(la)v[v', v']
where for isotropic turbulence, CT = 2/3.
The first term of the RHS of Eq.(7) is normally dominant
across most of the flow area, except for the nearwall
region. In the latter region, the magnitude of the second
term is significantly higher. Thus, the second term plays in
fact the role of a wall force, preventing small bubbles from
touching the wall.
For large bubbles, the turbulent dispersion force in the
radial direction can obtained from
MA = Cp,aVa (8)
The use of different formulations for different bubble sizes
(and also different values of the coefficient, C, ) is due to
the fact that the force on large bubbles is caused by
bubbleinduced turbulence rather than by the microscale
turbulence inside the liquidfield.
The lift force is used to account for the interfacial
momentum exchange between bubbles and the liquid field.
In the radial direction, this force can be written as
MW, = CQp, v(V, V,) x (Vx V, ) (9)
For turbulent twophase flows, this coefficient is typically
between 0.03 and 0.1. However, for large nonspherical
bubbles, the lift coefficient can not only reduce all the way
Paper No
to zero, but even change sign. The critical diameter at
which the lift coefficient changes sign has been estimated
at 4 mm by Kurul and Podowski (1988) and at 5.8 mm by
Tomiyama (1998).
Also, for large bubbles the alternate approach to model the
effect of wall force has been used in the form of a
lubrication force
M\ =C1pa n. (10)
where the wall force coefficient, C, is (Antal, 2005)
where the wall force coefficient, Ca,, is (Antal, 2005)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
=1l if i=j
0 otherwise
In Eq.(15), acr is the maximum concentration, or the
critical packing limit, of small bubbles (acr=0.74 for
spherical bubbles).
The relative probabilities of coalescence are normally
estimated based on experimental observations. In general
they decrease with the increasing bubble diameter (Prince
and Blanch, 1990; Luo and Svendsen, 1996). A typical
range for small bubbles can be estimated as
j 3) for y
Dfor y
for y> D
In the present simulations, it has been assumed that three
groups of bubbles of different diameters enter the flow
channel. As soon as the bubbles start coalescing, new,
larger than original, bubbles are formed. This in turn
leads to the formation of a large spectrum of bubbles of
different sizes. Eventually, the opposite process is also
likely to occur, i.e. large deformed bubbles start breaking
up, leaving behind small bubbles. Thus, it is very
important that a mechanistic and accurate, yet
computationally effective, mass transfer model be used to
model the interactions between bubbles of different sizes
(Kumbaro and Podowski, 2006). In the present model, the
total volumetric mass transfer for fieldk is defined as
m' = m,co mkbr (12)
i\ wh.l i"' c = mckco is the rate of mass increase per unit
volume of bubbles of groupk due to coalescence of
bubbles of all groups, and m"kb = ,'I is the mass
m
loss rate due to breakup of groupk bubbles into bubbles of
different sizes.
A common approach of the modeling of bubble
coalescence is based on a multiplegroup concept.
Specifically, if Db,, is the bubble diameter of bubble
groupi, the rate of coalescence between bubbles of groupi
and groupj is given in the form
myn =k Ka,a (13)
where K, is a measure of the probability of coalescence of
two bubbles. Several models have been proposed to date
for the coefficients, K,. The model used in current studies
is defined as
K,y = CykPf, (14)
The group coefficients, Cy are the measures of
probability that the coalescence of bubbles in group sizes i
and j produce bubble of sizek. They can be determined
using the modeling approach proposed by Podowski
(2009), which is based on directly using a discrete
combinatorial method.
The frequency of collisions is given by the following
expression (Antal et al, 2005):
IV, v d, +,Iv vd, IVV, I
fk ccd, +d 2) 2(c )3 (15)
where
0.4 < P( <0.8
(16a)
whereas for smalltolarge bubble coalescence we have
0.1
(16b)
The breakup of the large cap bubbles to form small
bubbles can be expressed as:
mk, = Kki k =Ck Pdfbr.k,k
It has been shown before that the volumetric mass transfer
rate for large bubbles to form small bubbles depends
mostly on the surface tension (Prince and Blanch, 1990;
Luo and Svendsen, 1996). It is also influenced by the local
turbulence level, and tends to increase with increasing
energy dissipation rate. Taking into the account both
factors, as well as the effects of the inertial and buoyancy
forces, the following expression can be used for the
frequency of large bubble breakup (Lehr et al, 2002)
27 0 4 67k, IVk V 0 2
fb0,kj 67J (18)
) Li1 P'?D [IDk1 4]
The group coefficients for break up model Ck, can be
evaluated in a manner similar to that for the bubble
coalescence model. The coefficient 77, is a measure of
the effect of departure of large bubbles' shape from
spherical.
The effective viscosity of the continuous liquid field of the
is given by
1V = + v + (21)
where v, is the molecular kinematic viscosity of the
liquid, C,k2 is the turbulent shearinduced kinematic
viscosity, and v,2 is the bubble bubbleinduced kinematic
viscosity.
The turbulent shearinduced kinematic viscosity of the
liquid component is modeled by a modified twophase
flow version of the High Reynolds Number ks model
(Antal et al. 2005).
The bubbleinduced kinematic viscosity is modeled using
the following expression (Sato and Sekoguchi, 1975)
= C,,Dk akVd,k V (22)
kl1
where Cb = 1.2.
C 1
0
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Experimental Facility
The experimental data used for the validation of the
present model have been taken at the thermaldynamic test
facility TOPFLOW (Prasser et al., 2007) at the Institute of
Safety Research of Forschungszentrum Rossendorf e.V.
The TOPFLOW facility has been equipped with a vertical
test section, which is a stainless steel pipe with an inner
diameter of 195.3 mm (DN200). The total height of the
section is 9 m. The flow rates in the test section DN200
(see Fig. 1) may be assigned over the following range: the
superficial gas velocity (f,) between 0.0025 and 7.772 m/s,
and the superficial liquid velocity (il) from 0.0405 to 1.61
m/s. Desalinated water was used in the experiments. The
metering system for the injected air flow supplied
volumetric flow rates related to standard conditions
(p=0.25MPa, T=300C). The test section DN200 has been
equipped with six gas injection locations which allowed
for injecting air or steam via orifices in the pipe wall.
This gas injection via wall orifices offers the advantage
that the twophase flow can raise smoothly to the
measurement plane, without being influenced by the feeder
within the tube at any other locations along the flow. The
inlet part of the test section has been connected to a gas
injection pipe and a compressed air system. The liquid
phase has been supplied from the bottom of the test
section. The measurement plane was always situated at the
upper end of the test section shown in Fig. 1.
Two wiremesh sensors were used. The lower sensor has
been used to obtain data on gas volume fraction profiles
and bubble size distributions, while the purpose of the
second sensor was to determine gas velocity by
crosscorrelation measurements between the two sensors.
Comparative measurements between the wiremesh sensor
and other research methods provided information on the
accuracy of the measurement technique and the evaluation
algorithms for the experimental determination of main
flow parameters. The accuracy of the gas volume
fraction averaged over the flow crosssection depends on
the twophase flow regime. Differences in the absolute
void fraction were determined for bubbly flows in the
range of + 1 %, and a systematic underestimation of about
 4 % was observed for slug flows. The database
established during the experiment contains data for the
evolution of the flow along the pipe, including radial
profiles of the void fraction and gas velocity, as well as
bubble size distributions.
Overview of the NPHASECMFD Computer Code
The NPHASECMFD code (Antal et al, 2000) is a
nominally pressurebased finite volume computational
multiphase fluid dynamics (CMFD) flow solver. The
individual conservation equations for mass, momentum
and energy, are solved for each field, together with the
equations for continuous field turbulence. The mixture and
field continuity equations are solved in coupled and
segregated (uncoupled) manner, using stationary
coefficient linearization. The code is fully unstructured
and can utilize secondorder accurate convection and
diffusion discretization. The technology used by the
NPHASECMFD code is an ensemble averaged multifield
model of twophase or multiphase flows.
water
injection
Figure 1. Vertical test section of the TOPFLOW facility
with variable gas injection system (DN 200).
Key features of NPHASECMFD code include the
following:
* Use of unstructured grids with arbitrary element types
* Capability to model an arbitrary number of fields
(fluid components and/or phases)
* Builtin or userdefined mechanistic modeling,
integrated with numerics
* Good robustness and numerical convergence
* Free surface modeling
* Parallel processing via MPI
Computational Grid and Boundary Conditions
The computational domain used in the NPHASECMFD
simulations was consistent with the geometry of the
TOPFLOW test section, and was shaped as a vertical
circular tube, L= 9 m long and D = 0.194 m in diameter. A
sample grid is shown in Figure 2. As it can be seen,
nonuniform nodalization schemes have been used in both
the radial and axial directions. Such an approach allowed
one to capture flow details near the inlet, as well as the
effect of high gradients of velocity and bubble
concentration near the wall of the tube. GRIDGEN was
used as a computational grid generator, to build the desired
mesh and specify the necessary boundary conditions.
Table 1 presents the initial conditions used in testing and
validation of the NPHASECMFD based model. These
conditions correspond to the TOPFLOW experimental runs
shown in the first column of Table 1.
Paper No
Paper No
Figure 2. Schematic of grid
NPHASECMFD simulations.
geometry
used in
Table 1. The initial conditions used in simulations
Sn Superficial liquid Superficial gas
velocity (m/s) velocity (m/s)
107 1.017 0.145
118 1.017 0.219
Four groups of bubbles of different sizes have been
identified in the TOPFLOW experiments. The majority
of NPHASECMFD simulations have been performed for
three bubble groups. This is because the measured
concentrations of the two smallest bubble groups were
very small (of the order of 3% or less), so that
differentiating between them would be below the
resolution level of the overall model.
The range of bubble sizes and the average bubble diameter
for each modeled group were as follows:
Group1: Size range from 0 to 5.8 mm; Dbl = 3 mm,
Group2: Size range from 5.8 mm to 7 mm; Db2= 6 mm,
Group3: Size range over 7 mm; Db3= 20 mm.
Parametric testing has shown that the differences between
the results obtained using the integer values of bubble
diameters given above and the arithmetic averages over
each range were negligible.
During the experiment the temperature of the twophase
mixture in the TOPFLOW experiments was kept constant
at 300C. However, since the height of the experimental
pipe was 9 m, the hydrostatic pressure drop between the
inlet and the exit was comparable to the atmospheric
pressure. This, in turn, has a significant effect on the
density of the air, which decreased by nearly a factor of
two compared to the initial density at a nearly atmospheric
pressure. Thus, it was deemed important that this factor be
accounted for in the model. Specifically, a
positiondependent gas density has been assumed in the
calculations, changing according to the ideal gas law
P(, \ / Pv,oUt
P (x) = (xo,, x) + Pv, ,, (23)
L
where P,,, =  is the air density at the inlet, pY,, t
RT
is the air density at the pressure boundary (conduit's
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
outlet), x,,, is an axial coordinate of pressure boundary, x
is the local axial coordinate, and L = xo, xn is the
total length of the conduit.
The inlet pressure was estimated based on the hydrostatic
head of air/water mixture corresponding to the average
void fraction observed in the experiments.
Results and discussion
Extensive testing of the proposed model has been
performed using as a reference the experimental conditions
corresponding to various data series of the TOPFLOW
experiments (Prasser et al., 2007).
Several parametric simulations have been performed, in
which the sensitivity of results to selected modeling
assumptions has been assessed. The results included in the
current paper show the parametric testing on the effects of
gas density changes and of the values of selected
interfacial force coefficients. The simulations have been
performed for the liquid and gas flow rate corresponding to
TOPFLOW Run 118.
To test the effect of gas density changes along the flow are
shown in Figure 3, where the calculated axial distribution
of gas density, along with the distributions of average void
fraction, are plotted for the cases of variable and constant
gas density. It can be seen that the average void fraction
increases considerably with decreasing gas density,
whereas it slightly decreases (due to gas acceleration) for a
fixed gas density case. Thus, the changing (localpressure
dependent) air density not only affects the magnitude of
void fraction, but may also overcome the effect of
buoyancydriven bubble velocity increase.
30 3
25
,

. 20
20
.t
15
S10
)
5
0
0 10 20 30
UD
2.5
2
1.5
 1
y 0.5
40 50
Figure 3. The effect of variable gas density on the average
void fraction distribution along the TOPFLOW test
section.
Figure 4 shows the effect of drag force on flow
distribution. The Wallis correlation, discussed above,
calculates the drag coefficient for small bubbles with a
good accuracy. However, in the case of large bubbles,
namely for Re>1000, its accuracy diminishes. The main
reason is that it does not properly account for the effect of
bubble deformation. As the bubble aspect ratio increases,
the drag coefficient also increases. The diameter of large
bubbles is still much smaller that the pipe diameter, so that
those bubbles will maintain a nearellipsoidal shape all the
way along the flow. Several values of the drag coefficient
 local average void fraction, Pin 0.25MPa
 local average void fraction, constant density
 gas density
Paper No
have been considered for test. They extend over the range
from spherical bubbles to flat circular disks facing the flow
(CD=1.2). The present study has been performed for a
twofield model of twophase flow, consisting of a
continuous liquid and large bubbles with Db = 20 mm.
As it can be seen in Figure 4, bubble velocity decreases
with increasing CD, as expected. At the same time, there is
an increase in the corresponding average void fraction.
Interestingly, the void fraction distribution flattens as CD
increases. At the same time, the liquid velocity has not
been affected by the value of drag coefficient. As a result
of the present parametric testing, CD=0.88 has been chosen
as representative value of the drag coefficient for large
bubbles for model validation purposes.
S 0.2
S0.15
0.1
0.05
0 0.02 0.04 0.06
Radius,m
2.5
2
1.5
1
0.5
0
0.08 0.1
C,
0 0.02 0.04 0.06 0.08
Radius,m
l Liquid velocity Cd=0.5  Liquid velocity Cd=0.6
 Liquid velocity Cd=0.8
 Gas velocity Cd 0.5
 Gas velocity Cd 0.8
o Gas velocity experiment
 Liquid velocity Cd 1.2
 Gas velocity Cd=0.6
 Gas velocity Cd 1.2
Liquid velocity experiment
Figure 4. The effect of drag coefficient on the radial
distributions of void fraction and phasic velocities. The
calculated parameters correspond to a distance, L/D = 40
from pipe inlet; Run 118.
The experimental studies reported in literature on the effect
of bubble deformation on the direction of lift force (Zun,
1980; Ervin & Tryggvason, 1997; Tomiyama, 1998)
clearly show that the direction of this force changes if a
substantial deformation of bubble shape occurs. Whereas
small spherical bubbles typically travel toward the wall,
large deformed bubbles behave in a different way. Namely,
asymmetric deformations, caused by the combined effects
of buoyancy and shear, lead to secondary flows around
a0
.^
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
such bubbles, which change both the interfacial shear and
pressure distribution around the bubble. This results in a
change of direction of the lift force, compared to that for
spherical bubbles. According to the observations of
Tomiyama (1998), the critical bubble diameter at which
lift force changes sign is approximately 5.8 mm for
air/water flows at ambient conditions.
The impact of both the sign and magnitude on the radial
distribution of large, 20 mm in diameter, bubbles and on
the corresponding velocity profiles are shown in Figure 5.
0.35
0.3
0.25
0.2
 CL=0.2
0.0.06 0.08 0.15
0.1 .CL=O0.1
0.1 ...... ...... CL=0
a CL 0.1
0.05
 CL=0.2
0 0.02 0.0 ,dis 0.06 0.08 0.1
2.1
0 0.02 0.04 0.06 0.08 0.1
Radius, m
Liquid velocity, CL=0.2 Gas velocity, CL=0.2
 Liquid velocity, CL 0.1 Gas velocity, CL=0.1
.. ...... Liquid velocity, CL 0 ...... ...... Gas velocity, CL=O
 Liquid velocity, CL 0.1 A Gas velocity, CL=0.1
Liquid velocity, CL 0.2 Gas velocity, CL=0.2
Figure 5. Radial void fraction and velocity profiles of
large bubbles for different lift force coefficients. The
calculated parameters correspond to a distance, L/D = 40
from pipe inlet; Run 118.
These results, which have been obtained for five different
values of the lift coefficient, from CL = 0.2 to CL = 0.2,
clearly demonstrate how the sign of lift coefficient affects
bubble concentration. As expected, as the sign changes
from positive to negative, the corresponding maximum
value of the volume fraction has moved from the near the
wall location to the center of the pipe.
The experimental data which have been used for the
validation of the proposed complete multifield model of
gas/liquid flows correspond to TOPFLOW Run 118 and
m~m~!!t
~nrrm~~ ~. ,,~~~
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Run 107. In each simulation, the local gas and liquid
velocity have been evaluated, along with the concentration
distribution for each bubble group and the total void
fraction. The radial profiles of those parameters have
been plotted at different streamwise locations and
compared against the experimental data.
Due the uncertainties associated with the effect of gas
injection condition on flow conditions there, the
measurements taken at two different distances (both short
compared to the total pipe height) downstream from the
injection zone have been used as inlet condition for
NPHASECMFD simulations.
Figures 6 and 7 show the results for Run 118 obtain using
the inlet conditions corresponding to the TOPFLOW data
collected at the axial location, L/D=13 along the pipe. As
it can be seen, the NPHASE results agree well with
experiment. The velocity distributions in Figure 6 are
consistent with data in the bulk region, while they are
slightly overpredicted near the wall. This is because there
is a swarm of small and intermediate bubbles near the wall
which affect the velocity of both phases in that region.
0
2.5
2
S1.5
" 1
0.5
0
0.05
E
0.04
0
o
0.03
20.0
> 0.02
0.03
0.025
 0.02
.o
0 0.015
U
0.01
0.005
0
(a)
L/D=13 nphase
L/D=22 nphase
L/D=40 nphase
L/D=13 experiment
X L/D=22 experiment
L/D=40 experiment
0 0.02 0.04 0.06 0.08 0.1
Radius
0 0.02 0.04 0.06 0.08 0.1
Radius
Figure 6. Evolution of a) liquid and b) gas velocity along
the vertical channel for the case of fourfield flow,
compared against experimental data for Run 118. Starting
point L/D=13.
0 0.02 0.04 0.06
Radius
0 0.02 0.04 0.06
Radius
0 0.02 0.04 0.06
Radius
0.08 0.1
0.08 0.1
0.08 0.1
0 0.02 0.04 0.06 0.08 0.1
Radius
Figure 7. Evolution of concentration along the vertical
channel for the case of fourfield flow: a) Db=3 mm; b)
Db=6 mm; c) Db=20 mm; d) total void fraction. Run 118,
starting point LD= 13.
Paper No
(b)
L/D=13 nphase
 L/D=22 nphase
 L/D=40nphase
L/D=13 experiment
L/D=22 experiment
L/D=40 experiment
2.5
1.5
0.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Although the lift force for this bubble group is large
enough to push the bubbles toward the wall, a small peak
in the predicted distribution can still be observed near the
center. On the other hand, the distribution of 3 mm
diameter bubbles closely follows the experimental trend.
It is interesting to notice that the velocity profiles change
only slightly along the flow between e, L/D=13 and
L D=40. The observed change in the gas concentration
profiles is more significant in magnitude, but there have
been no dramatic shape changes.
To better qualify the ability of the present model to predict
twophase flow parameter evolution in developing flows, a
series of simulations has been performed using the axial
location, L/D=7.7, as the starting point. The results of
simulations are shown in Figures 8 and 9. As it can be
seen, the volumetric concentration profile for small
bubbles (DB = 3 mm), shown in Figure 8(a) has been
captured quite well, including the nearwall peak. The
size of the intermediate group of bubbles (with
DB = 6 mm), belongs to the bubble diameter range for
which the lift force changes direction, so concentration
distribution in Figure 8(b) is nearly uniform, except for a
slight peak near the wall. As shown in Figure 8(c), the
agreement for large bubbles (DB = 20 mm) is also very
good, and the trend in the axial evolution of their radial
concentration distribution has been captured in a consistent
manner. Similar conclusions apply to the total void fraction
distributions, shown in Figure 8(d). The results in Figure 9
indicate that a very good agreement with the data has also
been obtained for the velocity distributions.
The evolution of the radial concentrations of differentsize
bubbles, as well as of the gas and liquid velocity
magnitude profiles, are shown in Figure 10 for the entire
computational domain. The left side of each color contour
represents the wall boundary. The results in Figure 10
correspond to the same conditions as those shown in
Figures 8 and 9. They allow one to better evaluate the
developing flow characteristics along the channel. It can
be seen that small bubbles stay near the wall all the way
along the channel, while bubbles with a diameter of 6 mm
are distributed uniformly along the pipe radius, with the
maximum value of volumetric concentration around 3%
near the outlet. At the same time, the large bubbles
experience highest concentrations near the center of the
pipe.
Since the previously discussed results have used a single
TOPFOW run as a reference, it was deemed important that
the new model be also validated against another set of flow
conditions. Run 107 has been selected for this purpose,
corresponding to the following flow conditions: jl= 1.017
m/s, j, = 0.14 m/s. These conditions also correspond to the
churturbulent flow regime. Again, two series of
simulations have been performed, one with the inlet
conditions corresponding to L/D=13, the other to L/D=7.7.
The results of simulations corresponding to the inlet
conditions at L/D=13 are shown in Figures 1112. The
NPHASECMFD results again agree quite well with the
experimental data.
Comparing the data for Run 107 (Fig. 11) with those for
Run 118 (see Fig. 6), it can be observed that the profiles of
void fraction for small and medium bubbles in the radial
direction in Run 107 are nearly uniform (see Fig. 11 (b)).
0
0.03
0.025
E
S0.02
o
S0.
S0.01
S0.01
0.005
0
0.35
0.3
E 0.25
S0.2
0.15
o 0.1
0 0.02 0.04 0.06
Radius
0.08 0.1
0.02 0.04 0.06 0.08 0.1
Radius
0 0.02 0.04 0.06
Radius
0.08 0.1
0.2 
> o.15
a  L/D=7.7 npha
H
> 0.15 L/D=13 nphase
L/D=22 nphaser
S0.1 L/D=40 nphase
SL/D=7.7 experiment
0.05 L/D=13 experiment
o L/D=22 experiment
L/D=40 experiment
0 0.02 0.04 0.06 0.08 0.1
Radius
Figure 8. Evolution of concentration along the vertical
channel for the case of fourfield flow: a) Db=3 mm; b)
Db=6 mm; c) Db=20 mm; d) total void fraction. Run 118,
starting point L/D=7.7.
Paper No
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Taking into account that the lift force coefficients for these
bubbles are the same as in the previous simulations, the
void fraction distributions of both small bubble groups has
been captured relatively well, although the nearwall
concentration of 3 mm diameter bubbles has been slightly
overpredicted.
1.5
0 0.02 0.04 0.06 0.08 0.1
Radius
L/D=7.7 nphase
L/D=13 nphase
 L/D=22 nphase
L/D=40 nphase
L/D=7.7 experiment
L/D=13 experiment
L/D=22 experiment
L/D=40 experiment
0.02 0.04Rad 0.06
Radius
0.05
S0.04
E
E
o 0.03
0.02
o
0
> ntr
0 0.02 0.04 0.06 0.08 0.1
Radius
0.03
(b)
E
 0.02
S0.01
0
0A
0 0.02 0.04 0.06
Radius
0.25
0.08 0.1
Figure 9. Evolution of flow velocity along the vertical
channel for the case of fourfield flow; Run 118, starting
point L/D=7.7.
 Db=3mm Db=6mm Db=20 mm
0644 0.0287 j.278
0.0600
o0.0400 ~0.0200 
0.00159 0.00443 0.000
Gas velocity
.12
I2.00
1.60
1.20
0.800
0.739
Figure 10. The side views of channel flow. The
distributions of void fraction for each dispersed field,
average gas and liquid velocities, Run 118.
0 0.02 0.04Radius0.06
Radius
0.08 0.1
0.08 0.1
0 0.02 0.04 0.06 0.08
Radius
L/D=13 nphase L/D=22 nphase
 L/D=40 nphase
X L/D=22 experiment
* L/D=13 experiment
* L/D=40 experiment
Figure 11. Evolution of concentration along the vertical
channel for the case of fourfield flow: a) Db=3 mm; b)
Db=6 mm; c) Db=20 mm; d) total void fraction. Run 107,
starting point L/D= 13.
Paper No
L/D=7.7 nphase
L/D=13 nphase
 L/D=22 nphase
 L/D=40 nphase
L/D=7.7 experiment
L/D=13 experiment
L/D=22 experiment
L/D=40 experiment
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
2.5
2 
E 1.5
c 1
0.5
0
0 0.02 0.04 0.06
Radius
0.08 0.1
0 0.02 0.04 Radiu 0.06 0.08 0.
Radius
 L/D=22 nphase gas
L/D=13 experiment gas
o L/D=40 experiment gas
 L/D=22 nphase liquid
L/D=13 experiment liquid
L/D=40 experiment liquid
Figure 12. The development of velocity of liquid and gas
along the channel. The comparison of radial distributions
obtained by NPHASE against the experimental data for the
Run 107, starting from L/D=13.
2
18 ( (a)
1.6 
1.4 
1.2
1
0 0.02 0.04 R s 0.06 0.08 0.1
Radius
1.5
_,, (b)
1.5 1 
11
0.9
0.7 ,
0 0.02 0.04 0.06 0.08 0.1
Radius
 L/D=7.7 NPHASE L/D=13 NPHASE
 L/D=22 NPHASE  L/D=40 NPHASE
L/D=7.7 experiment L/D=13 experiment
L/D=22 experiment L/D=40 experiment
Figure 13. The development of velocity of a) gas and b)
liquid along the channel. The comparison of radial
distributions obtained by NPHASE against the
experimental data for the Run 107. Starting point L/D=7.7.
0.02
0.015
0.01
0.005
n
0 0.02 0.04Radi 0.06
Radius
0 0.02 0.04 Radius.06
0.08 0.1
0.08 0.1
0 0.02 0.04Radu 0.06 0.08 0.1
Figure 14. Evolution of concentration along the vertical
channel for the case of fourfield flow: a) Db=3 mm; b)
Db=6 mm; c) Db=20 mm; d) total void fraction. Run 107,
starting point L/D=7.7.
Paper No
0.06
0.05
0.04
0.03
0.02
0.01
0
 L/D=7.7 NHPASE (a)
L/D=13 NPHASE
 L/D=22 NPHASE
 L/D=40 NPHASE
L/D=7.7 experiment rA
a L/D=13 experiment
X L/D=22 experiment
L/D=40 experiment 
'  ":5^
I
 L/D=13 nphase gas
 L/D=40 nphase gas
L/D=22 experiment gas
L/D=13 nphase liquid
 L/D=40 nphase liquid
L/D=22 experiment liquid
0.03
0.025
I~YI......"(b) 
~'1
A'Lu~

Paper No
The volume fraction distribution of large bubbles, as well
as the total void fraction (see Fig. 11 (d)), agree reasonably
well with the data. The agreement between the predicted
velocity profiles and the experimental measurements in
Figure 12 is also quite good.
Figures 13 and 14 show again the results for TOPFLOW
Run 107, but this time the experimental date at L/D=7.7
have been used at the inlet to the computational domain.
As before, the results of predictions have been directly
compared against the experimental data. Both the liquid
and gas radial velocity profiles along the initial section of
pipe, i.e. between L/D=7.7 and L/D=13, reach maximum
values relatively close to the wall and then decrease toward
the center of the pipe. Such velocity distributions affect
the predicted concentration distribution of large bubbles at
L/D=13, slowing down the transition to a centerpeaked
profile. The calculated concentration distributions for small
and intermediatesize bubbles, shown in Figure 14, agree
quite well with the data, although the center line
concentration of large bubble at L/D=13 is underpredicted.
In general, the results of NPHASECMFD predictions
agree with the experimental trends along the flow, and they
show good agreement at locations which are not too close
to the gas inlet zone. The effect of geometry of the gas
injection region on flow conditions and bubble distribution
there, augmented by the experimental uncertainties,
contribute to the modeling difficulties at distances close to
this region.
Conclusions
The present work has achieved several objectives. One of
them was to demonstrate the importance of proper physical
closure modeling of twophase flows in vertical conduits.
In particular, selected computational and modeling issues
have been investigated and resolved, associated with
multidimensional simulations of multiphase flows using a
multifield ensembleaveraged modeling framework. The
overall model has been implemented in the
NPHASECMFD solver and parametrically tested. The
other major objective of this work was to demonstrate the
ability of the NPHASECMFD code to predict the
evolution of adiabatic churturbulent gas/liquid flows in
vertical channels. The results of NPHASECMFDbased
computer simulations confirm both the modeling and
computational consistencies. The results of calculations
have been compared with the several experimental data
sets from the TOPFLOW test facility, and a good
agreement has been observed.
Needless to say, several unresolved modeling and
computational issues can still be identified, which require
further investigation and additional future work. Among
those are: the effect of interfacial forces on large deformed
bubbles, the mass transfer model between bubbles of
different sizes, and the mechanisms governing phase
distribution evolution in developing flows.
Acknowledgments
Funding for this research was supported by the U.S.
Department of Energy, Office of Nuclear Energy, under
DOE Idaho Operations Office Contract
DEAC0705ID14517. The experimental data were
obtained at the Institute of Safety Research of
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Forschungszentrum Rossendorf e.V in the frame of
research project funded by the German Federal Ministry of
Economics and Labor.
References
Antal, S.P, Ettorte, S.M., Kunz, R.F. and Podowski M.Z.,
Development of a next generation computer code for the
prediction of multicomponent multiphase flows, Proc. Int.
Conf. Trends in Num. and Phys. Modeling of Multiphase
Flow, (2000)
Antal S.P, Podowski, M.Z., Lahey, R.T. Jr., Barber, D. and
Delfino, C., Multidimensional Modeling of Developing
TwoPhase Flows in a Large Adiabatic Riser Channel,
Proceedings NURETH11, Avignon, France, (2005)
Drew, D.A., Analytical modeling of multiphase flows, In
Boiling Heat Transfer (Ed.: Lahey, R.T., Jr.), Elsevier, New
York, (1992)
Ervin, E.A. and Tryggvason, G., The rise of bubbles in a
vertical shear flow, Journal of Fluids Engineering, Vol. 119,
(1997)
Kumbaro, A. and Podowski, M.Z., The effect of bubble/
bubble interactions on local void distribution in twophase
flows, Proceedings 13th International Heat Transfer
Conference, Sydney, Australia, 2006.
Kurul, N. and Podowski, M.Z., Grid Generation for the
Analysis of Dispersed Phase Motion in TwoPhase Flows,
in Numerical Grid Generation in Computational Fluid
Mechanics, Pineridge Press, (1988)
Lehr, F, Millies, M. and Mewes, D., Bubble Size
Distributions and Flow Fields in Bubble Columns, AIChE
Journal, Vol.48, (2002)
Luo, F and Svendsen, H.F, Theoretical Model for Drop
and Bubble Breakup in Turbulent Dispersions, AIChE
Journal, Vol. 42, (1996)
Podowski, M.Z., On the consistency of mechanistic
multidimensional modeling of gas/liquid twophase flows,
Nuclear Engineering and Design, Vol.239, 933940,
(2009)
Prasser, H.M., Beyer, M., Carl, H., Manera, A., Pietruske,
H. and Schutz, P., Experiments on upwards gas/liquid flow
in vertical pipes, FZD482, (2007).
Prince, M.J. and Blanch, H.W., Bubble Coalescence and
Breakup in AirSparged Bubble Columns, AIChE Journal,
Vol.36, (1990)
Sato, Y. and Sekoguchi, K., Liquid velocity distribution in
twophase bubble flow, Int. J. Multiphase Flow, Vol.2,
(1975)
Tomiyama, A., Struggle with Computational Bubble
Dynamics, Proceedings International Conference on
Multiphase Flow, Lyon, France, (1998)
Wallis, G.B., One Dimensional Two Phase Flow, New
York: Mc GrawHill, (1969)
Zun, I., The transverse migration of bubbles influenced by
walls in vertical bubbly flow, International Journal of
Multiphase Flow, Vol.6, (1980)