Krepper2007-CFDMOdelingSubcooled Boiling.pdf

by user

Category: Documents





Krepper2007-CFDMOdelingSubcooled Boiling.pdf
Nuclear Engineering and Design 237 (2007) 716–731
CFD modelling of subcooled boiling—Concept, validation
and application to fuel assembly design
Eckhard Krepper a,∗ , Bostjan Končar b , Yury Egorov c
Forschungszentrum Rossendorf e.V. (FZR), Institute of Safety Research, P.O.Box 510119, 01314 Dresden, Germany
b Jožef Stefan Institute Jamova 39, 1000 Ljubljana, Slovenia
c ANSYS, Germany GmbH, CFX-Development, Staudenfeldweg 12, D-83624 Otterfing, Germany
Received 14 July 2006; received in revised form 19 October 2006; accepted 19 October 2006
The paper describes actual Computational Fluid Dynamics (CFD) approaches to subcooled boiling and investigates their capability to contribute
to fuel assembly design. In a prototype version of the CFD code CFX a wall-boiling model is implemented based on a wall heat flux partition
algorithm. It can be shown, that the wall boiling model is able to calculate the cross sectional averaged vapour volume fraction of vertical heated
tubes tests with good agreement to published experimental data. The most sensitive parameters of the model are identified. Needs for more detailed
experiments are established which are necessary to support further model development. The model is applied for investigation of the phenomena
inside a hot channel of a fuel assembly. Here the essential phenomenon is the critical heat flux. Although subcooled boiling represents only a
preliminary state towards the critical heat flux occurrence, essential parameters like swirl, cross flow between adjacent channels and concentration
regions of bubbles can be determined. By calculating the temperature of the rod surface the critical regions can be identified which may later on
lead to departure from nucleate boiling and possible damage of the fuel pin. The application of up-to-date CFD with a subcooled boiling model for
the simulation of a hot channel enables the comparison and the evaluation of different geometrical designs of the spacer grids of a fuel rod bundle.
© 2006 Elsevier B.V. All rights reserved.
1. Introduction
Subcooled flow boiling occurs in many industrial applications
and is characterized by large heat transfer coefficients. 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. The integrity of fuel rods of a fuel
assembly of a nuclear reactor core can be endangered by exceeding the critical heat flux. The critical heat flux depends on the
flow parameters and can be influenced by the geometrical design
of the fuel assemblies. Especially the spacer grids equipped with
mixing vanes play an important role to increase the permissible heat flux. The verification of design improvements and
their influence on the critical heat flux require expensive experiments. Therefore the supplementation or even the replacements
Corresponding author. Tel.: +49 351 260 2067; fax: +49 351 260 2383.
E-mail address: [email protected] (E. Krepper).
0029-5493/$ – see front matter © 2006 Elsevier B.V. All rights reserved.
of experiments by numerical analyses are of relevant interest in
fuel assembly design.
In the past lots of empirical correlations for critical heat flux
were developed fitted to tests and applied in purpose-developed
1D codes. But these correlations are valid only in the limited
specified region of fluid conditions and fluid properties and hold
for a defined geometry. Using lookup tables based on experiments, the fluid parameter validity range can be extended. But
also this method is limited only to a defined geometry. The
independence on the geometry can only be achieved by the application of CFD methods. The paper describes actual attempts
applying CFD to solve the described task and investigates their
capability to contribute to fuel assembly design.
2. Physical model
Currently the most conventional CFD approach to modelling two-phase flows with significant volume fractions of both
phases is the Eulerian two-fluid framework of interpenetrating
continua. Phase distribution results from solving the phase-
E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731
bubble influence factor
wall fraction cooled by single-phase convection
wall fraction cooled by quenching
interfacial area density
specific heat capacity of liquid
specific heat capacity of vapour
bubble departure diameter on the wall
mean bubble diameter in the bulk flow
bubble departure frequency
gravitational acceleration
interfacial heat transfer coefficient
evaporation enthalpy
coefficient for heat transfer by quenching
liquid heat conductivity
vapour heat conductivity
active nucleation site density
Nusselt number
Prandtl number of liquid
heat flux by evaporation
heat flux by single-phase convection
heat flux by quenching
total heat flux
Reynolds number
near-wall liquid temperature
staturation temperature
Tsub = Tsat − TL liquid subcooling
Tsup = TW − Tsat superheating
wall temperature
liquid characteristic near-wall temperature
waiting time
non-dimensional distance to the wall
vapour velocity
water velocity
bubble lift force
wall lubrication force
turbulent dispersion force
Greek symbols
volume fraction of bubbles
vapour density
liquid density
wall shear stress
specific continuity equations for volume fractions, and a separate
set of momentum equations is solved for each phase. For the
steam–water bubbly flow an energy equation is solved for liquid, while vapour is assumed to be saturated everywhere. The
exchange of mass, momentum and heat between phases are modelled using the correspondent source terms in the phase-specific
balance equations. For the dispersed bubbly flows the interfacial momentum transfer is modelled in terms of the drag force
due to the hydrodynamic resistance and the non-drag forces.
The consideration of the non-drag forces namely the lift, the
wall lubrication and the turbulent dispersion force showed good
agreement simulating air water bubbly flow (see Lucas et al.,
2004). Formulation of a particular two-fluid model therefore
consists of analytical or empirical correlations, used for calculating the interfacial forces, as well as the heat and mass fluxes,
as the functions of the average flow parameters. Since most of
these correlations are problem-specific, the range of their validity has to be kept in mind. The entire model has to be validated
against experiments.
A model of the internal bubbly flow with wall boiling, used
in this work, is outlined in the rest of this section.
2.1. Modelling of subcooled boiling at a heated wall
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 single-phase convective–conductive
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. A 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. Further
downstream the attached bubbles grow and then leave the wall
at certain critical size. This critical size may depend on the surface tension and on the flow regime of 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 modelling 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 modelling approach).
When steam bubbles move through the subcooled liquid, they
condense, releasing the latent heat.
The saturation temperature Tsat and the evaporation heat HLG
were specified for the system pressure level and kept constant in
all calculations, presented below. This assumption holds good
for high pressures typical for the nuclear reactor applications.
The model however can also operate with the locally calculated
Tsat and HLG . For water and steam, the phase-specific densities
ρL and ρG , the heat conductivity coefficients kL and kG , and the
heat capacities CpL and CpG are specified at saturation conditions for the purposes of this work, while the solver supports the
general temperature- and pressure-dependent properties as well.
Wall and liquid temperature values participate in several
constitutive model correlations in a relative form of the wall
superheat above saturation Tsup =TW − Tsat , and the liquid subcooling below saturation Tsub =Tsat − TL . Here TW is the local
wall temperature, and TL is the local liquid temperature of the
near-wall cell. Instead of local liquid temperature TL , a characteristic temperature of liquid TLW may be used to define the
magnitude of subcooling in some correlations, see later discussion.
E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731
In technical internal flows the controlled value is the heating intensity, i.e. the wall heat flux. The given external heat flux
Qtot , applied to the heated wall, is used for heating of the subcooled liquid up to and slightly above the saturation temperature
(sensible heat), and for vapour generation (latent heat). The sensible heat is commonly modelled as being transported by the
turbulent convection and also by the additional liquid mixing
mechanism due to the bubbles, periodically emerging from the
wall. The latter is termed quenching and represents the effects of
transient conduction through the patches of the fresh bulk liquid,
coming to the wall to replace each departing steam bubble. The
quenching heat flux is noticeable on the wall area fraction A2F ,
influenced by steam bubbles, whereas the convection heat flux
dominates on the remaining part of the wall A1F = 1 − A2F . The
overall heat balance at the wall, the so-called heat partitioning,
is usually written as a sum of the three parts:
The following model describes the nucleate wall boiling process
as the periodic release of vapour bubbles by each nucleation site:
Qtot = QF + QQ + QE
ṁW = ρG
where QF , QQ and QE are the heat flux components due to
single-phase turbulent convection, quenching, and evaporation,
respectively. The turbulent convection heat flux is calculated in
the CFX model version using the same temperature wall function as that for the pure liquid flow without boiling (Egorov and
Menter, 2004), but multiplied by the area fraction A1F :
QF = A1F
ρL CpL uW
(TW − TL ),
Ty++ (L)
is analytically calculated non-dimensional temperature (Kader, 1981; Wintterle, 2004) at the non-dimensional
distance from the near-wall cell y+ (L) and uW are the friction
velocity. Upgrading of the boiling model with wall functions,
specially customised for bubbly flows, will be the subject of
further research (Končar et al., 2005b).
To close the model, the remaining heat flux components QQ
and QE , as well as the area fraction A2F , have to be modelled
as functions of local flow parameters, including Tsup and Tsub .
Also a way of estimating the local characteristic temperature of
liquid TLW in Tsub must be defined. After that the Eq. (1) can
be solved iteratively for the local wall temperature TW , which
satisfies the wall heat flux balance.
Similar to the model by Kurul and Podowski (1991), the CFX
model represents QQ in terms of the quenching heat transfer
coefficient hQ
To obtain the grid independent solution for QQ and taking into
account the self-similarity of non-dimensional temperature profiles at different y+ , the characteristic temperature of liquid TLW
in Eq. (3) can be calculated as:
TLW = TW −
Ty++ (LW)
Ty++ (L)
(TW − TL ).
The evaporation heat flux QE is obtained via the evaporation
mass flux on the wall
QE = ṁW HLG .
Here dW is the bubble departure diameter, Na the surface density
of active nucleation sites, f the departure frequency, and a is the
bubble influence factor. A value of 2 is taken for the bubble
influence factor a, participating in the area fraction A2F (7). The
effect of area overlapping between the neighbouring sites is here
not modelled, therefore the evaporation rate ṁW is calculated
using the area-limited variant of the relation (6):
2 dW
A2F Na f
3 a2
In this work basically the same correlations for dW , Na , f, a, and
hQ are applied as those used by Kurul and Podowski (1991). The
active nucleation site density Na depends on the wall surface
roughness, on the wettability of the solid-liquid pair and above
all on the wall superheat. For water and steel Na is correlated
according to Lemmert and Chawla (1977) depending solely on
the wall superheat as
Na [m−2 ] = (185 Tsup [K])1.805
Ty++ (L)
QQ = A2F hQ (TW − TLW )
Na f
(adW )2
= min π
ṁW = ρG
For the bubble departure diameter dW , a following correlation
was applied (Tolubinsky and Kostanchuk, 1970):
, 1.4 [mm]
dW = min 0.6 [mm] exp
45 K
The model is derived for the high-pressure water boiling experimental data with the upper limit for the bubble departure
diameter (dBW = 1.4 mm). The correlation (10) depends solely
on characteristic liquid subcooling, which makes this correlation case dependent. In the literature, also mechanistic models
based on the microscopic phenomena may be found. The model
of (Ünal, 1976) presumes that the bubble is generated by evaporation of the superheated liquid micro-layer, while at the same
time the top of the bubble is subjected to condensation. The
bubble may slide or collapse on the heated surface during its
lifetime. After reaching the maximum size, the bubble departs
from the wall. Further investigation is necessary, since the model
of bubble departure diameter dW has a strong influence on the
share of the heat transfer components and also on the calculated
amount of steam production.
A simple estimation of the bubble departure frequency as
the terminal rise velocity over the bubble departure diameter is
adopted here (Cole, 1960):
4g(ρL − ρG )
f =
3dW ρL
The quenching heat transfer coefficient is calculated using the
analytical solution for one-dimensional transient conduction, as
E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731
suggested by Mikic and Rohsenow (1969):
2 hQ = √ f tW kL ρL CpL
where tW is waiting time between the bubble departure and the
appearance of a next bubble at the same nucleation site. A simple
assumption by Kurul and Podowski (1991) is adopted also here,
where the waiting time takes 80% of the bubble departure period:
tW =
Definition of the characteristic temperature TLW is a key point in
the CFD-implementation of the wall boiling model. Its importance arises from the spatially averaged representation of the
boiling phenomena by the mechanistic model, which must be
combined with the detailed local nature of the CFD method. The
constitutive correlations for the quenching heat flux (Eqs. (3)
and (12)) and for the bubble departure diameter (Eq. (10)) were
originally formulated for the one-dimensional thermo-hydraulic
models in terms of the mean-flow liquid temperature. Some of
the other popular correlations, not used in this work, also depend
on the core flow velocity. Formal implementation of such correlations as the CFD wall boundary conditions assumes replacing
the non-local core flow quantities by the near-wall local ones. It
can only be done for the extremely coarse grids, with the first
near-wall grid cell covering the whole boundary layer thick-
ness. With the near-wall grid, which is adequate for the CFD
turbulent wall treatment, such a straightforward approach can
significantly overestimate the vapour generation, since the CFD
solution attempts to resolve the internal superheated part of a
temperature boundary layer in the liquid phase (Menter, 2003).
In CFX, the characteristic temperature TLW is taken from the
reconstructed temperature profile using the temperature wall
function (4) and the given value of the non-dimensional distance
y+ to the wall.
In this work a constant value of 250 is used for y+ , although
more sophisticated choices are also possible, for example calculating y+ on the base of the local bubble departure diameter,
or on the base of the superheated liquid layer thickness. The
model, implemented in CFX, is found virtually invariant to the
near-wall grid cell size in a wide range covering the logarithmic
layer. Besides the improved heat and mass transfer representation, this implementation enables the use of the adequate grids
for the velocity wall function, which is a prerequisite for the
accurate calculation of the wall shear stress and the pressure
drop in internal flows.
Steam is assumed to be at saturation condition. Within the
subcooled liquid (TL < Tsat ) steam is condensing with the mass
transfer rate per unit volume:
ṁ = max
hLG (Tsat − TL )ALG
Fig. 1. Calculated distributions of void fraction (a) and temperature (b) in the channel; (c) measured vs. calculated averaged void fraction, averaged temperature,
axis temperature and wall temperature.
E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731
With superheated liquid, fluid is evaporating at the rate:
ṁ = max
hLG (TL − Tsat )ALG
2.2. Modelling of the momentum transfer
ALG is the interfacial area, and hLG is the interfacial heat transfer
coefficient, calculated according to Ranz and Marshall (1952):
hLG =
Nu =
(2 + 0.6Re1/2 Pr1/3 )
To close the phase transition model in the bulk bubbly flow with
a mean bubble diameter dB , Kurul and Podowski (1991) and also
Anglart et al. (1997) proposed to calculate the bubble diameter
dB locally as a linear function of liquid subcooling Tsub :
dB =
dB1 (Tsub − Tsub,2 ) + dB2 (Tsub,1 − Tsub )
Tsub,1 − Tsub,2
For typical nuclear energy applications these authors proposed
reference bubble diameters at the two reference subcooling
conditions: dB1 = 0.1 mm at Tsub,1 = 13.5 K and dB2 = 2 mm at
Tsub,2 = 5 K. 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.
For the liquid phase the Shear Stress Transport (SST) turbulence model by Menter (1994) was applied, which operates for
the considered here internal flows in long pipes and channels
similar to the k–ε model. The influence of the gas bubbles on
turbulence in liquid was modelled using Sato et al.’s (1981) eddy
viscosity model for bubble-induced turbulence.
Besides the drag forces, representing the flow resistance and
modelled here using the correlation by Ishii and Zuber (1979),
the non-drag forces also need to be modelled to predict the
non-homogeneous flow structure. Namely the turbulent dispersion, the lift force, and the wall lubrication force should
be considered. In the following expressions the forces acting
on gas bubbles are given. The non-drag forces were investigated analysing lots of experiments of adiabatic air/water bubbly
flow in a vertical tube for a large parameter range of superficial velocities performed in FZ-Rossendorf (see Lucas et al.,
The applied turbulent dispersion force is based on the Favre
average of the interfacial drag force (Burns et al., 2004) and is
calculated for a two phase flow according to
FDisp = −
3CD μt
grad α
(UG − UL )
4dB σt
Fig. 2. Calculated heat flux partition (a) and bubble size diameter at departure dW and in the bulk dB (b).
E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731
Fig. 3. Condensation rate (a) and cross sectional averaged condensation, bulk evaporation rate and interfacial area (b).
with CD as the bubble drag coefficient, μt the dynamic eddy
viscosity of liquid, and σ t the turbulent Schmidt number for the
volume fraction of the dispersed phase.
The lift force is defined in a standard way as:
using the gas and liquid velocities UG and UL and the liquid
density ρL . Tomiyama et al. (1995), Tomiyama (1998) has performed extensive investigations of the lift force coefficient Clift
as a function of the bubble size in the air–water two-phase flow.
He has found a sign change of Clift : when growing dB passes a
critical value, the lift coefficient becomes negative. For air water
flow at atmospheric pressure this value amounts to 5.8 mm. The
evaluation of Tomiyamas correlation for Clift (Tomiyama, 1998)
and experiments performed at TOPFLOW (Prasser et al., 2005)
show that for steam water flow at higher pressure the critical bubble size is shifted to smaller values: At 5 MPa to about 3.5 mm
and at 15 MPa to about 2 mm. Nevertheless, the bubble size in
Fig. 4. Measured (symbols) and calculated (solid lines) vapour volume fractions
dependent on the quality x.
Fig. 5. D = 12 mm, Q = 1.0 MW/m2 , G = 1000 kg/(m2 s).
FLift = −Clift ρL (UG − UL ) × rot UL
E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731
and another correlation by Tomiyama (1998)
FWall T
Fig. 6. D = 12 mm, Q = 0.8 MW/m2 , P = 7 MPa.
the here presented application cases, which amounts at departure
in the order of 0.5 mm and in the bulk in the order of 0.2–1.5 mm
will be lower than the critical value.
Investigations with air water flow indicate, that the system of
non-drag forces, which yields very good results with velocities
of about 1 m/s overestimates the gas fraction wall peak at high
water velocities >4 m/s (see Lucas et al., 2004).
The wall lubrication force pushes the bubbles away from the
wall. The available model correlations have been tested only for
adiabatic air water bubbly flows without bubble generation at
the wall, and therefore their use with the wall boiling model is
debatable. Known correlations are a correlation by Antal et al.
FWall A = −
ρW α
CW1 − CW2
Fig. 7. D = 12 mm, P = 7 MPa, G = 1000 kg/(m2 s).
dB α
= −CW
(D − y)2
ρL Urel
with n as the normal vector to the wall, dB as the bubble size
diameter, Urel as the velocity difference between the phases,
and α as the gas volume fraction. The force is dependent on
the distance to the wall y. Antal et al. (1991) proposed an
inverse proportionality to y, whereas Tomiyama (1998) proposed an inverse quadratic proportionality to y. The latter
model has been extensively validated for adiabatic air water
bubbly flows, whereas other investigations showed, that the
model of Antal seems to be more appropriate for wall boiling flows (see Končar et al., 2005a). A final statement on this
issue cannot yet be done and further investigations are necessary.
Although the implementation of the wall lubrication force
is necessary for the adiabatic two-phase flows, as it reproduces
the void fraction peak near the wall (Lucas et al., 2004) its use
at high-pressure wall boiling conditions may be questionable.
Pierre and Bankoff (1966) have found in their experiments no
void peak near the wall. Also, no void peak has been observed at
boiling of refrigerants, where the fluid properties are similar to
boiling of water at high pressures (Roy et al., 2002). At wall boiling the spherical-cap bubbles continuously grow and/or slide on
the heated wall. According to its formulation a zero void fraction at the wall is required. The Antal or Tomiyama formulation
of the wall force (Eqs. (20) and (21)) formally contradict with
the non-zero wall source of dispersed phase, since their integrals across the boundary layer are infinitely large. Still, the use
of wall lubrication force (Končar et al., 2004) was proved to be
advantageous for the conditions of low-pressure boiling of water,
where the void fraction peak was observed (Bartel et al., 2001).
Further research is necessary to improve the existing wall force
Fig. 8. D = 12 mm, P = 15 MPa, G = 2000 kg/(m2 s).
E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731
3. Validation of the model
An experimental arrangement suited for the comparison of
measured and calculated results is an upward water flow through
a heated tube. The boundary conditions can be clearly specified and the resulting cross sectional averaged steam volume
fraction can be measured by gamma densitometry with acceptable accuracy. Numerous experiments of this kind are reported
in the literature. Examples are the experiments published by
Bartolomej and Chanturiya (1967) and Bartolomej et al. (1982).
The boiling model, described in the previous chapter, was
applied to simulate these experiments to show the plausibility of
the calculated results. To identify the range of validity of applied
boiling model the calculated cross sectional averaged steam volume fractions were compared with the large set of measurements
at different operating conditions
3.1. Benchmark case
For the benchmark calculation the Bartolomej experiment
performed in 2 m long heated tube with the inner diameter
of 15.4 mm was selected. The heat flux was 5.7 × 105 W/m2
and the mass flow of the water at the pressure of 4.5 MPa
amounted to 900.0 kg/(s m2 ). In the calculations the inlet subcooling was set to 58.2 K so that the thermodynamic quality x = 0
was reached 1.75 m downstream of the tube inlet. For this particular experiment also the wall temperature along the tube was
The tube was modelled in 2D cylindrical geometry with 20
grid cells in radial and 150 grid cells in axial direction. The grid
refinement analysis performed in our previous work (Končar et
al., 2005a) showed that this grid is a good compromise between
the numerical accuracy and the computational effort. The following closure models were used for the reference calculation:
the SST turbulence model for the liquid phase, the Sato model
for the bubble induced turbulence, the Ishii–Zuber model of the
interfacial drag, the Tolubinsky model for the bubble departure
diameter (Eq. (10)) and the model for bulk bubble diameter size
according to Eq. (17). The Favre-averaged turbulent dispersion
force by Burns et al. (2004) was the only considered non-drag
force. Since for high-pressure boiling of water there is no available measurements of local two-phase flow parameters, only the
simplest closure models have been selected here. The decision
on the selected set of correlations is also based on good agreement with experiments over a wide range of flow conditions, as
shown in the next section.
Fig. 1a shows the calculated distribution of the steam volume
fraction in the bulk and Fig. 1b of the water temperature. The
presentation is stretched in radial direction. The tube is heated
from the right side and the central symmetry axis is on the left
Fig. 9. Influence of doubling the bubble diameter at departure on the wall heat flux partition (a) and on the condensation rate (b) (dotted lines—reference case).
E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731
Fig. 10. Influence of doubling the bubble diameter in the bulk on the condensation rate (a) and on the vapour volume fraction (b) (dotted lines—reference case).
side. Fig. 1c shows the comparison of measured (symbols) and
calculated (solid lines) cross sectional averaged values of steam
volume fraction and water temperature and also the wall temperature values along the axial tube length. Both, the point of
boiling inception and the increasing of the vapour volume fraction over the height show good agreement to the measured data.
Also the calculated temperatures correspond very well to the
To show the plausibility of the modelled results, further
parameters are extracted from the solution. The calculated wall
heat partitioning (Fig. 2a) shows, that at the beginning almost all
heat flux is transferred by the single-phase convection, whereas
Fig. 11. (a) Radial distribution considering all non-drag forces. The dotted line shows the reference case, where wall and lift bubble force were not considered. (b)
Influence of the non-drag forces on the cross sectional averaged vapour volume fraction (dotted lines—reference case).
E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731
Fig. 12. Vapour volume fraction distribution in the cross section at z = 0.475 m
without consideration of an inlet swirl.
with increasing tube length quenching and evaporation become
larger. The cross sectional averaged values of bubble departure
diameter and local bubble diameter in the bulk are compared in
Fig. 2b. According to Eq. (17) the bulk bubble diameter strongly
depends on the local water temperature in the bulk therefore it
increases significantly towards the end of the channel where
bulk evaporation occurs. Fig. 3a shows the condensation rate
calculated according to Eq. (14), which depends on the liquid
subcooling, the vapour volume fraction and the interfacial area.
Fig. 3b shows cross sectional averaged values of condensation
rate, bulk evaporation rate and interfacial area along the channel
presented model setup. The calculations were performed with
the same set of models as for the benchmark case. Fig. 4 shows
measured (symbols) and calculated (solid lines) cross sectional
averaged steam volume fractions dependent on the thermodynamic quality x for different pressures (see Bartolomej and
Chanturiya, 1967). The first case at 1.5 MPa corresponds to
the test presented in the previous chapter. The tube diameter
was 15.4 mm, the heat flux 3.8 × 105 W/m2 and the mass flow
amounted to 900.0 kg/(s m2 ). The point of boiling inception is
very well predicted and also the trends of the volume fraction
development with increased thermodynamic quality show good
agreement with experiments.
Figs. 5–8 show the results of further test series performed in
a tube with inner diameter of 12 mm (Bartolomej et al., 1982).
The pressure effect is presented in Fig. 5. A reasonable agreement with experiments was obtained at pressures between 3
and 11 MPa, whereas the averaged void fraction is somewhat
underpredicted at 15 MPa. The influence of mass flux at constant pressure and heat flux is presented in Fig. 6. A significant
under-estimation of void fraction was found at the lowest mass
flux of G = 400 kg/(m2 s), whereas a good agreement with experiments may be observed at higher mass fluxes. The influence of
heat flux at constant values pressure p = 7 MPa and mass flux
G = 1000 kg/(m2 s) is presented in Fig. 7. A very good prediction
may be observed at moderate heat fluxes, whereas some discrepancy occurs at higher heat flux values of 1.7 and 2 MW/m2 .
3.2. Comparison with experiments—range of validity
To investigate the range of validity of the subcooled boiling
model at different pressure, mass flow and heat flux conditions,
numerous experiments performed in different circular tubes and
published by Bartolomej (1967, 1982) were compared to the
Fig. 13. Velocity field at the inlet: axial velocity 5 m/s, velocity at the circumference of the swirl 3 m/s (case d).
Fig. 14. Axial profiles of the swirl (see Eq. (21)) for the investigated cases
(vz = 5m/s).
E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731
Fig. 15. Stream line presentation of the simulated heat channel.
The heat flux effect at higher pressure p = 15 MPa and mass
flux G = 2000 kg/(m2 s) may be observed in Fig. 8. Here the
void fraction is over-predicted at most of the heat flux values,
whereas the discrepancy is the smallest at the highest heat flux
(2.2 MW/m2 ). Discrepancies at high pressures (15 MPa) may be
caused due the use of simple model for bubble departure diameter (Tolubinsky and Kostanchuk, 1970), which is not pressure
Nevertheless, in all presented cases, the qualitative trends are
reasonably predicted. Summarizing this comparison, the model
presented in the previous chapter is best suited for pressures
from 3 to 11 MPa, for heat fluxes up to 1.2 MW/m2 and for mass
flow rates at about 1000 kg/(m2 s).
the cross sectional averaged void fraction however is partially
compensated by a higher condensation rate (see Fig. 9b).
The bubble size in the bulk has via the interfacial area an
influence on the condensation rate (see Fig. 10a). The effect on
the vapour volume fraction (see Fig. 10b) is partially compensated by the radial temperature profile (see Fig. 11). Most of
the vapour is found in areas of temperatures near the saturation
temperature and small condensation rates.
Although the influence of the bubble sizes on the vapour
volume fraction is limited, the calculations show the influence
3.3. Identification of most sensitive model parameters
In the actual calculations a quite simple correlation for the
bubble size diameter at departure dW fitted to experiments was
used (see Eq. (10)). In the literature correlations based on mechanistic models of the microscopic phenomena are known (e.g.
Ünal, 1976) and their application may extend the parameter
range of validity. To demonstrate the sensitivity of the calculated
vapour volume fraction on the modeling of dW , the bubble bubble departure diameter was doubled. The bubbles grow longer
at the wall before they leave into the bulk. The calculation with
the doubled bubble size in the Figs. 9 and 10 is represented
by the solid lines, whereas the dotted lines show the reference
case. In the case of doubled diameter a higher fraction of the
heat flux contributes to evaporation (see Fig. 9a). The effect on
Fig. 16. Radial profiles of the non-drag forces for the rod bundle geometry at
VR = 3m/s.
E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731
on the relation of the microscopic mechanisms. The bubble size
at departure has a large influence on the wall heat partition. The
simulated bubble size in the bulk influences the condensation
rate mainly via the interfacial area density. Unfortunately the
measurement of the bubble sizes by actual techniques is very
difficult. Difficulties arise due to very narrow channels of about
10 mm, where the size of the vapour bubbles is expected in the
range between 0.5 and 2 mm. Furthermore, the pressure and
temperature of the working fluid typical for nuclear reactor is
very high. The design of experiments using model fluids enables
the enlargement of the geometry and easier determination of
cross sectional vapour distribution. Recent attempts applying
tomographic methods are very promising.
3.4. Influence of non-drag forces
To show the influence of the lift and the wall force, corresponding simulations were performed. First test calculations
revealed, that the assumed lift coefficient Clift = 0.288 like in air
water adiabatic flow is too high. Caused by the shifted water
velocity profile by the generated steam all generated vapour
remained at the rod surface.
Therefore, for the calculations presented in this study, a
smaller lift force coefficient of Clift = 0.06 was assumed. A wall
force according to Eq. (20) (Antal et al., 1991) with the coefficients CW1 = −0.025 and CW2 = 0.075 was considered.
The non-drag forces have influence on the radial vapour volume fraction profiles. In Fig. 11a radial profiles of the liquid
velocity (blue line, left axis), gas velocity (green line, left axis),
gas volume fraction (red line, right axis) and liquid temperature
(black line, left temperature axis) are presented. At z = 0.5 m the
flow is still single phase and a typical profile of a liquid flow in a
tube can be found. At 1.75 m steam is present. The vapour velocity flow profile is shown by the green line, the liquid profile is
seen accelerated by the vapour. The consideration of the bubble
forces in the calculation caused that the maximum value of the
vapour volume fraction decreases and is shifted somewhat away
from the wall. The maximum liquid temperature (black solid
line) is calculated just at the wall. The influence on the cross
sectional averaged vapour volume fraction is almost negligible
(see Fig. 11b).
Fig. 17. Vapour volume fraction distribution in the cross section at z = 0.475 m
considering a swirl (case d).
of 12.6 mm. The thermal hydraulic and transport water properties were set for a pressure of 15.7 MPa, typical for PWR
conditions. The heat flux at the rod surface was assumed to
be 1.0 × 106 W/m2 and the subcooling at the inlet was set to
12 K expecting the generation of vapour in the simulated section. The axial water velocity was set to Uz = 5 m/s. The faces at
the low and high x respective at low and high y were simulated
as periodic boundary conditions.
4. Simulation of a hot channel of a fuel assembly
In this chapter CFX calculations of a hot channel between
four rods are presented. Although the CFD modelling of critical
heat flux is not yet possible, the simulations shall demonstrate
the capability of CFD supporting the fuel assembly design. In
the calculations only subcooled boiling is simulated, which is
here considered as a preliminary phenomenon towards departure
of nucleate boiling (DNB). DNB might occur at the thermal
hydraulic conditions of a PWR. A situation was investigated,
when at full power and full pressure the inlet temperature rise
caused undesired boiling in the channel.
A section between two spacer grids having a length of
z = 0.5 m was simulated. The grid represents a subchannel
between 4 rods having a diameter of 9 mm and a rod distance
Fig. 18. Influence of the swirl on the cross sectional averaged vapour volume
fraction and water temperature for case d).
E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731
The same model setup for mass, momentum and heat transfer
as in the benchmark case of Section 3.1 was applied. The result
is a vapour volume fraction distribution in the channel cross
section at z = 0.475 m (see Fig. 12). The figure shows maximum
vapour generation at the smallest distance between the rods.
To investigate the effect of a swirl caused by mixing vanes,
additional horizontal velocity components UX and UY were
given at the inlet according to Eq. (22).
= −UR
Ux = UR cos α +
Uy = UR sin α +
= UR ,
for 1 mm < r < 3.8 mm
Besides the reference case (a) with zero transversal velocity (UR = 0 m/s), the following cases with different transversal
velocities at the circumference of the circular swirl VR were
simulated: (b) 1 m/s, (c) 2 m/s and (d) 3 m/s to simulate examples for different geometrical design of mixing vanes. Fig. 13
shows the resulting inlet velocity field for the case (d) with
UR = 3 m/s.
Fig. 14 presents the relation of the horizontal velocity components to the overall value of the velocity corresponding to Eq.
(23) along the vertical axis z.
Ux2 + Uy2
swirl = Ux2 + Uy2 + Uz2
Due to the centrifugal force, the heavier fluid component – the
water – is pushed outwards, whereas a large amount of the lighter
component – the vapour – is accumulated in the centre of the
channel. The streamline presentations of Fig. 15a for water and
Fig. 15b for vapour clearly show this phenomenon.
The centrifugal force Fcent attracting at the gaseous phase can
be calculated according to
F cent = −ρL U
=∇ ×U
and the vorticity ω.
with the velocity vector U
Fig. 16 shows the radial profiles
bubble forces
projected to a normal vector (1/ 2) (1; 1; 0) pointing from the
centre of the channel to the centre of a rod. FDisp is calculated
according to Eq. (18) and Flift according to Eq. (19). The figure
shows that in this geometry and under these flow conditions
the flow is mainly determined by the centrifugal force and the
turbulent dispersion force. The latter force acts mainly in the
Fig. 19. Cross sectional averaged heat flux components (a) and superheating temperature at the rod surface (b) for the case d) (solid lines) and the reference case
without swirl (dotted lines).
E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731
centre of the channel. Here a large vapour volume fraction and a
large gradient of the vapour volume fraction can be found. The
lift force plays only a limited role.
Fig. 17 shows the calculated vapour volume fraction distribution at z = 0.475 m considering the counter clockwise swirl.
In the shadow of the rods bubble accumulations are calculated.
In the calculations like in the basic case only a monodispersed
bubble size distribution in the bulk dependent on the liquid temperature according Eq. (17) was modelled. Further work has to
consider however for this high vapour volume fraction in the
centre bubble coalescence and the occurrence of larger bubbles.
Hence in the actual calculation the interfacial area and consequently the condensation in the bulk are overestimated. In the
reality, higher vapour volume fractions might be expected.
In Fig. 18 the influence of the swirl on the axial profile of
the cross sectional averaged vapour volume fraction for the case
(d) is presented. Indeed the cross sectional averaged vapour volume fraction is reduced by the swirl as desired—in this example
Fig. 20. Distribution of the wall superheating temperature at the surface of rod 1 without swirl (a) and considering a swirl (cases b–d) (for the definition of the angle
Φ, here in radians, see Fig. 11).
E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731
from about 8% to about 6.7%, whereas the averaged water temperature in the channel is almost not influenced. The tangential
velocity near the wall is increased by the swirl. As a consequence, higher portions of the wall heat flux are transformed
by single-phase convection whereas the share of evaporation is
decreased (see Fig. 19a). Furthermore the cross sectional averaged wall superheating temperature is reduced (see Fig. 19b).
In the presented example the reduction at z = 0.1 m amounts to
about 1.5 K.
Very informative is the distribution of the superheating Tsup
at the rod surface of rod 1 shown in Fig. 20 for the four cases
(a–d). Without swirl, a symmetric distribution is found showing
maximum values at locations with the smallest distance of the
adjacent rods.
In contrast, the assumed swirl generates a non-equally distributed wall superheating. Hot spots (e.g. at z = 0.15, Φ = 0.3)
might be the first locations, where later on the critical heat flux
could occur.
5. Summary and conclusions
Capabilities of the up-to-date two-fluid CFD method of simulating steam–water bubbly flows with phase transition are
demonstrated. To simulate boiling at a heated wall, the phenomena on the micro-scale have been modelled by appropriate
closure relations. The given heat is removed from the wall by
different mechanisms. The size of bubbles, emerging from the
wall, has been proven to have a large influence on the heat distribution between the different heat flux components. The size of
bubbles in the bulk is correlated to the local liquid subcooling,
since it determines the interfacial area and the condensation rate
of the vapour bubbles.
Further improvements of the microscopic models are necessary. The model for the bubble size at departure has to be replaced
by a mechanistic model. The simulation of the bubble size in the
bulk has to consider bubble coalescence and fragmentation.
The experimental design needs to be optimised to deliver
the relevant information for the CFD model improvement.
The determination of the spatial vapour volume distribution is
necessary to support further model developments, especially
modelling of the momentum exchange between gas bubbles
and the liquid. The determination of the bubble size would
be desirable for correct simulation of the bubble behaviour at
the heated surface and in the bulk. Nevertheless, despite its
current developmental state, the used CFD model yields to correct quantitative results and predicts right trends in parametric
Although only a preliminary state towards the critical heat
flux is simulated, essential parameters like the swirl, the cross
flow between adjacent channels and concentration regions of
bubbles can be determined. By predicting the temperature at the
rod surface, critical regions can be identified, which might later
on lead to departure from nucleate boiling and possible damage
of the fuel pin. The avoiding of uneven temperature distributions
at the rod surface can be considered as a design criterion for the
spacer grids. The application of actual CFD with a subcooled
boiling model for the simulation of a hot channel enables at
least the comparison and the evaluation of different geometrical
designs of the spacer grids.
Anglart, H., Nylund, O., Kurul, N., Podowski, M.Z., 1997. CFD prediction of
flow and phase distribution in fuel assemblies with spacers. NURETH-7,
Saratoga Springs, New York, 1995. Nucl. Eng. Des. 177, 215–228.
Antal, S.P., Lahey, R.T., Flaherty, J.E., 1991. Analysis of phase distribution in
fully developed laminar bubbly two-phase flow. Int. J. Multiphase Flow 7,
Bartel, M.D., Ishii, M., Masukawa, T., Mi, Y., Situ, R., 2001. Interfacial area
measurements in subccoled flow boiling. Nucl. Eng. Des. 210, 135–155.
Bartolomej, G.G., Chanturiya, V.M., 1967. Experimental study of true void fraction when boiling subcooled water in vertical tubes. Thermal Engineering,
vol. 14, pp. 123–128 (translated from Teploenergetika, no. 2, vol. 14, 1967,
pp. 80–83).
Bartolomej, G.G., et al., 1982. An experimental investigation of true volumetric
vapour content with subcooled boiling in tubes. Thermal Engineering, vol.
29, pp. 132–135 (translated from Teploenergetika, no. 3, vol. 29, 1982, pp.
Burns, A.D., Frank, T., Hamill, I., Shi, J.-M., 2004. The favre averaged drag
model for turbulence dispersion in Eulerian multi-phase flows. In: 5th Int.
Conf. on Multiphase Flow, ICMF’2004. Yokohama, Japan.
Cole, R., 1960. A photographic study of pool boiling in the region of the critical
heat flux. AIChE J. 6, 533–542.
Egorov, Y., Menter, F., 2004. Experimental Implementation of the RPI Wall
Boiling Model in CFX-5.6. Technical Report ANSYS/TR-04-10. ANSYS
Germany GmbH.
Ishii, M., Zuber, N., 1979. Drag coefficient and relative velocity in bubbly,
droplet or particulate flows. AIChE J. 25, 843–855.
Kader, B.A., 1981. Temperature and concentration profiles in fully turbulent
boundary layers. Int. J. Heat Mass Transfer 24, 1541–1544.
Končar, B., Kljenak, I., Mavko, B., 2004. Modelling of local two-phase flow
parameters in upward subcooled flow boiling at low pressure. Int. J. Heat
Mass Transfer 47, 1499–1513.
Končar, B., Krepper, E., Egorov, Y., 2005a. CFD Modelling of Subcooled
Flow Boiling for Nuclear Engineering Applications International Conference Nuclear Energy for New Europe, Bled, Slovenia, September 5–8, 2005
(Paper 140).
Končar, B., Mavko, B., Hassan, Y.A., 2005b. Two-phase wall function for
modeling of turbulent boundary layer in subcooled boiling flow. In: 11th
International Topical Meeting on Nuclear Reactor Thermal-Hydraulics
(NURETH-11), Avignon, France, October 2005 (Paper 443).
Kurul, N., Podowski, M.Z., 1991. On the modeling of multidimensional effects
in boiling channels. In: ANS Proceedings of 27th National Heat Transfer
Conference, Minneapolis, MN.
Lemmert, M., Chawla, J.M., 1977. Influence of flow velocity on surface boiling
heat transfer coefficient. In: Hahne, E., Grigull, U. (Eds.), Heat Transfer
in Boiling. Academic Press and Hemisphere, ISBN 0-12-314450-7, pp.
Lucas, D., Shi, J.-M., Krepper, E., Prasser, H.-M., 2004. Models for the forces
acting on bubbles in comparison with experimental data for vertical pipe
flow. In: 3rd International Symposium on Two-Phase Flow Modelling and
Experimentation, Pisa, Italy (Paper ha04).
Lucas, D., Krepper, E., Prasser, H.-M., 2005. Development of co-current airwater flow in a vertical pipe. Int. J. Multiphase Flow 31, 1304–1328.
Menter, F.R., 2003. Personal communication.
Menter, F.R., 1994. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA J. 32, 8.
Mikic, B.B., Rohsenow, W.M., 1969. A new correlation of pool-boiling data
including the fact of heating surface characteristics. ASME J. Heat Trans.
91, 245–250.
Pierre, C.C., Bankoff, S.G., 1966. Vapor volume profiles in developing twophase flow. Int. J. Heat Mass Transfer 10, 237–249.
Prasser, H.-M., Beyer, M., Carl, H., Gregor, S., Lucas, D., Pietruske, H., Schütz,
P., Weiss, F.-P., (Paper 399) 2005. Evolution of the structure of a gas-
E. Krepper et al. / Nuclear Engineering and Design 237 (2007) 716–731
liquid two-phase flow in a large vertical pipe. In: NURETH-11, Avignon,
Ranz, W.E., Marshall, W.R., 1952. Evaporation from drops. Chem. Eng. Progr.
48 (3), 141–146.
Roy, R.P., Kang, S., Zarate, J.A., Laporta, A., 2002. Turbulent subcooled boiling
flow—experiments and simulations. J. Heat Trans. 124, 73–93.
Sato, Y., Sadatomi, M., Sekoguchi, K., 1981. Momentum and heat transfer in
two-phase bubble flow. Int. J. Multiphase Flow 7, 167–177.
Tolubinsky, V.I., Kostanchuk, D.M., 1970. Vapour bubbles groth rate and heat
transfer intensity at subcooled water boiling; Heat Transfer 1970, Preprints
of papers presented at the 4th International Heat Transfer Conference, vol.
5, Paris (Paper No. B-2.8).
Tomiyama, A., et al., 1995. Effects of Eötvös number and dimensionless liquid
volumetric flux on lateral motion of a bubble in a laminar duct flow. Adv.
Multiphase Flow, 3–15.
Tomiyama, A., 1998. Struggle with computational bubble dynamics. In: Third
International Conference on Multiphase Flow, ICMF’98, Lyon, France, June
8–12, 1998.
Ünal, H.C., 1976. Maximum bubble diameter, maximum bubble growth time
and bubble growth rate. Int. J. Heat Mass Trans. 19, 643–649.
Wintterle, T. 2004. Development of a numerical boundary condition for the
simulation of nucleate boiling at heated walls. Diploma Thesis University.
Stuttgart (IKE-8-D-014).
Fly UP