Science and Technology of Nuclear Installations
Volume 2009 (2009), Article ID 953527, 8 pages
doi:10.1155/2009/953527
Research Article

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

Christophe Morel1 and Jérôme M. Laviéville2

1Commissariat à l'Energie Atomique, DEN/DER/SSTH/LMDL, 17 Rue des Martyrs, 38054 Grenoble Cedex 9, France
2R&D Division, Electricité de France, 6 Quai Watier, 78400 Chatou, France

Received 19 December 2007; Accepted 9 June 2008

Academic Editor: Fabio Moretti

Copyright © 2009 Christophe Morel and Jérôme M. Laviéville. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 𝑘 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 𝜋 𝜎 𝑑 − { l n ( 𝑑 / 𝑑 0 0 ) } 2 / 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 𝑑 0 0 ( 𝑥 , 𝑡 ) 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: 𝑛 = 𝑆 0 , 𝑑 1 0 𝑆 = 1 𝑛 , 𝑎 𝑖 = 𝜋 𝑆 2 , 𝛼 = 𝜋 𝑆 3 6 , ( 5 ) where 𝑑 1 0 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 𝑑 1 0 is one of them, the Sauter mean diameter 𝑑 3 2 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 𝑆 1 and 𝑆 2 and of the void fraction 𝛼 (see [5]): 𝜎 = l n 6 𝛼 𝑆 1 𝜋 𝑆 2 2 , 𝑑 0 0 = 6 𝛼 𝜋 𝑆 2 𝑒 − 5 𝜎 2 / 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 0 0 𝑒 𝜎 2 / 2 ( 𝛾 2 − 9 ) . ( 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 𝑆 1 and 𝑆 2 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 + 𝜑 𝛾 , ( 1 0 ) which is used for 𝛾 = 1 , 2 . In this equation, 𝑣 𝛾 , 𝐺 𝛾 , and 𝜑 𝛾 are defined by the following relations: 𝑣 𝛾 ∫ = ∞ 0 𝑣 ( 𝑑 ) 𝑑 𝛾 𝑃 ( 𝑑 ) 𝛿 𝑑 ∫ ∞ 0 𝑑 𝛾 , 𝐺 𝑃 ( 𝑑 ) 𝛿 𝑑 𝛾 ∫ = ∞ 0 𝐺 ( 𝑑 ) 𝑑 𝛾 𝑃 ( 𝑑 ) 𝛿 𝑑 ∫ ∞ 0 𝑑 𝛾 , 𝜑 𝑃 ( 𝑑 ) 𝛿 𝑑 𝛾 = ∞ 0 𝜑 ( 𝑑 ) 𝑑 𝛾 𝛿 𝑑 . ( 1 1 ) 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 𝑆 1 and 𝑆 2 read (see [5]) 𝜕 𝑆 1 𝜕 𝑡 + ∇ ⋅ ( 𝑆 1 𝑣 𝐺 ) = 8 𝜋 3 6 𝛼 𝜋 2 / 3 𝐶 𝑡 𝜀 𝐿 1 / 3 √ 1 . 6 1 ( 2 1 / 3 − 2 ) 𝑆 1 4 / 3 𝑓 ( 1 , 𝜎 , 𝑃 0 0 ) , 𝜕 𝑆 2 𝜕 𝑡 + ∇ ⋅ ( 𝑆 2 𝑣 𝐺 ) = 8 𝜋 3 6 𝛼 𝜋 1 / 3 𝐶 𝑡 𝜀 𝐿 1 / 3 √ 1 . 6 1 ( 2 2 / 3 − 2 ) 𝑆 2 5 / 3 𝑓 ( 2 , 𝜎 , 𝑃 0 0 ) , ( 1 2 ) where the two velocities 𝑣 1 and 𝑣 2 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 , 𝜎 , 𝑃 0 0 ) and 𝑓 ( 2 , 𝜎 , 𝑃 0 0 ) are numerical functions of the width parameter 𝜎 and the coalescence probability (based on the molecular chaos assumption) 𝑃 0 0 of two bubbles having a diameter 𝑑 0 0 (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 𝜌 𝜕 𝜌 𝜕 𝑡 + 𝑣 ⋅ ∇ 𝜌 , ( 1 3 ) 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 𝜌 𝜕 𝜌 𝜕 𝑡 + 𝑣 𝑑 ⋅ ∇ 𝜌 ∣ 𝐷 = 𝑑 ≅ − 3 1 𝜌 𝐺 𝜕 𝜌 𝐺 𝜕 𝑡 + 𝑣 𝐺 ⋅ ∇ 𝜌 𝐺 . ( 1 4 ) 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: 𝜕 𝑆 1 𝜕 𝑡 + ∇ ⋅ ( 𝑆 1 𝑣 𝐺 ) c o m p r e s s i b i l i t y = 𝐺 0 𝑆 0 𝑆 = − 1 3 𝜌 𝐺 𝜕 𝜌 𝐺 𝜕 𝑡 + 𝑣 𝐺 ⋅ ∇ 𝜌 𝐺 , 𝜕 𝑆 2 𝜕 𝑡 + ∇ ⋅ ( 𝑆 2 𝑣 𝐺 ) c o m p r e s s i b i l i t y = 2 𝐺 1 𝑆 1 2 = − 3 𝑆 2 𝜌 𝐺 𝜕 𝜌 𝐺 𝜕 𝑡 + 𝑣 𝐺 ⋅ ∇ 𝜌 𝐺 . ( 1 5 )

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 ̇ 𝑚 𝜌 𝐺 , ( 1 6 ) 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 ∣ 𝐷 = 𝑑 𝜌 𝐺 ⟨ ̇ 𝑚 ∣ 𝐷 = 𝑑 ⟩ . ( 1 7 ) 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: 𝜕 𝑆 1 𝜕 𝑡 + ∇ ⋅ ( 𝑆 1 𝑣 𝐺 ) p h a s e c h a n g e = 𝐺 0 𝑆 0 = 2 𝜌 𝐺 Γ 𝐺 − Γ 𝑁 𝑎 𝐼 𝑆 0 , 𝜕 𝑆 2 𝜕 𝑡 + ∇ ⋅ ( 𝑆 2 𝑣 𝐺 ) p h a s e c h a n g e = 2 𝐺 1 𝑆 1 = 4 𝜌 𝐺 Γ 𝐺 − Γ 𝑁 𝑎 𝐼 𝑆 1 , ( 1 8 ) where 𝑆 0 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 𝜌 𝐺 𝜋 𝑑 3 𝑑 / 6 , the volumetric number of bubbles source term is obtained. The corresponding source terms in the 𝑆 1 and 𝑆 2 equations are simply obtained by multiplying this volumetric number of bubbles source term by the diameter 𝑑 𝑑 or by its square 𝑑 2 𝑑 , 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 . ( 1 9 ) 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 N u 𝐿 to characterize the liquid-to-interface heat transfer, (19) becomes 𝑞 ̇ 𝑚 ≅ 𝐿 ⋅ 𝑛 𝐿 ℎ 𝐺 𝐼 − ℎ 𝐿 𝐼 = N u 𝐿 𝜆 𝐿 ( ℎ 𝐺 𝐼 − ℎ 𝐿 𝐼 ) 𝐷 ( 𝜃 𝐿 − 𝜃 s a t ) , ( 2 0 ) where 𝜆 𝐿 is the liquid thermal conductivity and 𝜃 𝐿 and 𝜃 s a t represent the liquid and saturation temperatures, respectively. The Nusselt number is given by the Ranz-Marschall law [12] N u 𝐿 = 2 + 0 . 5 5 R e 1 / 2 𝑃 𝑟 𝐿 1 / 3 , ( 2 1 ) 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 𝑞 𝐿 𝐼 𝑎 𝐼 = Γ 𝐺 ( ℎ 𝐺 𝐼 − ℎ 𝐿 𝐼 Γ ) w i t h 𝐺 = ⟨ ̇ 𝑚 ∣ 𝐷 = 𝑑 ⟩ 𝜋 𝑑 2 𝑞 𝑛 𝑃 ( 𝑑 ) 𝛿 𝑑 , 𝐿 𝐼 𝑎 𝐼 = ⟨ 𝑞 𝐿 ⋅ 𝑛 𝐿 ∣ 𝐷 = 𝑑 ⟩ 𝜋 𝑑 2 𝑛 𝑃 ( 𝑑 ) 𝛿 𝑑 , ( 2 2 ) 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 𝑞 𝐿 𝐼 𝑎 𝐼 = N u 𝐿 𝜆 𝐿 𝐷 ( 𝜃 𝐿 − 𝜃 s a t ) ∣ 𝐷 = 𝑑 𝜋 𝑑 2 𝑛 𝑃 ( 𝑑 ) 𝛿 𝑑 . ( 2 3 ) 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 𝑞 𝐿 𝐼 𝑎 𝐼 ≅ 𝜋 𝑑 N u 𝐿 𝜆 𝐿 ( 𝜃 𝐿 − 𝜃 s a t ) 𝑛 𝑃 ( 𝑑 ) 𝛿 𝑑 ≅ 𝑎 𝐼 N u 𝐿 ( 𝑑 3 2 ) 𝑑 3 2 𝜆 𝐿 ( Θ 𝐿 − Θ s a t ) . ( 2 4 ) In this expression, it is implicitly assumed that all the bubbles have the same diameter which is given by the Sauter mean diameter 𝑑 3 2 = 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 𝜋 𝜆 𝐿 𝑆 1 + 0 . 5 5 𝜋 𝜆 𝐿 𝑃 𝑟 𝐿 1 / 3 𝑣 𝐿 1 / 2 | 𝑣 𝑟 | 1 / 2 𝑛 𝑑 3 / 2 0 0 9 e x p 8 𝜎 2 × ( Θ 𝐿 − Θ s a t ) . ( 2 5 ) 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 | 𝑣 𝑟 | ≈ 𝑣 𝑇 = 4 3 Δ 𝜌 𝑔 𝑑 𝜌 𝐿 𝐶 𝐷 , ( 2 6 ) 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]) 𝐶 𝐷 = 2 4 R e ( 1 + 0 . 1 R e 3 / 4 ) . ( 2 7 ) Assuming that the bubble Reynolds number is sufficiently high in order to approximate the preceding relation by 𝐶 𝐷 = 2 . 4 / R e 1 / 4 , we obtain the following relation for the square root of the terminal velocity: 𝑣 𝑇 1 / 2 = Δ 𝜌 𝑔 1 . 8 𝜌 𝑐 𝑣 𝑐 1 / 4 2 / 7 𝑑 5 / 1 4 . ( 2 8 ) Reevaluating the liquid-to-interface heat transfer from (4), (21), (23), and (28) gives a third model: 𝑞 𝐿 𝐼 𝑎 𝐼 = 2 𝜋 𝜆 𝐿 𝑆 1 + 0 . 5 5 𝜋 𝜆 𝐿 𝑃 𝑟 𝐿 1 / 3 𝑣 𝐿 1 / 2 Δ 𝜌 𝑔 1 . 8 𝜌 𝐿 𝑣 𝐿 1 / 4 2 / 7 × 𝑛 𝑑 1 3 / 7 0 0 e x p 1 6 9 9 8 𝜎 2 ( Θ 𝐿 − Θ s a t ) . ( 2 9 ) 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 𝑆 1 and 𝑆 2 , namely the terms 𝐺 0 and 𝐺 1 with 𝐺 given from (17). The results are the following: 𝜕 𝑆 1 𝜕 𝑡 + ∇ ⋅ ( 𝑆 1 𝑣 𝐺 ) p h a s e c h a n g e = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 4 𝜌 𝐺 𝜆 𝐿 ℎ 𝐺 𝐼 − ℎ 𝐿 𝐼 𝑛 𝑑 0 0 e x p 𝜎 2 2 + 1 . 1 𝜌 𝐺 𝜆 𝐿 ℎ 𝐺 𝐼 − ℎ 𝐿 𝐼 𝑃 𝑟 𝐿 1 / 3 𝑣 𝐿 1 / 2 Δ 𝜌 𝑔 1 . 8 𝜌 𝐿 𝑣 𝐿 1 / 4 2 / 7 𝑛 𝑑 1 / 7 0 0 e x p 𝜎 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎦ 9 8 × ( Θ 𝐿 − Θ s a t ) , 𝜕 𝑆 2 𝜕 𝑡 + ∇ ⋅ ( 𝑆 2 𝑣 𝐺 ) p h a s e c h a n g e = 8 𝜌 𝐺 𝜆 𝐿 ℎ 𝐺 𝐼 − ℎ 𝐿 𝐼 𝑛 ( Θ 𝐿 − Θ s a t ) + 2 . 2 𝜌 𝐺 𝜆 𝐿 ℎ 𝐺 𝐼 − ℎ 𝐿 𝐼 𝑃 𝑟 𝐿 1 / 3 𝑣 𝐿 1 / 2 Δ 𝜌 𝑔 1 . 8 𝜌 𝐿 𝑣 𝐿 1 / 4 2 / 7 𝑛 𝑑 6 / 7 0 0 e x p 1 8 4 9 𝜎 2 × ( Θ 𝐿 − Θ s a t ) . ( 3 0 ) 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 𝑑 0 0 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): 2 0 2 7 k g / m 2 s .(iii) Liquid inlet temperature: 2 8 . 5 2 ∘ C .(iv) Wall to fluid heat flux density: 7 6 2 0 0 W / m 2 .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 𝑑 3 2 (Figure 5), and the mean diameter 𝑑 2 0 defined by (6) (Figure 6). All the radial profiles (experimental and calculated) have been taken at the end of the heated length.

53527.fig.001
Figure 1: comparison of the different models on the void fraction profile.
53527.fig.002
Figure 2: comparison of the different models on the vapor velocity profile.
53527.fig.003
Figure 3: comparison of the different models on the liquid temperature profile.
53527.fig.004
Figure 4: comparison of the different models on the interfacial area profile.
53527.fig.005
Figure 5: comparison of the different models on the Sauter mean diameter profile.
53527.fig.006
Figure 6: comparison of the different models on the mean diameter 𝑑 2 0 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 𝑑 2 0 . 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.

References

  1. C. Colin, X. Riou, and J. Fabre, “Turbulence and shear-induced coalescence in gas-liquid pipe flows,” in Proceedings of the 5th International Conference on Multiphase Flow (ICMF '04), Yokohama, Japan, May-June 2004, paper no. 425.
  2. J. M. Delhaye, Thermohydraulics of Two-Phase Systems for Industrial Design and Nuclear Engineering, McGraw-Hill, New York, NY, USA, 1981.
  3. J. Garnier, E. Manon, and G. Cubizolles, “Local measurements on flow boiling of refrigerant 12 in a vertical tube,” Multiphase Science and Technology, vol. 13, no. 1-2, pp. 1–111, 2001.
  4. M. Ishii, “Two-fluid model for two-phase flow,” Multiphase Science and Technology, vol. 5, no. 1–4, pp. 1–63, 1990.
  5. M. Ishii and T. Hibiki, Thermo-Fluid Dynamics of Two-Phase Flow, Springer, New York, NY, USA, 2006.
  6. A. M. Kamp, Ecoulements turbulents à bulles dans une conduite en micropesanteur [Ph.D. thesis], Institut National Polytechnique de Toulouse, Toulouse, France, 1996.
  7. A. M. Kamp, A. K. Chesters, C. Colin, and J. Fabre, “Bubble coalescence in turbulent flows: a mechanistic model for turbulence-induced coalescene applied to microgravity bubbly pipe flow,” International Journal of Multiphase Flow, vol. 27, no. 8, pp. 1363–1396, 2001. View at Publisher · View at Google Scholar
  8. N. Kurul and M. Z. Podowski, “Multidimensional effects in forced convection subcooled boiling,” in Proceedings of the 9th International Heat Transfer Conference, vol. 1, pp. 21–26, Jerusalem, Israel, August 1990, paper no.-04.
  9. C. Morel, J. Pouvreau, J. M. Laviéville, and M. Boucker, “Numerical simulations of a bubbly flow in a sudden expansion with the NEPTUNE code,” in Proceedings of the 3rd International Symposium on Two-Phase Flow Modeling and Experimentation, Pisa, Italy, September 2004.
  10. C. Morel, S. Mimouni, J. M. Laviéville, and M. Boucker, “R113 boiling bubbly flow in an annular geometry simulated with the NEPTUNE code,” in Proceedings of the 11th International Topical Meeting on Nuclear Reactor Thermalhydraulics (NURETH-11), Avignon, France, October 2005, paper: 248.
  11. W. E. Ranz and W. R. Marschall, “Evaporation from drops,” Chemical Engineering Progress, vol. 48, pp. 173–180, 1952.
  12. H. C. Ünal, “Maximum bubble diameter, maximum bubble-growth time and bubble-growth rate during the subcooled nucleate flow boiling of water up to 17.7 MN/m2,” International Journal of Heat and Mass Transfer, vol. 19, no. 6, pp. 643–649, 1976. View at Publisher · View at Google Scholar
  13. W. Yao and C. Morel, “Volumetric interfacial area prediction in upward bubbly two-phase flow,” International Journal of Heat and Mass Transfer, vol. 47, no. 2, pp. 307–328, 2004. View at Publisher · View at Google Scholar