# Modeling of Multisize Bubbly Flow and Application to the Simulation of Boiling Flows with the Neptune_CFD Code

Authors: Christophe Morel; Jérôme M. Laviéville
Journal: Science and Technology of Nuclear Installations
 (2009)

---

## Abstract

This paper describes the modeling of boiling multisize bubbly flows and its application to the simulation of the DEBORA experiment. We follow the method proposed originally by Kamp, assuming a given mathematical expression for the bubble diameter pdf. The original model is completed by the addition of some new terms for vapor compressibility and phase change. The liquid-to-interface heat transfer term, which essentially determines the bubbles condensation rate in the DEBORA experiment, is also modeled with care. First numerical results realized with the Neptune_CFD code are presented and discussed.


---

## Body

## 1. Introduction

This paper
describes the modeling of
boiling multisize bubbly flows and its application to
the simulation of the DEBORA experiment (see [1]). First of all, the two-fluid model we use
for the calculation of boiling bubbly flows is presented (see
[2, 3]). In its first version, this two-fluid model assumed
that the bubbles can be characterized by a single bubble diameter which is the
classical Sauter mean diameter. Of course, this diameter is not a constant, but
is determined as a function of the bubbles interfacial area concentration and
the void fraction. As a consequence, the mean bubble diameter evolves according
to several physical phenomena: bubble coalescence and breakup, but also phase
change phenomena, including bubble nucleation and vapor compressibility inside
the bubbles. All these phenomena are taken into account as different source
terms in the interfacial area concentration balance equation.In real flows, a
bubble diameter spectrum is often observed, rather than a single bubble
diameter. In order to take this diameter spectrum into account, Kamp et al. (see
[4–6]) introduced a log-normal law
in order to mathematically describe the bubble probability density function
(pdf). As a consequence, the bubble diameter pdf is a function of only two
parameters which can be analytically expressed as functions of two particular
moments of the bubble diameter distribution function. The interfacial area concentration
is one of these two moments. Therefore, only one additional moment transport
equation is necessary in order to be able to reconstruct the bubble diameter
distribution function, hence extending our previous models to the multisize
effects appearing in boiling bubbly flows. The original work of Kamp was
devoted to the study of bubble coalescence in adiabatic flows under
microgravity. Here, we extend his model ability to boiling bubbly flows by
taking into account the vapor compressibility as well as phase change terms in
the different balance equations.The outline of
the paper is the following. The classical two-fluid model for boiling bubbly
flows is summarized in Section2. In Section
3, the original model of Kamp is
briefly recalled. This model is extended in Section 4 in order to take into
account the vapor compressibility as well as the phase change effects,
including wall nucleation, into the different balance equations. Section 5 is
devoted to the modeling of the liquid-to-interface heat transfer term including
the multisize effect of the bubbles. Numerical simulations of one DEBORA
experimental test realized with the Neptune_CFD code are illustrated in Section
6, before concluding and making some perspectives for future work in
Section 7.
This activity was realized in the context of the European NURESIM project.


## 2. Two-Fluid Model for Boiling Bubbly Flow

The two-fluid model we use for our boiling bubbly flow
calculations is constituted of the following six balance equations
(e.g., [3]):(i)
two mass balance equations:(1)∂αkρk∂t+∇⋅(αkρkV¯k)=Γkk=L,G, where t is the time, αk,ρk,V¯k denote the volumetric fraction of phase k, its averaged
density and velocity and Γk is the interfacial mass transfer per unit volume and unit time; the phase index
k takes the values L for the liquid phase and G for the gas phase;(ii)
two momentum balance equations:(2)∂αkρkV¯k∂t+∇⋅(αkρkV¯kV¯k)=−αk∇¯p+M¯k+αkρkg¯+∇⋅[αk(τ¯¯k+τ¯¯kT)]k=L,G,where p is the pressure, g is the gravity
acceleration, M¯k is the interfacial momentum transfer per
unit volume and unit time, and τ¯¯k and τ¯¯kT denote the molecular and turbulent stress
tensors, the latter being also called the Reynolds stress tensor;(iii)
two total enthalpy balance equations:(3)∂∂t[αkρk(hk+Vk22)]+∇⋅(αkρk(hk+Vk22)V¯k)=αk∂p∂t+αkρkg¯⋅V¯k+Γk(hki+Vk22)+qki′′ai+qwk′′′−∇⋅[αk(q¯k+q¯kT)]k=L,G, where hk is the phase-averaged enthalpy for phase k and hki is the
interfacial-averaged enthalpy. We have assumed that the two phases are governed
by the same averaged pressure field p and we make no distinction between the
pressures in the two phases or between the bulk pressure and the interface
pressure for simplicity. The three terms Γk, Mk, and qki′′ai denote the interfacial transfer terms of mass,
momentum, and heat, the quantity ai being the interfacial area
concentration. The terms qwk′′′ denote the wall-to-fluid heat transfer per
unit volume and unit time for each phase. The heterogeneous (wall) nucleation
of bubbles is included as a part of the term Γk (the bubble heterogeneous
nucleation corresponds to an interfacial transfer of mass between phases, even
if it is located onto the heated wall surface). The two terms q¯k and q¯kT denote the molecular and turbulent heat fluxes
inside phase k. Except for the liquid-to-interface heat transfer term qLi′′ai (see
Section 5), the modeling of the different averaged transfer terms
(interfacial, turbulent, and wall transfer terms) is largely discussed in our
previous papers (see [2, 3, 7]). Here, we
only present what has been changed through the introduction of the bubble multisize
modeling.


## 3. Kamp’s Model for the Coalescence of Multisize Bubbles

Kamp et al. (see
[4–6]) propose to model the bubble
size distribution function by a log-normal law. The bubbles are assumed to
remain spherical and are characterized by their diameter d. The γ’s moment of the bubble distribution
function and the log-normal law for the bubble diameter pdf are given by(4)Sγ=^n∫P(d)dγδd,P(d)=12πσ^de−{ln(d/d00)}2/2σ^2,where P(d)
denotes the bubble diameter pdf and the function f=nP is the bubble diameter
distribution function, n being the total number of bubbles per unit volume of
the two-phase mixture. It should be noted that n and P(d) are local
instantaneous quantities (defined in the sense of the ensemble average). In
fact, n and P(d) should be noted as n(x̲,t) and P(d;x̲,t) but we
have omitted their dependence in the position x̲ and time t for
simplicity. The diameter pdf P depends on (x̲,t) through the dependence
of the two parameters d00(x̲,t) and σ^(x̲,t). It should be noted that d, x̲,
and t are independent variables. The first four remarkable moments of the
bubble diameter distribution function are related to more usual quantities by
the following relations:(5)n=^S0,d10=^S1n,ai=^πS2,α=^πS36, where d10 is a mean bubble diameter, ai is the interfacial area concentration,
and α is the void fraction (vapor volumetric
fraction). An infinity of mean diameters can be constructed by using the
following definition:(6)dγδ=^(SγSδ)1/(γ−δ).The diameter d10 is one of them, the Sauter mean diameter d32 constructed with the second
and third moments is another one. The two parameters of the diameter pdf are
analytically expressed by the following functions of the particular moments S1 and S2 and of the void fraction α (see
[5]):(7)σ^=ln(6αS1πS22),d00=6απS2e−5σ^2/2.Then, with these
two parameters being determined, any moment defined by the first relation (4)
can be expressed analytically by the following relation:(8)Sγ=6απd00γ−3eσ^2/2(γ2−9).With the void
fraction field α being part of the solution of the equations of
the two-fluid model (1)–(3), it is
sufficient then to use two additional balance equations for the two particular
moments S1 and S2 to completely close the geometrical
quantities (provided that these two equations are themselves closed). The departure
point is the so-called Liouville equation for the bubble diameter distribution
function which can be derived in a very general manner as(9)∂nP(d)∂t+∇⋅[nP(d)v¯(d)]=−∂∂d[GnP(d)]+φ(d,x¯,t),where v̲(d)
is the velocity of a bubble having a diameter d, the quantity G(d;x̲,t)
is the Lagrangian derivative of the bubble diameter measured along the bubble
path. It represents the continuous bubble size variation due to vapor
compressibility and phase change. The last term φ(d,x¯,t) represents all phenomena which are responsible
for discontinuous bubble diameter changes, namely the bubble coalescence and
breakup as well as bubble nucleation and collapse. Multiplying (9) by dγ and integrating the resulting
equation with respect to d from zero to infinity, the following balance
equation for the γ’s moment is obtained:(10)∂Sγ∂t+∇⋅(Sγv¯γ)=γGγ−1Sγ−1+φγ,which is used
for γ=1,2. In this equation, v̲γ, Gγ,
and φγ are defined by the following relations:(11)v¯γ=^∫0∞v¯(d)dγP(d)δd∫0∞dγP(d)δd,Gγ=^∫0∞G(d)dγP(d)δd∫0∞dγP(d)δd,φγ=^∫0∞φ(d)dγδd.In the original
model of Kamp, the term proportional to Gγ−1 is neglected and the
coalescence is the only effect taken into account in the modeling of φγ.
At the end, the modeled equations for S1 and S2  read (see [5])(12)∂S1∂t+∇⋅(S1v¯G)=8π3(6απ)2/3CtεL1/31.61(21/3−2)S14/3f(1,σ^,P00),∂S2∂t+∇⋅(S2v¯G)=8π3(6απ)1/3CtεL1/31.61(22/3−2)S25/3f(2,σ^,P00), where the two
velocities v̲1 and v̲2 have been assumed to
be equal to the averaged gas velocity v̲G. In the right-hand
side of (12), εL is the dissipation rate of
liquid turbulence, Ct is the ratio of the dispersed phase turbulent
fluctuations by the continuous phase turbulent fluctuations, and f(1,σ^,P00) and f(2,σ^,P00) are numerical functions of the width parameter σ^ and the coalescence probability (based on the
molecular chaos assumption) P00 of two bubbles having a diameter d00 (see the original work of Kamp for more information).


## 4. Extension of Kamp’s Model to Vapor Compressibility and Phase Change

In the original
work of Kamp, the termG in (9) was neglected. Here we calculate the
corresponding terms Gγ−1 appearing in (10) when G is due
to the vapor compressibility (Section 4.1) and to phase change
(Section 4.2).


### 4.1. Vapor Compressibility

LetD(x̲,t)
be the bubble diameter in physical space (the notation d stands for the bubble
diameter in phase space). In the absence of phase change, the bubble mass is
conserved along its trajectory. As a consequence, the bubble diameter variation
and the vapor density variation measured along the bubble path are related
through(13)DGDDt=^∂D∂t+v¯⋅∇D=−D3ρ(∂ρ∂t+v¯⋅∇ρ),where ρ denotes the vapor density and v¯ the
bubble velocity. The quantity G appearing in (9) is the conditional expectation
of the Lagrangian derivative given by (13) conditioned by the equality D(x¯,t)=d:(14)G=〈DGDDt∣D=d〉=−d3〈1ρ(∂ρ∂t+v¯⋅∇ρ)∣D=d〉≅−d31ρG(∂ρG∂t+v¯G⋅∇ρG).If we assume
that the vapor density ρ and the bubble velocity v̲ do not depend
on the bubble diameter D, the conditional average appearing in the second
expression of (14) can be replaced by the unconditional one. The last equality
in (14) simply assumes that we neglect nonlinear effect when we average the
term under the brackets, therefore keeping only first order effects. This is
clearly an approximation where the microscopic quantities ρ and v̲ characterizing the bubble are
replaced by the macroscopic (averaged) ones ρG and v¯G.Substituting the
expression (14) into (10) and (11) gives the following expressions
for the effect of the vapor compressibility on the two particular moments
considered:(15)[∂S1∂t+∇⋅(S1v¯G)]compressibility=G0S0=−S13ρG(∂ρG∂t+v¯G⋅∇ρG),[∂S2∂t+∇⋅(S2v¯G)]compressibility=2G1S1=−23S2ρG(∂ρG∂t+v¯G⋅∇ρG).


### 4.2. Phase Change

In the previous
paragraph, we made the assumption of no phase change in order to derive the
compressibility terms in the equations forS1,2. Here we assume a
constant vapor density ρ=ρG in order to derive the
phase change terms in the same equations. At the end, the two effects are taken
into account simultaneously by assuming that the corresponding terms can be
linearly superimposed into the S1,2 balance equations. For a
constant vapor density, the bubble diameter variation in physical space is
simply given by(16)DGDDt=^∂D∂t+v¯.∇D=2m˙ρG,where m˙ denotes the bubble mass gain per unit surface
per unit time due to vaporization or condensation through the bubble surface.
Substituting (16) into (14) gives(17)G=〈DGDDt∣D=d〉=〈2m˙ρG∣D=d〉=2ρG〈m˙∣D=d〉.The closure
issue consists to propose an expression for the conditional expectation 〈m˙∣D=d〉.
The simplest choice is to assume that 〈m˙∣D=d〉 is simply given by the ratio (ΓG−ΓN)/aI, where ΓG is the vapor source per
unit volume per unit time and ΓN is the part of this source
due to newly nucleated bubbles, therefore ΓG−ΓN is the part of the phase
change through the surfaces of the already existing bubbles. Making this simple
choice, we obtain the following expressions for the effect of the phase change
on the two particular moments considered:(18)[∂S1∂t+∇⋅(S1v¯G)]phasechange=G0S0=2ρGΓG−ΓNaIS0,[∂S2∂t+∇⋅(S2v¯G)]phasechange=2G1S1=4ρGΓG−ΓNaIS1, where S0 is the volumetric number of bubbles n obtained by making γ=0 into the relation (8). A more sophisticated model taking into account the
dependence of 〈m˙∣D=d〉 on the diameter d will be presented in the
next section.


### 4.3. A Few Words on Nucleation

In subcooled
boiling flows, as this is the case for the DEBORA experiment, vapor bubbles are
essentially nucleated onto the heated wall surface (heterogeneous nucleation)
rather than in the liquid bulk, because the liquid is subcooled in the major
part of the flow. Therefore, we can neglect a possible homogeneous nucleation
in comparison to the heterogeneous nucleation. The newly nucleated bubbles are
supposed to be generated with a unique size: the so-called detachment diameterdd. In our model, this detachment diameter is calculated with the
help of Ünal’s correlation (see [8]), used together with the model of
Kurul and Podowski [9] for the active nucleation site density. These models
are used to estimate the vapor produced in the meshes adjacent to the wall due
to nucleation ΓN. Dividing ΓN by the mass of a nucleated
bubble ρGπdd3/6, the
volumetric number of bubbles source term is obtained. The corresponding source
terms in the S1 and S2 equations are simply obtained by
multiplying this volumetric number of bubbles source term by the diameter dd or by its square dd2, according to the order of the
moment considered.


## 4.1. Vapor Compressibility

LetD(x̲,t)
be the bubble diameter in physical space (the notation d stands for the bubble
diameter in phase space). In the absence of phase change, the bubble mass is
conserved along its trajectory. As a consequence, the bubble diameter variation
and the vapor density variation measured along the bubble path are related
through(13)DGDDt=^∂D∂t+v¯⋅∇D=−D3ρ(∂ρ∂t+v¯⋅∇ρ),where ρ denotes the vapor density and v¯ the
bubble velocity. The quantity G appearing in (9) is the conditional expectation
of the Lagrangian derivative given by (13) conditioned by the equality D(x¯,t)=d:(14)G=〈DGDDt∣D=d〉=−d3〈1ρ(∂ρ∂t+v¯⋅∇ρ)∣D=d〉≅−d31ρG(∂ρG∂t+v¯G⋅∇ρG).If we assume
that the vapor density ρ and the bubble velocity v̲ do not depend
on the bubble diameter D, the conditional average appearing in the second
expression of (14) can be replaced by the unconditional one. The last equality
in (14) simply assumes that we neglect nonlinear effect when we average the
term under the brackets, therefore keeping only first order effects. This is
clearly an approximation where the microscopic quantities ρ and v̲ characterizing the bubble are
replaced by the macroscopic (averaged) ones ρG and v¯G.Substituting the
expression (14) into (10) and (11) gives the following expressions
for the effect of the vapor compressibility on the two particular moments
considered:(15)[∂S1∂t+∇⋅(S1v¯G)]compressibility=G0S0=−S13ρG(∂ρG∂t+v¯G⋅∇ρG),[∂S2∂t+∇⋅(S2v¯G)]compressibility=2G1S1=−23S2ρG(∂ρG∂t+v¯G⋅∇ρG).


## 4.2. Phase Change

In the previous
paragraph, we made the assumption of no phase change in order to derive the
compressibility terms in the equations forS1,2. Here we assume a
constant vapor density ρ=ρG in order to derive the
phase change terms in the same equations. At the end, the two effects are taken
into account simultaneously by assuming that the corresponding terms can be
linearly superimposed into the S1,2 balance equations. For a
constant vapor density, the bubble diameter variation in physical space is
simply given by(16)DGDDt=^∂D∂t+v¯.∇D=2m˙ρG,where m˙ denotes the bubble mass gain per unit surface
per unit time due to vaporization or condensation through the bubble surface.
Substituting (16) into (14) gives(17)G=〈DGDDt∣D=d〉=〈2m˙ρG∣D=d〉=2ρG〈m˙∣D=d〉.The closure
issue consists to propose an expression for the conditional expectation 〈m˙∣D=d〉.
The simplest choice is to assume that 〈m˙∣D=d〉 is simply given by the ratio (ΓG−ΓN)/aI, where ΓG is the vapor source per
unit volume per unit time and ΓN is the part of this source
due to newly nucleated bubbles, therefore ΓG−ΓN is the part of the phase
change through the surfaces of the already existing bubbles. Making this simple
choice, we obtain the following expressions for the effect of the phase change
on the two particular moments considered:(18)[∂S1∂t+∇⋅(S1v¯G)]phasechange=G0S0=2ρGΓG−ΓNaIS0,[∂S2∂t+∇⋅(S2v¯G)]phasechange=2G1S1=4ρGΓG−ΓNaIS1, where S0 is the volumetric number of bubbles n obtained by making γ=0 into the relation (8). A more sophisticated model taking into account the
dependence of 〈m˙∣D=d〉 on the diameter d will be presented in the
next section.


## 4.3. A Few Words on Nucleation

In subcooled
boiling flows, as this is the case for the DEBORA experiment, vapor bubbles are
essentially nucleated onto the heated wall surface (heterogeneous nucleation)
rather than in the liquid bulk, because the liquid is subcooled in the major
part of the flow. Therefore, we can neglect a possible homogeneous nucleation
in comparison to the heterogeneous nucleation. The newly nucleated bubbles are
supposed to be generated with a unique size: the so-called detachment diameterdd. In our model, this detachment diameter is calculated with the
help of Ünal’s correlation (see [8]), used together with the model of
Kurul and Podowski [9] for the active nucleation site density. These models
are used to estimate the vapor produced in the meshes adjacent to the wall due
to nucleation ΓN. Dividing ΓN by the mass of a nucleated
bubble ρGπdd3/6, the
volumetric number of bubbles source term is obtained. The corresponding source
terms in the S1 and S2 equations are simply obtained by
multiplying this volumetric number of bubbles source term by the diameter dd or by its square dd2, according to the order of the
moment considered.


## 5. Heat and Mass Transfers through the Surface of the Existing Bubbles

For a relatively
low speed flow, the local instantaneous interfacial balance of total energy (see
[10, 11]) can be simplified:(19)(q¯G−q¯L)⋅n¯G−m˙(hGI−hLI)≅0.In this
relation, q̲G⋅n̲G and q̲L⋅n̲G are the heat flux densities on the two sides of the interface, n̲G being the unit vector normal to the interface directed towards the exterior of
the bubble. The relation (19) is valid at each point of the bubble interface.
The quantities hGI and hLI are the vapor and liquid
enthalpies at the bubble interface, they are generally assumed to be equal to
the saturation enthalpies. The relation (19) gives the link between the rate of
phase change m˙ and the two heat flux densities between the
two phases and the interface. In a bubbly flow with relatively small bubbles,
the temperature of the vapor inside the bubbles is generally close to the
saturation temperature. In subcooled liquid, the rate of bubble condensation is
essentially determined by the heat exchange on the liquid side, rather than the
one on the vapor side. Neglecting q̲G⋅n̲G in
comparison to q̲L⋅n̲G in the previous
equation and introducing a Nusselt number NuL to characterize
the liquid-to-interface heat transfer, (19) becomes(20)m˙≅q¯L⋅n¯LhGI−hLI=NuLλL(hGI−hLI)D(θL−θsat),where λL is the liquid thermal
conductivity and θL and θsat  represent the liquid and
saturation temperatures, respectively. The Nusselt number is given by the Ranz-Marschall law  [12] (21)NuL=2+0.55Re1/2PrL1/3,where Re and PrL are the bubble Reynolds number and the liquid Prandtl number, respectively.
Neglecting the fluctuations of the enthalpies hGI and hLI,
the averaging of the first relation (20) gives(22)qLI′′aI=^ΓG(hGI−hLI)withΓG=^∫〈m˙∣D=d〉πd2nP(d)δd,qLI′′aI=^∫〈q¯L⋅n¯L∣D=d〉πd2nP(d)δd, where ΓG and  qLI′′aI  are the average interfacial mass and heat
transfer terms appearing in (1) and (3). The term ΓG is related to the vapor
phase because m˙ is the bubble mass gain and qLI′′aI is the liquid-to-interface heat transfer per
unit volume per unit time, the heat transfer density qLI′′ being expressed per unit area per unit time, aI being the interfacial area concentration. Equation (22) is very useful because
it directly relates the mass transfer rate to the liquid-to-interface heat transfer
rate (under the assumption of a nearly saturated vapor inside the bubbles). The
averaging of the second relation (20) gives(23)qLI′′aI=^∫〈NuLλLD(θL−θsat)∣D=d〉πd2nP(d)δd.The relations
(21) and (23) are the basic tools to calculate the average interfacial liquid-to-interface
heat transfer. In what follows, we give three different approximated
expressions for this term. The first one is a single-size approximation of the
liquid-to-interface heat flux, it is given by(24)qLI′′aI≅∫πdNuLλL(θL−θsat)nP(d)δd≅aINuL(d32)d32λL(ΘL−Θsat).In this expression,
it is implicitly assumed that all the bubbles have the same diameter which is
given by the Sauter mean diameter d32=6α/aI (single-size
approximation). The relation (21) is directly calculated as a function of the
Sauter mean diameter and the average relative velocity of the bubbles into the
Reynolds number. The second (more precise) approximation consists in using the
log-normal law for the bubble diameter pdf (4). The result is(25)qLI′′aI=[2πλLS1+0.55πλLPrL1/3vL1/2|vr|1/2nd003/2exp(98σ^2)]×(ΘL−Θsat).The dependence
of the bubble Reynolds number Re on the bubble diameter d has been taken into
account in the derivation of (25), but not the variation of the relative
velocity |vr|, also appearing in the Reynolds number, as a function
of the bubble size. When using (25), the norm of the relative velocity is
assumed to be given by the averaged fields. In order to alleviate this
approximation, a third model has been constructed by assuming that the norm of
the relative velocity can be accurately given by the terminal velocity(26)|vr|≈vT=43ΔρgdρLCD,where Δρ is the absolute value of the density
difference |ρL−ρG| and CD is the
drag coefficient, which can be given for spherical bubbles by (e.g.,
[13])(27)CD=24Re(1+0.1Re3/4).Assuming that
the bubble Reynolds number is sufficiently high in order to approximate the
preceding relation by CD=2.4/Re1/4, we obtain the
following relation for the square root of the terminal velocity:(28)vT1/2=(Δρg1.8ρcvc1/4)2/7d5/14.Reevaluating the
liquid-to-interface heat transfer from (4),
(21), (23), and
(28) gives a third
model:(29)qLI′′aI=[2πλLS1+0.55πλLPrL1/3vL1/2(Δρg1.8ρLvL1/4)2/7×nd0013/7exp(16998σ^2)](ΘL−Θsat).Now, having the
model given by (20)-(21) for the bubble mass gain m˙,
we can evaluate more precisely the corresponding terms in the equations for the
two moments S1 and S2, namely the terms G0 and G1 with G given from (17). The results are the following:(30)[∂S1∂t+∇⋅(S1v¯G)]phasechange=[4ρGλLhGI−hLInd00exp(σ^22)+1.1ρGλLhGI−hLIPrL1/3vL1/2(Δρg1.8ρLvL1/4)2/7nd001/7exp(σ^298)]×(ΘL−Θsat),[∂S2∂t+∇⋅(S2v¯G)]phasechange=8ρGλLhGI−hLIn(ΘL−Θsat)+2.2ρGλLhGI−hLIPrL1/3vL1/2(Δρg1.8ρLvL1/4)2/7nd006/7exp(1849σ^2)×(ΘL−Θsat). To conclude the modeling sections, we can summarize the
complete modeling of boiling bubbly flows with multisize bubbles in its present
state. The basic equations are the mass, momentum and total enthalpy balance
equations written for the two phases, they are given by
(1)–(3). A multisize
model is used in conjunction with the two-fluid model to evaluate the evolution
of the bubble diameter spectrum under several physical effects: coalescence,
vapor compressibility, and phase-change. It should be noted that a bubble
breakup term is removed in the present state of the model. A log-normal law is
used to describe the bubble diameter pdf. This law involves two parameters d00 and σ^ which are analytically determined from two
particular moments of the bubble diameter pdf (7). Therefore, a set of two
additional transport equations are written for these two moments. These moments
balance equations have the general form (10) with their right-hand side including
a coalescence model (12), a vapor compressibility model
(15), and two possible
phase change models (18) or (30). The model also includes wall nucleation
source terms which are not presented here (see [2]). At the
end, we have proposed several expressions for the liquid-to-interface heat
transfer term appearing in the liquid enthalpy balance (3). This term is
responsible for the bubble condensation in the core of the flow. The first one
(24) is the classical expression used in usual (single size) bubbly flow models.
The two others (25) and
(29) take into account the bubble diameter spectrum as
modeled by the log-normal law. The last one (29) also takes into account the
dependence of the bubble relative velocity on the bubble diameter. All these
models should be tested. In the following section, we give the results of a
comparison of Neptune_CFD simulations with an experimental test of the DEBORA
database.


## 6. Numerical Simulations of a Boiling Bubbly Flow

The DEBORA experiment, carried out at the French Commissariat
à l’Energie Atomique (see [1]), has been chosen to evaluate
our model. In this experiment, the R-12 has been adopted as the working fluid
to simulate the PWR (pressurized water reactor) conditions under low pressure.
Some liquid R-12 flows upwardly inside a vertical pipe having an internal
diameter equal to 19.2 mm.
A part of the tube wall of 3.5 m long is electrically heated. Numerous vapor bubbles
are nucleated onto the heated wall surface, and condense in the core of the
flow, where the liquid is subcooled. At the end of the heated section, the
radial profiles of the void fraction, bubble diameter, and vertical bubble
velocity are measured by means of an optical probe and the liquid temperature
is measured by means of a thermocouple.Here we evaluate the different models presented above on a
particular DEBORA experimental test. The controlling parameters of this test
are the following.(i)
Test pressure: 14.59 bar.(ii)
Inlet mass flux density (liquid only):2027kg/m2s.(iii)
Liquid inlet temperature:28.52°C.(iv)
Wall to fluid heat flux density:76200W/m2.All the models have been implemented in the multiphase
Neptune_CFD code. With the flow being axi-symmetric, we use a two-dimensional
axi-symmetric calculation grid. This grid is characterized by 80 axial meshes
and 10 radial meshes. We assume that this relatively coarse grid is sufficient
for a first test of the different models. The models tested are the three
expressions of the liquid-to-interface heat flux (24),
(25), and (29) and the multisize
model for bubbly flows with the two different expressions (18) and (30) for the
phase change terms. The compared quantities are the void fraction
(Figure 1),
the vapor mean velocity (Figure 2), the liquid temperature
(Figure 3), the
interfacial area concentration (Figure 4),
the Sauter mean diameter d32
(Figure 5), and the mean diameter d20 defined by
(6)
(Figure 6). All
the radial profiles (experimental and calculated) have been taken at the end of
the heated length.Figure 1
comparison of the different models on the void fraction profile.Figure 2
comparison of the different models on the vapor velocity profile.Figure 3
comparison of the different models on the liquid temperature profile.Figure 4
comparison of the different models on the interfacial area profile.Figure 5
comparison of the different models on the Sauter mean diameter profile.Figure 6
comparison of the different models on the  mean diameterd20 profile.Each figure compares the results of five calculations
together with the experimental values.According to these comparisons, the model giving better
results on the void fraction, the vapor mean velocity, the liquid temperature,
and the interfacial area concentration seems to be the one using
(29) and (30)
with the interfacial averaged enthalpies being equal to the phase averaged
ones. It can also be seen that the void fraction is strongly overestimated for
all the models tested, but the interfacial area concentration is overestimated
too, therefore giving approximately the right values for the bubble Sauter mean
diameter and the mean diameter d20. The overestimation of the void
fraction could be attributed to the inadequacy of the Ranz and Marschall law
(21) to describe the condensation of bubbles in a subcooled liquid. This issue
will be studied in future work.


## 7. Conclusions

A multisize model for boiling bubbly flows has been presented
in detail. This model is not completely achieved (e.g., a bubble breakup term
is missing in the equations). Nevertheless, we made first calculations to
evaluate the capabilities of this model in its present state. Five different
calculations have been done and compared to a single experimental test of the
DEBORA experiment. It has been shown that the results are sensitive to the
expressions used for the liquid-to-interface heat transfer (24),
(25), or (29)
and to the phase change terms in the geometrical equations (18) or (30). The
model using the particular expressions (29) and (30), with the interfacial averaged
enthalpies being given by the phase averaged ones, gives better results. Unfortunately,
the void fraction and the interfacial area density are strongly overestimated. Further
investigations are needed to clarify these questions.Several issues can be raised for future developments:(i)
a bubble breakup model should be added for completeness of
the model;(ii)
a bubble collapse model is missing too;(iii)
the physical model given by Ranz and Marschall [12] (21)
could not be adapted for condensing bubbles in subcooled liquid. The validity
of this model should be evaluated first (perhaps by using DNS calculations).



---

*Source: 953527-2008-10-14.xml*