7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Assessment of the Heat Transfer Model and Turbulent Wall Functions for CFD
Simulations of Subcooled and Saturated Boiling
Deoras Prabhudharwadkar't, Martin Lopez de Bertodano1'*, John Buchanan Jr.2
'Purdue University, School of Nuclear Engineering, 400 Central Drive, West Lafayette, IN47906, USA
tdeoras@purdue.edu, *bertodan @purdue.edu
2Bechtel Marine Propulsion Corporation, West Mifflin PA15122, USA
buchanan@bettis.gov
Keywords: TwoFluid Model, Subcooled boiling, Bubbly Flow, Turbulence Wall Functions
Abstract
This paper describes the details of validation of the twofluid model for subcooled boiling simulations using the CFD code
CFX. Since subcooled boiling occurs at a superheated wall placed in a subcooled liquid, most of the heat and mass transfer
takes place close to the wall. Therefore, this study was focused on the assessment of two important aspects of near wall
predictions, viz., the wall heat transfer model and the turbulent wall functions. Validation was performed using the stateof
theart multidimensional experimental data available in the literature.
The wall heat transfer model is based on splitting the wall heat flux into three components, viz., the single phase
convection, the evaporation at the bubble surface and the quenching effect at the heated wall after the bubble departure. The
evaporation and quenching components use the following three closure parameters: 1. Bubble Nucleation Site Density, 2.
Bubble Departure Diameter, and 3. Bubble Departure Frequency. It is difficult to test three separate uncertainties
simultaneously and hence only those data sets where the bubble diameter has been measured have been selected.
The turbulent wall functions are used to prescribe the near wall boundary conditions for momentum and turbulence
quantities and hence influence the near wall velocity and turbulence distribution. The stateoftheart twofluid model uses a
twophase ke turbulence model (Lopez de Bertodano et al., 1994) with the standard wall function (Launder and Spalding,
1974).
Although the lawofthewall is universal for single phase flows, the same is not true when the boundary layer contains a
twophase bubbly mixture. Previous experiments (Marie et al., 1997) have shown that the velocity profile still follows a
logarithmic profile but the slope and the intercept constant vary with the concentration of dispersed fluid and the relative
magnitude of the shear and buoyant forces. A modification of the law of the wall coefficients has been suggested by Marie et
al. (1997) based on the experimental data obtained from adiabatic airwater bubbly flow over a vertical flat plate. The wall
law coefficients of CFX were modified using the twophase wall function model and significant improvements were noticed
in the turbulent parameter predictions.
Introduction
The prediction of boiling flows depends on accurate
representation of the heat and mass transfer processes at
the wall where the vapor is formed. The accuracy of the
model prediction depends on how well the following
three independent parameters are estimated: 1. Bubble
Nucleation Site Density, 2. Bubble Departure Diameter,
and 3. Bubble Departure Frequency. These parameters
depend on complex physical phenomena that are not
completely understood. Different closure relations exist
in the literature for each of the above parameters which
are discussed in detail in the next section.
One of the goals of the present study is to perform an
independent assessment of each of these parameters from
the perspective of CFD implementation. However, it is
difficult to test three separate uncertainties simultaneously
and hence only those data sets have been selected where the
bubble diameter has been measured. A review of the
departure frequency models was done to select the one best
suited for the present data sets. Then, using the known
departure bubble size and the chosen departure frequency
model, two well known models for the nucleation site
density were implemented and the results were compared.
The heat transfer at the wall also depends on the turbulence
intensity near the wall and this should be correctly
predicted as well. The wall function approach that is
usually used with the Standard kepsilon model to obtain
the near wall turbulence quantities using coarser meshes
needs a correction that accounts for the effect of presence
of a twophase mixture. Such an approach has been
implemented successfully for adiabatic flows (Marie et
al., 1997) and this has been extended for boiling flows
here.
The twofluid model used here was previously tested for
adiabatic flows where the interfacial momentum transfer
terms were validated (Prabhudharwadkar et al., 2009).
The energy equation along with the interfacial heat and
mass transfer correlations were incorporated into the
model. The next section describes in detail the twofluid
model used in this study.
Nomenclature
ai Interfacial Area Concentration
CD Drag Coefficient
CVM Coefficient of Virtual Mass (0.5)
Db Bubble Diameter
Ddep Bubble departure diameter
ek Phase enthalpy
fdep Bubble departure frequency
g Gravitational Acceleration
hi, Latent heat of evaporation
jL Superficial Velocity of Liquid
jG Superficial Velocity of Gas
k Turbulence Kinetic Energy
Ma Interfacial Momentum Source
Amk Mass Source
nste Nucleation site density
q" Heat Flux
vk Velocity (vector)
vR Relative Velocity (vector)
T Temperature
Greek Symbols
ka Volume Fraction of kh phase
a Void Fraction
F Dissipation of Turbulence Kinetic Energy
Pk Density
S Turbulent (Reynolds) Stress
S Dynamic viscosity
A Thermal Conductivity
Superscript
D Drag
TD Turbulent Diffusion
L Lift
W Wall
Subscript
k Phase identifier
1, 1, f, L Continuous (Liquid) Phase
2, g, G Dispersed (Gas) Phase
dep Departure
t Turbulent
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Twofluid Model with Heat and Mass Transfer
For diabetic flows, the ensemble averaged twofluid model
equations (Ishii and Hibiki, 2006) governing the motion of
each phase has the following form:
Mass Conservation:
ak Pk + V.OkPkV = Amk
d t
Momentum conservation:
dv = akvPk + RY)
,Pk + V ak Pk k k kVP + k +Z=k
at
+ apk g_ + Mt + ,,, ,v
(2)
Energy conservation:
z, z, Re
d ape + Va ,U,e,= V.k c+ qR)
aL H (3)
qkk ktk V, (4)
q(k = kVTk qk kVek
In the present calculations of subcooled boiling, the vapor
generated was assumed to be at saturated temperature and
hence was treated as isothermal phase.
Interfacial momentum transfer
The interfacial momentum transfer force comprises of force
terms due to drag, turbulent diffusion, lift and a wall force:
M =MD L TD W
Mki Mi +M +M +M
mki  ki ki ki
The drag force of the bubbles is given by,
D 3 C R R(6)
M Cap1 D VR (6)
4 D,
where the drag coefficient, CD, is given by the IshiiZuber
correlation (1979). Db is the Sauter mean diameter of the
bubble field and the relative velocity is given by
VR =V2 )1
The wall force model accounts for the effect that keeps the
centers of the bubbles no closer than approximately one
bubble radius from the wall. This force is important when
the lift force is present. It is given by (Antal et al., 1991):
M2 =C ,C (. . n,
Ca= min 0 c ,, (7)
Db Ywall / ]
where n and ywal are the normal vector and the
distance from the wall respectively. The CFX values for
the coefficients are: Cw = 0.01 and Cw2= 0.05. In a
previous study (Prabhudharwadkar et al., 2009) it was
found that the Antal's CFX Default coefficients work
well for adiabatic bubbly flows through vertical pipe.
The turbulent diffusion force is closed using the Lopez de
Bertodano Model (Lopez de Bertodano, 1991):
TD
M = CTDPk P Va (8)
The turbulent diffusion coefficient used for the present
calculations is CTD = 0.25.
The lift force is given by (Auton (1987)):
M = CL, pla 2 v)x(vxv,) (9)
The lift coefficient used for the calculations is CL = 0.1
(Lopez de Bertodano, 1991). This value is of the same
order as Tomiyama's experimental value (2002) for a
bubble in Couette flow, CL = 0.288.
Interfacial heat transfer
The interface to liquid heat transfer (k=l) is expressed as,
qj, = h ,a,(Ta,, T,), (10)
where, h1i is the heat transfer coefficient between the
liquid and the interface, closed as follows,
h Nu (11)
d
In the above equation, Al is the liquid thermal
conductivity and d, is a length scale which is assumed
equal to the bubble diameter. The Nusselt number (Nu)
in the above closure equation is given by following
expression (Ranz and Marshall, 1952):
Nu = 2 + 0.6Re5 Pr033 (12)
Reynolds number (Reb) used in the above closure
defined as,
Reb =
vl
The interfacial mass transfer is given as,
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
h,(,T,)a h, (Tl T,)a,
cond,L max 0 eap, = in , 0
(14)
where, Hig is the latent heat of phase change.
Wall heat transfer
Wall heat transfer is the most important aspect of subcooled
boiling as it provides the sources for energy and phase mass
balance equations. Most wall heat transfer models are based
on partitioning the wall heat flux into three components.
These were implemented in the multidimensional CFD
calculations by Kurul and Podowski (1990). The three
components of the wall heat flux are 
a) Evaporation heat flux (qe) The part of heat flux
utilized in formation of vapor at the wall,
b) Forced convection to the liquid (qc),
c) Wall quenching by liquid phase (transient conduction)
(qq) This accounts for the heat transfer to the
subcooled liquid that replaces the detached bubble at
the wall.
The wall surface is assumed to be split into two parts (A1,
A2) each under the influence of one phase. Fraction A2 is
influenced by the vapour bubbles formed on the wall and
participates in the evaporation and quenching heat transfer.
Fraction A1 is the remaining part of the wall surface,
(AI=1A2) and participates in the convective heat transfer
to the liquid. A1 and A2 are related to the nucleation site
density per unit wall area (n,i,) and to the influence area of
a single bubble forming at the wall nucleation site. The
KurulPodowski model assumes, that the diameter of the
bubble influence zone is twice as big as the bubble
departure diameter (a = 2)
A2 = (a Dde)2 nse A = 1 A2
4
Evaporation Heat Flux:
The evaporation heat flux is obtained from the rate of vapor
generation at the wall which is given as a product of mass
of a bubble at detachment, the detachment (departure)
frequency and the nucleation site density,
SD3
, de PG fdep ns1te
q, = fm A
es is
The specification of evaporation heat flux using above
equation requires closure for bubble departure diameter,
(13) bubble departure frequency and the site density. The default
available model in CFX uses following closure relations:
The bubble departure diameter is given by the empirical
correlation of Tolubinsky and Kostanchuk (1970), the
Nucleation Site Density is obtained using the correlations of
Lemmart and Chawla (1977) and the departure frequency is
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
closed using the simplest available expression (Cole,
1960), which uses a characteristic bubble velocity
(terminal velocity of bubble rise) and a characteristic
bubble size (departure diameter).
The frequency of bubble departure obtained by Cole
(1960) was derived under the assumption of bubble
moving with terminal velocity once detached from the
surface. This is a fair assumption for pool boiling
scenario at low heat fluxes where bubble detachment is
hydrodynamically governed, but may not be accurate for
flow boiling at high heat fluxes where the
thermodynamics governs the bubble growth and
detachment. The CFX default closure relations are given
below:
Ddep = Dteexp( (T Tiq )/Tt),
Dr = 0.6 mm, Tef =45K (18)
nsite = n/ ((LTwa Ts ) / ATrs )' 805
n, = 7.9384x105 m2, ATre =10K (19)
= 4gAp ..
d 3DdepPliq
Quenching Heat Flux:
As mentioned previously, the quenching heat flux is the
component of wall heat flux utilized to heat the cold
liquid replacing the detached bubble adjacent to the
heated wall. In order to evaluate this component, Mikic
and Rohsenow (1969) used an analytical approach
starting with transient conduction in semiinfinite
medium with heated wall being the only boundary where
temperature is specified.
DT(x,t) la2T(x,t)
t )=a ; T(x,o)=TF,T(o,t)=T,T(oo,t)=T
dt dx
(21)
Kurul and Podowski (1991) made an assumption that
quenching occurs between detachment of one bubble and
appearance of next bubble (nucleation), and this time
was assumed to be 80% of the detachment period. The
final form of quenching heat transfer closure of Kurul
Podowski model is as follows:
S2A f t=0.8 (22)
rzra f
q = A2 h(T T ) (23)
Single Phase Convection Heat Flux:
This is evaluated under the standard assumption of a
logarithmic temperature profile across the turbulent
boundary layer (Kader (1981), ANSYS CFX Solver
Theory Guide ,"2 II .,1).
qw,lph = + w w Thql
T+
T+ =Pr, y er +(2.121n(y*)+fp)el/F,
,=f(Pr),F =/f(y),
u =C kY2 y p u* Ay/p,
(25)
where, Ay is the distance of the wall adjacent node from the
wall, ki is the liquid turbulence kinetic energy. In the above
equation fl is a function of liquid Prandtl number and F is a
function of y (ANSYS CFX Solver Theory Guide i'i I
The convective heat flux component through the liquid area
fraction is thus given as,
qc =Alfhc( T,q ), h = C
T+
It should be noted here that all the wall heat flux model
correlations above use a liquid temperature (Tliq). The
single phase heat transfer correlation above uses the wall
adjacent node temperature. However, in the quenching heat
transfer correlation which was based on a onedimensional
model, the liquid temperature refers to the bulk mean
temperature. As a good estimate, CFX approximates this
temperature with the temperature at a fixed y (250). This is
done in order to have a mesh size independent evaluation of
wall heat flux partitions.
Review of the Bubble Departure Frequency and
Nucleation Site Density models:
Situ et al. 21"') recently reviewed all the available
departure frequency models and proposed a correlation for
bubble departure frequency based on experimental data of
subcooled boiling flow for wide ranges of pressure and heat
flux (Basu et al. ,'1ll" ), Thorncroft et al. (1998) and Situ
2I'i1n,4 The departure frequency was correlated to the
boiling heat flux as follows:
N = 10.7N .634 (27)
Nfd is a nondimensional representation of frequency:
Nfd fdep ,
Nfd
aL,q
where, aLiq is the liquid thermal diffusivity.
NqNB is the dimensionless nucleate boiling heat flux
obtained from Chen's correlation.
qNB Ddep
NqNB NB D,
aLiqPg hf
q NB hNB (T Tsa,,, )
S(0 79 045 049
hN, = S(0.00122) o1029h24~n
(29)
(30)
T T )024Ap 075
(31)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Ap = p(T) p(Ts,) (32)
S is a suppression factor which is described in detail in
the next section on site density (equation (41)).
However, when the frequency data of Thorncroft et al.
(1998), Basu et al. (2005) and Situ iI,1 i), was compared
with the prediction of Situ et al.'s correlation and Cole's
model, the accuracy of the models was found comparable
as shown in the figure below (Fig. 1 (a) and (b)).
(a)
Nfd (experiment)
Nfd (experiment)
Figure 1: Comparison of the bubble departure frequency
models
The data scatter has been contained by the dotted lines in
the figures above which are spaced by an order of
magnitude from the correlation. The data scatter about
the prediction is comparable and hence selection of the
Cole's model (default CFX model) was reasonable for
the present problem.
A correlation for the nucleation site density has been
proposed by Hibiki and Ishii (2003) which is valid over
1198 bar and for most practical combinations of fluid
and surface material (e.g., WaterStainless steel, Water
Copper, WaterZr4, R113Nichrome etc.). This
correlation is essentially an improvement of the previously
established correlation of Kocamustaffaogullari and Ishii
(1983) validated against water data for 050 bar. The
KocamustaffaogullariIshii (KI) correlation relates the site
density to the wall superheat as follows:
n,,e = R44 fR ), (33)
where,
nsie=nte Ddep (34)
R = Rc/(Dde /2) (35)
R, is the critical radius of the surface cavity which
represents a minimum cavity size which can be activated at
a given wall superheat. R, is given as,
2 a(1+ (p, p)) /P
Rc t.1(, AT /(RTgsa, ))1)
Under the following conditions:
Pg << p, hg, (T, Ts,) (RTTa,)<<1,
R, can be simplified as,
Rc = 2oTst,, (, hf ATs,,,).
The density dependent parameter in equation (33) is given
as,
f(p*)= 2.157x10 7p*32(1+0.0049p*)413
P = (P, )/P,
The effective wall superheat in equation (36) is usually less
than the actual wall superheat. The reason for this is as
follows. A bubble nucleated at the wall grows through a
liquid film adjacent to the wall where a considerably high
temperature gradient exists. So, in reality, it experiences a
lower mean superheat than the wall superheat. In case of
pool boiling, the difference is not significant and the
superheat based on the wall temperature can be taken as the
effective superheat. In case of forced convective boiling,
there exists a steeper temperature gradient in the near wall
field due to thermal boundary layer. Hence, the effective
superheat experienced by the bubble is smaller than the
actual wall superheat and it is given as,
The multiplier S in the above equation is the superheat
suppression factor and is given as (Chen, 1966),
S=1/( +1.5x105 Re1,)
The twophase Reynolds number is given as,
Re, = [G (1 x)dh /, ]F125
Situ et al.(2008) Correlation Prediction
(Rel
AT,,,, = SAT,, where, AT, = T7 T1
'" sate or"a at w at
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
RESULTS
Validation of the Heat and Mass Transfer Models
The Martinelli parameter can be approximated as,
X=,,(= 9 'g 5, f U (44)
Sx A ) p [g)
Thus, equations (33)(44) represent the complete model
of KocamustafaogullariIshii (1983), which has been
incorporated in CFX12 for the present study.
Turbulence Transport
The closure for the Reynolds stresses in equation (2) is
based on the twophase ke model developed by Lopez de
Bertodano et al. (1994). This model assumes that the
shear induced and bubble induced turbulent stresses are
added together:
Re Re) +(Re) (45)
=k =k SI =k BI
The resulting expression for the diffusivity of
momentum (effective viscosity) of the liquid phase is:
k2 
Va C +CDBaDb R
(46)
where, the first term on the RHS corresponds to the ke
model for the shear induced turbulence viscosity and the
second term corresponds to Sato's model (1981) for the
bubble induced turbulence viscosity. The coefficient C, =
0.09 is the standard value according to the ke model. A
value of 0.6 is recommended for CDB. It is important to
note that the kE model transport equations for the liquid
phase are solved together with the continuity and
momentum equations (Eqs. (1) and (2)). The standard
coefficients of the ks model are left untouched.
Standard Wall Functions (Launder and Spalding (1974))
are used to model turbulence quantities near the wall
which avoids need to resolve the complete boundary
layer. Turbulence kinetic energy dissipation rate
(epsilon) is specified using an algebraic closure
assuming turbulence production due to shear is balanced
by dissipation. Shear stress at the wall is connected to the
velocity in the fully turbulent boundary layer
(computational node next to the wall) using a
logarithmic law of wall which is described in detail later.
This completes the description of the twofluid model
and all the closure relations used in the present study.
The results of the comparison of two previously
mentioned nucleation site density model are discussed in
the subsequent section. The results reported here are
restricted to the prediction of vapor void fraction
distribution which of primary interest in this study.
1. Subcooled boiling of R12 at 1527 bar (Morel et al
(2003))
Morel et al. (2003) performed subcooled boiling
experiments in a vertical pipe of 19.2 mm internal diameter
and a length of 5 m. The pipe had a heated section of 3.5 m
preceded by and followed by unheated lengths of 1 m and
0.5 m respectively. The input parameters for the four tests
reported were as summarized in Table 1. The liquidvapor
density ratio in these experiments corresponds to that for
watersteam at 95150 bar. Dual sensor fiber optics probes
were used to measure void fraction
Table 1: Simulation input parameters for tests tpl and tp6
Parameter Test1 Test2 Test3 Test4
(Deb5) (Deb6) (Debl3) (DeblO)
Mass Flux 1986 1984.9 2980.9 2027.0
(kg/m2)
Pressure 26.15 26.15 26.17 14.59
(bar)
Inlet 18.12 16.11 18.12 23.24
Subcooling
(C)
Heat flux 73.89 73.89 109.42 76.24
(W/m2) __
A 2" wedge of the pipe is used as the domain. The mesh
had 20 radial and 250 axial uniformly spaced elements.
Fig. 2 shows the domain and crosssection of the mesh
used.
Figure 2: Domain and Radial Mesh for R12 boiling
simulation case
All the data sets include radial profiles of bubble diameter.
Note that with Freon at these test pressures, the surface
tension is of the order of 0.001 N/m which results into very
F=1.0 ,forX,, >10
=2.35(0.213+1/X,,)0736,forX,, <10
(
small bubble sizes (i.e., 0.5 mm). For the CFD
calculations the mean of the radial profile value is used
for the bubble diameter in the bulk, whereas, the value at
the wall is used as the departure bubble size (Fig. 3). It
was found during all the cases that a better match with
experimental data is obtained with the coefficient of the
turbulent diffusion force CTD = 0.5 instead of the value
obtained for adiabatic air water flows, CTD = 0.25. This
change may be attributed to bubble diameters in these
cases which are one order of magnitude smaller than that
in atmospheric airwater systems (Prabhudharwadkar et
al., 2009). Therefore there are more interactions of the
bubbles with smaller turbulent eddies.
Fig. 4 shows the prediction of radial variation of void
fraction obtained using two different site density models.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
The LemmartChawla (1977) correlation which is used in
the KurulPodowski model was found to underestimate the
void fraction in all the cases. KocamustafaogullariIshii
model (termed as KI Model in the figures below) predicted
the data well in all the cases. The prediction however
deteriorated at low pressure data (P = 14.59 bar) and this
needs further investigation.
P = 26.15 bar, q" = 73.89 kW/m2, G = 1986 kg/m2 s,
dTsub = 18.12 C
Experiment (Morel et a] (2003))
Mean
r ^ ^ ^ urr^ ^ U1ur
2 E04
0 E+00
0 0002 0004 0006 0008 001
r (m)
P = 26.17 bar, q" = 109.42 kW/m2, G = 2980.9 kg/m2 S
dTsub= 18.12 C
Experiment (Morel et al (2003))
Mean
4 E04 f at
2 E04
0 E+00
0 0002 0004 0006 0008
r (m)
P = 26.15 bar, q" = 73.89 kW/m2, G = 1984.9 kg/m2 s,
dTsub = 16.11 C
Experiment (Morel et al (2003))
Mean
0 E+00
0 0 002 0 004 0 006 0 008 0 01
r (m)
P= 14 59 bar, q" = 76 24 kW/m2, G = 2027 kg/m2 s,
dTsub= 2324 C
1 5E 03
Experiment(Morel et al (2003))
Mean
1 OE 03 (m*
**%
5 OF 04 
0 OE+00
0 0 002 0004 0006 0008 001
r (m)
Figure 3: Bubble diameter approximation for the simulations
1 E03
8 E04
6 E04
4 E04
2 E04
1 E03
8 E04
6 E04
~*,
~55
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
P = 26.15 bar, q" = 73.89 kW/m2, G = 1986 kg/m2 s, P = 2.15 bar, q" 73.89 kW/m2, G =1984.9 kg/m2 s,
dTsub = 18.12 OC, D= 0.452 mm, De,= 0.235 mm dTsub =16.11 C, D = 0.498 mm, Ddp= 0.25 mm
P = 26.17 bar, q" = 109.42 kW/m2, G = 2980 kg/m2 s,
dTsub 17.44 C,D= 0.384 mm, Dd= 0.146 mm
 SiteDensiy LemmarChawla(1977)
^^ ^
Figure 4: Void Fraction prediction using the two Nucleation Site Density models
2. Subcooled Boiling of R113 at 2.69 bar in an
Annulus (Roy et al (2002))
Roy et al. (2002) performed subcooled boiling
experiments in a vertical annulus of 15.78 mm internal
diameter and 38.02 mm outer diameter. The total
length of the channel was 3.66 m of which the initial
0.91 m was unheated (adiabatic). The pipe had a heated
section of 2.75 m.
Six experiments were reported (designated as tpltp6)
and two of them (having minimum and maximum
mass flux, tpl and tp6) have been selected for
simulation purpose here. Radial Profiles of Void
Fraction were measured with a dualsensor fiberoptic
probe. The input data for the two representative cases is
stated in Table 3:
Table 3: Simulation input parameters for tests tpl and
tp6_
Parameter tp tp6
Mass Flux (kg/m2) 568 784
Pressure (bar) 2.69 2.69
Inlet Subcooling (C) 42.7 50.2
Inner wall heat flux (W/m2) 95000 116000
Similar to previous problem, a 2 section of the annulus
is used as the domain. The mesh had 25 radial and 266
axial uniformly spaced elements. Fig. 5 shows the
domain and crosssection of the mesh used.
Figure 5: Mesh crosssection for R113 simulations
These experiments had a characteristic near wall
boundary layer of vapor (boiling layer) and the bulk of
the pipe (more than 50%) had only liquid. Bubble
diameter was reported for the test tp6. As the pressure
was same in all the experiments and the range of other
parameters was also narrow, same values were used for
the other simulated case (tpl). The bubble diameter
profile is shown below in Fig. 6.
P = 14 59 bar, q" = 76 24 kW/m2, G = 2027 8 kg/m2 s,
dTb = 23 24 oC, D=0 916 mm, Ddp=0 37mm
SExperment(Morel et a (2003))
Site Density K (198)
SiteDensty Le a Chawa(1977)
b.. *..*44W0
 MipnMe ditmeern
S  Sutr mc dinamccr
2 : n Ie r
S A,
assumed to be 1 mm and the departure bubble diameter
0.6 mm arrived at from the "most probable" diameter
values reported in the experiments (Fig. 6).
Simulations were also done with the Sauter mean
diameter values in the bulk and at the wall; however
the predictions with most probable values matched the
data well. As stated earlier, the simulations were
performed with the KI site density model only since
this data set is in the pressure range for which the KI
model is already found to be best suited through
previous benchmark problems. The void fraction
prediction for the two cases is shown below in Fig. 7.
The width of the boiling layer was well predicted
with the most probable bubble sizes; however the near
wall void fraction was underpredicted in both the
cases.
Twophase Wall Function Implementation
The velocity profile in the turbulent boundary layer
next to the wall follows a logarithmic profile as given
below:
U =1log(y+)+c
ic
where, U+=U y+ = uy = (48)
w, Y
U, V p
K is the von Karman constant (= 0.41) and C is a log
layer constant (= 5.2). U is the velocity component
tangential to the wall and y is the distance normal to
the wall.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
0.6
Experiment(Roy etal. 2002)(tpl)
0.5 Db = 1mm, Ddep =0.6 mm
SDb= 1.7 mm, Ddep 1.2 mm
0.4
0.2
0.1
0.2 0.4 0.6
(rri)l/rori)
0.8 1
0.2 0.4 0.6 0.8
(rri)/(rori)
Figure 7: Void fraction distribution prediction for tpl
and tp6 with most probable and mean diameter values
The above definitions are in context of single phase flow.
Marie et al. (1997) performed experiments where a flat
plate was placed parallel to bubbly airwater flow. Void
fractions and velocity profiles were measured in the
boundary layer next to the wall for different inlet
concentrations of the air. It was found that the velocity
profile still follows a logarithmic lawofthewall but
with a modified slope and log layer constant. When the
streamwise momentum equation for liquid phase was
integrated through the fully developed twophase
boundary layer, following equation was obtained:
(ia(vL UL = L u g(a )BL (49)
where, (ULVL) is the turbulent shear stress (Reynolds
stress), a is the average void fraction in the boundary
layer and aO is the average void fraction in the free
stream. YBL is the boundary layer thickness. The second
term on the right side of the above equation represents
the buoyancy due to presence of dispersed phase.
Neglecting the viscous stresses, the above equation can
be simplified for dilute dispersed flows as,
(ULVL) "L g(a a )BL = 2 (50)
The twophase log layer equation proposed by Marie et
al. (1997) uses u, as the velocity scale instead of u, in
equations (4748). When this velocity scale was used
for the law of the wall it was found that the log layer
profiles for all the inlet void fraction were parallel and
the log layer coefficient C increased with the inlet void
fraction. Fig. 8 below shows the results obtained by
Marie et al (1997) for different inlet void fractions after
using a modified velocity scale.
35
30
25 7
15
S0 a 0 a = o%
10 A a a = 2% 
V a a =3.5%
5 0 0 a 6..
o10 10' 12 /1 103
Figure 8: Log layer velocity profiles measured by
Marie et al. (1997) for different air void fractions
The log layer constant for twophase mixture, C*, was
related to the single phase C with the following
relation.
C*= C+yo+ (1/l 1)(1/Ic)log1 (51)
y +represents the edge of the laminar sublayer which
is usually about 11 and
/=U (52)
U
Based on the above model, the log layer coefficients
were modified in the CFD code CFX and turbulence
quantities were compared with single phase and two
phase law of the wall. The above mentioned
experimental benchmark of Roy et al (2002) included
measurements for turbulence kinetic energy and
Reynolds stress. A single phase case was first tried to
validate the single phase wall function included in the
CFD code CFX. Fig. 9 shows turbulence kinetic energy
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
and Fig. 10 shows the turbulent shear stress profile for a
single phase experiment reported by Roy et al (test spl:
Mass Flux = 568 kg/m2, Pressure = 2.69 bar, Inlet
Temperature = 315.85 K, Wall Heat Flux = 16 kW/m2).
The results were satisfactory confirming adequacy of the
mesh size and model accuracy for the prediction of near
wall turbulence quantities.
0005
Experiment Royetal (2002)
0004 CFXResults
0003
0002
0001
0 02 04 06
(rri)/(rori)
08 1
Figure 9: Turbulence kinetic energy prediction for
single phase case
0000o
0000
0 ooo0005
0 001
Figure 10: Reynolds stress prediction for single phase
case
The two phase log law coefficients were then
implemented in CFX to compare the effective turbulence
kinetic energy and turbulent stress. The effective
turbulence kinetic energy superposes the shear induced
and bubble induced components as follows:
kL, = kL,SI + kL,BI; kBI = 4 avR (53)
Significant improvement was noticed in the prediction
of the turbulence kinetic energy (Fig. 11) for both the
previously reported cases especially near the boiling
boundary layer edge where a single phase wall function
showed a characteristic dip.
The Reynolds stress (obtained using equations (4546))
was found to be more accurate with the modified wall
law coefficients (Fig. 12).
Experiment (Royetal (2002))(spl)
CFX Results
4 06 0
S 04 06 08
(rri)/(rori)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
S 0.2 0.4 0. 0 O 1 0 0.2 0.4 0.6 08 1
(rri)/(rori) (rri)/(rori)
Figure 11: Turbulence kinetic energy prediction for wallboiling cases tpl and tp6 of Roy et al. (2002)
0002 0.002
0.001 0.001
0 0
0.2 0.4A 0.6 0.8 0.2 04 0.6 0.8
S 1 0.001
A A
S p Experiment(Roy etal. (2002))(tp5)
o002 0002
Experiment(Roy et al. (2002))(tpl) 0.002
Originalwall function ginWall Function
003 Twophae wall function 0.003 Twophase wall function
0004 0.004
(rri/lirori) .00(rri/(rori)
Figure 12: Reynolds Stress prediction for wallboiling cases tpl and tp6 of Roy et al. (2002)
assemblies with spacer. Nuclear Engineering and
Conclusions Design, Vol. 177, 215228 (1997)
1. The heat flux splitting model for wallboiling
flows was validated using data sets where bubble
diameters were known.
2. The nucleation site density model by
Kocamustaffaogullari and Ishii (1983) in
combination with the Bubble Departure Frequency
model of Cole (1960) is found to predict the data
well.
3. Predictions were found to deteriorate at lower
pressures of the R12 data; however the results were
still satisfactory.
4. Use of a modified log law coefficients based on the
approach of Marie et al. (1997) improved
prediction of the turbulence quantities near the
wall.
References
Anglart, H, Nylund, O., Kurul, N., Podowski, M.Z.
CFD prediction of flow and phase distribution in fuel
ANSYS CFX Solver Theory Guide, ANSYS CFX
Release 11.0, (2006)
Antal, S.P., Lahey, R.T. Jr., Flaherty, J.E. Analysis of
phase distribution in fully developed laminar bubbly
twophase flow. International Journal of Multiphase
Flow, Vol. 17, No. 5, 635652 (1991)
Auton T. R., Hunt J. C. R. and Prud'homme M. The
force exterted on a body in inviscid unsteady non
uniform rotational flow. Journal of Fluid Mechanics,
Vol. 197, 24157 (1987)
Basu, N., Warrier G. R., Dhir, V.K. Wall heat flux
partitioning during subcooled flow boiling: Part II
Model validation. Journal of Heat Transfer, Vol. 127, 73
93 (2005)
Chen, J.C. Correlation for boiling heat transfer to
saturated fluids in convective flow. I& EC Proc. Des.
Dev., Vol. 5, 322329 (1966).
Experiment(Roy etal. (2002))(tpl)
 OriginalWall Function
 Twophase Wall Function
* *^ e"V e  .
E
0.004
0002
Cole, R. A Photographic Study of Pool Boiling in the
Region of The Critical Heat Flux. AIChE Journal, Vol.
6, 533 534 (1960)
Egorov, Y, Menter, E Experimental Implementation of
the RPI Wall Boiling Model In CFX11, Ansys CFX,
Staudenfeldweg 12, 83624 Otterfing, Germany (2004)
Hibiki, T. & Ishii, M. Active nucleation site density in
boiling systems. International Journal of Heat & Mass
Transfer, Vol. 46, 25872601 (2003)
Ishii, M. & Hibiki, T., Thermofluid dynamics of two
phase flow, Springer (2006)
Ishii, M. & Zuber N., Drag coefficient and relative
velocity in bubbly, droplet or particulate flows. AIChE
Journal, Vol. 25, No.5, 843855 (1979)
Kader, B.A. Temperature and concentration profiles in
fully turbulent boundary layers. International Journal of
Heat and Mass Transfer, Vol. 24, Issue 9, 15411544
(1981)
Kocamustaffaogullari G. & Ishii M. Interfacial area and
Nucleation Site Density in boiling systems.
International Journal of Heat and Mass Transfer, Vol.
26, 13771397 (1983).
Kurul, N. & Podowski, M. Z., Multidimensionnal
effects in forced convection subcooled boiling,
International. Heat Transfer Conference, Jerusalem,
Vol. 1, 2126 (1990)
Lahey, R.T., Jr., Lopez de Bertodano M. and Jones
O.C., Jr., Phase distribution in complex geometry
conduits, Nuclear Engineering and Design, Vol. 141,
pp. 177201, 1993.
Launder B.E., and Spalding D.B., The numerical
computation of turbulent flows, Computer Methods in
Applied Mechanics and Engineering, Vol.3, pp. 269
289, 1974.
Lemmert, M., Chawla, J.M. Influence of flow velocity
on surface boiling heat transfer coefficient, In: Hahne,
E., Grigull, U. (Eds.), Heat Transfer in Boiling,
Academic Press and Hemisphere, 237247 (1977)
Lopez de Bertodano, M. Turbulent Bubbly Flow in a
Triangular Duct, Ph. D. Thesis, Rensselaer Polytechnic
Institute, Troy New York, 1991.
Lopez de Bertodano, M. A., Lahey, R. T., Jr., and
Jones, O. C., Jr. Development of a ke model for bubbly
twophase flow. Journal of Fluids Engineering, Vol.
116, 128134 (1994)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Mikic, B.B. & Rohsenow, W.M. A new correlation of
poolboiling data including the fact of heating surface
characteristics, Journal of Heat Transfer, 91, 245250,
(1969)
Morel C., Yao W. & Bestion D. Three dimensional
modeling of boiling flow for the NEPTUNE code,
Proceedings of the 10th International Topical Meeting on
Nuclear Reactor Thermal Hydraulics, Korea (2003)
Prabhudharwadkar, D.M., Bailey, C.A., Lopez de
Bertodano, M.A. & Buchanan Jr., J.R. TwoFluid CFD
model of adiabatic airwater upward bubbly flow through
a vertical pipe with a OneGroup Interfacial Area
Transport Equation, Paper FEDSM200978306,
Proceedings of the ASME 2009 Fluids Engineering
Division Summer Meeting, Vail, Colorado, USA (2009)
Ranz, W.E. & Marshall, W.R. Evaporation from drops:
Parts I&II. Chem. Eng. Prog., Vol. 48, Issue 3, 141146
(1952)
Roy, R.P, Kang, S., Zarate, J.A., Laporta, A., Turbulent
subcooled boiling flow: Experiments and Simulations.
Journal of Heat Transfer, Vol. 124, 7393 (2002)
Sato, Y, Sadatomi, M. & Sekoguchi, K. Momentum and
heat Transfer in Twophase Bubble Flow. Int. J.
Multiphase Flow, Vol. 7, 179190 (1981).
Situ, R., Experimental and theoretical investigation of
adiabatic bubbly flow and subcooled boling flow in an
annulus, Ph.D. Thesis, Purdue University, USA (2004).
Situ R., Ishii M., Hibiki, T., Tu, J.Y, Yeoh, G.H., Mori,
M. Bubble departure frequency in forced convective
subcooled boiling flow. International J. Heat and Mass
Transfer, Vol. 51, 62686282 (2008)
Thorncroft, GE., Klausner, J.E, Mei, R. An
experimental investigation of bubble growth and
detachment in vertical upflow and downflow boiling.
International Journal of Heat and Mass Transfer, Vol.
41, 38573871 (1998)
Tolubinsky, V.I. & Kostanchuk, D.M. Vapour bubbles
growth rate and heat transfer intensity at subcooled water
boiling. Proceedings of the 4th International Heat
Transfer Conference, Vol. 5, Paris (Paper No. B2.8),
(1970)
Tomiyama, A., Tamai, H., Zun, I., Hosokawa, S.
Transverse migration of single bubbles in simple shear
flows. Chemical Engineering Science, Vol. 57, 1849
1858 (2002)
