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

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.

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 Section 2. 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) where 𝑡 is the time, 𝛼_𝑘,𝜌_𝑘,𝑉_𝑘 denote the volumetric fraction of phase 𝑘, its averaged density and velocity and Γ_𝑘 is the interfacial mass transfer per unit volume and unit time; the phase index 𝑘 takes the values 𝐿 for the liquid phase and 𝐺 for the gas phase;

(ii) two momentum balance equations:𝜕𝛼_𝑘𝜌_𝑘𝑉_𝑘𝜕𝑡+∇⋅(𝛼_𝑘𝜌_𝑘𝑉_𝑘𝑉_𝑘)=−𝛼_𝑘∇𝑝+𝑀_𝑘+𝛼_𝑘𝜌_𝑘𝑔+∇⋅[𝛼_𝑘(𝜏_𝑘+𝜏^𝑇_𝑘)]𝑘=𝐿,𝐺, (2)where 𝑝 is the pressure, 𝑔 is the gravity acceleration, 𝑀_𝑘 is the interfacial momentum transfer per unit volume and unit time, and 𝜏_𝑘 and 𝜏^𝑇_𝑘 denote the molecular and turbulent stress tensors, the latter being also called the Reynolds stress tensor;

(iii) two total enthalpy balance equations:𝜕𝜕𝑡[𝛼_𝑘𝜌_𝑘(ℎ_𝑘+𝑉²_𝑘2)]+∇⋅(𝛼_𝑘𝜌_𝑘(ℎ_𝑘+𝑉²_𝑘2)𝑉_𝑘)=𝛼_𝑘𝜕𝑝𝜕𝑡+𝛼_𝑘𝜌_𝑘𝑔⋅𝑉_𝑘+Γ_𝑘(ℎ_𝑘𝑖+𝑉²_𝑘2)+𝑞^′′_𝑘𝑖𝑎_𝑖+𝑞^′′′_𝑤𝑘−∇⋅[𝛼_𝑘(𝑞_𝑘+𝑞^𝑇_𝑘)]𝑘=𝐿,𝐺, (3) where ℎ_𝑘 is the phase-averaged enthalpy for phase 𝑘 and ℎ_𝑘𝑖 is the interfacial-averaged enthalpy. We have assumed that the two phases are governed by the same averaged pressure field 𝑝 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 Γ_𝑘, 𝑀_𝑘, and 𝑞^′′_𝑘𝑖𝑎_𝑖 denote the interfacial transfer terms of mass, momentum, and heat, the quantity 𝑎_𝑖 being the interfacial area concentration. The terms 𝑞^′′′_𝑤𝑘 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 Γ_𝑘 (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 𝑞_𝑘 and 𝑞^𝑇_𝑘 denote the molecular and turbulent heat fluxes inside phase 𝑘. Except for the liquid-to-interface heat transfer term 𝑞^′′_𝐿𝑖𝑎_𝑖 (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 𝑑. The 𝛾’s moment of the bubble distribution function and the log-normal law for the bubble diameter pdf are given by𝑆_𝛾=𝑛𝑃(𝑑)𝑑^𝛾1𝛿𝑑,𝑃(𝑑)=√𝑒2𝜋𝜎𝑑^{−{ln(𝑑/𝑑00)}2/2}𝜎², (4)where 𝑃(𝑑) denotes the bubble diameter pdf and the function 𝑓=𝑛𝑃 is the bubble diameter distribution function, 𝑛 being the total number of bubbles per unit volume of the two-phase mixture. It should be noted that 𝑛 and 𝑃(𝑑) are local instantaneous quantities (defined in the sense of the ensemble average). In fact, 𝑛 and 𝑃(𝑑) should be noted as 𝑛(𝑥,𝑡) and 𝑃(𝑑;𝑥,𝑡) but we have omitted their dependence in the position 𝑥 and time 𝑡 for simplicity. The diameter pdf 𝑃 depends on (𝑥,𝑡) through the dependence of the two parameters 𝑑₀₀(𝑥,𝑡) and 𝜎(𝑥,𝑡). It should be noted that 𝑑, 𝑥, and 𝑡 are independent variables. The first four remarkable moments of the bubble diameter distribution function are related to more usual quantities by the following relations:𝑛=𝑆₀,𝑑₁₀𝑆=¹𝑛,𝑎_𝑖=𝜋𝑆₂,𝛼=𝜋𝑆₃6, (5) where 𝑑₁₀ is a mean bubble diameter, 𝑎_𝑖 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:𝑑_𝛾𝛿𝑆=(^𝛾𝑆_𝛿)^{1/(𝛾−𝛿)}. (6)The diameter 𝑑₁₀ is one of them, the Sauter mean diameter 𝑑₃₂ 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 𝑆₁ and 𝑆₂ and of the void fraction 𝛼 (see [5]):𝜎=ln(6𝛼𝑆₁𝜋𝑆²₂),𝑑₀₀=6𝛼𝜋𝑆₂𝑒⁻⁵𝜎²_/2. (7)Then, with these two parameters being determined, any moment defined by the first relation (4) can be expressed analytically by the following relation:𝑆_𝛾=6𝛼𝜋𝑑^𝛾−3₀₀𝑒𝜎²_/2(𝛾²₋₉₎. (8)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 𝑆₁ and 𝑆₂ 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)where 𝑣(𝑑) is the velocity of a bubble having a diameter 𝑑, the quantity 𝐺(𝑑;𝑥,𝑡) 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 𝜑(𝑑,𝑥,𝑡) 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 𝑑^𝛾 and integrating the resulting equation with respect to 𝑑 from zero to infinity, the following balance equation for the 𝛾’s moment is obtained:𝜕𝑆_𝛾𝜕𝑡+∇⋅(𝑆_𝛾𝑣_𝛾)=𝛾𝐺_𝛾−1𝑆_𝛾−1+𝜑_𝛾, (10)which is used for 𝛾=1,2. In this equation, 𝑣_𝛾, 𝐺_𝛾, and 𝜑_𝛾 are defined by the following relations:𝑣_𝛾∫=^∞0𝑣(𝑑)𝑑^𝛾𝑃(𝑑)𝛿𝑑∫^∞₀𝑑^𝛾,𝐺𝑃(𝑑)𝛿𝑑_𝛾∫=^∞0𝐺(𝑑)𝑑^𝛾𝑃(𝑑)𝛿𝑑∫^∞₀𝑑^𝛾,𝜑𝑃(𝑑)𝛿𝑑_𝛾=^∞₀𝜑(𝑑)𝑑^𝛾𝛿𝑑. (11)In the original model of Kamp, the term proportional to 𝐺_𝛾−1 is neglected and the coalescence is the only effect taken into account in the modeling of 𝜑_𝛾. At the end, the modeled equations for 𝑆₁ and 𝑆₂ read (see [5])𝜕𝑆₁𝜕𝑡+∇⋅(𝑆₁𝑣_𝐺)=8𝜋3(6𝛼𝜋)^2/3𝐶_𝑡𝜀_𝐿^1/3√1.61(2^1/3−2)𝑆₁^4/3𝑓(1,𝜎,𝑃₀₀),𝜕𝑆₂𝜕𝑡+∇⋅(𝑆₂𝑣_𝐺)=8𝜋3(6𝛼𝜋)^1/3𝐶_𝑡𝜀_𝐿^1/3√1.61(2^2/3−2)𝑆₂^5/3𝑓(2,𝜎,𝑃₀₀), (12) where the two velocities 𝑣₁ and 𝑣₂ have been assumed to be equal to the averaged gas velocity 𝑣_𝐺. In the right-hand side of (12), 𝜀_𝐿 is the dissipation rate of liquid turbulence, 𝐶_𝑡 is the ratio of the dispersed phase turbulent fluctuations by the continuous phase turbulent fluctuations, and 𝑓(1,𝜎,𝑃₀₀) and 𝑓(2,𝜎,𝑃₀₀) are numerical functions of the width parameter 𝜎 and the coalescence probability (based on the molecular chaos assumption) 𝑃₀₀ of two bubbles having a diameter 𝑑₀₀ (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 term 𝐺 in (9) was neglected. Here we calculate the corresponding terms 𝐺_𝛾−1 appearing in (10) when 𝐺 is due to the vapor compressibility (Section 4.1) and to phase change (Section 4.2).

4.1. Vapor Compressibility

Let 𝐷(𝑥,𝑡) be the bubble diameter in physical space (the notation 𝑑 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𝐷_𝐺𝐷𝐷𝑡=𝜕𝐷𝜕𝑡+𝑣𝐷⋅∇𝐷=−(3𝜌𝜕𝜌𝜕𝑡+𝑣⋅∇𝜌), (13)where 𝜌 denotes the vapor density and 𝑣 the bubble velocity. The quantity 𝐺 appearing in (9) is the conditional expectation of the Lagrangian derivative given by (13) conditioned by the equality 𝐷(𝑥,𝑡)=𝑑:𝐷𝐺=⟨^𝐺𝐷𝑑𝐷𝑡∣𝐷=𝑑⟩=−3⟨1𝜌(𝜕𝜌𝜕𝑡+𝑣𝑑⋅∇𝜌)∣𝐷=𝑑⟩≅−31𝜌_𝐺(𝜕𝜌_𝐺𝜕𝑡+𝑣_𝐺⋅∇𝜌_𝐺). (14)If we assume that the vapor density 𝜌 and the bubble velocity 𝑣 do not depend on the bubble diameter 𝐷, 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 𝑣 characterizing the bubble are replaced by the macroscopic (averaged) ones 𝜌_𝐺 and 𝑣_𝐺.

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:[𝜕𝑆₁𝜕𝑡+∇⋅(𝑆₁𝑣_𝐺)]_{compressibility}=𝐺₀𝑆₀𝑆=−¹3𝜌_𝐺(𝜕𝜌_𝐺𝜕𝑡+𝑣_𝐺⋅∇𝜌_𝐺[),𝜕𝑆₂𝜕𝑡+∇⋅(𝑆₂𝑣_𝐺)]_{compressibility}=2𝐺₁𝑆₁2=−3𝑆₂𝜌_𝐺(𝜕𝜌_𝐺𝜕𝑡+𝑣_𝐺⋅∇𝜌_𝐺). (15)

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 for 𝑆_1,2. Here we assume a constant vapor density 𝜌=𝜌_𝐺 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 𝑆_1,2 balance equations. For a constant vapor density, the bubble diameter variation in physical space is simply given by𝐷_𝐺𝐷𝐷𝑡=𝜕𝐷𝜕𝑡+𝑣2𝑚̇.∇𝐷=𝜌_𝐺, (16)where 𝑚̇ 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𝐷𝐺=⟨^𝐺𝐷2𝑚̇𝐷𝑡∣𝐷=𝑑⟩=⟨𝜌_𝐺2∣𝐷=𝑑⟩=𝜌_𝐺⟨𝑚̇∣𝐷=𝑑⟩. (17)The closure issue consists to propose an expression for the conditional expectation ⟨𝑚̇∣𝐷=𝑑⟩. The simplest choice is to assume that ⟨𝑚̇∣𝐷=𝑑⟩ is simply given by the ratio (Γ_𝐺−Γ_𝑁)/𝑎_𝐼, where Γ_𝐺 is the vapor source per unit volume per unit time and Γ_𝑁 is the part of this source due to newly nucleated bubbles, therefore Γ_𝐺−Γ_𝑁 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:[𝜕𝑆₁𝜕𝑡+∇⋅(𝑆₁𝑣_𝐺)]_phasechange=𝐺₀𝑆₀=2𝜌_𝐺Γ_𝐺−Γ_𝑁𝑎_𝐼𝑆₀,[𝜕𝑆₂𝜕𝑡+∇⋅(𝑆₂𝑣_𝐺)]_phasechange=2𝐺₁𝑆₁=4𝜌_𝐺Γ_𝐺−Γ_𝑁𝑎_𝐼𝑆₁, (18) where 𝑆₀ is the volumetric number of bubbles 𝑛 obtained by making 𝛾=0 into the relation (8). A more sophisticated model taking into account the dependence of ⟨𝑚̇∣𝐷=𝑑⟩ on the diameter 𝑑 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 diameter 𝑑_𝑑. 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 Γ_𝑁. Dividing Γ_𝑁 by the mass of a nucleated bubble 𝜌_𝐺𝜋𝑑³_𝑑/6, the volumetric number of bubbles source term is obtained. The corresponding source terms in the 𝑆₁ and 𝑆₂ equations are simply obtained by multiplying this volumetric number of bubbles source term by the diameter 𝑑_𝑑 or by its square 𝑑²_𝑑, 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:(𝑞_𝐺−𝑞_𝐿)⋅𝑛_𝐺−𝑚̇(ℎ_𝐺𝐼−ℎ_𝐿𝐼)≅0. (19)In this relation, 𝑞_𝐺⋅𝑛_𝐺 and 𝑞_𝐿⋅𝑛_𝐺 are the heat flux densities on the two sides of the interface, 𝑛_𝐺 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 ℎ_𝐺𝐼 and ℎ_𝐿𝐼 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 𝑚̇ 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 𝑞_𝐺⋅𝑛_𝐺 in comparison to 𝑞_𝐿⋅𝑛_𝐺 in the previous equation and introducing a Nusselt number Nu_𝐿 to characterize the liquid-to-interface heat transfer, (19) becomeṡ𝑞𝑚≅^𝐿⋅𝑛_𝐿ℎ_𝐺𝐼−ℎ_𝐿𝐼=Nu_𝐿𝜆_𝐿(ℎ_𝐺𝐼−ℎ_𝐿𝐼)𝐷(𝜃_𝐿−𝜃_sat), (20)where 𝜆_𝐿 is the liquid thermal conductivity and 𝜃_𝐿 and 𝜃_sat represent the liquid and saturation temperatures, respectively. The Nusselt number is given by the Ranz-Marschall law [12] Nu_𝐿=2+0.55Re^1/2𝑃𝑟_𝐿^1/3, (21)where Re and 𝑃𝑟_𝐿 are the bubble Reynolds number and the liquid Prandtl number, respectively. Neglecting the fluctuations of the enthalpies ℎ_𝐺𝐼 and ℎ_𝐿𝐼, the averaging of the first relation (20) gives𝑞^′′_𝐿𝐼𝑎_𝐼=Γ_𝐺(ℎ_𝐺𝐼−ℎ_𝐿𝐼Γ)with_𝐺=⟨𝑚̇∣𝐷=𝑑⟩𝜋𝑑²𝑞𝑛𝑃(𝑑)𝛿𝑑,_{′′𝐿𝐼}𝑎_𝐼=⟨𝑞_𝐿⋅𝑛_𝐿∣𝐷=𝑑⟩𝜋𝑑²𝑛𝑃(𝑑)𝛿𝑑, (22) where Γ_𝐺 and 𝑞^′′_𝐿𝐼𝑎_𝐼 are the average interfacial mass and heat transfer terms appearing in (1) and (3). The term Γ_𝐺 is related to the vapor phase because 𝑚̇ is the bubble mass gain and 𝑞^′′_𝐿𝐼𝑎_𝐼 is the liquid-to-interface heat transfer per unit volume per unit time, the heat transfer density 𝑞^′′_𝐿𝐼 being expressed per unit area per unit time, 𝑎_𝐼 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𝑞^′′_𝐿𝐼𝑎_𝐼=⟨Nu_𝐿𝜆_𝐿𝐷(𝜃_𝐿−𝜃_sat)∣𝐷=𝑑⟩𝜋𝑑²𝑛𝑃(𝑑)𝛿𝑑. (23)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𝑞^′′_𝐿𝐼𝑎_𝐼≅𝜋𝑑Nu_𝐿𝜆_𝐿(𝜃_𝐿−𝜃_sat)𝑛𝑃(𝑑)𝛿𝑑≅𝑎_𝐼Nu_𝐿(𝑑₃₂)𝑑₃₂𝜆_𝐿(Θ_𝐿−Θ_sat). (24)In this expression, it is implicitly assumed that all the bubbles have the same diameter which is given by the Sauter mean diameter 𝑑₃₂=6𝛼/𝑎_𝐼 (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𝑞^′′_𝐿𝐼𝑎_𝐼=[2𝜋𝜆_𝐿𝑆₁+0.55𝜋𝜆_𝐿𝑃𝑟_𝐿^1/3𝑣_𝐿^1/2|𝑣_𝑟|^1/2𝑛𝑑^3/2₀₀9exp(8𝜎²)]×(Θ_𝐿−Θ_sat). (25)The dependence of the bubble Reynolds number Re on the bubble diameter 𝑑 has been taken into account in the derivation of (25), but not the variation of the relative velocity |𝑣_𝑟|, 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|𝑣_𝑟|≈𝑣_𝑇=43Δ𝜌𝑔𝑑𝜌_𝐿𝐶_𝐷, (26)where Δ𝜌 is the absolute value of the density difference |𝜌_𝐿−𝜌_𝐺| and 𝐶_𝐷 is the drag coefficient, which can be given for spherical bubbles by (e.g., [13])𝐶_𝐷=24Re(1+0.1Re^3/4). (27)Assuming that the bubble Reynolds number is sufficiently high in order to approximate the preceding relation by 𝐶_𝐷=2.4/Re^1/4, we obtain the following relation for the square root of the terminal velocity:𝑣_𝑇^1/2=(Δ𝜌𝑔1.8𝜌_𝑐𝑣_𝑐^1/4)^2/7𝑑^5/14. (28)Reevaluating the liquid-to-interface heat transfer from (4), (21), (23), and (28) gives a third model:𝑞^′′_𝐿𝐼𝑎_𝐼=[2𝜋𝜆_𝐿𝑆₁+0.55𝜋𝜆_𝐿𝑃𝑟_𝐿^1/3𝑣_𝐿^1/2(Δ𝜌𝑔1.8𝜌_𝐿𝑣_𝐿^1/4)^2/7×𝑛𝑑^13/7₀₀exp(169𝜎98²)](Θ_𝐿−Θ_sat). (29)Now, having the model given by (20)-(21) for the bubble mass gain 𝑚̇, we can evaluate more precisely the corresponding terms in the equations for the two moments 𝑆₁ and 𝑆₂, namely the terms 𝐺₀ and 𝐺₁ with 𝐺 given from (17). The results are the following:[𝜕𝑆₁𝜕𝑡+∇⋅(𝑆₁𝑣_𝐺)]_phasechange=4𝜌_𝐺𝜆_𝐿ℎ_𝐺𝐼−ℎ_𝐿𝐼𝑛𝑑₀₀𝜎exp(²2)+1.1𝜌_𝐺𝜆_𝐿ℎ_𝐺𝐼−ℎ_𝐿𝐼𝑃𝑟_𝐿^1/3𝑣_𝐿^1/2(Δ𝜌𝑔1.8𝜌_𝐿𝑣_𝐿^1/4)^2/7𝑛𝑑^1/7₀₀𝜎exp(²)98×(Θ_𝐿−Θ_sat[),𝜕𝑆₂𝜕𝑡+∇⋅(𝑆₂𝑣_𝐺)]_phasechange=8𝜌_𝐺𝜆_𝐿ℎ_𝐺𝐼−ℎ_𝐿𝐼𝑛(Θ_𝐿−Θ_sat)+2.2𝜌_𝐺𝜆_𝐿ℎ_𝐺𝐼−ℎ_𝐿𝐼𝑃𝑟_𝐿^1/3𝑣_𝐿^1/2(Δ𝜌𝑔1.8𝜌_𝐿𝑣_𝐿^1/4)^2/7𝑛𝑑^6/7₀₀exp(18𝜎49²)×(Θ_𝐿−Θ_sat). (30) 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 𝑑₀₀ 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/m²s.(iii) Liquid inlet temperature: 28.52^∘C.(iv) Wall to fluid heat flux density: 76200W/m².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 𝑑₃₂ (Figure 5), and the mean diameter 𝑑₂₀ 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 diameter $d_{20}$ 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 𝑑₂₀. 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).

Acknowledgments

Neptune_CFD is a three-dimensional two-fluid code developed more especially for nuclear reactor applications within the framework of the Neptune project, financially supported by CEA (Commissariat à l'Energie Atomique), EDF (Electricité de France), IRSN (Institut de Radioprotection et de Sûreté Nucléaire), and AREVA-NP. This activity was realized in the context of the NURESIM project financed by the European Union.