# Mechanistic Multidimensional Modeling of Forced Convection Boiling Heat Transfer

Authors: Michael Z. Podowski; Raf M. Podowski
Journal: Science and Technology of Nuclear Installations
 (2009)

---

## Abstract

Due to the importance of boiling heat transfer in general, and boiling crisis in particular, for the analysis of operation and safety of both nuclear reactors and conventional thermal power systems, extensive efforts have been made in the past to develop a variety of methods and tools to evaluate the boiling heat transfer coefficient and to assess the onset of temperature excursion and critical heat flux (CHF) at various operating conditions of boiling channels. The objective of this paper is to present mathematical modeling concepts behind the development of mechanistic multidimensional models of low-quality forced convection boiling, including the mechanisms leading to temperature excursion and the onset of CHF.


---

## Body

## 1. Introduction

Because
of the complexity of phenomena governing boiling heat transfer in general, and
subcooled boiling in particular, the predictions of CHF have traditionally been
based on correlating data obtained from numerous experimental measurements. At
the same time, our understanding of the local physical
mechanisms governing the near-wall phase change and transport has been steadily
improving. Given the progress made in the computational fluid dynamics methods,
the possibility of using complete multidimensional models to predict boiling
heat transfer before, and up to, the onset of CHF becomes an attractive option
complementing the traditional phenomenological approach used in the past.The
objective of this paper is to discuss various physical and mathematical
modeling concepts for local heat transfer phenomena in boiling systems, and to
show that the proposed approach can be combined with mechanistic
multidimensional models of two-phase flow and used to predict various
parameters characterizing low-quality forced-convection boiling, including the
mechanisms leading to temperature excursion and the onset of CHF.


## 2. Multidimensional Multifield Model of Forced-Convection Boiling

The multifield modeling concept is based on coupling a complete mechanistic
multidimensional model of two-phase flow with the models governing local heat
transfer and phase change phenomena in heated channels. The major components of
the overall model are the following:(a)
conservation equations for the individual fields,(b)
the model of turbulence,(c)
closure laws for the interfacial mass transfer,(d)
closure laws for the interfacial momentum transfer,(e)
closure laws for the interfacial energy transfer,(f)
kinematic boundary conditions for phasic velocities and turbulence,(g)
thermal boundary conditions for the near-wall heat transfer.


### 2.1. Multifield Conservation Equations

The multifield model of two- and multiphase flows assumes that each phase
can be represented by either one or several fields. In the latter case, the
individual fields represent topologically different flow structures within a
given phase (such as liquid film and liquid droplets or groups of gas bubbles
of different sizes). Typical ensemble-averaged conservation equations for a
multifield model can be written as follows:Mass(1)∂(αkρk)∂t+∇⋅(αkρkuk)=Γk,Momentum(2)∂(αkρkuk)∂t+∇⋅(αkρkukuk)=−∇(αkpk)+∇⋅(αkτktot)+αkρkg+M^ki,Energy(3)∂(αkρkek)∂t+∇⋅(αkρkukek)=−∇⋅(αkqktot)+∇⋅[αk(−pkI+τktot)⋅uk]+αkρkg⋅uk+E^ki,
where αk, ρk, uk,
and ek are the volume fraction, density, velocity,
and specific internal energy, respectively, of field-k,  Γk is the net mass transfer rate (this and all
other terms are per unit volume) of field-k, τktot is the total shear stress, M^ki is the total interfacial force, qktot is the local heat flux, and E^ki is the interfacial heat transfer rate, all for
field-k.Details concerning the
various terms in (1)–(3) can be found
in [1].In the
modeling of low-quality nucleate boiling flows, a two-field model can normally
be used, which accounts for the continuous liquid (l) and dispersed
vapor (v) fields. This is shown in Figure 1(a), where upon departing
from the heated wall, bubbles move into the subcooled liquid region and
condense.An illustration of the multifield model of two-phase flow in a boiling channel.
(a)
Two-field model(b)
Three-field modelFor high
wall heat fluxes, vapor concentration near the heated wall increases, so that
the individual departing bubbles may coalesce and form elongated bubbles, as
shown in Figure1(b). Since the
mechanisms governing the motion and heat transfer of elongated bubbles may be
quite different from those for small spherical bubbles (field-v1),
it is appropriate to use a separate field for the former, denoted here as
field-v2.The
multidimensional multifield conservation equations, (1)–(3), must be
complemented by appropriate boundary conditions and closure laws. The boundary conditions of particular
interest to this work are associated with the effect of heated channel wall on
local flow and phase change in the near-wall region. The necessary closure laws
and models are those for turbulence, interfacial forces, interfacial mass
transfer, and heat transfer. A brief description of the major closure laws and
boundary conditions pertinent to low-quality two-phase flows in boiling
channels is given below.


### 2.2. Multifield Model of Turbulence Model

Two-phase flow turbulence is normally modeled using theκ-ε model [2] for
the continuous liquid field (k=l), modified to include the effect of
bubble-induced turbulence. The turbulent kinetic energy, κl, and the energy
dissipation, εl, of the liquid field are, respectively, given
by
(4)∂(αlκl)∂t+∇⋅(αlulκl)=∇⋅(αlμltρlσκ∇κl)+Slκ,∂(αlεl)∂t+∇⋅(αlulεl)=∇⋅(αlμltρlσε∇εl)+Slε.For givenκlεl, the turbulent
viscosity, μlt,
can be expressed as
(5)μlt=Cμρlκleεl+μBt,
where the
bubble-induced turbulence can be written as [3]
(6)μBt=CμBαlρl∑k=12(uvk−ul).


### 2.3. Interfacial Mass Transfer

Since in nucleate boiling, evaporation occurs only at (or
near) the heated wall; the interfacial mass transfer between the liquid and
vapor fields is practically only due to vapor condensation in contact with
subcooled liquid. In addition, after
elongated bubbles start forming, another possible interfield mass exchange
mechanism is due to small bubble coalescence.
Hence, the volumetric mass transfer terms can be expressed as(7)Γl=−Γevap+Γcond=−(Γv1+Γv2),(8)Γv1=Γv1evap−Γv1cond−Γv1clsc,(9)Γv2=Γv2evap−Γv2cond+Γv1clsc.In (7), Γevap=Γv1evap+Γv2evap and Γcond=Γv1cond+Γv2cond are the total evaporation and condensation
rates per unit volume, respectively, and Γv1clsc is the coalescence rate of small bubbles per
unit volume (which becomes a source term for the elongated bubble field).


### 2.4. Interfacial Momentum Transfer

In general, the total interfacial force on
phase-k can be expressed as a superposition of several component forces:
(10)M^ki=pki∇αk−τki⋅∇αk+Γkuki+∑jMk,ji.In
dispersed bubbly flows, the major interfacial forces are the following (for
details and source references, see [1]):Drag Force(11)Ml,vkD=−Mvk,lD=−18CD,kAvk′′′ρl(uvk−ul)|uvk−ul|,Virtual Mass Force(12)Ml,vkVM=−Mvk,lVM=−CVM,kαvkρl(DvkuvkDt−DlulDt),Lift Force(13)Ml,vkL=−Mvk,lL=−CL,kαvkρl(uvk−ul)×∇×ul,Turbulent Dispersion Force(14)Ml,vkTD=−Mvk,lTD=CTD,kαvkρl∇(αlκl)for k=1,2, where αl=1−(αv1+αv2).As it can be seen from (11)–(14), direct
momentum exchange occurs between the continuous liquid and each vapor field,
whereas bubble-bubble interactions are due to the coalescence of small bubbles
to form elongated bubbles (see Section 2.5 below) and they mainly involve mass
and energy transfers.


### 2.5. Interfacial Energy Transfer

Given a
very low-vapor superheat, the interfacial volumetric condensation rates can be
expressed as(15)Γkcond=qk′′ihfg(k=v1,v2),
where the interfacial energy transport from vapor field-k to the continuous liquid field is given by
(16)qk′′i=Hk−li(Tsat−Tl)Ak′′′(k=v1,v2).In (16), Ak′′′ is the interfacial area density for vapor
field-k, and Hk−li is the local interfacial heat transfer
coefficient that can be obtained from
(17)Hk−li=kldk[2+0.6(ρldkvr,kμl)0.5(cplμlkl)0.33],(k=v1,v2).Equations (16) and (17) imply
(18)ql′′i=qv1′′i+qv2′′i=[Hv1−liAv1′′′+Hv2−liAv2′′′](Tsat−Tl).The rate of bubble coalescence can be obtained from(19)Γv1clsc=bv1−v1αv12+bv1−v2αv1αv2,where bv1−v1 and bv1−v2 are the corresponding rate coefficients.


## 2.1. Multifield Conservation Equations

The multifield model of two- and multiphase flows assumes that each phase
can be represented by either one or several fields. In the latter case, the
individual fields represent topologically different flow structures within a
given phase (such as liquid film and liquid droplets or groups of gas bubbles
of different sizes). Typical ensemble-averaged conservation equations for a
multifield model can be written as follows:Mass(1)∂(αkρk)∂t+∇⋅(αkρkuk)=Γk,Momentum(2)∂(αkρkuk)∂t+∇⋅(αkρkukuk)=−∇(αkpk)+∇⋅(αkτktot)+αkρkg+M^ki,Energy(3)∂(αkρkek)∂t+∇⋅(αkρkukek)=−∇⋅(αkqktot)+∇⋅[αk(−pkI+τktot)⋅uk]+αkρkg⋅uk+E^ki,
where αk, ρk, uk,
and ek are the volume fraction, density, velocity,
and specific internal energy, respectively, of field-k,  Γk is the net mass transfer rate (this and all
other terms are per unit volume) of field-k, τktot is the total shear stress, M^ki is the total interfacial force, qktot is the local heat flux, and E^ki is the interfacial heat transfer rate, all for
field-k.Details concerning the
various terms in (1)–(3) can be found
in [1].In the
modeling of low-quality nucleate boiling flows, a two-field model can normally
be used, which accounts for the continuous liquid (l) and dispersed
vapor (v) fields. This is shown in Figure 1(a), where upon departing
from the heated wall, bubbles move into the subcooled liquid region and
condense.An illustration of the multifield model of two-phase flow in a boiling channel.
(a)
Two-field model(b)
Three-field modelFor high
wall heat fluxes, vapor concentration near the heated wall increases, so that
the individual departing bubbles may coalesce and form elongated bubbles, as
shown in Figure1(b). Since the
mechanisms governing the motion and heat transfer of elongated bubbles may be
quite different from those for small spherical bubbles (field-v1),
it is appropriate to use a separate field for the former, denoted here as
field-v2.The
multidimensional multifield conservation equations, (1)–(3), must be
complemented by appropriate boundary conditions and closure laws. The boundary conditions of particular
interest to this work are associated with the effect of heated channel wall on
local flow and phase change in the near-wall region. The necessary closure laws
and models are those for turbulence, interfacial forces, interfacial mass
transfer, and heat transfer. A brief description of the major closure laws and
boundary conditions pertinent to low-quality two-phase flows in boiling
channels is given below.


## 2.2. Multifield Model of Turbulence Model

Two-phase flow turbulence is normally modeled using theκ-ε model [2] for
the continuous liquid field (k=l), modified to include the effect of
bubble-induced turbulence. The turbulent kinetic energy, κl, and the energy
dissipation, εl, of the liquid field are, respectively, given
by
(4)∂(αlκl)∂t+∇⋅(αlulκl)=∇⋅(αlμltρlσκ∇κl)+Slκ,∂(αlεl)∂t+∇⋅(αlulεl)=∇⋅(αlμltρlσε∇εl)+Slε.For givenκlεl, the turbulent
viscosity, μlt,
can be expressed as
(5)μlt=Cμρlκleεl+μBt,
where the
bubble-induced turbulence can be written as [3]
(6)μBt=CμBαlρl∑k=12(uvk−ul).


## 2.3. Interfacial Mass Transfer

Since in nucleate boiling, evaporation occurs only at (or
near) the heated wall; the interfacial mass transfer between the liquid and
vapor fields is practically only due to vapor condensation in contact with
subcooled liquid. In addition, after
elongated bubbles start forming, another possible interfield mass exchange
mechanism is due to small bubble coalescence.
Hence, the volumetric mass transfer terms can be expressed as(7)Γl=−Γevap+Γcond=−(Γv1+Γv2),(8)Γv1=Γv1evap−Γv1cond−Γv1clsc,(9)Γv2=Γv2evap−Γv2cond+Γv1clsc.In (7), Γevap=Γv1evap+Γv2evap and Γcond=Γv1cond+Γv2cond are the total evaporation and condensation
rates per unit volume, respectively, and Γv1clsc is the coalescence rate of small bubbles per
unit volume (which becomes a source term for the elongated bubble field).


## 2.4. Interfacial Momentum Transfer

In general, the total interfacial force on
phase-k can be expressed as a superposition of several component forces:
(10)M^ki=pki∇αk−τki⋅∇αk+Γkuki+∑jMk,ji.In
dispersed bubbly flows, the major interfacial forces are the following (for
details and source references, see [1]):Drag Force(11)Ml,vkD=−Mvk,lD=−18CD,kAvk′′′ρl(uvk−ul)|uvk−ul|,Virtual Mass Force(12)Ml,vkVM=−Mvk,lVM=−CVM,kαvkρl(DvkuvkDt−DlulDt),Lift Force(13)Ml,vkL=−Mvk,lL=−CL,kαvkρl(uvk−ul)×∇×ul,Turbulent Dispersion Force(14)Ml,vkTD=−Mvk,lTD=CTD,kαvkρl∇(αlκl)for k=1,2, where αl=1−(αv1+αv2).As it can be seen from (11)–(14), direct
momentum exchange occurs between the continuous liquid and each vapor field,
whereas bubble-bubble interactions are due to the coalescence of small bubbles
to form elongated bubbles (see Section 2.5 below) and they mainly involve mass
and energy transfers.


## 2.5. Interfacial Energy Transfer

Given a
very low-vapor superheat, the interfacial volumetric condensation rates can be
expressed as(15)Γkcond=qk′′ihfg(k=v1,v2),
where the interfacial energy transport from vapor field-k to the continuous liquid field is given by
(16)qk′′i=Hk−li(Tsat−Tl)Ak′′′(k=v1,v2).In (16), Ak′′′ is the interfacial area density for vapor
field-k, and Hk−li is the local interfacial heat transfer
coefficient that can be obtained from
(17)Hk−li=kldk[2+0.6(ρldkvr,kμl)0.5(cplμlkl)0.33],(k=v1,v2).Equations (16) and (17) imply
(18)ql′′i=qv1′′i+qv2′′i=[Hv1−liAv1′′′+Hv2−liAv2′′′](Tsat−Tl).The rate of bubble coalescence can be obtained from(19)Γv1clsc=bv1−v1αv12+bv1−v2αv1αv2,where bv1−v1 and bv1−v2 are the corresponding rate coefficients.


## 3. Near-Wall Heat Transfer in Nucleate Boiling

In the
nucleate subcooled boiling in a heated channel, the wall heat is partially used
to form bubbles and the remaining portion is transferred to the liquid. The
heat transfer from the wall in the vicinity of a nucleation site occurs during
two distinct periods: the bubble growth time and the waiting time. The total convective heat flux from the wall
is the sum of three models [4]:
(20)qNB′′=q1ϕ′′+qe′′+qQ′′,
where q1ϕ′′ is the single-phase convective heat flux, qe′′ is the heat flux associated with phase change
(evaporation), and qQ′′ is the so-called quenching heat flux, which is
transferred to the liquid phase during the waiting time.Outside
of the influence area of the bubbles, the heat transfer from wall to the liquid
can be determined from(21)q1ϕ′′=A1ϕ′′StPρlcplul,P(Tw−Tl,P),
where A1ϕ′′ is the fraction of the wall unaffected by the
nucleation sites, StP is the Stanton number calculated from a heat
transfer correlation in terms of the local liquid velocity and Prandtl number, Tw, is the wall temperature,
and Tl,P is the
local liquid temperature near the heated wall. In numerical calculations, the
local fluid properties are normally determined at the center of the near-wall
computational cell. As long as a consistent turbulence model is used (see
Section 2.2), the calculated single-phase heat flux component is independent of
the distance between point-P and the wall, or, in other words, is
grid-independent.The evaporation heat flux is given by(22)qe′′=π6ddet3ρvfdetN′′hfg,
where ddet is the critical bubble diameter at detachment, fdet is the frequency of nucleation, and N′′ is the number of nucleation sites per unit
area (nucleation site density).The
quenching heat flux has been analytically calculated by Del Valle and Kenning
[5] as
(23)qQ′′=ΘDfdetA2ϕ′′2kl(Tw−Tl)πΘDkl/(ρlcpl),
where ΘD is the waiting (or dwell) time elapsed between
the detachment of a bubble and the nucleation of a subsequent one. The term A2ϕ′′ is the fraction of the wall area participating
in the quenching heat flux.Additional relationships needed to close the model include the following [6].(i) Bubble diameter at departure [7]
(24)ddet=0.0014⋅Exp(−ΔTsub45),
where ΔTsub is in C°, and ddet is in m. It is known that the liquid velocity
has a significant effect on the bubble diameter at detachment.(ii) Nucleation site density [8]
(25)N′′=185(Tw−Tl)1.805.Two very important terms in the model described above are the bubble waiting time and the nucleation frequency. A mechanistic approach to determine both parameters is described in Section4 based on the analysis of the bubble ebullition cycle.


## 4. Bubble Ebullition Cycle

The mechanisms of nucleation at a cavity on a heated wall
and the subsequent bubble growth have been investigated very extensively
before. However, most existing results are normally given in the form of
relationships between selected parameters, rather than as complete analytical
models. It turns out that using a rigorous analytical approach to combine
transient heat transfer solutions for the heated wall and for the liquid filling
the space vacated by departing bubbles, a consistent model can be derived for
bubble ebullition in forced convection-subcooled boiling. A schematic of the
ebullition cycle is shown in Figure2. The cycle consists of two major stages:
the quenching period and the bubble growth period.Figure 2
Surface temperature oscillations in the wall during the ebullition cycle.The
quenching (or dwell) period is initiated when the subcooled liquid fills the
space near the heated wall vacated by the departing bubble. During that time, the space inside a given
cavity is gradually filled with vapor which eventually forms a hemispherical
bubble on the top of the cavity. The quenching period is followed by a bubble
growth period during which the bubble quickly expands forming a thin liquid
sublayer separating the bubble from the heated surface. Since the liquid
sublayer is very thin, the temperature drop across this layer is very small, so
that the local wall temperature quickly drops to the saturated vapor temperature
which itself only exceeds slightly the saturation temperature corresponding the
system pressure.Accounting
to the periodic nature of the process, a closed-form solution can be obtained
for the parameters characterizing the timing of bubble ebullition. The model is
based on using coupled solutions to the transient heat conduction equations for
the heated wall and the fluid laminar sublayer near the wall:(26)∂T∂t=a∇2T,
where a=k/ρcp is the thermal diffusivity of the
respective material.Since
the characteristic time of surface temperature fluctuations during nucleation
is very short and, thus, the distance across the wall affected by a change in
the surface temperature is small compared to the size of the surface area
exposed to quenching by cold water, the wall heat conduction can be
approximated by a one-dimensional model.
Using a steady-state as a reference, the time-dependent temperature
distribution across the heated wall during the quenching period (i.e., between
0 andΘD in Figure 2) becomes
[9]
(27)Tw(y,t)=To+qH′′kwy+aw2π∫0ty[aw(t−t′)]1.5Exp[−y24aw(t−t′)]×[Ti(t′)−To]dt′,
where qH′′>0 is the constant heating rate per unit surface
area of the heated wall, To is the average (constant) temperature of the
heated wall surface in contact with the fluid, and y is the distance
across the heater away from the wetted surface.Using (27),
the following expression can be derived for the time-dependent surface heat
flux during the quenching (dwell) period
(28)qw′′(t)=−qH′′+kwπaw{∫0tdTidt′dt′t−t′+[Ti(0+)−Tot]},
where Ti(0+) is the wall temperature at the beginning of
the dwell period (the “+” superscript is used to indicate that a temperature
discontinuity occurs at t=0, as shown in Figure 2).Since
the near-wall space vacated by the departing bubble is immediately filled by
subcooled liquid, the liquid temperature starts quickly increasing in contact
with the heated wall. The time- and position-dependent liquid temperatures
across the laminar boundary layer (where the effect of flow-driven convection
is negligible) can be determined by solving the transient heat conduction
equation, (26). Assuming that the liquid temperature outside the boundary layer
is constant, Tbl=constant,
whereas the liquid temperature in contact with the heated wall is Ti(t),
yields the following expression:
(29)tl(y,t)=tbl+al2π∫0ty[al(t−t′)]1.5exp[−y24al(t−t′)][ti(t′)−tbl]dt′In a
manner similar to that used for the wall, the following expression is obtained
from (29) for the instantaneous time-dependent surface heat flux into the fluid:
(30)qw′′(t)=−kw∂Tw∂y|y=0=−klπal{∫0tdTidt′dt′t−t′+Ti(0+)−Tblt}.Combining
(28) and (30) and taking the limit as t→0+ yields the surface temperature at the
beginning of the quenching (or dwell) period:
(31)Ti(0+)=(kwToaw+klTblal)(kwaw+klal)−1.Substituting(31) back into the combined (28) and (30), one obtains
(32)∫0tdTi(t′)dt′dt′t−t′=qH′′kw/πaw+kl/πal.Solving (32) for Ti(t) with the initial condition Ti(t)=Ti(0+) yields
(33)Ti(t)=Ti(0+)+2qH′′π(kw/aw+kl/al)t.Substituting
(33) into (29) and rearranging yield the following expression for the
instantaneous surface heat flux during the dwell period:
(34)qw′′(t)=−kl∂Tl∂y|y=0=−klTi(0+)−Tblπalt−qH′′1+(kw/kl)al/aw.


## 4.1. Dwell TimeΘD

The time
of transition from the dwell period to the growth period is reached when the bubble
radius becomes equal to the wall cavity radius. At this time instant, the steam
temperature inside the hemispherical bubble outside the cavity must be equal to
the temperature of the surrounding liquid at a distance from the wall equal to
the cavity radius. Using the Clausius-Clapeyron equation to determine the steam
superheat, and combining the resultant expression with (34), one arrives at the
following relationship:
(35)2qH′′/πkw/aw+kl/alΘD−[Ti(0+)−Tbl]rcπalΘD+Ti(0+)−Tsat−qH′′rc/kl1+(kw/kl)al/aw−2σrcTsatυfghfg=0,
where rc is the cavity radius and σis the surface tension.As can
be seen, for given fluid and wall properties and thermal conditions, (35)
becomes an algebraic (quadratic) equation for the bubble dwell time, ΘD.Interestingly,
if the dwell time is known, (35) can be used to determine the average surface
heat flux during the quenching period. Specifically, integrating (34) from 0 to ΘD yields
(36)|q¯D′′|=2kl[Ti(0+)−Tbl]πalΘD+qH′′1+(kw/kl)al/aw.


## 4.2. Growth TimeΘG

As soon
as the bubble starts growing outside the cavity, the vapor is produced mainly
from the evaporating liquid sublayer between the wall and the bubble. Since the
temperature drop across the thin sublayer is very small, the local surface
temperature during this period remains close to the saturation temperature. The
energy balance for the growing bubble can be written as(37)ρghfg2d(dB)dt=|qw′′(t)|.The
time-dependent wall heat flux during the bubble growth time can be obtained
from (28) by noticing that during this period, that is, for ΘD<t<ΘD+ΘG,
the wall temperature remains approximately constant and equal to the saturation
temperature:
(38)qw′′(t)=−qH′′−kwTi(ΘD)−Tsatπaw(t−ΘD).Substituting
(38) into (37) and integrating between ΘD and ΘD+ΘG yield the following expression for the maximum
bubble diameter at detachment:
(39)ddet=2ρvhfg∫ΘDΘD+ΘG|qw′′(t)|dt=2ρvhfg[qH′′ΘG+2kwTi(ΘD)−TsatπawΘG],
where (see (33))
(40)Ti(ΘD)=Ti(0+)+2qH′′/πkw/aw+kl/alΘD.The
bubble diameter at detachment can be evaluated using one of several different
models that have been developed to date. In particular, using the force balance
for a single bubble, Staub [10] developed the following expression:
(41)ddet=0.5{[3τw2(ρl−ρv)g]2+3σF(βc)(ρl−ρv)g−3τw2(ρl−ρv)g},
where τw is the wall shear stress and F(βc) is an experimentally determined function of
the contact angle.For
given bubble diameter at detachment and bubble dwell time, (39) can be readily
solved for the bubble growth time ΘG.Interestingly,
if the bubble growth time is known, (39) can be used again, this time to
determine the average wall heat flux during the bubble growth period:
(42)|q¯G′′|=1ΘG∫ΘDΘD+ΘG|qw′′(t)|dt=qH′′+2kwTi(ΘD)−TsatπawΘG.Finally,
the bubble frequency of detachment and the total detachment time can be
obtained from(43)fdet=1Θdet=1ΘD+ΘG.Naturally,
summing up (34) and (42) yields
(44)|q¯D′′|+|q¯G′′|=qH′′,
where qH′′ is the surface heat flux of the heater (see
(27)).


## 5. Wall Temperature Excursion (CHF)

### 5.1. Physical Concept

The wall
temperature excursion (or CHF) in low-quality boiling is associated with the
ability of the bubbles formed at the nucleation sites to depart from the wall,
so that the vacated space can be filled with fresh liquid, so that the
quenching and bubble growth processes can continue [9, 11]. Thus, two conditions must
be satisfied simultaneously to avoid wall temperature excursion, one concerned
with dispersed bubble concentration in the near-wall region, the other with the
velocity at which the bubble departs from the wall to make room for the next
generation bubble formed at the same nucleation site. Since the size of bubbles
departing from the nucleation sites is normally very small (of the order of 1 mm or less), their shape is nearly spherical. Thus, the two conditions
mentioned above can be formulated as follows.(1) The bubble maximum void fraction of dispersed bubbles in the near-wall
region cannot exceed the maximum packing factor for spherical particles, the
theoretical value of which (obtained from geometrical considerations) is(45)α≤αmax=π32≈0.74.(2) Taking into account that to form a new undeformed bubble at a given
cavity at the time of detachment, the distance between the previously formed
bubble at the same cavity and the heated wall at the same time instant must be
at least equal to the bubble diameter, the maximum evaporation heat flux to
avoid temperature excursion must satisfy the following condition:(46)qe′′hfg=π6ddet3ρvfdetN′′≤ρvudet,minα,
where
(47)udet,min=mddetfdet.The
coefficient,m≥1,
reflects the fact that in reality the axial motion of bubbles may require a
distance from the wall larger than one bubble diameter to avoid the coalescence with other
bubbles formed in the adjacent upstream cavities.It is interesting to notice that (46) can be rewritten as
(48)π6ddet2N′′≤mα.


### 5.2. Near-Wall Heat Transfer in the Presence of Elongated Bubbles

A schematic illustrating the near-wall conditions in the presence of elongated
bubbles is shown in Figure3. It follows from the previous discussion that at
any axial location along the channel the time-dependent wall temperature
fluctuates periodically, increasing in the poor heat transfer (elongated bubble)
region and decreasing, where heat transfer is more efficient (in the nucleate
boiling region).Figure 3
Phase distribution in the near-wall region.According
to Figure3, averaging the local wall heat flux over both regions always yields
the given constant heating rate of the heater per unit surface area. In
particular, the overall average heat flux can be partitioned into the following
terms:
(49)qH′′=qv1′′+qv2′′.In (49), qv1′′ is the nucleate boiling component
corresponding to small dispersed bubbles:
(50)qv1′′=qNB′′Av1′′,
where qNB′′ is given by (20), and Av1′′ is the area density (wall area fraction) for
the small bubble region.The heat
flux component corresponding to elongated bubbles,qv2′′, is given by[9, 11]
(51)qv2′′=qls′′ALB′′LlsLLB+qdry′′ALB′′(1−LlsLLB),
where ALB′′ is the wall area density for elongated
bubbles, Lls is the length of the liquid sublayer
underneath an elongated bubble, qls′′ is the corresponding average heat flux across
the liquid layer, and qdry′′ is  the heat transfer rate per unit wall in the dry region.The
length of large bubbles has been measured by Gersey and Mudawar [12] and
found to agree very well with the critical wavelength of the Helmholtz
instability vapor/liquid interface. This
length can be calculated from
(52)LLB=2πσρl(vv−vl)2(1+ρlρv).Another
parameter characterizing elongated bubbles is their distance from the wall. Assuming that the initial distance (at the
tip of the bubbles) corresponds to the viscous sublayer thickness [13], we obtain
(53)LLB=10vlρlτw,
where τw is the wall shear stress.


### 5.3. Wall Temperature Excursion Criterion

It can
be readily noticed that as the total wall heat flux increases, the elongated
bubble area density will also increase and so will the dry region under the
elongated bubbles. Consequently, the
local heat flux in the elongated bubble region,qv2′′,
will go down and that in the dispersed bubble region, qv1′′,
will go up. With the corresponding area density decreasing, nucleate boiling in
the region between elongated bubbles will intensify, eventually leading to a
situation when replenishment of the near-wall region with liquid will be no
longer possible. As a result, a sudden wall temperature excursion will occur. The critical condition beyond which small
bubbles cannot be removed away from the wall and replaced by fresh liquid can
be written as
(54)π6ddet2N′′≤mαNB=mαv1Av1′′,
where αNB is the volumetric fraction of dispersed
bubbles in the nucleate boiling region and αv1 is the local near-wall volume fraction (per
unit volume of the channel) of the dispersed bubbles.Naturally,
the local near-wall void fraction is the sum of the partial vapor volume
fractions of the dispersed bubbles and the elongated bubbles:(55)α=αv1+αv2.


## 5.1. Physical Concept

The wall
temperature excursion (or CHF) in low-quality boiling is associated with the
ability of the bubbles formed at the nucleation sites to depart from the wall,
so that the vacated space can be filled with fresh liquid, so that the
quenching and bubble growth processes can continue [9, 11]. Thus, two conditions must
be satisfied simultaneously to avoid wall temperature excursion, one concerned
with dispersed bubble concentration in the near-wall region, the other with the
velocity at which the bubble departs from the wall to make room for the next
generation bubble formed at the same nucleation site. Since the size of bubbles
departing from the nucleation sites is normally very small (of the order of 1 mm or less), their shape is nearly spherical. Thus, the two conditions
mentioned above can be formulated as follows.(1) The bubble maximum void fraction of dispersed bubbles in the near-wall
region cannot exceed the maximum packing factor for spherical particles, the
theoretical value of which (obtained from geometrical considerations) is(45)α≤αmax=π32≈0.74.(2) Taking into account that to form a new undeformed bubble at a given
cavity at the time of detachment, the distance between the previously formed
bubble at the same cavity and the heated wall at the same time instant must be
at least equal to the bubble diameter, the maximum evaporation heat flux to
avoid temperature excursion must satisfy the following condition:(46)qe′′hfg=π6ddet3ρvfdetN′′≤ρvudet,minα,
where
(47)udet,min=mddetfdet.The
coefficient,m≥1,
reflects the fact that in reality the axial motion of bubbles may require a
distance from the wall larger than one bubble diameter to avoid the coalescence with other
bubbles formed in the adjacent upstream cavities.It is interesting to notice that (46) can be rewritten as
(48)π6ddet2N′′≤mα.


## 5.2. Near-Wall Heat Transfer in the Presence of Elongated Bubbles

A schematic illustrating the near-wall conditions in the presence of elongated
bubbles is shown in Figure3. It follows from the previous discussion that at
any axial location along the channel the time-dependent wall temperature
fluctuates periodically, increasing in the poor heat transfer (elongated bubble)
region and decreasing, where heat transfer is more efficient (in the nucleate
boiling region).Figure 3
Phase distribution in the near-wall region.According
to Figure3, averaging the local wall heat flux over both regions always yields
the given constant heating rate of the heater per unit surface area. In
particular, the overall average heat flux can be partitioned into the following
terms:
(49)qH′′=qv1′′+qv2′′.In (49), qv1′′ is the nucleate boiling component
corresponding to small dispersed bubbles:
(50)qv1′′=qNB′′Av1′′,
where qNB′′ is given by (20), and Av1′′ is the area density (wall area fraction) for
the small bubble region.The heat
flux component corresponding to elongated bubbles,qv2′′, is given by[9, 11]
(51)qv2′′=qls′′ALB′′LlsLLB+qdry′′ALB′′(1−LlsLLB),
where ALB′′ is the wall area density for elongated
bubbles, Lls is the length of the liquid sublayer
underneath an elongated bubble, qls′′ is the corresponding average heat flux across
the liquid layer, and qdry′′ is  the heat transfer rate per unit wall in the dry region.The
length of large bubbles has been measured by Gersey and Mudawar [12] and
found to agree very well with the critical wavelength of the Helmholtz
instability vapor/liquid interface. This
length can be calculated from
(52)LLB=2πσρl(vv−vl)2(1+ρlρv).Another
parameter characterizing elongated bubbles is their distance from the wall. Assuming that the initial distance (at the
tip of the bubbles) corresponds to the viscous sublayer thickness [13], we obtain
(53)LLB=10vlρlτw,
where τw is the wall shear stress.


## 5.3. Wall Temperature Excursion Criterion

It can
be readily noticed that as the total wall heat flux increases, the elongated
bubble area density will also increase and so will the dry region under the
elongated bubbles. Consequently, the
local heat flux in the elongated bubble region,qv2′′,
will go down and that in the dispersed bubble region, qv1′′,
will go up. With the corresponding area density decreasing, nucleate boiling in
the region between elongated bubbles will intensify, eventually leading to a
situation when replenishment of the near-wall region with liquid will be no
longer possible. As a result, a sudden wall temperature excursion will occur. The critical condition beyond which small
bubbles cannot be removed away from the wall and replaced by fresh liquid can
be written as
(54)π6ddet2N′′≤mαNB=mαv1Av1′′,
where αNB is the volumetric fraction of dispersed
bubbles in the nucleate boiling region and αv1 is the local near-wall volume fraction (per
unit volume of the channel) of the dispersed bubbles.Naturally,
the local near-wall void fraction is the sum of the partial vapor volume
fractions of the dispersed bubbles and the elongated bubbles:(55)α=αv1+αv2.


## 6. Model Testing and Validation

### 6.1. Nucleate Boiling

The
theoretical model of bubble ebullition discussed in Sections3 and 4 was used
to evaluate various parameters characterizing nucleate boiling, and compare the
result of predictions against the results of measurements. Typical results are
shown in Tables 1 and 2 and in Figure 4.Table 1
Predicted bubble diameter at detachment.
ddet [mm] (23)ddet [mm]
Current model at p=5MPaddet [mm]
Current model at p=7MPaddet [mm]
Current model at p=10MPa1.11.161.020.68Table 2
Predicted nucleation frequency.
Pressure [MPa]fdetddet [m1/2/s]
Ceumern-Lindenstjerna [12]fdetddet [m1/2/s] Current model53.673.6273.543.24103.472.67Figure 4
A comparison between the predicted and experimental evaporation heat flux [14].The
predictions for the bubble diameter at departure and the nucleation frequency
are shown in Tables1 and 2. First, the predicted bubble diameter has been
compared against the expression given by (24), obtained for a liquid velocity
of 0.2 m/s. The calculations using the
present models have been performed for a wall heat flux of qw′′=4.73⋅103kW/m2, a subcooling of 10°C,
and three different pressures. The results are shown in Table 1 [9].As can
be seen, whereas the predicted bubble diameter at detachment compares well
against the simple correlation of Tolubinsky and Kostanchuk [7], the
current model also allows for quantifying the effect of system pressure (and
other parameters, such as velocity) which is ignored in the expression given by
(24).Next, the calculated
bubble departure frequency was compared against the correlation of Ceumern-Lindenstjerna [15]. The calculations were performed for qw′′=1.4⋅103kW/m2. The results for three different pressures are
shown in Table 2.The predicted evaporation heat flux in a boiling channel
used by Bartolomei and Chanturia [13] is shown in Figure 4. In this
comparison, the nucleation sites density was obtained from the expression
proposed by Lemmert and Chawla [8] (see (25)).
The measured wall superheat was about 9.5°C but changed
slightly along the channel. In order to
account for the uncertainties in the wall temperature measurements, the results
for two fixed values are shown in Figure 4. As can be seen, a very good
agreement has been obtained, especially for liquid subcooling between 5°C
and 15°C.It can
be noticed that the current model is capable of quantifying the effect of
various physical parameters which are normally not accounted for when a
phenomenological approach, based solely on correlating experimental data, is
used.


### 6.2. Temperature Excursion

To
perform predictions of the temperature excursion, the three-dimensional
three-field model of two-phase flow discussed in Section2 has been combined
with both models discussed in Sections 3–5. The
calculations have been performed for the experimental conditions of Hino and Ueda [16], in which R113 was used at a pressure of 147 kPa. The heated test
section was 0.357 m long and 0.018 m ID tube, with a centrally located heated
rod, 0.008 m diameter. The outer tube wall was insulated, and there was an
unheated section installed upstream of the annulus-shaped heated section,
allowing the flow to reach fully developed conditions at the entrance to the
heater.Typical
radial and axial distributions of various local flow parameters are shown in
Figure5. The axial distributions of the near-wall concentrations of both small
bubbles and elongated bubbles are shown in Figure 6. As it can be seen, the
rate of bubble concentration increase in the dispersed-bubble region
accelerates as the elongated bubbles start being formed. This is due to a
dramatic reduction in the heat transfer removal rate in the elongated bubble
region, which in
turn increases the local heat flux in the nucleate-boiling region. As it can be
seen in Figure 6, whereas the average wall heat flux is fixed at a level of 241kW/m2, the nucleate boiling heat flux experiences a dramatic growth
as the concentration of elongated bubbles stars increasing. Indeed, due to the
combined effects of increasing evaporation rate and increasing fraction of the
wall area occupied by elongated bubbles, the local void fraction approached the
value, α=0.75,
which already exceeded the limit given by (44). This, in turn, stopped the wall
replenishment by liquid phase and caused a sudden temperature excursion. Thus,
one concludes that the assumed wall heat flux was just beyond the onset of CHF.
Converting the difference between the actual volume fraction value of 0.75 and
the critical value of 0.74 into the corresponding power level difference yields
the critical heat flux (CHF) of about 238kW/m2. Since the uniform
heat flux of the heater used in the calculations (241kW/m2)
corresponded to the experimental onset of temperature excursion, the estimated
prediction error in the present case was less than 2%. Similar calculations performed for other
conditions produced errors of the order of ±15%.Figure 5
Color contours of the
calculated phasic velocities, temperature, and void fraction in subcooled
boiling:massflux=512kg/m2s, heatflux=241kW/m2,
inlet subcooling = 30 K.Axial
distributions of near-wall channel parameters for the conditions shown in
Figure5: (a) the concentrations of small bubbles and elongated bubbles, (b) a
comparison between the heat transfer removal rate in the elongated bubble
region, average wall heat flux, and critical heat flux.
(a)(b)Figures7 and 8 show the calculated CHF for various flow conditions, in comparison with
typical data based on experimental correlations. The predicted critical heat
flux is within the range of the measured values for heated channels operating
at subcooled boiling or low-quality boiling conditions. What is particularly
important is that
the predicted trends in the critical heat flux agree well with the existing
experimental evidence. Specifically, qc′′ gradually decreases with increasing velocity
(and, thus, flow rate) and decreases with increasing vapor concentration (flow
quality).Figure 7
Typical calculated values of critical heat flux, compared against experimental data for subcooled and low-quality boiling of water at 4.5 MPa [17].Figure 8
Typical calculated values of critical heat flux, compared against experimental data for subcooled and low-quality boiling of water at 4.5 MPa [17].


## 6.1. Nucleate Boiling

The
theoretical model of bubble ebullition discussed in Sections3 and 4 was used
to evaluate various parameters characterizing nucleate boiling, and compare the
result of predictions against the results of measurements. Typical results are
shown in Tables 1 and 2 and in Figure 4.Table 1
Predicted bubble diameter at detachment.
ddet [mm] (23)ddet [mm]
Current model at p=5MPaddet [mm]
Current model at p=7MPaddet [mm]
Current model at p=10MPa1.11.161.020.68Table 2
Predicted nucleation frequency.
Pressure [MPa]fdetddet [m1/2/s]
Ceumern-Lindenstjerna [12]fdetddet [m1/2/s] Current model53.673.6273.543.24103.472.67Figure 4
A comparison between the predicted and experimental evaporation heat flux [14].The
predictions for the bubble diameter at departure and the nucleation frequency
are shown in Tables1 and 2. First, the predicted bubble diameter has been
compared against the expression given by (24), obtained for a liquid velocity
of 0.2 m/s. The calculations using the
present models have been performed for a wall heat flux of qw′′=4.73⋅103kW/m2, a subcooling of 10°C,
and three different pressures. The results are shown in Table 1 [9].As can
be seen, whereas the predicted bubble diameter at detachment compares well
against the simple correlation of Tolubinsky and Kostanchuk [7], the
current model also allows for quantifying the effect of system pressure (and
other parameters, such as velocity) which is ignored in the expression given by
(24).Next, the calculated
bubble departure frequency was compared against the correlation of Ceumern-Lindenstjerna [15]. The calculations were performed for qw′′=1.4⋅103kW/m2. The results for three different pressures are
shown in Table 2.The predicted evaporation heat flux in a boiling channel
used by Bartolomei and Chanturia [13] is shown in Figure 4. In this
comparison, the nucleation sites density was obtained from the expression
proposed by Lemmert and Chawla [8] (see (25)).
The measured wall superheat was about 9.5°C but changed
slightly along the channel. In order to
account for the uncertainties in the wall temperature measurements, the results
for two fixed values are shown in Figure 4. As can be seen, a very good
agreement has been obtained, especially for liquid subcooling between 5°C
and 15°C.It can
be noticed that the current model is capable of quantifying the effect of
various physical parameters which are normally not accounted for when a
phenomenological approach, based solely on correlating experimental data, is
used.


## 6.2. Temperature Excursion

To
perform predictions of the temperature excursion, the three-dimensional
three-field model of two-phase flow discussed in Section2 has been combined
with both models discussed in Sections 3–5. The
calculations have been performed for the experimental conditions of Hino and Ueda [16], in which R113 was used at a pressure of 147 kPa. The heated test
section was 0.357 m long and 0.018 m ID tube, with a centrally located heated
rod, 0.008 m diameter. The outer tube wall was insulated, and there was an
unheated section installed upstream of the annulus-shaped heated section,
allowing the flow to reach fully developed conditions at the entrance to the
heater.Typical
radial and axial distributions of various local flow parameters are shown in
Figure5. The axial distributions of the near-wall concentrations of both small
bubbles and elongated bubbles are shown in Figure 6. As it can be seen, the
rate of bubble concentration increase in the dispersed-bubble region
accelerates as the elongated bubbles start being formed. This is due to a
dramatic reduction in the heat transfer removal rate in the elongated bubble
region, which in
turn increases the local heat flux in the nucleate-boiling region. As it can be
seen in Figure 6, whereas the average wall heat flux is fixed at a level of 241kW/m2, the nucleate boiling heat flux experiences a dramatic growth
as the concentration of elongated bubbles stars increasing. Indeed, due to the
combined effects of increasing evaporation rate and increasing fraction of the
wall area occupied by elongated bubbles, the local void fraction approached the
value, α=0.75,
which already exceeded the limit given by (44). This, in turn, stopped the wall
replenishment by liquid phase and caused a sudden temperature excursion. Thus,
one concludes that the assumed wall heat flux was just beyond the onset of CHF.
Converting the difference between the actual volume fraction value of 0.75 and
the critical value of 0.74 into the corresponding power level difference yields
the critical heat flux (CHF) of about 238kW/m2. Since the uniform
heat flux of the heater used in the calculations (241kW/m2)
corresponded to the experimental onset of temperature excursion, the estimated
prediction error in the present case was less than 2%. Similar calculations performed for other
conditions produced errors of the order of ±15%.Figure 5
Color contours of the
calculated phasic velocities, temperature, and void fraction in subcooled
boiling:massflux=512kg/m2s, heatflux=241kW/m2,
inlet subcooling = 30 K.Axial
distributions of near-wall channel parameters for the conditions shown in
Figure5: (a) the concentrations of small bubbles and elongated bubbles, (b) a
comparison between the heat transfer removal rate in the elongated bubble
region, average wall heat flux, and critical heat flux.
(a)(b)Figures7 and 8 show the calculated CHF for various flow conditions, in comparison with
typical data based on experimental correlations. The predicted critical heat
flux is within the range of the measured values for heated channels operating
at subcooled boiling or low-quality boiling conditions. What is particularly
important is that
the predicted trends in the critical heat flux agree well with the existing
experimental evidence. Specifically, qc′′ gradually decreases with increasing velocity
(and, thus, flow rate) and decreases with increasing vapor concentration (flow
quality).Figure 7
Typical calculated values of critical heat flux, compared against experimental data for subcooled and low-quality boiling of water at 4.5 MPa [17].Figure 8
Typical calculated values of critical heat flux, compared against experimental data for subcooled and low-quality boiling of water at 4.5 MPa [17].


## 7. Conclusions

Several aspects of
mechanistic multidimensional modeling and computer simulations of two-phase
flows and boiling heat transfer have been discussed. The specific models
included the mechanisms of local-subcooled boiling heat transfer in forced convection
flows, a mechanistic approach to bubble ebullition cycle, and criteria for
temperature excursion (CHF) in low-quality flows.The results of model testing indicate that the proposed
physical mechanisms are consistent with the existing experimental evidence
regarding the phase and temperature distributions, and wall temperature
excursion. The current approach is particularly suitable for implementation in
general CFD computational models using a multifield concept of two-phase flow.



---

*Source: 387020-2008-12-22.xml*