7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
A Comprehensive Mlechanistic Mlodel for the Prediction of
Vertical Subcooled Boiling Flows
G.H. Yeoha~b, S.C.P. Cheunge, J.Y. TuC and M.K.M Hoa
"Australian Nuclear Science Technology Organisation (ANSTO), PMB 1, Menai, NSW 2234, Australia
bSchool of Mechanical and Manufacturing Engineering, University of NSW, NSW 2052, Australia
"School of Aerospace, Mechanical and Manufacturing Engineering, RMIT University, Victoria 3083, Australia
guan.yeoh@anto.gov. au, chipok.cheung@rmit. edu.au, j iyuan. tu@rmit. edu.au and mark, ho @ansto.gov. au
Keywords: Subcooled boiling flow, fractal model, wall heat flux partitioning, twofluid model, population balance
Abstract
Based on our previous investigation in Yeoh et al. (2008), an improved wall heat flux partitioning model at the heated surface has been
developed based on the consideration of: (i) accounting sliding as well as stationary bubbles, (ii) bubble frequency which accounts for the
complete ebullition description of the bubble growth, sliding and liftoff and (iii) a force balance model acting on departing bubbles. This
model coupled with a threedimensional twofluid model and population balance equations based on the class method MUltiple SIze Group
(MUSIG) model have led to satisfactory agreement being achieved between the model predictions and experimental measurements.
Nevertheless, one of the shortcomings is the reliance of empirical correlations for the active nucleation site density in the wall heat flux
partitioning model. This discrepancy brings about the uncertainties within the heat flux partitioning model in appropriately evaluating the
vapor generation rate, which subsequently influences the model prediction particularly on the axial and radial void fraction profiles within the
bulk fluid flow. Combining the fractal model with the aforementioned subcooled boiling flow model in the absence of empirical correlations
for the active nucleation site density, a comprehensive mechanistic model to predict vertically oriented subcooled boiling flows is developed.
The proposed model is assessed against the experimental data of radial measurements of Yun et al. (1997) and Lee et al. (2002) and axial
measurements of Zeitoun and Shoukri (1996) for upward subcooled boiling flows within annular channels. Improved model predictions are
obtained. Discussions on the agreement of other twophase flow parameters are also presented.
Introduction
Subcooled boiling flow is widely applied in many industrial
processes. It is inherently complex, which still makes the
analysis of such flow extremely difficult. Of particular
concern is the development of suitable models for the
prediction of lowpressure subcooled boiling flow which
remains as challenging as ever.
Modeling subcooled boiling flow at low pressures can
be broadly classified according to two categories: (i) Heat
transfer and wall heat flux partitioning during subcooled
boiling flow at the heated wall and (ii) Twophase flow and
bubble behaviors in the bulk subcooled flow away from the
heated wall. For the second category, it has been found that
the use of twofluid model coupled with the population
balance approach the MUSIG model can aptly predict the
local and axial distributions of the bubble size, void fraction,
interfacial area concentration and velocity profiles for
upward vertical bubbly subcooled boiling flows (Yeoh & Tu,
2004, 2005a).
The first category generally involves the development
of models to predict the heat transfer rate during subcooled
nucleate boiling flow by appropriately partitioning the wall
heat flux. Mechanistic models based on relevant heat transfer
mechanisms occurring during the boiling process can be
readily employed for both the prediction of the wall heat flux
as well as the partitioning of the wall heat flux between the
liquid and vapor phases. Most computational fluid dynamics
(CFD) investigations on subcooled boiling flow are mainly
based on the use of such models. The wall heat flux can be
partitioned into the respective heat flux components:
singlephase, transient conduction and evaporation.
Singlephase heat transfer is assumed to occur in areas of the
heated surface unaffected by the bubbles. Transient
conduction (quenching heat flux component) occurs over the
area of the heater surface under the influence of bubbles. The
model consists of three unknown parameters, which should
be modeled accurately. They are: active nucleation site
density (N,), departing bubble diameter (Db) and bubble
frequency (f).
In order to remove the application uncertainty of
correlations empirically determined for the departing bubble
diameter, a force balance model has been developed in Yeoh
& Tu (2005b). It has been demonstrated that the bubble
growth and its maximum detached bubble size contributed
significantly towards the void growth and heat transfer
within the bulk liquid. Also, in order to eliminate the
uncertainty of evaluating the bubble frequency, it is proposed
that the fundamental theory based on the description of an
ebullition cycle in nucleate boiling, which are: (1) waiting
period tz, (transient conduction of heat to liquid); (2) growth
period t,: (a) bubble growth rate, (b) evaporation process, (c)
agitation of liquid around the bubble and (d) termination of
bubble whether by departure or collapse, is employed instead,
viz.,
1
/ t+
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Fdu
F,
Fs
F,
F,
F,'''
F,
Turbulent chspersion
F;,
g
G
Gs
h
H
k
K
hfg
1,
n1,
N,
N,~,to
Qw,
Qc
Qe
Or,
Qtcst
r
uHSteady drag force due to asymmetrical
growth of the bubble (N)
surface tension force (N)
shear lift force (N)
forces along the xdirection (N)
forces along the ydirection (N)
8CliOH Of interf8Ci81 foTCCS from vapor on
liquid (N)
action of interfacial forces from liquid on
vapor (N)
drag force (N)
lift force (N)
wall lubrication force (N)
turbulent dispersion force (N)
gravitational constant (m s )
mass flux (kg m"2 s ')
dimensionless shear rate
interfacial heat transfer coefficient
enthalpy (J kg )
Jakob number
thermal conductivity (W m2 K ) or
turbulent kinetic energy (m s 2)
proj ected area of bubble
latent heat of vaporization (J kg ')
sliding distance (m)
number density of the discrete bubble ith
class (m )
active nucleation site density (m 2)
total number of nucleation sites per unit
area (m )>
wall heat flux (W m )
heat transfer due to forced convection (W
m )
heat transfer due to evaporation (W m )
heat transfer (transient conduction) due to
stationary bubble (W m )
heat transfer (transient conduction) due to
sliding bubble (W m )
bubble radius at heated wall (m) or flow
spacing within annular channel (m)
curvature radius of the bubble at heated
surface (m)
cavity radius at heated surface (m)
bubble Reynolds number
ratio of the actual number of bubbles
lifting off to the number of active
nucleation sites
radius of inner heated wall (m)
radius of outer unheated wall (m)
spacing between nucleation sites (m)
additional source terms due to coalescence
and breakage (kg m" s ')
Stanton number
time (s)
bubble growth period (s)
bubble sliding period (s)
bubble waiting period (s)
temperature (K)
difference in temperature (K)
pressure (N m )
Velocity (m s )
It has been indicated in Basu et al. (' l l 2005b) and
Sateesh et al. (2005) that transient conduction due to sliding
bubbles can become the dominant mode of heat transfer. In
Yeoh et al. (2I III, part of the improvements made to the wall
heat partition model has been to incorporate the sliding
period as well as mechanistically determined the bubble
frequency through equation (1). Nevertheless, the wall heat
partition model remains lacking especially in determining the
active nucleation site density which still requires the use of
empirical relationships.
In the quest of reducing the application of empirical
correlations for the active nucleation site density, Xiao and
Yu (2 II) have proposed a fractal model for subcooled
boiling flow. Essentially, the model considers the fractal
distribution of sizes of active cavities on the heated surface.
This mechanistically derived model results in a dependence
on the wall superheat, liquid subcooling, bulk velocity of
fluid, fractal dimension, minimum and maximum active
cavity sizes contact and physical properties of fluid.
The objectives of this present study are thus: (i) to
formulate and incorporate a fractal model of determining the
active nucleation site density making the wall heat partition
model comprehensively mechanistic within the framework
of twofluid and modified MUSIG models for lowpressure
subcooled boiling flows in the generic computational fluid
dynamics (CFD) computer code CFX4.4; and (ii) to evaluate
the improved wall heat partition model through validation
against experimental measurements. Comparisons of
computational predictions for a range of different mass and
heat fluxes and inlet subcoolings are performed against axial
measurements of Zeitoun & Shoukri (1996) and local
measurements of Yun et al. (1997) and Lee et al. (I II == '
Nomenclature
alf interfacial area concentration (m )
4, fraction of heater area occupied by bubbles
C1, C2 constants defined in equation (10)
CD drag coefficient
CL shear lift coefficient
C, specific heat (J kg' K')
C, constant defined in equation (18)
C, acceleration coefficient
d vapor bubble diameter at heated surface
(m)
df fractal dimension
d,, surface/bubble contact diameter (m)
D average bubble diameter (m)
Db departing bubble diameter (m)
D, active cavity diameter (m)
Dc,n,,, minimum diameter of active cavity (m)
Dc,n,, maximum diameter of active cavity (m)
Dd bubble departure diameter (m)
Di bubble liftoff diameter (m)
D, bubble Sauter diameter (m)
f bubble frequency (Hz)
Af scalar fraction related to the number
density of the discrete bubble classes
F degree of surface cavity flooding
Fb buoyancy force (N)
Fi, force due to the hydrodynamic pressure
(N)
friction velocity (m s )
specific volume of discrete bubble ith class
(m3 kg')
Cartesian coordinate along x
nondimensional normal distance from
heated wall
Cartesian coordinate along y
advancing angle (rad)
vapor void fraction
liquid void fraction
receding angle (rad)
Thermal boundary layer thickness (m)
turbulent dissipation rate (m s )
thermal diffusivity (m s ')
effective viscosity (Pas)
viscosity (Pas)
bubble contact angle (rad)
inclination angle (rad)
Density (kg m')
surface tension (N m )
interfacial mass transfer from vapor to
liquid (kg m' s )
interfacial mass transfer from liquid to
vapor (kg m s')
Greek letters
a
ag
at
ri
ine
1,,
sa
Subsit
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
nucleation sites is derived from Yu and Cheng ('I III) as
2 Dr.
Dc~r
The term Dr.,, in the above equation represents the
averaged value over all maximum active cavities which can
be determined according to Xiao and Yu (2I II) as
1
Dr,,, r.,
where n =(Tw T,, ()/STw and &T, is designated a constant
with a value of 0. 10C in which T,,., is calculated according to
Tw~ = T,r +i(STw,) with i = 1, 2,..., n .
The evaluation of minimum and maximum active cavity
diameters are given by Hsu (1962) for nucleation site
distribution, vi .,
3, (T,,. Ter)
D = 2
axial distribution
vapor
channel entrance
liquid
local distribution
surface heater
saturation
subcooled
3,(T, Ea)(Tw Ea,) 4 (T ,
0 (T,, T,)(Tw T,) S(T,, T,)
(7)
where T,, is the temperature of the heater surface, Ti is the
temperature of the liquid (~= 2aTsar/Pg hg
C= (1+cos 9)lsin I and C=1+ cos 0. The angle a
represents the bubble contact angle while the thickness of the
thermal boundary layer 4 can usually be expressed as
Mathematical Formulation
Fractal Characteristics of Nucleation Sites on
Boiling Surfaces
Based on Xiao and Yu (2 II), the active cavities that are
formed on the heated surface are analogous to pores in
porous media. Thus, the cumulative number of active cavities
with diameters greater than and equal to particular active
cavity diameter, D,, can be described by
k,=
In the above equation, the singlephase heat transfer
coefficient for forced convection hi can be evaluated from the
local Stanton number St for turbulent convection.
Fractal Wall Heat Partition Model Incorporating
Bubble
A fractal model for the wall heat partition model is derived
based on the idea that the nucleation site size distribution
follows the fractal power law as described by equation (2).
Differentiating equation (2), the number of active cavities of
sizes lying between D, and D, + do, can be obtained as
D < D
cann c c~mm
where Dc,,c, and D, ,,,, are respectively the maximum and
minimum diameters of active cavity. The total number of
nucleation sites per unit area from the minimum to maximum
active cavity can thus be obtained from equation (2) as
dNa = d Dc D dDe
For a boiling system, the fractal dimension df of
(7, E,)45(~
ct. )
cr.
D
ator D
<,,,,,,
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
where do, > 0 and dN, > 0. A fractal model fo1
enhancement in heat transfer during forced convel
boiling can be attributed by the presence of both sliding
stationary bubbles can be derived as follows. There
essentially two mechanisms: (1) the latent heat transfer d
microlayer evaporation and (2) transient conduction a
disrupted thermal boundary layer reforms during the wa
period (i.e. incipience of the next bubble at the a
nucleation site).
Transient conduction occurs in regions at the poil
inception and in regions being swept by sliding bubbles. 1
stationary bubble, a fractal model for the heat flux fron
minimum site Dc~a to the maximum site Dc~a is given
Q~c = T 2 pl (T T) oK t f(da
ir(t, +tst) 4
J lji(t,+tst) 4
where Dd is the bubble departure diameter. The a
equation indicates that some fraction of the nucleation
will undergo transient conduction while the remaining w
in the growth period. For a sliding bubble, a fractal mode
the heat flux due to transient conduction that takes I
during the sliding phase and the area occupied by the sli
bubble at any instant of time from the minimum site Dc,>
the maximum site Dc,ma is given by
Qtcsr = 2 kPC (T Ti)R 1,Kotf(d, fa)
D,msn ti (w sl)
(Tw T,) x
Tb Tw
where x is the normal distance from the wall. Based on the
criterion of the incipience of boiling from a bubble site inside
the thermal boundary layer, the bubble internal temperature
for a nucleus site (cavity) with radius R, (= D, / 2) is
Tb = Tsa + 2as'at x = CR
b t Rehf,p
where C2 = /Sin0 By substituting equation (14) into
equation (13), the waiting time tw can be obtained as
(10) tw = 1 (Tw T,) C~r
ira (T T,,,) 2oTs,/lC2 g hr,,r
The reduction factor Rf appearing in equations (10) and
(11) depicts the ratio of the actual number of bubbles lifting
off per unit area of the heater surface to the number of active
nucleation sites per unit area, viz., Rf = 1/(1,/s) where lis i
the sliding distance and s is the spacing between nucleation
sites. In the present study, it shall be assumed that the
nucleation sites are distributed in a square grid and that the
bubbles slide only in the direction of the fluid flow (Basu et
al. 2005). The spacing between nucleation sites can thus be
approximated as s = 1/ The factor R, is obtained
alongside with the sliding distance evaluated from the force
balance model (to be described later). The significance of this
) factor provides the information whereby the bubble
departing from its site of origin merges with other nucleating
bubbles at adjacent sites. It is noted that for the case where
the sliding distance 1, is less than the spacing s, Rf = 1
The heat flux attributed to vapor generation is given by
the energy carried away by the bubbles lifting off from the
heated surface. It also represents the energy of vaporization
whereby the bubble size of DIis produced. Afractal model for
this heat flux from the minimum site D mnto the maximum
site Dc,ma is can be expressed by c
(11
dN,)
where the average bubble diameter D is given by
D = (Dd + D, )/2 and DI is the bubble liftoff diameter,
The bubble leaving the heater surface draws in liquid
from an area called the area of influence that is some factor K
times the projected area of the bubble. As suggested by
Kenning and Victor (1981), a value of K= 4 has been applied.
Equations (10) and (11) require the evaluation of the
bubble frequency fand waiting time t,. The bubble frequency
is giving by equation (1) while the bubble waiting time can be
obtained based on the analytical expression derived below.
Assuming that the heat capacity of the heater wall
p,C,,S, is very small, the conduction process can be modeled
by considering onedimensional transient heat conduction
into a semiinfinite medium with the liquid at a temperature Ti
and the heater surface at a temperature T,. The wall heat flux
can be approximated by
Based on the characteristics of fractal media, an
expression which relates the pore volume fraction to fractal
dimension, minimum and maximum pore size in porous
media:
k,(T, T,)
e,=
~J~J~J~J~J~J~J~J~J~J~J~J~J~J~J~J~~
~iat
where d = 2 in a twodimensional space. A fractal model for
the forced convection, which will always prevail at all times
in areas of the heater surface that are not influenced by the
stationary and sliding bubbles, can be obtained as
If the temperature profile inside this layer is taken to be linear
(according to Hsu & Graham, 1976), it can thus be expressed
as
+ 2 k;rP IC;(T T)R fsF 2 w
Q, = ..RISm f io : p~h (d~)
:dd,
D
where ut is the adj acent liquid velocity.
The total ivall heat flux Q,, is thus the combination of the
following heat flux components: Q, .= Qe + Qt + Qrei+ Q,. It
can be seen from equations (10), (11) and (16) that the bubble
frequency and waiting time are related to the sizes of active
cavities. Numerical integration via the Simpson's rule is
performed to obtain the respective heat flux components
during the calculations.
Bubble LiftOff and Bubble Departure Model
The development of the force balance model concentrates on
the various forces that influence the growth of a bubble
during flow conditions in the directions parallel and normal to
a vertical heating surface. These forces are formulated
according to the studies performed by Klausner et al. (1993)
and Zeng et al. (1993). Figure 1 illustrates the forces acting on
the bubble in the xdirection and vdirection; they are
respectively,
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
e, = StpICpilll (1 K/) (Tw, Ti)
Heated
Wall
IFx = F, + Fdzw + FsL + Fi, + Fep
, I~Lift Off
I.~ Sliding
Bubble Growth
Sliding
Distance
IF, = F,. + Fdlit + Fqs + Fb
where F, is the surface tension force, Fdit iS the unsteady drag
due to asymmetrical growth of the bubble and the dynamic
effect of the unsteady liquid such as the history force and the
added mass force, FsL is the shear lift force, Fi, is the force due
to the hydrodynamic pressure, Fep iS the contact pressure
force accounting for the bubble being in contact with a solid
rather than being surrounded by liquid, F,, is the quasi
steadydrag in the flow direction, and Fb iS the buoyancy
force. In addition, g indicates the gravitational acceleration; a,
p and 8, are the advancing, receding and inclination angles
respectively; d,,.is the surface/bubble contact diameter; and d
is the vapor bubble diameter at the wall.
The forces acting in the xdirection can be estimated
from:
Fsx = d oG[cos p cos a]; F U= P,, cosO,;
e. P
Figure 1: Schematic drawings illustrating the forces acting
on a growing vapor bubble (above) and a bubble departing,
sliding and lifting off from a vertical heated surface (below).
The drag coefficient CD and shear lift coefficient CL
appearing in the drag and lift forces are determined according
to the relationships proposed by Klausner et al. (1993), viz.,
1 9 rd
', = C,p,AC ~irr; F, = piA ; F
2 4 4
ird2 2&
4 r,
In the ydirection, they are:
Fsy = dwo. [sina+sinp]; F,,,.
(aP)
F, sin8, ;
P', = 6CD / r~~r F, = r(P Pg
From the various forces described along the xdirection and
vdirection, r is the bubble radius, AC is the relative velocity
between the bubble centre of mass and liquid, CD and CL are
the respective drag and shear lift coefficients and r, is the
curvature radius of the bubble at the reference point on the
surface x = 0, which is r, ~ 5r (Klausner et al., 1993).
where n = 0.65 and Re = piUld/p, is the bubble Reynolds
number. The dimensionless shear rate G, is:
(dP/dx)(r/AC) The gradient dP/dx can be determined
through the universal velocity profile for turbulent flow:
v,
2 12
CD = ++076
3 Re
C, =~n' 3.87G +0014'
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
departure or liftoff. Here, a correlation based on the
experimental data of Maity (2I III) as a function of the bubble
contact angle a is employed to determine the ratio of the
bubble base diameter d, to the bubble diameter at departure
or liftoff, which is given as C= 1exp(2006)
Experimental observations by Klausner et al. (1993) and
Bibeau and Salcudean (1994) have indicated that the
advancing angle a and receding angle P varied quite
substantially during the sliding phase. Considering the
complexity of the bubble departure and bubble liftoff, and
the difficulty in obtaining the measurements, the advancing
and receding angles can be reasonably evaluated through the
bubble contact angle 8as a= 8+ 8' andp= 8 8'. Klausner
et al. (1993) have employed an angle of 4.50 in their
theoretical analysis while Bibeau and Salcudean (1994) have
reported a value of 2.50. According to Winterton (1984), the
angle 8' has nonetheless been postulated to be as high as 100.
In the present study, an angle of 50 is adopted. For the
inclination angle a ,, a value of 100 that gave the best fit to the
data by Klausner et al. (1993) is employed for the current
force balance model. Based on the experimental observations
of Zeitoun and Shoukri (1996) and Lee et al. :L2 2I), the
bubble contact angles have been taken to be at 350 and 450 for
the present investigation.
TwoFluid Model
The twofluid model based on the EulerianEulerian
framework is adopted in the present study of which the
governing equations are based on the ensembleaveraged of
mass, momentum and energy transport equations for each
phase. With the liquid phase (at) taken as the continuum
phase while the vapor (bubbles) as disperse phase (ag), these
equations can be written as:
Continuity Equation ofLiquid Phase
Continuity Equation of Vapour Phase
= 2.51n(9.8x ) (23)
where u, is the friction velocity and x' = plu,x/4u is the
nondimensional normal distance from the heated wall. The
velocity profile in (23) is assumed to be applicable for the
timeaveraged velocity distribution in the vicinity of the
heated wall. Adjacent velocities, determined through the
twofluid model, are used to obtain the varying local friction
velocities through equation (23). These friction velocities are
subsequently used to evaluate the gradients dU/dx along
the heated wall to determine the shear rate G,.
The growth force Fdu iS modeled by considering a
hemispherical bubble expanding in an inviscid liquid, which
is given by Zeng et al. (1993) as
F, = r2 s2 r (24)
where ( ) indicates differentiation with respect to time. The
constant C, is taken to be 20/3 according to Zeng et al.
(1993). In estimating the growth force, additional information
on the bubble growth rate is required. As in Zeng et al. (1993),
a diffusion controlled bubble growth solution by Zuber
(1961) is adopted:
rT 2b o~;a IO;CplTsa, k~
r () = Ja ; J =; p= (25)
where Ja is the Jakob number, rl is the liquid thermal
diffusivity and b is an empirical constant that is intended to
account for the asphericity of the bubble. For the range of
heat fluxes investigated in this investigation, b is taken to be
0.21 based on a similar subcooled boiling study performed by
Helfried et al. (2005), which has been experimentally verified
through their inhouse measurements with water as the
working fluid.
While a vapor bubble remains attached to the heated wall,
the sum of the parallel and normal forces must satisfy the
following conditions: IFx = 0 and IF, = 0. For a sliding
bubble case, the former establishes the bubble departure
diameter (Dd) while the latter yields the bubble liftoff
diameter (DI). The growth period tz appearing in equation (1)
can be readily evaluated based on the availability of the
bubble size at departure from its nucleation site through
equation (25). The lift off period tl can also be similarly
calculated based on the bubble liftoff diameter. The
difference between the bubble liftoff and bubble growth
periods provides the period for the sliding bubble; the sliding
distance 1, can subsequently be determined (see Figure 1). An
estimation on this sliding distance can be determined
according to the experimental correlation of Maity (2I III) as
1, = (2/3) C~t 2 Where ts is the sliding time (tt tg) and C, is
an acceleration coefficient correlated in terms of the
tangential liquid velocity (ut) adjacent to the heated surface:
Cy = 3.2ut +1. This coefficient reflects the increase in bubble
velocity with time after it begins to slide away from a
nucleation site.
In reality, the surface/bubble contact diameter d,
evolves from the point of inception until the point of
ap,a
at .(pta~zil) = Flg
at
@ap, J+V p~zl J S
M2lomentum Equations of Liquid and Vapour Phases
atm+ 9. (psa~i m)= am? + ampm
+Va,,,p Vl,, +(Vl,,)lj]+(FmnimEnds)+Fm
Liquid Phase : m=1,n= g
Vapor Phase: m=g,n=1
Enerev Equations of Liquid and Vapour Phases
For the local measurements performed by Yun et al. (1997)
and Lee et al. :al III), the experimental setup consisted of a
vertical concentric annulus with an inner heating rod of 19
mm outer diameter. The heated section was a 1.67 m long
Inconel 625 tube with 1.5 mm wall thickness and filled with
magnesium oxide powder insulation. The rod was uniformly
heated by a 54 kW DC power supply. The outer wall
comprised of two stainless steel tubes with 37.5 mm inner
diameter. The plane for measuring the radial distribution was
located at 1.61 m downstream of the beginning of the heated
section. Demineralised water was used as the working fluid.
Local gas phase parameters such as local void fraction,
bubble frequency and bubble velocity were measured by a
twoconductivity probe method. The bubble Sauter diameters
(assuming spherical bubbles) were determined through the
IAC, calculated using the measured bubble velocity spectrum
and bubble frequency.
For the axial measurements performed by Zeitoun and
Shoukri (1996), the test section was a vertical concentric
annular test section. The inner tube, which had a 12.7 mm
outside diameter, was a 30.6 cm long, thickwalled
stainlesssteel tube (0.25 mm thick) that was electrically
heated. The entire inner tube was connected to a 55 kW DC
pOwer supply. The outer tube was a 25.4 mm inner diameter
plexiglass tube that permitted visual observation.
Distilleddegassed water was used as the working fluid. A
digital image processing technique was used to analyze the
highspeed video information and to measure bubble size
distributions along the subcooled boiling region. A single
beam gamma densitometer was used for the void fraction
measurements. Experimental conditions for the local and
axial data that have been used for comparison with the
simulated results are presented in Table 1.
Table 1: Experimental conditions for local (L1, L2, L3)
and axial (Al, A2, A3) cases.
Pmler T,niet Tsub at. G
Case [MPa] [oC] (inlet) [kW/m ] [kg/m s]
[oC]
L1 0. 143 96.9 13.4 152.9 474.0
L2 0. 137 94.9 13.8 197.2 714.4
L3 0. 143 92.1 17.9 251.5 1059.2
Al 0. 137 91.9 14.9 286.7 156.2
A2 0. 150 94.6 16.6 508.0 264.3
A3 0. 150 88.9 22.5 705.0 411.7
Numerical Details
The conservation equations for mass, momentum and energy
of each phase were discretised using the finite volume
technique. A total number 15 bubble classes were prescribed
for the dispersed phases, which are illustrated in Table 2 for
the local and axial cases. This represents an additional set of
15 transport equations of which they were progressively
solved and coupled with the flow equations during the
simulations.
The advantage of an annular geometrical shape was
utilized by modeling only one quarter of the annulus as the
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
bubbleinduced turbulent viscosity.
ap,,a,,H
+(T,,,H,,r,H,,,H,
Liquid Phase :m=1,n=g
Vapor Phase :m=g,n=1 2)
On the right hand side of equation (27), S, represents the
additional source terms due to coalescence and breakage for
the range of bubble classes that can exist within the vapor
phase. Here, the coalescence rate of Prince and Blanch (1990)
and breakage rate of Luo and Svendsen (1996) are applied. For
bubble coalescence, turbulent collision in the inertial subrange
of isotropic turbulence is considered while for bubble breakage,
the behavior of the socalled daughter bubble size
fragmentation caused by the interaction of the bubbles with
turbulent eddies is adopted. Detailed expressions of these rates
can be found in Yeoh and Tu (2005a). The total interfacial
force F,,,, appearing in equations (28) are formulated through
appropriate consideration of different subforces affecting the
interface between each phase. For the liquid phase, the
interfacial forces Fig ar
F = "'"# +Fif "+F~arallhbruation+F""swerbuent ,,wm (30)
The terms appearing on the right hand side of equation (30)
are: drag force, lift force, wall lubrication force and turbulent
dispersion force. More detail descriptions of these interfacial
forces can be found in Anglart and Nylund (1996). Note that
for the gas phase, Fgi = Fig. The interfacial mass transfer rate
due to condensation in the bulk subcooled liquid in equations
(29) can be expressed by the following expression:
r, = ha,, (Tsar fT )lh, (31)
Here, h indicates the interphase heat transfer coefficient. The
wall vapor generation rate is modeled in a mechanistic
manner by considering the total mass of bubbles detaching
from the heated surface as:
gl fg at /hf + Cpl(T,, T,)
For subcooled boiling flows, the wall nucleation rate is
treated as a specified boundary condition for equation (29).
The boundary mass flux is achieved by apportioning the
generation rate based on the size of bubbles at liftoff on the
heated surface. The term frig in equatiOn (21) represents the
mass transfer due to condensation redistributed for each of
the discrete bubble classes within the MUSIG model. The gas
void fraction ag along with the scalar fraction/f are related to
the number density of the discrete bubble ith class n,
(similarly to the jth class n,) and specific volume v, as ag/ =
n,v,. A kE turbulence model is employed for the continuous
liquid and dispersed gaseous phases. The effective viscosity
1,", in the momentum and energy equations is taken as the
sum of the molecular viscosity and turbulent viscosity. The
turbulent viscosity is considered as the total of the
shearinduced turbulent viscosity and Sato et al. (1981)
Experimental Facility
Bubble class diocal (mm) d,,ai (mm)
1 0.45 0.34
2 0.94 0.70
3 1.47 1.10
4 2 02 1.51
5 2.58 1.93
6 3.14 235
7 3.71 277
8 4.27 3.19
9 4.83 3.61
10 5.40 4.03
11 5.96 4.46
12 6.53 4.88
13 7.10 5.30
14 7.66 5.72
15 8.23 6.15
Results and Discussion
The improved wall heat partition model is assessed against
subcooled flow boiling of local radial and axial experimental
data covering a wide range of different mass and heat fluxes,
For the local cases of L1, L2 and L3, the measured and
predicted radial profiles of the bubble Sauter diameter (D,)'
vapour void fraction and interfacial area concentration (alf=
6 a/D,) located at the measuring plane 1 .61 m downstream of
the beginning of the heated section are shown in Figures 24.
Note that in all the figures, the dimensionless parameter
(rR,)/(RoR,) = 1 indicates the inner surface of the unheated
flow channel wall while (rR,)/(RoR,) = 0 indicates the
surface of the heating rod in the channel.
Peak local void fractions observed near the heated
surface are generally typical of subcooled boiling flows. The
presence of larger bubbles away from the heated wall
confirms the occurrence of bubble coalescence. Further away
from the heated wall, the reduction of the bubble sizes
indicates the bubbles condensing due to the subcooled bulk
liquid. As seen from Figures 24, good agreement was
achieved between the predicted parameters and experimental
data. At the measuring location, the bubble departure
diameters for L1, L2 and L3 were predicted to be about 0.58
mm, 0.65 mm and 0.68 mm while the bubble liftoff
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
diameters were approximately attained as 1.45 mm 1.32 mm
and 1.21 mm respectively. From these bubble sizes, the ratio
of the bubble liftoff diameter to the bubble departure
diameter achieved an approximate value of two, which
showed good agreement against experimental observations
of Basu et al. (2005).
At the point of bubble lifting off into the bulk liquid, the
shear lift force was seen as the dominant force in overcoming
the surface tension forces acting on the vapor bubble. The
calculated shear lift and surface tension forces are respectively
for Ll: F, ~ 5.06 x 10' N and F, ~ 1.78 x 104 N; L2: F, ~
5.06 x 10' N and F, ~ 1.80 x 10 N; and L3: F, ~ 5.06 x
10' N and F, ~ 1.79 x 104 N.
For the axial Case 3 and 4, the measured and predicted
profiles of the void fraction, bubble Sauter diameter
nrormalized by thre length scale o/grp, p and
Interfacial area concentration along the heated section are
presented in Figures 57.
The increasing bubble sizes along the axial length clearly
showed the bubble coalescence intensifying as the bubbles
traveled downstream. Predicted bubble Sauter diameter
profiles are seen to agree well against experimental data. The
void fraction clearly showed a typical plateau representing
initial stages of boiling being established near the channel
entrance. Meanwhile, a sharp increase of the void fraction
profiles existed towards the channel exit. In this region,
reduction in the relative importance of condensation at the
bubble interface due to decreasing liquid subcooling
encouraged the likelihood of substantial bubble coalescence.
Conclusions
Formulation of a fractal model for the wall heat flux
partitioning model that accounted stationary and sliding
bubbles alongside the calculation of the bubble frequency
through the fundamental consideration of bubble frequency
theory during lowpressure subcooled flow boiling has been
presented. This model when coupled with a force balance
model to predict forces acting on a vapour bubble growing
under subcooled flow boiling demonstrated the capability to
accommodate more complex analyses of bubble growth,
bubble departure and bubble liftoff over a wide range of wall
heat fluxes and flow conditions. The models were assessed
against local and axial subcooled boiling flow measurements.
Comparison of the predicted results against local
measurements of Yun et al. (1997) and Lee et al. (~II)
showed good agreement for the local radial profiles of the
void fraction, bubble Sauter diameter and interfacial area
concentration. Predicted growth and waiting times of the
bubble frequency at particular local wall superheat and
subcooling temperatures were found to be in reasonable
agreement with experimental data of Basu et al. (2005).
Comparison against measurements of Zeitoun and Shoukri
(1996) also revealed good agreement when comparing the
axial profiles of the void fraction, bubble Sauter diameter and
interfacial area concentration. The model clearly represented
a plateau occurring at the initial boiling stages, typical of
subcooled flow boiling at low pressures, followed by the
significant rise of the void fraction downstream of the
annular channel.
domain for both the local and axial cases since uniform wall
heat flux was applied. A bodyfitted conformal system was
employed to generate the threedimensional mesh within the
annular channel resulting in a total of 13 (radial) x 30 (height)
x 3 (circumference) control volumes for the local case while a
total of 8 (radial) x 20 (height) x 3 (circumference) control
volumes were used for the axial case. Since wall function was
used in the present study to bridge the wall and the fully
turbulent region away from heater surface, the normal
distance between the wall and the first node in the bulk liquid
should be such that the corresponding x+ was greater than 30.
Grid independence was examined. In the mean parameters
considered, further grid refinement did not reveal significant
changes to the twophase flow parameters.
Table 2: Evaluated diameters of
classes for the local and axial cases.
each discrete bubble
0 0.5 1
(a)
Case L2
e Measurement
Prediction
1 0.5 1
(b)
Case L3
Measurement
Prediction
****,
Case L3
 e eMeasurement
Prediction
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
0.5
0.4
S0.3
0.1
Case L1
*Measurement
. Prediction
*
0.01
E
S0.0075
0.005
i
22 0.0025
1 0
Case L1
e M~easur"emet
 Pedcton
b I
0.5
(a)
0.5
0.4
t; 0.3
S0.2
0.1
0
0.01
S0.0075
0.005
a,0.0025
0
0.5 1
(b)
0.01
S0.0075
0.005
a, 0.0025
0
0L .2
0.1
0
(c)
Figure 2: Local void fraction profiles at the measuring
plane.
(c)
Figure 3: Local bubble Sauter diameter profiles at the
measuring plane.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
,800
S700
S400
S300
S200
Case Al
: easurementi
: Prediction
0 0.1 0.2 0.
Height (m)
Case L1
I Measurement
0 0.2 0.4 0.6 0.8 1
S800
700
S600
a,500
0 400
m 300
200
S100
0
0.5
0.4
S0.3
S0.2
0.1
0
0 0.2 0.4 0.6 0.8 1
0 0.1 0.2
Height (m)
700
: Case L3
S600
E Measurement
c 500
8  Prediction
0 400
2 300 *
~ 0
0 0.2 0.4 0.6 08
0.5
Case A3
0.4
eMeasuremenl
& 0.3 Prediction
0.2
0.1 ~~
1 0 0.1 0.2
Height (m)
(c)
n profiles at Figure 5: Axial void fraction profiles.
t
(c)
Figure 4: Local interfacial area concentration
the measuring plane.
eMeasurement
S Prediction
Case Al
0 0.1 0.2 0
Height (m)
Pedcin M e a s u r e m e n t
Case A2
0 0.1 0.2 0.3
Height (m)
(b)
CaMeeasurement
Prediction
0 0.2 0. .
Height (m)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
10000
1000
:00
10
1
Case Al
2L eMeasurement 
0 0.1 0.2 0.3
Height (m)
10000
1000
100
10
1
0 0.1 0.2
Height (m)
10000
1000
100
10
1
E
0n
Case A3
eMeasurement
Prediction
0 0.1 0.2 0
Height (m)
(c)
Figure 6: Axial dimensionless bubble Sauter diameter void
fraction profiles.
(c)
Figure 7: Axial interfacial area concentration profiles.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Yeoh, G. H. & Tu, J. Y., Population balance modelling for
bubbly flows with heat and mass transfer, Chem. Eng. Sci.,
Vol. 59, 31253139 :a II 4).
Yeoh, G. H. & Tu, J. Y., Thermalhydrodynamic modelling
of bubbly flows with heat and mass transfer, AIChE J., Vol.
Yeoh, G. H. & Tu, J. Y., A unified model considering force
balances for departing vapour bubbles and population
balance in subcooled boiling, Nuc. Des. Eng., Vol. 235,
12511265 (2005b).
Yeoh, G. H., Cheung, S. C. P., Tu, J. Y. & Ho, M. K. M.,
Fundamental consideration of wall heat partition of vertical
subcooled boiling flows, Int. J. Heat Mass Transfer, Vol. 51,
38403853 (21 II).
Zeitoun, O. & Shoukri, M., Bubble behavior and mean
diameter in subcooled flow boiling, ASME J. Heat Transfer,
Vol. 118, 110116 (1996).
Zeng, L. Z., Klausner, J. F., Bemnhard, J. F. & Mei, R., A
unified model for the prediction of bubble detachment
diameters in boiling systems II. Flow boiling, Int. J. Heat
Mass Transfer, Vol. 36, 22712279 (1993).
Zuber, N., The dynamics of vapor bubbles in nonuniform
temperature fields, Int. J. Heat Mass Transfer, Vol. 2, 8398
(1961).
Zuber, N., Nucleate boiling the region of isolated bubbles
and the similarity with natural convection, Int. J. Heat Mass
Transfer, Vol. 6, 5378 (1963).
References
Anglart. H & Nylund O., CFD application to prediction of
void distribution in twophase bubbly flows in rod bundles,
Nuc. Sci. Eng., Vol. 163, 8198 (1996).
Basu, N., Warrier, G. R. & Dhir, V. K., Wall heat flux
partitioning during subcooled flow boiling: Part I model
development, ASME J. Heat Transfer, Vol. 127, 131140
I ( II l~ a
Basu, N., Warrier, G. R. & Dhir, V. K., Wall heat flux
partitioning during subcooled flow boiling: Part II model
validation, ASME J. Heat Transfer, Vol. 127, 141148
(2005b).
Hsu, Y. Y. & Graham, R. W., Transport process in boiling
and twophase systems, Washington: Hemisphere (1976).
Klausner, J. F., Mei, R., Bemnhard, D. M. & Zeng, L. Z.,
Vapor bubble departure in forced convection boiling, Int. J.
Heat Mass Transfer, Vol. 36, 651662 (1993).
Prince, M. J. & Blanch, H. W., Bubble coalescence and
breakup in airsparged bubble column, AIChE J., Vol. 36,
14851499 (1990).
Lee, T. H., Park, G.C & Lee, D. J., Local flow
characteristics of subcooled boiling flow of water in a
vertical annulus, Int. J, Multiphase Flow, Vol. 28, 13511368
Luo, H. and Svendsen, H., Theoretical model for drop and
bubble breakup in turbulent dispersions. AIChE J., Vol. 42,
12251233 (1996).
Maity, S., Effect of velocity and gravity on bubble dynamics,
MS thesis, University of Califomnia, Los Angeles (I II n n,.
Prodanovic, V., Fraser, D. & Salcudean, M., Bubble
behaviour in subcooled flow boiling of water at low
pressures and low flow rates, Int. J. Multiphase Flow, Vol.
Sato, Y., Sadatomi, M. & Sekoguchi, K., Momentum and
heat transfer in twophase bubbly flow I, Int. J. Multiphase
Flow, Vol. 7, 167178 (1981).
Sateesh, G., Das, S. K. & Balakrishnan, A. R., Analysis of
pool boiling heat transfer: effect of bubbles sliding on the
heating surface, Int. J. Heat Mass Transfer, Vol. 48,
15431553, (2005).
Xiao, B. & Yu, B., A fractal analysis of subcooled flow
boiling heat transfer, Int. J. Multiphase Flow, Vol. 33,
11261139(~ I n,.
Yun, B. J., Park, G.C., Song, C. H. & Chung, M. K.,
Measurements of local twophase flow parameters in a
boiling flow channel, Proceedings of the OECD/CSNI
Specialist Meeting on Advanced Instrumentation and
Measurement Techniques (1997).
