7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Prediction of Subcooled Wall Boiling in a Heated Annulus Including Conjugate Heat
Transfer in the Central Heated Rod
Th. Frank*, P. Beckstein*, C. Lifante*, A.D. Burnst
ANSYS Germany GmbH, Staudenfeldweg 12, Otterfing, D83624, Germany
T School of Process Material and Environmental Engineering, CFD Centre, University of Leeds, LS2 9JT, UK
Email: Thomas.Frankg(ansvs.com
Keywords: CFD, watersteam flow, subcooled boiling, wall boiling, CHT, phase change
Abstract
The aim of this paper is to present the validation of a newly implemented wall boiling model for the prediction of subcooled
nucleate boiling under pressurized and normal conditions, so e.g. in rod bundles and fuel assemblies of PWR and BWR, first made
available with full GUI support in ANSYS CFX 12.0 (ANSYS Inc., 2009). The model formulation is based on the socalled RPI wall
boiling model published by Kurul & Podowski (1991) but has been extended in ANSYS 12.0 by a special approach regarding the prediction
of the liquid subcooling temperature in order to make the model grid independent in the context of a CFD simulation. Furthermore the wall
boiling heat partitioning algorithm has been coupled to the prediction of conjugate heat transfer in the solid material of the heated boundary
of the fluid domain, allowing for the definition of volumetric energy sources in e.g. the rods of nuclear fuel assemblies heated by nuclear
fission reaction. The Paper discusses the theory of the implemented model as well as the application of the model to boiling flow in
vertically directed circular annulus with a centralized heated rod, as investigated by Lee et al. (2008). CFD results are compared to
experiments of Lee et al. applying a hierarchy of consistently refined meshes and different CFD setup configurations including the
prediction of heat transfer in the solid material of the heated rod by CHT and by prescribing a volumetric thermal energy release in the solid
material instead of a prescribed wall heat flux or wall temperature at the heated boundary of the fluid flow domain.
1. Introduction
Subcooled flow boiling occurs in many industrial
applications and is characterized by large heat transfer
coefficients. In particular the regime of subcooled flow
boiling is of interest for flows through nuclear reactor fuel
assemblies of PWR and BWR. However, this efficient heat
transfer mechanism is limited by the critical heat flux where
the heat transfer coefficient decreases leading to a rapid
heater temperature excursion potentially leading to heater
melting and destruction. Therefore the occurrence of the
critical heat flux flow regime has to be safely avoided.
Furthermore the flow in the regime of nucleate subcooled
boiling and the transition to critical heat flux is affected by
many flow parameters and can be influenced by the
geometrical design of the fuel assemblies. Multiphase flow
morphology and bubble/slug dynamics play an important
role in the radial fluid temperature and vapour distribution,
substantially influencing the boiling process. From the point
of view of geometry design especially the spacer grids
equipped with mixing vanes are of importance in increasing
the permissible heat flux.
The study of fuel assembly designs and factors of influence
on the boiling process usually requires extraordinary
expensive experiments. Therefore the supplementation or
even the replacement of experiments by numerical analyses
are of relevant interest in the study of boiling flows. Since
flow patterns in boiling flows in e.g. fuel assemblies are
complex 3d with lateral crossflows and resulting steam
redistribution on heater surfaces, it is highly desirable to
derive CFD methods for the accurate prediction of boiling
processes. The paper describes the actual implementation of
a wall boiling model in ANSYS CFX 12.0 and its
application to a validation testcase under normal pressure
conditions similar to conditions in fuel assemblies of BWR.
Special attention has been given to the raise of the constant
wall heat flux assumption by facilitating the coupled heat
transfer prediction in adjacent solid and fluid domains of the
flow configuration for wall boiling processes and thereby
predicting the wall heat flux directly.
2. Nomenclature
bubble influence factor
wall fraction cooled by singlephase convection
wall fraction cooled by quenching
interfacial area density
specific heat capacity of liquid
specific heat capacity of vapour
bubble diameter in the bulk flow
bubble departure diameter on the wall
bubble departure frequency
gravitational acceleration
mass flow rate
interfacial heat transfer coefficient
evaporation enthalpy
coefficient for heat transfer by quenching
liquid heat conductivity
height of the measurement cross section
length of the experimental test section
active nucleation site density
Nusselt number
Prandtl number of liquid
wall heat flux
heat flux due to singlephase convection
heat flux due to evaporation
heat flux due to quenching
total wall heat flux
radius
inner radius of outer tube
outer radius of inner tube/heating rod
dimensionless width of the circular annulus
Reynolds number
Volume fraction of gaseous phase
nearwall liquid temperature
saturation temperature
=TsatTL; liquid subcooling temperature
=TwTsat; wall superheating
wall temperature
liquid characteristic nearwall temperature
bubble waiting time
nondimensional distance to the wall
vapour velocity
liquid velocity
Greek letters
PG density of the gaseous phase
PL density of the liquid phase
Tw wall shear stress
Subsripts
G gaseous phase
L liquid phase
in inlet properties
max maximum
ref reference quantity
sub subcooling
sup superheating
tot total
3. The Experimental Testcase
The present investigation of CFD model development and
validation for nucleate subcooled boiling under almost
normal pressure conditions has been based on a recently
published series of experiments by Lee et al. (2008). Further
experimental investigations on the same experimental test
facility have been published in recent years (Lee et al., 2002
& 2009).
The test loop utilized in the cited experiments is aimed to
the investigation of subcooled nucleate boiling in a vertical
circular annulus with a centralized heater rod (see Fig. 1),
where subcooled boiling occurs on the heated surface. The
test section is a 2376mm long vertical concentric annulus
with a heated inner tube. The inner tube is of Ro=9.5mm
outer radius and electrically heated. The length of the heated
test section is LT=1670 mm consisting of an Inconel 625
tube with a 1.5 mm thickness. The tube is uniformly heated
by a 54 kW DC power supply. The outer tube is comprised
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
of two stainless steel tubes of R=18.75mm inner radius,
which are connected, below the measuring plane, by a
transparent glass tube of 50 mm length, therefore facilitating
optical measurements of radial distributions of flow
parameters. The circular annulus of the test section has a
gap width of 9.25mm and therefore a hydraulic diameter of
18.5mm. The heated section begins at the same elevation as
the test section inlet, and the measuring plane is located at
LM=1610mm downstream of the beginning of the heated
section inlet.
t Outlet
LM LT
t Inlet
Measuring
Plane
[for experimental
and numerical
Results)
Heated Wall
AdabaticWall
InnerTube
(Heaing Rod)
Annulus
OuterTube
Figure 1: Schematic of the experimental test section of the
vertical circular annulus in the experiments of Lee et al.
(2008)
Experiments were carried out for the various levels of heat
flux, mass flux, inlet liquid temperature and inlet pressure. A
total of 32 data sets were investigated, where detailed
measurement data for radial distributions of local void
fraction, the vapour and the liquid velocities, Sauter mean
vapour bubble diameter and the local interfacial area
concentration are available for only 12 out of the 32
experimental conditions from (Lee et al., 2008). The local
vapour phase parameters were measured by a doublesensor
conductivity probe method traversed over 13 points in the
measurement cross section of the annulus with the
dimensionless width of the annulus:
R, =(rRo)/(RRo) (1)
between 0.11 and 0.90 to obtain the radial profiles of the
local flow parameters. Here, r, R and Ro are the radial
location measured from the symmetry axis, inner radius of
Axis
outer tube and outer radius of inner tube, respectively. The
local liquid velocity was measured by the Pitot tube method.
Local bubble diameter and interfacial concentration has
been reconstructed from these measured data in a
postprocessing step.
Set No.
25
16
q
[kW m2]
220.0
320.4
G
[kg m2sl
1057.2
718.8
Tin
90.1
83.8
Pin
rkPa]
134.4
121.1
Table 1: Investigated operating validation testcase
conditions from the experiments of Lee et al. (2008)
Two distinct operating points of the experimental test
facility, characterized by the applied wall heat flux q", the
liquid mass flux at the test section inlet G, the inlet liquid
temperature Tin and the inlet pressure level Pin, have been
selected for CFD investigation. The data of these operating
points are given in Table 1, where the data set numbers
correspond to naming convention in the paper of Lee et al.
(2008). Operating conditions had been selected such, that
the Set25 represents the case, where the least of all vapor is
produced, while Setl6 represents the case with the largest
amount of vapor.
4. The Physical Model
The flow under investigation is described in the framework
of the currently most conventional CFD approach to
modeling twophase flows with significant volume fractions
of both phases the Eulerian twofluid model derived under
the assumption of interpenetrating continue. Material
properties for both vapor and liquid had been specified by
defining material properties based on IAPWSIF97
water/water steam tables defined for the given range of
temperature and pressure of the testcases. Phase distribution
results from solving the phasespecific continuity equations
for volume fractions, and separate sets of momentum
equations are solved for each phase, where buoyancy and
interfacial momentum transfer is taken into account.
Momentum transport equations are supplemented by
turbulence model equations, where the shear stress transport
model (SST) has been applied to the continuous phase and a
zeroequation disperse phase turbulence model together
with the Sato bubble enhanced turbulence model have been
used to describe the turbulence effects arising from the
bubbly phase (for details see ANSYS Inc., 2009).
For the steamwater bubbly flow an energy equation is
solved for liquid, while for the description of the nucleate
subcooled boiling processes under consideration the vapour
is assumed to be saturated at all time. The exchange of
mass, momentum and heat between phases are modeled
using the correspondent source terms in the phasespecific
balance equations. For the dispersed bubbly flow assumed
for the nucleate subcooled boiling processes the interfacial
momentum transfer is modeled in terms of the Grace drag
force due to the hydrodynamic resistance and the nondrag
forces.
Regarding the consideration of the nondrag forces the
current framework in ANSYS CFX (ANSYS Inc., 2009)
allows for the inclusion of the lift, the wall lubrication, the
virtual mass and the turbulent dispersion force. Further
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
forces can be added by userdefined source terms. Previous
investigations have shown good agreement for simulating
adiabatic airwater bubbly flows (see Frank et al., 2008),
polydispersed airwater and steamwater flows (see Krepper
et al., 2008) as well as recondensing steamwater flows (see
Lifante et al., 2009) using the formulations of Tomiyama
(1995, 1998) for the lift force, the generalized wall
lubrication force formulation by Frank (2005) and the Favre
averaged drag turbulent dispersion force formulation by
Bums et al. (2004). In the present investigation nondrag
forces with the exclusion of the wall lubrication and virtual
mass forces have been applied.
The modeling approach for wall boiling will be described in
a separate section. Once the steam is produced at the wall, it
will be assumed, that the steam is at local pressure
dependent saturation temperature at all time. Further the
steam condensates in the bulk subcooled liquid (TL < Tsat)
with the mass transfer rate per unit volume:
z = max h (Ta T (2)
S HLG
With superheated liquid, fluid is evaporating at the rate:
S=max hLG(TL Tsat)ALGO (3)
HLG
ALG is the interfacial area, and hLG is the interfacial heat
transfer coefficient, calculated according to Ranz and
Marshall (1952):
hLG =kNu =kL(2+0.6Rel2Prl/3) (4)
d, d
This relationship is valid for mass transfer at the interface of
rather small bubbles with diameters well below 0.5mm. So
recently Tomiyama (2009) proposed a slightly changed
relationship for evaporation/condensation at the phase
interface of larger bubbles in the range of dB1,...,3mm:
hLG k Nu k (2+0.15Re Pro5) (5)
dB dB
To close the phase transition model in the bulk bubbly flow,
a phasic characteristic length scale for the mean bubble
diameter dB has to be provided. This can be obtained from
applying a population balance model like
homogeneous/inhomogeneous MUSIG model or a
DQMOM model. Here we follow the simplified approach of
providing a local mean bubble diameter as proposed by
Kurul and Podowski (1991) as well as Anglart et al. (1997),
where both proposed to calculate the local bubble diameter
dB as a linear function of liquid subcooling Tsub:
d dBl (T* b Tub,2) +dB2 (T ub,l ub
B (6)
sub,1 sub,2
For typical nuclear energy applications these authors
proposed for subcooled nucleate boiling under PWR
conditions (so, high pressure conditions) reference bubble
diameters at the two reference subcooling conditions: dB1 =
0.1mm at Tsub,1 = 13.5K and dB2 =2mm at Tsub,2= 5K. The
bubble size in the bulk has a direct influence on the
interfacial area density and on the condensation respective
evaporation rate in the bulk. It is clear, that these
assumptions cannot be applied without reconsideration to
the nucleate subcooled boiling under normal pressure
conditions, since already the experimental observations
show, that Sauter mean bubble diameters of up to 4mm
occurred in the experiments under some of the operating
conditions. Therefore this relationship (5) was a first
candidate for necessary model modifications for bulk and
wall boiling under low pressure conditions, as will be
discussed in a later section.
4.1. Modeling Nucleate Subcooled Boiling at Heated Walls
The current implementation and exposure of the wall
boiling model in the graphical CFD preprocessor of ANSYS
CFX 12.0 had predecessors in earlier versions of ANSYS
CFX as beta model capabilities. Therefore more detailed
descriptions of the wall boiling modeling approach exist
from earlier publications, referring to Egorov et al. (2004)
and Krepper et al. (2007). All this model development
follow the general outline of the socalled wall heat flux
partitioning algorithm developed by Kurul & Podowski
(1991). Since this initial model development was aimed
more on 1d thermohydraulic modeling of the phenomenon,
model enhancements and adjustments were necessary in
various places of the model algorithm formulation in order
to accommodate for the specific requirements of an
implementation into a general 3d CFD solver.
Subcooled boiling is observed at heated surfaces, when the
heat flux applied to the wall is too high to be transferred to
the core flow of liquid by the singlephase convec
tiveconductive mechanisms. The term "subcooled" means,
that the saturation temperature is exceeded only in a local
vicinity of the wall, whereas the average temperature in the
bulk is still below saturation.
The point, where the local wall temperature reaches the
saturation temperature, is considered as the onset of
nucleate boiling. Steam bubbles are generated at the heated
surface at nucleation sites, with the surface density of these
sites depending on different factors, including the wall
superheat. With increasing wall superheat ATsup=TwTsat the
attached bubbles grow and then leave the wall at certain
critical size. This critical size, called bubble departure
diameter, may depend on the surface tension and on the
forces acting on the bubbles from the surrounding fluid.
Heat transfer from the wall is then described as being
carried by turbulent convection of liquid, by transient
conduction due to the departing bubbles and by evaporation.
Distribution of the entire wall heat flux between these
mechanisms (wall heat partitioning) can be calculated by
modeling each mechanism in terms of the nucleation site
density, the size of departing bubbles, their detachment
frequency, and waiting time until the next bubble appears
on the same site (mechanistic modeling approach). This
mechanistic modeling approach of the wall boiling process
is required in the framework of the CFD code, since for
technical applications it is mostly impossible to fully
resolve the microphenomenon of steam bubble formation at
the heated wall on the underlying numerical mesh and with
the applied time scale of integration. Instead the resulting
steam production and enhanced heat transfer to the liquid is
taken into account by the mechanistic model of wall boiling
based on the wall heat flux partitioning algorithm. Once the
steam bubbles are released from the nucleation sites, they
move through the subcooled liquid and condensate,
releasing the latent heat again in correspondence to eq. (1).
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Following the modeling approach of wall heat flux
partitioning, the applied wall heat flux on the heated surface
is split into 3 parts: Qc, the turbulent convective heating of
the liquid, QQ, the quenching heat flux and QE, the
evaporative heat flux:
Q, = Q +QQ QE (7)
As already mentioned, in this model vapour is assumed to
be saturated everywhere, and no part of the wall heat flux is
arranged for superheating of the vapour phase. The heat
partitioning model considers the whole heated wall surface
as being separated into two fractions: a) fraction A2
influenced by the vapour bubbles, formed on the wall and b)
fraction A1 being the remaining wall surface area with
Ai=1A2. The wall area fraction Al represents the part of
the wall surface that is not affected by the growing steam
bubbles. Therefore the wall heat flux for this part of the
surface is modelled in a similar way as for the singlephase
convective heat transfer into pure liquid, by using the
turbulent wall function procedure as outlined in Egorov et al.
(2004). Given that, the convective heat flux can be written
as:
Qc Ah= (T T) (8)
where he is the turbulent heat transfer coefficient, which
depends on the velocity field and is modelled using the
turbulent temperature wall function (see Egorov et al., 2004).
The wall area fraction A2 represents the remaining part of
the surface, which exchanges heat with both phases. The
already mentioned evaporative heat flux QE is consumed for
evaporation of the initially subcooled liquid:
E =m(hG,,, h,) (9)
with:
j3
mh = p 6 a
2
= PG F df (10)
3
=pG w a ,5 dw Naf
S4 }3
resulting in:
QEP = Pcd min r a 5 N fh (11)
3 4 a G
where m is the evaporation mass transfer rate per unit
wall area, A'2F is the nonlimited wall area influenced by
vapour bubble formation, hQsat and hL are the specific
enthalpies of the saturated vapour and subcooled liquid
respectively, dw is the bubble departure diameter, Na is the
nucleation site density and f is the bubble detachment
frequency. The quenching heat flux due to transient vapour
bubble departure and cooling of the wall area A2 by
substituting fresh subcooled liquid is modelled as:
Q =42h(T, T (12)
where hq is the quenching heat transfer coefficient. In the
above relationships the area A2 influenced by the growing
vapour bubbles is related to the nucleation site density and
the bubble departure diameter:
2 2
A2 =min(r a NT,1) (13)
4
where dw is the bubble departure diameter, Na is the
nucleation site density and a is a influence factor introduced
by Kurul & Podowski (1991) and is assumed to be a=2.
In order to arrive at a closed model formulation for the
above wall heat flux partitioning scheme, a larger number of
closure models have to be provided. These are required for
the following model parameters:
Na, wall nucleation site density
bw, bubble departure diameter
f, bubble detachment frequency
hQ, quenching heat transfer coefficient
bubble waiting time
The required closure relationships are provided from
correlations, following in most cases the used correlations in
the original model formulation of Kurul & Podowski (1991),
but providing alternatives or the possibility for the model
user to introduce his own model correlation as a
userdefined relationship instead. For more details on the
different submodels please refer to (ANSYS Inc., 2009).
One particular and rather important correlation used in this
model closure is introduced for the bubble departure
diameter. Here Kurul & Podowski (1991) adopted a
correlation established by Tolubinski & Kostanchuk (1970):
d, = min dref exp 
I r <
A,b d
* ,4.o
The parameters of the original model are dimensional
Grid 1 Grid 2 Grid 3
No. of nodes 6.342 24.682 97.362
No. of elements 20x150 40x300 80x600
y (Set25) 88 45 25
Max. aspect ratio 24 24 24
Max. cell volume 1.01 1.01 1.01
ratio
(dmax=1.4mm, dref=0.6mm, ATref=45K) and ATsub refers to
the local liquid subcooling. These model data are specific
for the model application to nucleate subcooled boiling
under pressurized conditions and need to be revised in case
of model application to different operating conditions.
4.2. Boundary Conditions for the Wall Boiling Model
The implementation of the wall boiling model for nucleate
subcooled boiling in ANSYS CFX 12.0 supports the
specification of either a prescribed wall heat flux or a
prescribed wall temperature at the surface of the heated wall.
Eq. (6) provides in both cases the relationship to predict
either the resulting wall temperature in dependence on the
prescribed wall heat flux or vice versa.
Another, and in practical cases even more interesting
possibility, is the specification of a volumetric energy
source in the solid material of the heater and the prediction
of the heat transfer due to conduction in the solid material
using conjugate heat transfer (CHT) prediction. In this case
both the wall heat flux and the wall temperature are part of
the solution from a coupled simulation of CHT in the solid
material and multiphase flow CFD in the fluid domain of
the application. ANSYS CFX 12.0 supports this type of
simulation with both 1:1 and nonconforming meshes at the
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
fluidsolid interface.
5. Model Validation based on CFD Simulation of the
Testcase of Lee et al. (2008)
Formerly the wall boiling model, described in the previous
section, was already thoroughly validated (e.g. Krepper et
al., 2007) for nucleate subcooled boiling under pressurized
conditions using the data from vertical channel or pipe
experiments with heated outer walls for comparison, e.g.
experiments published by Bartolomej and Chanturiya
(1967) and Bartolomej et al. (1982).
In the present investigation the wall boiling model is
applied to nonpressurized conditions, using the testcase
configuration from the experimental test facility and
experimental data of Lee et al. (2008), as described in
section 3. To identify the range of validity of applied boiling
model and undertaken model parameter changes the
calculated steam volume fraction and steam/water velocity
profiles were compared with the large set of measurements
at different operating conditions.
5.1. CFD Simulations with Prescribed Wall Heat Flux
5.1.1. CFD Geometry and Numerical Meshes
Due to the radial symmetry of the circular annulus of the
experimental test facility of Lee et al. (2008) the CFD
simulations have been carried out twodimensional on a 1
degree symmetry sector of the geometry. A hierarchy of
consistently refined hexahedral meshes have been generated.
In all cases the meshes had 1 mesh cell in circumferential
direction and mesh nodes in radial and axial direction of the
circular annulus were uniformly distributed. Main meshing
parameters are listed in Table 2, where the dimensionless
wall distance y+ is evaluated at the heater wall of the
configuration.
Table 2: Mesh hierarchy for the CFD investigation.
5.1.2. Boundary conditions
Fig. 2 shows the schematic of the 1 degree slice of the
circular annulus representing the fluid domain for the CFD
simulation of the testcases of Lee et al. (2008). As shown,
symmetry boundary conditions are applied to both sides of
the symmetry segment, the outer wall of the circular annulus
was set to adiabatic and a constant wall heat flux boundary
condition has been applied to the surface of the centralized
heater rod.
For the inlet boundary conditions it was assumed, that the
test section of Lee et al. allows for fully developed turbulent
inflow conditions. In order to account for that, single phase
isothermal water flow at specified pressure, temperature and
mass flow rate through the same 1670mm long test section
was previously calculated, applying the SST turbulence
model. The boiling flow simulation was afterwards
initialized with the profile data of all velocity components,
turbulent kinetic energy and turbulent eddy dissipation from
the outlet cross section of this predecessor singlephase flow
calculation. The same data where prescribed as inlet
boundary conditions of the liquid phase for the boiling flow
simulation. For the steam phase a very small volume
fraction of 1015 and the saturation temperature
corresponding to the prescribed pressure level including
hydrostatic pressure in the circular annulus of the Lee et al.
testcase geometry was prescribed at the inlet cross section.
The reference pressure of the individual calculations was set
to the corresponding pressure as specified by Lee et al. for
the individual datasets and an average static pressure outlet
boundary condition was applied.
Figure 2: Geometry configuration and boundary con
ditions of CFD simulations with prescribed wall heat flux.
5.1.3. Wall Boiling Submodels and Parameter Modifications
As already mentioned, the wall boiling model depends on a
larger number of submodels and model parameters, where
most of them had been derived for wall boiling processes
under pressurized conditions. In the present investigation
the following settings for the submodels of the wall boiling
model have been used:
Wall nucleation site density: Lemmert & Chawla
model
Liquid quenching heat transfer coefficient: Del
Valle Kenning model
Bubble detachment frequency: Cole model, derived
from terminal bubble rise velocity over the bubble
departure diameter
Bubble waiting time: Tolubinski & Kostanchuk
model, which sets the bubble waiting time to 80%
of the time between bubble departures
Bubble diameter influence factor: default value of
2.0
Fixed y+ for evaluation of liquid subcooling
temperature from turbulent wall functions: default
value of 250
For details of the named submodels please refer to ANSYS
Inc. (2009). Special attention has been directed to the
modelling of the bubble diameter in the bulk liquid dB, the
bubble departure diameter dw and the maximum area
fraction of bubble influence A2,max.
The default for the modelling of the bubble diameter in the
bulk liquid under pressurized conditions is the relationship
introduced by Kurul & Podowski (1991) in accordance with
eq. (5), relating the bubble diameter in the bulk liquid in a
piecewise linear relationship to the local liquid subcooling
temperature with default model parameters resulting in a
maximum bubble diameter of dB<1.48mm (see Fig. 4).
From the experimental data from Lee et al. (2008) it can be
seen, that e.g. under the experimental conditions of the
dataset Set25 maximum bubble diameters of about
dBmx=2.8mm have been measured. Under different
operating conditions even larger bubble diameters of up to
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
dB,max~4.2mm have occurred in the experiments under
normal pressure conditions. Partially the large bubble
diameters in the bulk liquid can result from bubble
coalescence due to the higher steam volume fraction in the
vicinity of the heater wall. Unfortunately without applying a
population balance model in combination with the wall
boiling model we cannot account for these bubble
coalescence effects. Instead it was tried to take into account
the larger observed bubble diameters, which are especially
of importance for the rate of steam recondensation in the
subcooled liquid, carrying out a parameter study with the
settings from Table 3. Following a suggestion of Egorov et
al. (2004), the original piecewise linear relationship of
Kurul & Podowski (1991) was additionally smoothed by
applying tanh functions as smoothing function. The resul
ting bubble diameters in the bulk liquid in dependence on
the local liquid subcooling can be seen in Fig. 4 for the ori
ginal Kurul & Podowski law and the made modification 3.
Modification dB1 Tsub dB2 Tsub,2
Kurul & Po 0.1mm 13.5K 2.0mm 5K
dowski (1991)
dB, Mod. 1 0.15mm 13.5K 4.0mm 5K
dB, Mod. 2 0.15mm 13.5K 4.0mm 5K
dB, Mod. 3 0.15mm 25.0K 4.0mm 5K
Table 3: Parameter modifications regarding the modeling of
the bubble diameter in the bulk liquid.
Modification dw
Tolubinski & Kostanchuk (1970) from eq. (13)
dw, Mod. 1 dw=lmm
dw, Mod. 2 dw=2mm
dw, Mod. 3 dw=3mm
dw, Mod. 4 from eq. (13)
multiplied by 4.0
Table 4: Parameter modifications regarding the modelling
of the bubble departure diameter.
A direct consequence of the modifications made for the
modelling of the bulk bubble diameter is a resulting change
of the bubble departure diameter, since both bubble
diameter laws should give consistent values close to the
surface of the heater. Furthermore it is even more important,
that the bubble departure diameter is of main influence on
the area of bubble influence on the heated wall, bubble
departure frequency and thereby finally on the evaporation
rate and steam production resulting from the wall boiling
model. For the given operating conditions of Set25 the
standard relationship of Tolubinski & Kostanchuk (1970)
results in a bubble departure diameter of about 0.5mm,
which seems to be much too small in comparison to the near
wall bubble diameter measurements. Therefore a parametric
study has been carried out for the bubble departure diameter
dw with the different model settings as listed in Table 4.
Finally, in a number of previous investigations the wall area
influenced by the vapour bubbles A2 was limited by a
maximum area fraction of bubble influence A2 ,mx=0.5 by
default. However it was found, that for the present testcase
with the increased bubble departure diameters this artificial
limitation of A2 results in a rather early limitation of the
evaporative heat flux and thereby the amount of produced
steam volume fraction from the wall boiling model was
substantially underpredicted in comparison to the
experiments. Therefore this limitation was raised by setting
A2,max 1.0 throughout the carried out CFD simulations.
5.1.4. Results and Comparison to Experimental Data
First CFD investigations based on the specification of the
Set25 dataset from the Lee et al. (2008) experiments were
related to investigation of numerical parameters and the
adherence to CFD best practice guidelines. Simulations had
been carried out in steadystate, applying the ANSYS CFX
highresolution scheme (2nd order TVD scheme) as
advection scheme for the hydrodynamic system of equations
and 1st order upwind scheme for the turbulence model
equations. For convergence control a RMS residual criteria
of 104 has been applied. Additionally global mass,
momentum and energy as well as temperature and volume
fraction data in 60 axially and radially distributed
monitoring point locations have been monitored and were
found to finally result in stationary values. Furthermore it
was investigated and found, that the steadystate solution
algorithm requires for the desired level of convergence
integration time scales in the range of At0.2,...,100ms,
which is case and mesh resolution dependent. For the mesh
level 3 and Set25 operating conditions an integration time
scale of 0.4ms was required.
In a first series of CFD simulations the influence of the
parametric change for the bubble diameter in the bulk liquid
has been investigated. Fig. 5 shows the resulting radial
distributions for the vapour volume fraction and bubble
diameter at the measurement location in comparison to the
data of Set25 (the local bubble diameter being a direct result
from the Kurul & Podowski law in accordance to eq. (3)
and the local liquid subcooling temperature Tsub=TsatTL). It
can be seen from Fig. 5a), that the changes in the modelling
parameters for the 3 modifications of the bulk bubble
diameter have led to a decrease in the maximum value of
the steam volume fraction directly at the heater surface and
to an increase in the steam volume fraction further apart
from the heated wall. The slope of the steam volume
fraction profiles are in better agreement with the
experimental data for higher bulk bubble diameters closer to
the experimentally measured values, as can be expected.
Nevertheless the total amount of produced steam is still
underpredicted and it can be seen from the last 2
modifications that a further increase in bubble bulk diameter
seems not to have a further influence. With the 3rd
modification of the bulk bubble diameter the results for the
predicted bubble diameters are now in the same order as the
measured values (see Fig. 5b).
Setting dB=2mm and A2 max1.0 a second series of
parameter investigation has been carried out for Set25
operating conditions. Fig. 6 a)c) show the CFDexperiment
comparison of profiles of vapour volume fraction, axial
water and vapour velocities at the measurement cross
section. With the 3rd modification of the bubble departure
bubble diameter in accordance to Table 4 and Fig. 4 a
reasonable good agreement between the numerically
predicted and measured vapour volume fraction profiles
could be established (Fig. 6a). Water velocities are in good
agreement to measured data in a certain distance from the
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
heated wall (Fig. 6b). But close to the heated wall the water
velocities from the CFD simulation clearly show an axial
acceleration due to buoyancy and the applied free slip wall
boundary condition of the vapour phase, resulting in a larger
difference to the measured water velocities. A
corresponding difference between predicted and measured
vapour velocities can be observed from Fig. 6c) in the
vicinity of the wall, while the comparison to data in a
certain distance from the wall is quite satisfactory. An
explanation of the observed differences is, that in the reality
of a wall boiling process on a vertical heated surface the
steam bubbles initially grow on the heated surface without
moving. Later a sliding motion of the growing steam
bubbles has been observed in experiments for bubble
diameters close to the departure diameter and if the
buoyancy force on the growing steam bubble becomes large
enough. Nevertheless the sliding velocity of these steam
bubbles is substantially slower than the rising velocity of
steam bubbles in free stream motion, which can explain the
observed differences in vapour velocity. A changed wall
boundary condition based on an estimated bubble sliding
velocity could contribute to improvement of this
comparison. On the other hand side some doubt regarding
the measured vapour velocities is remaining, where no
buoyancy effect from the approx. 30% vapour volume
fraction close to the wall is remarkable in the data.
Finally Fig. 7 shows the result of a grid convergence study,
where the final model parameter settings for Set25 have
been applied to the series of even 4 subsequently refined
meshes. The radial volume fraction profiles show a
tendency to a slight further increase in the predicted vapour
volume fractions, also the difference between mesh 3 and 4
becomes rather small and is in fair agreement to the
experimental data from lee et al. (2008).
Figure 3: Geometry configuration and boundary con
ditions of CFD simulations with prescribed volumetric
energy source in the solid core material of the heating rod.
5.2. CFD Simulations with Prescribed Volumetric Energy
Source in the Solid Core Material of the Heating Rod
An alternative to the specification of the wall heat flux or a
given wall temperature is the specification of a volumetric
energy source in the solid material of the heater. In fact in
many cases this is the more desirable simulation approach,
since e.g. the thermal energy output of the heater is known
whereas the wall temperature and the wall heat flux are
rather sought as part of the solution. In ANSYS CFX 12.0
the prediction of conjugate heat transfer (CHT) in solid
materials has been interfaced with the heat partitioning
algorithm of the wall boiling model. The combination of
both physical models is supported on fluidsolid interfaces
with 1:1 and nonconforming meshes (GGI General Grid
Interface).
5.2.1. Numerical Meshes and Boundary Conditions
In order to demonstrate the capability of ANSYS CFX 12.0
to combine CHT and wall boiling model, a corresponding
setup for the heating rod with a volumetric energy source
term has been set up. For the operating conditions of the
dataset Set25 the prescribed wall heat flux of q"=220 kW
m2 corresponds to a volumetric energy source of E=82.3
MW m3. Hereby the volumetric energy source has been
specified in an inner core of the solid material of the heating
rod with a radius Rcore=0.75 Ro=7.125mm, as can be seen
from Fig. 3. The outer nonheated cladding has been
specified in this case as 0.75 Ro
thereby a solidsolid interface between the homogeneously
heated core and the nonheated cladding material and a
fluidsolid interface between the cladding and the fluid
domain, where the wall boiling will occur. Both core and
cladding solid materials have been set to steel properties.
Symmetry boundary conditions apply to the boundary
patches in circumferential direction (front and rear), the top
boundaries of the core and cladding domains have been set
to adiabatic, while for the bottom boundaries of both solid
domains an isothermal boundary condition with the fluid
inlet temperature of the dataset under consideration has
been applied.
Investigations and comparisons to existing CFD solutions
from section 5.1 and the experimental data have been
carried out on all 3 mesh levels using conforming meshes
and on mesh level 3 only for nonconforming mesh
interfaces. In this case the following fluidsolid domain
meshing combinations have been investigated:
Configuration 1 (mesh levels 13):
Fluid domain mesh from Table 2;
Core & Cladding solid hexahedral meshes with
same mesh resolution as the fluid domain;
all mesh interfaces are treated as 1+1 conforming
meshes with direct data transfer
Configuration 2:
same mesh as grid level 3 of Configuration 1;
all mesh interfaces are treated as GGI interfaces'
Configuration 3:
Fluid domain 80x603 cell hexahedral mesh;
Core solid 80x600 cell hexahedral mesh;
Cladding solid 80x600 cell hexahedral mesh;
all mesh interfaces are treated as GGI interfaces
5.2.2. Results and Comparison to Experimental Data
Fig. 8 shows the resulting solid/fluid temperature and steam
volume fraction distribution in the vertical plane for the
CFD simulation with a prescribed volumetric energy source
in the heated core of the heater rod. From the left hand side
image it can be seen, that the temperature distribution in the
solid material is quite homogeneous over the height of the
heating rod, while in the fluid domain the fluid temperature
is continuously increasing from the inlet to the outlet cross
section. The latter corresponds to a continuously increasing
amount of steam being produced at the heated surface of the
1 GGI General Grid Interface
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
heater rod (right image).
For the measurement cross section at LM=1610mm the solid
and fluid temperatures have been compared for
Configuration 1 on three consistently refined meshes (see
Fig. 9). From this figure it can be seen, that the resulting
solid temperatures are almost identical on all 3 meshes,
while the fluid temperatures slightly differ. Additionally in
the fluid domain the fluid temperatures are compared to the
CFD result from the simulation with the prescribed constant
wall heat flux, where only a rather minor difference in the
results could be observed. Fig. 10 shows the same
comparison for results on 3 different meshes of
Configuration 1 for the radial steam volume fraction
distribution at the measurement cross section with the
additional comparison to experimental data. Here the
resulting steam volume fraction on mesh level 3 for the
CHT Configuration 1 leads to even a slightly higher steam
volume fraction in a certain distance from the wall, what
corresponds well to the slightly higher liquid temperatures
from Fig. 9 (and therefore a slightly lower steam
recondensation rate at this radial location).
Finally CFD results from the 3 different meshing
configurations defined in section 5.2.1 using two different
approaches for the data transfer between meshes on both
sides of a meshing interface (1+1 direct data exchange and
GGI) have been compared to each other in Fig. la)b). The
diagrams in Fig. 11 show that identical results for radial
temperature and steam volume fraction distributions are
obtained with all three meshing & data exchange methods,
proofing the correctness of the corresponding algorithmic
implementations for the combination of CHT and the wall
boiling model inANSYS CFX 12.1.
Conclusions
A wall boiling model for the CFD prediction of nucleate
subcooled boiling has been implemented in ANSYS CFX
12.0 together with the necessary underlying submodels for
model closure and GUI support in the ANSYS CFX
physical preprocessor. In earlier publications (see Egorov et
al. (2004), Krepper et al. (2007, 2008)) the model had been
successfully validated for nucleate subcooled boiling under
pressurized conditions. In the present work, the model has
been applied to the validation experiment of Lee et al.
(2008), which is aimed to nucleate subcooled boiling under
normal pressure conditions in flow geometry with
dimensions very similar to the arrangement of heated rods,
which can be found in nuclear reactor fuel assemblies. In
this paper the model parameters have been identified, which
have to be changed for a successful CFD prediction of wall
boiling and steam production under the operating conditions
of the selected testcase. Special model modification with
regard to the prediction of the local liquid subcooling
temperature in the nearest wall mesh cell lead to a
stabilization of the model and thereby to almost grid
independent results of the CFD prediction of the wall
boiling process. The comparison to the experimental data of
Lee et al. (2008) led to satisfactory agreement, taking into
account the measurement accuracy in the experiments and
the remaining uncertainties in the correlations used in wall
boiling model closure.
Furthermore the wall boiling heat partitioning algorithm has
been coupled to the prediction of conjugate heat transfer
(CHT) in the solid material of the heated boundary of the
fluid domain, allowing for the definition of volumetric
energy sources in e.g. the rods of nuclear fuel assemblies
heated by nuclear fission reaction instead of fluid domain
thermal boundary conditions. The model implementation
supports both conforming and nonconforming meshes at
domain interfaces. In the validation exercise it has been
proven that the different meshing and data exchange
approaches lead to identical results of the CFD simulation.
Acknowledgements
This research has been supported by the German Ministry of
Economy (BMWi) under contract number 150 1328 in the
framework of the German CFD Network on Nuclear
Reactor Safety Research and Alliance for Competence in
Nuclear Technology, Germany.
References
Anglart, H., Nylund, O., Kurul, N., Podowski, M.Z.: "CFD
prediction of flow and phase distribution in fuel assemblies
with spacers", NURETH7, Saratoga Springs, New York,
1995. Nuclear Eng. & Design (NED), Vol. 177 (1997), pp.
215228.
ANSYS Inc., ANSYS CFX 12.0: Users Manual, 2009
Bartolomej, G.G, Chanturiya,V.M.: "Experimental study of
true void fraction when boiling subcooled water in vertical
tubes", Thermal Engineering, Vol. 14 (1967), pp. 123128
(translated from Teploenergetika, no. 2, Vol. 14 (1967), pp.
8083).
Bartolomej, G.G., et al.: ,,An experimental investigation of
true volumetric vapour content with subcooled boiling in
tubes", Thermal Engineering, Vol. 29 (1982), pp. 132135
(translated from Teploenergetika, no. 3, vol. 29, 1982,
pp.2023).
Bums A.D., Frank Th., Hamill I., Shi J.M.: "The Favre
Averaged Drag Model for Turbulent Dispersion in Eulerian
MultiPhase Flows", 5th International Conference on
Multiphase Flows, ICMF'2004, Yokohama, Japan, May 31 
June 3, 2004, paper No. 392, pp. 117.
Egorov, Y, Menter, F.: "Experimental Implementation of the
RPI Wall Boiling Model in CFX5.6", Technical Report
ANSYS/TR0410, ANSYS Germany GmbH (2004).
Frank Th.: "Progress in the numerical simulation (CFD) of
3dimensional gasliquid multiphase flows", 2. NAFEMS
CFDSeminar "Simulation komplexer Stromungsvorgfinge
(CFD) Anwendungen und Entwicklungstendenzen",
Wiesbaden, Germany, 25.26. April 2005, pp. 118.
Frank Th., Zwart PJ., Krepper E., Prasser H.M., Lucas D.:
"Validation of CFD models for mono and polydisperse
airwater twophase flows in pipes", J. Nuclear Engineering
& Design (NED), Vol. 238, pp. 647659, March 2008.
Krepper E., Koncar B., Egorov Y.: "CFD modeling of
subcooled boiling Concept, validation and application to
fuel assembly design", Nuclear Engineering and Design
(NED), Vol. 237, pp. 716731, 2007.
Krepper E., Lucas D., Frank Th., Prasser H.M., Zwart PJ.:
"The inhomogeneous MUSIG model for the simulation of
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
polydispersed flows", Nuclear Engineering and Design
(NED), Vol. 238, pp. 16901702, 2008.
Kurul N., Podowski M.Z., On the Modelling of
Multidimensional Effects in Boiling Channels, ANS Proc.
27th National Heat Transfer Conference, Minneapolis, July
2831, 1991.
Lee T.H., Park G.C., Lee D.J.: ,,Local flow characteristics
of subcooled boiling flow in a vertical concentric annulus",
International J. Multiphase Flows, Vol. 28 (2002), pp.
13511386.
Lee T.H., Yun B.J., Park G.C., Hibiki T., Kim S.O.:
"Local flow structure of subcooled boiling flow of water in
a heated annulus", Proceedings of the 16th International
Conference on Nuclear Engineering (ICONE16), May
1115, 2008, Orlando, Florida, USA, Paper No.
ICONE1648170, pp. 12.
Lee T.H., Situ R., Hibiki T., Park H.S., Ishii M., Mori M.:
"Axial developments of interfacial area and void
concentration profiles in subcooled boiling flow of water",
Int. J. Heat and Mass Transfer, Vol. 52 (2009), pp. 473487.
Lifante C., Frank Th., Bums A.D.: "Prediction of Poly
disperse Steam Bubble Condensation in Subcooled Water
using the Inhomogeneous MUSIG Model", ANSYS Confe
rence & 27. CADFEM Users Meeting, Congress Center
Leipzig (CCL), 18.20. November 2009, pp. 130.
Ranz, W.E., Marshall, W.R.: "Evaporation from drops",
Chem. Eng. Progr., Vol. 48 (3), pp. 141146 (1952).
Tolubinsky, VI., Kostanchuk, D.M.: "Vapor bubbles growth
rate and heat transfer intensity at subcooled water boiling",
Heat Transfer, Preprints of papers presented at the 4th
International Heat Transfer Conference, Vol. 5, Paris (Paper
No. B2.8), 1970.
Tomiyama A., Sou A., Zun I., Kanami N., Sakaguchi T.:
"Effect of Eotvos number and dimensionless liquid
volumetric flux on lateral motion of a bubble in a laminar
duct flow", Advances in Multiphase Flow, p. 35, Elsevier
Science (1995).
Tomiyama A.: "Struggle with computational bubble
dynamics", ICMF'98, 3rd Int. Conf. Multiphase Flow
(ICMF), Lyon, France, 1998.
Tomiyama A.: Private communication (2009).
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
dI
_ *Kurul & Podowski (orig.)
.' Kurul & Podowski (smoothed)
N : . .db, Mod. 3
db. Mod. 3 (smoothed)
ld
______ j
0 Tsub2 10
20 Tsub, 30
Tb [K]
Figure 4: Bubble diameter in the bulk fluid in dependence on the local liquid subcooling temperature: Kurul & Podowski
(1991) vs. 3rd modification of the bubble bulk diameter relationship.
0.
E 0.
E u
b .
S
S0.
E 0.
10
E 0
0.
30 I
Kimd & PodaMkll
25  =Mod. 1
EF=Moa. 2
20
 =Mod. 3
15  Expernental Data
15 
10 D
05
105 "T___      
00 I I
0 0.1 0.2 0.3 0.4 0.5 (.6 0.7 0.8 0.9 1
Radial Position Rp [
0.004 Kun & PadowkI  b=Mod. 1
0.0035 = ............ ...... db=Mod. 2 l=Mod. 3
0.003 OExpuunfllDaa
0.002
nnno.;  
0.001
0.0005
0
S
0.1 0.2 0.3 0.4 0.5 0.6 0.7
Radial Position Rp []
0.8 0. i
Figure 5: a) Steam volume fraction and b) bubble bulk diameter radial profiles for Set25 operating conditions on mesh
level 3 for the different parameter variations in the bubble bulk diameter dependence on local liquid subcooling temperature.
I I i' , .
~Zh
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
a)
0
o
L.
E
45
IK &Podo Packi
AMud. Z dWl nm
S&=Mod. 2 dWl=2mm
0 0 =Mod.5 Zd. 2; 0=3m
 ExIpai.t la
0 01 0.2 03 04 05 06 07 08 09 1
Radial Position Rp []
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
Radial Position R, []
1.4
I1.2
0.8
? 0.6
E
0.4
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7
Radial Position Rp []
0.8 0.9 1
Figure 6: a) Steam volume fraction, b) water axial velocity and c) steam axial velocity radial profiles at measurement
location for Set25 operating conditions on mesh level 3 for the different parameter variations in the bubble bulk diameter dB
and bubble departure diameter dw.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
WAW
0.35
C 0.30
0
0.25
1Gridl; amMod. 2; dV&2nm
L 0.20
U1 Grid2; dBMod. 2; d=2mnm
E 0.15
.2 Grid3; c=Mod. 2; dVWnmm
0.10 Grid:; ldMod. 2; dW=2mm
E 0.05 ExperimentalData
S o.oo ^ :  
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Radial Position Rp []
Figure 7: Investigation of grid independence of solution for steam volume fraction at measurement location for Set25
testcase and given model parameter changes for bubble bulk diameter and bubble departure diameter.
0 0.01 (m) 0 0.01 (m)
a 
Figure 8: Solid/fluid temperature (left) and steam volume fraction (right) distribution in a vertical cross section of the
heating rod and circular annulus of the Lee et al. testcase configuration.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
II
I
20.00
Radius [mm]
Figure 9: Comparison of temperature distributions in the solid materials of the heated core, the cladding material and water
temperatures in the adjacent fluid domain at the measurement cross section of the Lee et al. testcase, Set25. Comparison to
CFD results from simulation with specified wall heat flux at the heated wall (blue line).
040
0 35 Set 25, Exp Data
3 Grid3(wall heatflux)
,2 030 Grid1 volumetricc heat source)
\  Grid2 volumetricc heat source)
025 
S 0  Grid3 volumetricc heat source)
I 020
015
E
010
0 05
000
900 11 50 14 00 16 50 1900
Radius [mm]
Figure 10: Comparison of experimental data for radial steam volume fraction distribution to CFD results for the Lee et al.
testcase under Set25 operating conditions with specified wall heat flux at the heated wall and with specified volumetric
energy source in the solid material of the heater.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
415.0
410.0
435.0
430.0
335.0
3300
375.0
370.0
335.0
aiona
i FluiiTemperature (walllhestflw)
 FluidTemrperature (conformingmesh, 1 1)
+ FluidTemperature (conforming mesh, 001)
FluicTemperanure (noncon'ormingmenl, G01)
 SoliTenperature (conforming mesh,1 1)
SolldTenperaure (conforming mesh, GGI)
X SolidTenperatures(nonconforming mesh,O0)
5.00 10.OD 15.00
Radius [mm]
20.00
Figure 11: Comparison of a) radial temperature distributions in the solid materials of the heated core, the cladding material
and water temperatures in the adjacent fluid domain and b) radial steam volume fraction distribution at measurement location
for the three different configurations of conforming and nonconforming meshes with 1+1 and GGI data connection at the
meshing interfaces of adjacent domains.
= &.
L"
I
 Heated Solid
Core
040
Set25, Exp Data
035 
SSteamVF (wall heat flux)
C 0 SteamVF (conforming mesh,1 1)
S 030
0 <* SteamVF (conformingmesh, GGI)
0 25 i SteamVF (nonconformingmesh, GGI)
U
E 020
0 15
E
(f
.2 010
0 05
 10   
0 00
000 020 040 060 080 1 00
Radius [mm]
I I
I I
1 1
I a
