# CFD for Subcooled Flow Boiling: Parametric Variations

Authors: Roland Rzehak; Eckhard Krepper
Journal: Science and Technology of Nuclear Installations
 (2013)

---

## Abstract

We investigate the present capabilities of CFD for wall boiling. The computational model used combines the Euler/Euler two-phase flow description with heat flux partitioning. Very similar modeling was previously applied to boiling water under high pressure conditions relevant to nuclear power systems. Similar conditions in terms of the relevant nondimensional numbers have been realized in the DEBORA tests using dichlorodifluoromethane (R12) as the working fluid. This facilitated measurements of radial profiles for gas volume fraction, gas velocity, liquid temperature, and bubble size. Robust predictive capabilities of the modeling require that it is validated for a wide range of parameters. It is known that a careful calibration of correlations used in the wall boiling model is necessary to obtain agreement with the measured data. We here consider tests under a variety of conditions concerning liquid subcooling, flow rate, and heat flux. It is investigated to which extent a set of calibrated model parameters suffices to cover at least a certain parameter range.


---

## Body

## 1. Introduction

Subcooled flow boiling occurs in many industrial applications where large heat transfer coefficients are required. However, the efficient heat transfer mechanism provided by vapor generation is limited at a point where liquid is expelled from the surface over a significant area. This occurs at the critical heat flux where the heat transfer coefficient begins to decrease with increasing temperature leading to an unstable situation. In this event a rapid heater temperature excursion occurs which potentially leads to heater melting and destruction. For a given working fluid, the critical heat flux depends on the flow parameters as well as the geometry of the flow domain. The verification of design improvements and their influence on the critical heat flux requires expensive experiments. Therefore, the supplementation or even the replacement of experiments by numerical analyses is of high interest in industrial applications.In the past, many different empirical correlations for critical heat flux were developed and fitted to data obtained from experimental tests. These have been implemented mainly in purpose-specific 1D codes and applied for engineering design calculations. However, these correlations are valid only in the limited region of fluid properties, working conditions, and geometry corresponding to the tests to which they were fitted. Using large look-up tables based on a great number of experiments, a significant range of fluid properties and working conditions can be covered. But this method is still limited to only that specific geometry for which they were developed. Independence of the geometry can only be achieved by the application of CFD methods. Existing CFD models, however, are not yet able to describe critical heat flux reliably. A precondition would be the complete understanding and simulation of boiling as a preliminary state towards critical heat flux.For engineering calculations, currently the most widely used CFD approach to model two-phase flows with significant volume fractions of both phases is the Eulerian two-fluid framework of interpenetrating continua (see, e.g., [1–3]). In this approach, balance equations for mass, momentum, and energy are written for each phase, that is, gas and liquid, separately and weighted by the so-called volume fraction which represents the ensemble averaged probability of occurrence for each phase at a certain point in time and space. Exchange terms between the phases appear as source/sink terms in the balance equations. These exchange terms consist of analytical or empirical correlations, expressing the interfacial forces, as well as heat and mass fluxes, as functions of the average flow parameters. Since most of these correlations are highly problem-specific, their range of validity has to be carefully considered and the entire model has to be validated against experiments.For the case of boiling flows, where heat is transferred into the fluid from a heated wall at such high rates that vapor is generated, additional source terms describing the physics of these processes at the heated wall have to be included. A CFD wall boiling model following the lines of Kurul and Podowski [4, 5] was calibrated and validated by several authors, for example, Krepper et al. [6], against the experimental results of Bartolomej and Chanturiya [7].  In these tests, subcooled flow boiling of water at high pressure flowing upwards in a vertical pipe heated from the outside was investigated and measurements of the axial development of gas volume fraction, wall temperature, and cross-sectionally averaged liquid temperature were provided.An application of continuing high interest is the thermal hydraulic flow in a nuclear reactor. However, typical flow conditions encountered in this application do not particularly lend themselves to experimental investigation. High pressure, high temperature, narrow channels, and small expected sizes of steam bubbles represent significant challenges for measurements. The use of refrigerants can greatly relieve this burden. In the French DEBORA tests [8, 9], dichlorodifluoromethane (R12) was used as the working medium. Advantages are that this allows a choice of test parameters that are more convenient for the measurement compared to the water/steam system at high pressure. The same vapor/liquid density ratio can be achieved at a much lower system pressure and the same Reynolds number can be achieved at a larger diameter of the heated pipe. This enables a measurement of radial profiles for gas volume fraction, temperature, liquid and gas velocities, and bubble sizes which allows a stringent validation of CFD models.The applicability of CFD models to the DEBORA tests was recently studied in Krepper and Rzehak [10, 11] to demonstrate their general validity and identify specific weak points. By coupling a population balance approach to the wall boiling model, it was shown that a reasonable prediction of the measured radial profiles of void fraction, bubble size, liquid velocity, and liquid temperature could be achieved. Based on a review of the theoretical and experimental basis of correlations used in the wall boiling model, a careful assessment of the necessary recalibrations to describe the DEBORA tests was given. Most important are the bubble size at detachment and the nucleation site density. Within a limited range of conditions, different tests could be simulated with a single set of model parameters.In the present work, the previous methodology is applied to a number of further test cases for varied conditions of liquid subcooling, flow rate, and heat flux. In this way, the robustness of the model formulation can be assessed and the (in)dependence of model parameters on the experimental conditions be investigated.The paper is organized as follows. A brief summary of the DEBORA test facility and selected data is given in Section2. In Section 3 the details of the models used for the simulations are described. The presentation is adapted from our previous work [10, 11] to make the paper self-contained. The setup of the simulations together with necessary calibrations of the models is described in Section 4. The results presented in Sections 5 and 6 show that the coupling of wall boiling with the population balance model enables the simulation both of the measured bubble size profiles and the gas volume fraction profiles provided the models are suitably calibrated. A discussion of the achievements so far and needs for further research are finally offered in Section 7.


## 2. The Investigated Experiments

### 2.1. Test Cases and Parameters

A detailed description of the DEBORA test facility can be found in Manon [8] and Garnier et al. [9]. In a vertical heated pipe having an inner diameter of 19.2 mm, dichlorodifluoromethane (R12) is heated over a pipe length of 3.5 m as sketched in Figure 1. The facility is operated at various pressures, flow rates, heat fluxes, and subcoolings. The radial profiles for gas volume fraction and gas velocity at the end of the heated length are measured by means of an optical probe. Furthermore, profiles of bubble size at this position are available. In addition, radial liquid temperature profiles as well as axial profiles of the wall temperature are measured by thermocouples. Unfortunately, temperature- and phase-related measurements are not both available for all experimental conditions.Figure 1
Sketch of the DEBORA test geometry.Table1 lists the parameters of the DEBORA tests selected for the present investigation. The first two cases, P26-G2-Q74-T16 and P26-G2-Q74-T18, have been considered previously, for example, in Yao and Morel [12, 13], Boucker et al. [14], and Krepper and Rzehak [10, 11]. The next three cases extend the range of inlet subcooling at the same flow rate and heat flux. The second series of four cases covers the same range of inlet subcooling but at different values of flow rate and heat flux. In the last group, the heat flux is varied at constant values of inlet subcooling and flow rate. For all investigated tests, the system pressure is kept at a fixed value of P=2.62 MPa. The question is to what extent the model is capable to reproduce the data under all of these conditions with as little recalibration as possible.Table 1
System parameters for the selected test cases,P=2.62 MPa.
Mass flow rate(kg m−2 s−1)
Wall heat flux(kW m−2)
Inlet subcooling(K)
Inlet temperature(°C)
Outlet equilibrium vapor quality
P26-G2-Q74-T16
15.7
70.7
0.0848
P26-G2-Q74-T18
17.9
68.5
0.058
P26-G2-Q74-T20
2000
74
19.8
66.6
0.0324
P26-G2-Q74-T26
25.7
60.8
−0.0205
P26-G2-Q74-T28
27.6
58.3
−0.0719
P26-G3-Q118-T16
16.0
70.7
0.1056
P26-G3-Q118-T20
3000
118
20.2
66.2
0.0483
P26-G3-Q118-T24
24.2
62.2
−0.0027
P26-G3-Q118-T28
27.6
58.3
−0.0523
P26-G3-Q129-T29
129
28.6
57.9
−0.0334
P26-G3-Q139-T28
3000
139
27.8
58.6
0.0055
P26-G3-Q148-T28
148
28.3
58.2
0.0259


### 2.2. Fluid Properties and Dimensionless Groups

To compare results obtained for the DEBORA tests to other typical experiments and applications, where the working medium is water at different pressure levels, the values of the relevant dimensionless groups have to be considered. For boiling phenomena, the pipe Reynolds number, the liquid-to-gas density ratio, the Jakob number, and the Boiling number play a role. Focusing on the bubble dynamics furthermore the bubble Reynolds number, the Eötvös number, and the Morton number are important. Conditions at which the DEBORA tests were performed have been chosen to match values of these dimensionless parameters to those for water at conditions typical for pressurized water reactors. The replacement of water by R12 allows measurements at more convenient pressures and temperatures. Also advantageous is the possible increase of pipe diameter which enables the measurement of radial profiles.For R12 liquid and vapor, the relevant material properties were taken from the National Institute of Standards and Technology (NIST) Standard Reference Database on Thermophysical Properties of Fluid Systems (http://webbook.nist.gov/chemistry/fluid). The saturation properties for the value of pressure considered, that is, P=2.62 MPa, are listed in Table 2. Based on the NIST data, fluid property tables were generated and used in the calculations.Table 2
Material properties of R12 and water atP=2.62 MPa.
Material
R12
Pressure condition (MPa)
2.62
Saturation temperature (K)
360
Surface tension (N m−1)
0.00180
Enthalpy of vaporization (J kg−1)
293 × 103
Vapor density (kg m−3)
0.172 × 103
Liquid density (kg m−3)
1.02 × 103
Liquid specific heat capacity (J kg−1 K−1)
1.42 × 103
Liquid viscosity (kg m−1 s−1)
89.5 × 10−6
Liquid thermal conductivity (W m−1 K−1)
0.0457


## 2.1. Test Cases and Parameters

A detailed description of the DEBORA test facility can be found in Manon [8] and Garnier et al. [9]. In a vertical heated pipe having an inner diameter of 19.2 mm, dichlorodifluoromethane (R12) is heated over a pipe length of 3.5 m as sketched in Figure 1. The facility is operated at various pressures, flow rates, heat fluxes, and subcoolings. The radial profiles for gas volume fraction and gas velocity at the end of the heated length are measured by means of an optical probe. Furthermore, profiles of bubble size at this position are available. In addition, radial liquid temperature profiles as well as axial profiles of the wall temperature are measured by thermocouples. Unfortunately, temperature- and phase-related measurements are not both available for all experimental conditions.Figure 1
Sketch of the DEBORA test geometry.Table1 lists the parameters of the DEBORA tests selected for the present investigation. The first two cases, P26-G2-Q74-T16 and P26-G2-Q74-T18, have been considered previously, for example, in Yao and Morel [12, 13], Boucker et al. [14], and Krepper and Rzehak [10, 11]. The next three cases extend the range of inlet subcooling at the same flow rate and heat flux. The second series of four cases covers the same range of inlet subcooling but at different values of flow rate and heat flux. In the last group, the heat flux is varied at constant values of inlet subcooling and flow rate. For all investigated tests, the system pressure is kept at a fixed value of P=2.62 MPa. The question is to what extent the model is capable to reproduce the data under all of these conditions with as little recalibration as possible.Table 1
System parameters for the selected test cases,P=2.62 MPa.
Mass flow rate(kg m−2 s−1)
Wall heat flux(kW m−2)
Inlet subcooling(K)
Inlet temperature(°C)
Outlet equilibrium vapor quality
P26-G2-Q74-T16
15.7
70.7
0.0848
P26-G2-Q74-T18
17.9
68.5
0.058
P26-G2-Q74-T20
2000
74
19.8
66.6
0.0324
P26-G2-Q74-T26
25.7
60.8
−0.0205
P26-G2-Q74-T28
27.6
58.3
−0.0719
P26-G3-Q118-T16
16.0
70.7
0.1056
P26-G3-Q118-T20
3000
118
20.2
66.2
0.0483
P26-G3-Q118-T24
24.2
62.2
−0.0027
P26-G3-Q118-T28
27.6
58.3
−0.0523
P26-G3-Q129-T29
129
28.6
57.9
−0.0334
P26-G3-Q139-T28
3000
139
27.8
58.6
0.0055
P26-G3-Q148-T28
148
28.3
58.2
0.0259


## 2.2. Fluid Properties and Dimensionless Groups

To compare results obtained for the DEBORA tests to other typical experiments and applications, where the working medium is water at different pressure levels, the values of the relevant dimensionless groups have to be considered. For boiling phenomena, the pipe Reynolds number, the liquid-to-gas density ratio, the Jakob number, and the Boiling number play a role. Focusing on the bubble dynamics furthermore the bubble Reynolds number, the Eötvös number, and the Morton number are important. Conditions at which the DEBORA tests were performed have been chosen to match values of these dimensionless parameters to those for water at conditions typical for pressurized water reactors. The replacement of water by R12 allows measurements at more convenient pressures and temperatures. Also advantageous is the possible increase of pipe diameter which enables the measurement of radial profiles.For R12 liquid and vapor, the relevant material properties were taken from the National Institute of Standards and Technology (NIST) Standard Reference Database on Thermophysical Properties of Fluid Systems (http://webbook.nist.gov/chemistry/fluid). The saturation properties for the value of pressure considered, that is, P=2.62 MPa, are listed in Table 2. Based on the NIST data, fluid property tables were generated and used in the calculations.Table 2
Material properties of R12 and water atP=2.62 MPa.
Material
R12
Pressure condition (MPa)
2.62
Saturation temperature (K)
360
Surface tension (N m−1)
0.00180
Enthalpy of vaporization (J kg−1)
293 × 103
Vapor density (kg m−3)
0.172 × 103
Liquid density (kg m−3)
1.02 × 103
Liquid specific heat capacity (J kg−1 K−1)
1.42 × 103
Liquid viscosity (kg m−1 s−1)
89.5 × 10−6
Liquid thermal conductivity (W m−1 K−1)
0.0457


## 3. The Models

The general equations for diabatic two-phase flow in the Euler/Euler framework of interpenetrating continua have been reviewed in many places before (e.g., [1–3]). Therefore, we here focus on those issues that are particularly relevant for the present investigation as a starting point serves the subcooled boiling model implemented in ANSYS CFX which has been used in a large number of previous studies (e.g., Krepper et al. [6] and references therein).The first major block is the wall boiling model describing vapor generation at the wall and transfer of sensible heat to the liquid. Here, the CFX model closely follows the heat flux partitioning approach of Kurul and Podowski [4, 5]. It has been emphasized in Krepper and Rzehak [10] that some of the commonly used correlations are not universally applicable but have to be recalibrated carefully to the specific conditions under investigation. The basis on which such a recalibration can be carried out for varying conditions is discussed. The final outcome of the calibration will be presented in Section 4.2.A second issue is the modeling of interfacial area that determines the exchange of mass, momentum, and energy between the phases. For the bubbly flow regime, it is convenient to use an equivalent Sauter diameter and work with the bubble size. Obviously, in boiling flows the bubble size may change due to both condensation/evaporation and bubble coalescence/breakup. The importance of taking into account the latter processes has been shown in Krepper et al. [11]. This may be achieved by a population balance model with a source at the wall, the size of the generated bubbles being given by the bubble detachment diameter as calculated from the wall boiling model. Since the rates for bubble coalescence and breakup depend sensitively on the turbulence, bubble-induced contributions to the turbulence have to be included as well. To simplify matters, the vapor bubbles are assumed to be at saturation temperature everywhere which is a rather good approximation except close to the critical heat flux.Turbulent fluctuations are modeled by a shear stress turbulence (SST) model according to Menter [15, 16] applied to the liquid phase. This corresponds to a k-ω model near the walls and a k-ε model far from the walls. The frequently used prescription of Sato et al. [17] for the bubble-induced turbulence has been replaced by including source terms in the turbulence equations following Politano et al. [18]. In addition, a wall function for boiling flows based on analogy to a rough wall is employed, which could be shown in Krepper and Rzehak [10] to give an improved prediction of the velocity profiles.For momentum exchange between the phases, finally, lift and turbulent dispersion forces are included in the model in addition to the ubiquitous drag force. For adiabatic flows, there is in addition also a force that pushes a bubble translating in an otherwise quiescent liquid parallel to a wall in close proximity away from this wall. Different models for this so-called wall force were investigated, but the influence turned out to be small. Therefore, in the present study, this force was neglected. In general the applicability of a wall force to flow boiling should be investigated further.


### 3.1. Modeling of Boiling at a Heated Wall

In boiling, heat is transported from the hot wall to the fluid by several different mechanisms. On parts of the wall, where no bubbles reside, heat flows directly to the subcooled liquid in the same way as that in single-phase flow. On parts of the wall where bubbles grow, heat is consumed by the vapor generation which occurs at the so-called nucleation sites. Moreover, there is a liquid mixing mechanism due to the bubbles which leave the wall. As a consequence of the recirculation around the detaching bubble, a cold liquid from the bulk of the flow is brought into contact with the hot wall which leads to additional cooling. This mechanism, which is obviously not present in single-phase flows, is termed quenching.Accordingly, the given external heat fluxQtot, applied to the heated wall, is written as a sum of three parts:
(1)Qtot=QC+QQ+QE,
where QC, QQ, and QE denote the heat flux components due to single-phase turbulent convection, quenching, and evaporation, respectively. The individual components in this heat flux partitioning are then modeled as functions of the wall temperature and other local flow parameters. Once this is accomplished, (1) can be solved iteratively for the local wall temperature TW, which satisfies the wall heat flux balance. Denoting the fraction of area influenced by the bubbles as AW, the heat flux components are expressed as discussed in the following.The turbulent convection heat flux is calculated in the CFX model version (see Wintterle [19]) in much the same way as for a pure liquid flow without boiling but multiplied by the fraction of area unaffected by the bubbles; that is,
(2)QC=(1-AW)hC(TW-TL).
Here hC is the heat transfer coefficient which is written using the temperature wall function T+(y+) known from Kader [20] as
(3)hC=ρCPuτT+,
where nondimensional variables (indicated by superscript “+”) and the friction velocity uτ are defined as usual. Note that (2)-(3) may be evaluated at any location y, provided it is used consistently whenever a variable depends on position.Q
Q is represented in terms of the quenching heat transfer coefficient hQ:
(4)QQ=AWhQ(TW-TL).
A grid independent solution for QQ is obtained by evaluating the nondimensional temperature profile at a fixed value of y+.The evaporation heat fluxQE is obtained via the evaporation mass flux at the wall:
(5)QE=m˙WHLG,
where HLG is the latent heat of evaporation and the generated vapor mass m˙W is expressed in terms of the bubble diameter at detachment dW, bubble generation frequency f, and nucleation site density N as
(6)m˙W=ρGπ6dW3fN.
The bubble size at detachment is known to depend on a large number of variables, namely, the liquid subcooling, the flow rate, the heat flux, the wall superheat, the fluid properties at the system pressure, and the contact angle between the fluid and the wall. Although several differing proposals have been made (e.g., Ünal [21], Winterton [22], Zeng [23], Kandlikar et al. [24], Klausner et al. [25], Situ et al. [26], Basu et al. [27], and references therein), a generally valid empirical or theoretical correlation covering all of these dependencies is not available at present.Some more readily usable results are possible for restricted conditions, for example, for water at atmospheric pressure. However, for the present application to R12, unfortunately even such more specialized results are not available. Yet more complex is the dependence on the contact angle, which varies significantly not only with the material combination of working fluid and heating surface but also depending strongly on ill-characterized factors like surface roughness and contaminants (e.g., Griffith and Wallis [28] and Hong et al. [29]). For the present application, the value of the contact angle is not available.With respect to the other variables, it may be noted that for each individual test, the values of flow rate and heat flux are constant along the pipe. The available data on axial dependence of the wall superheat (cf., Section5) suggest that this is also the case for the wall superheat except probably close to the beginning of the heated section. Thus further simplification may be sought by fixing these conditions and focusing on the dependence with respect to the liquid subcooling.Common practice in engineering simulations of flow boiling largely relies on an experimental investigation of the bubble size at detachment for water at atmospheric pressure by Tolubinsky and Kostanchuk [30], where the observed dependence on the liquid subcooling can be fitted to a correlation
(7)dW=drefe-(Tsat-TL)/ΔTrefd.
Here again TL is obtained by evaluating the nondimensional temperature profile of Kader [20] at a fixed value of y+. The values of the fit parameters dref and ΔTrefd then must be expected to depend on all the other variables. By suitable calibration of these fit parameters, it has been shown for a number of systems that reasonable agreement with data can be obtained [6, 10].The approach taken in the present investigation is to allow a different calibration for different values of flow rate and heat flux. Eventually it will turn out that within certain limits the same calibration can be used even for different conditions. Detailed values will be given in Section4.2.1.Concerning the nucleation site density, most available correlations are expressed in the form of power laws depending on the wall superheat as(8)N=Nref(TW-TsatΔTrefN)p.
Notably, the only variable that appears is the wall superheat, although a recent compilation [31] shows that vastly different parameter values are required to match different data sets. This may be rationalized based on the microscopic theory of nonclassical heterogeneous nucleation [32]. Accordingly, bubble growth is initiated from preexisting nuclei trapped in surface corrugations of micrometer size as sketched schematically in Figure 2.Figure 2
Sketch of bubble nucleus.An important variable characterizing the thermodynamic stability of such nuclei is the local equilibrium temperature at the curved liquid vapor interface [33] which may be roughly identified with the wall superheat; namely,
(9)TW-Tsat=2RnσTsatρGHLG.
Due to the small spatial scales involved, this process may be expected to be independent of other parameters like liquid subcooling, flow rate, and heat flux. This provides some rationale for the appearance of the wall superheat as the only variable in correlations for the nucleation site density.The diversity of specific expressions for these correlations is likely related to the geometry of corrugations on a given heating surface which depends strongly on the processes that were used to finish the surface. These processes are very diverse and in most boiling experiments not specifically controlled. Moreover, a characterization of the resulting surface topography is rarely available.Parameter variations showed that the assumed nucleation site density has almost no influence on the calculated liquid temperature, a small influence on the calculated gas volume fraction, but a strong influence on the calculated wall superheatTW-Tsat. Due to the common case of missing information on the nucleation site density, the approach taken in the present investigation is to make use of the measured wall temperature to calibrate the correlation equation (8). Since this calibration depends only on properties of the heating surface, it can then be used for all different conditions, as will be corroborated by the results. Detailed values will be given in Section 4.2.1.In terms of bubble detachment diameter and nucleation site density given by (7) and (8), the wall area fraction influenced by vapor bubbles AW is given by
(10)AW=π(adW2)2N.
Here, a is the so-called bubble influence factor, for which a value of 2 is commonly used [4, 5]. Direct experimental evidence concerning this quantity is rather scarce. Probably the most relevant source is Han and Griffith [34] who in some “rough” experiments determined the hydrodynamic disturbance caused by lifting a spherical particle from a horizontal surface and found that it has a range of twice the size of the particle. A similar size has been claimed by Cieslinski et al. [35] from PIV measurement of the flow field around departing bubbles, although the quality of the images presented is rather poor.SinceAW=1 corresponds to the case where the whole surface is under the influence of bubbles, AW as calculated by (10) has to be limited to values smaller than this limit. Moreover, it should be kept in mind that already as AW approaches 1, the assumptions of the model are not really satisfied anymore.The bubble detachment frequencyf is given according to Cole [36] as a function of the detachment size dW as
(11)f=4g(ρL-ρG)3CDdWρL.
Relations of this type have been reviewed critically by Ivey [37] and Ceumern-Lindenstjerna [38].The quenching heat transfer coefficient is calculated using the analytical solution for one-dimensional transient conduction, as suggested by Mikic and Rohsenow [39]:
(12)hQ=2πftwaitkLρLCP,L,
where twait is the waiting time between the departure of a bubble and the appearance of the next one at the same nucleation site. A simple assumption by Kurul and Podowski [4, 5] is adopted here, too, where the waiting time takes a fixed fraction of the bubble departure period:
(13)twait=0.8f.
Data supporting this simple estimate come from the work of del Valle and Kenning [40] which, however, is limited to the regime where the heat flux is 75% of the critical heat flux and larger.


### 3.2. Population Balance Approach to Bubble Size Distribution

To describe polydispersed flows within a purely Eulerian approach, a number of different (multiple) bubble size groupsi=1⋯M is considered, each representing bubbles of typical size di. The fraction of gas volume contained in each bubble size group is denoted as αi so that the total gas volume fraction is given by
(14)αG=∑i=1Mαi.
It is useful to also define occupation numbers fi=αi/αG giving the contribution of each size group to the total gas volume fraction. Obviously we then have ∑ifi=1.For each size group, the equation of mass conservation assumes the form(15)∂∂t(αiρG)+∇(αiρGuG)=Γitopo+Γiphase,
where the right-hand side gives the net source of mass for group i which results from topological changes due to coalescence and breakup as well as phase change due to condensation and evaporation. Models for these two contributions will be discussed in the following subsections.For the homogeneous MUSIG model only one momentum and energy equation for the total amount of vapor is considered as well as the conservation equations of the liquid of course. In these equations, the total gas volume fractionαG is calculated according to (14) above. In addition, also the bubble size dB appears which is taken in the Sauter sense representing the interfacial area AI=6αG/dB. In order to preserve this interpretation, dB is calculated from the occupation number and bubble size for each group as
(16)dB=(∑i=1Mfidi)-1.
The advantage of the homogeneous MUSIG approach is that a large number of bubble size groups can be considered while keeping the computational effort within reasonable bounds. On the other hand, profound effects of bubble size are missed entirely like, for example, the change in the sign of the lift force as discussed in Section 3.4.2. To capture such phenomena, provision has to be made that bubbles of different size may move with different velocities. To overcome these limitations, the inhomogeneous MUSIG model [41, 42] was developed and applied to boiling phenomena in Krepper et al. [11]. For the present investigation, however, it will be seen that this more elaborate treatment is not needed.


#### 3.2.1. Coalescence and Breakup

The net mass source for size groupi due to bubble coalescence and breakup can be expressed as the sum of bubble birth rates due to the breakup of larger bubbles from groups j>i to group i and coalescence of smaller bubbles from size groups j<i to group i as well as bubble death rates due to breakup of bubbles from size group i to smaller bubbles in groups j<i and the coalescence of bubbles from size group i with bubbles from any other group to even larger ones which belong to groups j>i. That is,
(17)Γitopo=Bibreak-Dibreak+Bicoal-Dicoal.
The birth and death rates in turn are commonly expressed in terms of the coalescence and break-up kernels such that
(18)Bibreak=ρGαG∑j>ib(mj,mi)fjmimjΔmi,Dibreak=12ρGαG∑j<ib(mi,mj)fiΔmj,Bicoal=12ρG2αG2∑j≤i∑k≤ic(mj,mk)Xjkifjfkmj+mkmjmk,Dicoal=ρG2αG2∑jc(mi,mj)fifj1mj.
Here Xjki is an approximation to the delta-function δ(mj+mk-mi). For the break-up and coalescence kernel functions b and c, the commonly used break-up models according to Luo and Svendsen [43] and the coalescence models of Prince and Blanch [44] are applied in the present work but were adjusted by factors to match the measured bubble sizes. In this way the applicability of the general framework is demonstrated, but of course further developments will be necessary to improve the physical models and overcome such tuning procedures.


#### 3.2.2. Condensation and Evaporation Including Boiling at the Wall

When condensation or evaporation occurs, the volume fraction in size groupi changes for two reasons: (i) mass is transferred directly between the bubbles and the liquid and (ii) since due to this direct mass transfer, the bubbles are shrinking or growing, they may subsequently belong to a different size group.Written as a source term for size groupi the direct mass transfer to the liquid is given by
(19)Γ~i=-AI,iHLGhL,i(TL-Tsat),
where similar to previous work [10] the assumption has been made that the gas is at saturation temperature. The total source terms for size class i including also the ensuing change of bubble size, that is, Γiphase in (15), have been derived recently by Lucas et al. [45] as
(20)Γiphase={mimi-mi-1Γ~i-mimi+1-miΓ~i+1forΓ~i<0,thatis,condensation,mimi-mi-1Γ~i-1-mimi+1-miΓ~iforΓ~i>0,thatis,evaporation,
where mi=ρπd3/6 is the mass of each bubble in size group i. Basing the calculation on bubble mass rather than size for compressible flows has the advantage that since mass is conserved, no extra terms arise in the equations. Conversion to the corresponding bubble size which depends on the local density can be done straight forwardly as needed. For incompressible flows, no differences between mass- and size-based groups arise.In principleΓ~i should be evaluated with the group size di, but for practical reasons, an approximation is used where a direct transfer term Γ~j is calculated according to (21) only for each velocity group j using the Sauter mean diameter of (16) and the area based fraction thereof is used for Γ~i; that is,
(21)Γ~i=AI,iAI,jΓ~j.
In this way, the size dependency of the factor AI,i in (21) is treated exactly, but that of the factor hI,i is not. The liquid side heat transfer coefficient, finally, is calculated according to Ranz and Marshall [46] as
(22)hL=kLdBNu=kLdB(2+0.6Re1/2Pr1/3).
In addition to the source terms for the continuity equations for the bubble size groups, there is also a mass source for the liquid phase continuity equation which is given by
(23)ΓL=-∑iΓ~i.
Moreover, corresponding secondary sources appear in the momentum and energy equations.A validation of the above procedure against experimental data has been given by Krepper et al. [47].To include the generation of vapor bubbles at the wall, an additional source term,Srpi, is included in (19) following [48]. This source term applies only to the equation corresponding to the size group whose diameter is the closest to the bubble detachment diameter dW. It is given by the evaporation mass flux computed in the wall heat partitioning distributed evenly throughout the grid cells adjacent to the heated wall; that is,
(24)Srpi=m˙WSV,
where m˙W is given by (6) and S and V are wall surface area and volume of the corresponding grid cell, respectively.


### 3.3. Modeling of Turbulence in a Two-Phase Flow

Liquid-phase turbulence is modeled by a shear stress transport (SST) model [15, 16] with additional source terms describing the effects of the bubbles as detailed in Section 3.3.1. The reverse effect exerted by the turbulent fluctuations on the bubbles is taken into account by the turbulent dispersion force as will be described in Section 3.4.3. In addition, a turbulent wall function for boiling flows was employed [10] which is summarized in Section 3.3.2. Turbulence in the gas phase is neglected based on the small density as argued in Troshko and Hassan [49].


#### 3.3.1. Bubble-Induced Turbulence

The effects of the dispersed gas bubbles on the turbulence in the continuous liquid phase are modeled by introducing appropriate source terms in thek- and ε- or ω-equations.The bubble-induced sourceSLk in the k-equation describes additional generation of turbulent kinetic energy due to the bubbles. Assuming that all energy lost by the bubble due to drag is converted into turbulent kinetic energy in the wake of the bubble, this source can be calculated as
(25)SLk=|Fdrag·(uG-uL)|,
where Fdrag is the drag force given by (31) and uG-uL is the relative velocity.The bubble-induced sourceSLε in the ε-equation describes additional dissipation of turbulent kinetic energy due to the bubbles. Following the same dimensional reasoning used in deriving the single phase ε-equation, this term is postulated proportional to the source in the k-equation divided by a suitable time scale τ; that is,
(26)SLε=Cε3SLkτ.
Unfortunately no theoretical derivation is available for the time scale and a number of different expressions have been proposed. For the present simulation, we follow Lee et al. [50], Pfleger and Becker [51], and Politano et al. [18] who assumed that the usual turbulence time scale may be used; that is,
(27)τ=kε.
The coefficient Cε3=1.92 was found to give satisfactory results.An equivalent source term for theω-equation is derived from SLk and SLε as
(28)SLω=1CμkSLε-ωkSLk,
where Cμ=0.09, a standard parameter in the k-ε turbulence model. In the SST model, this term is used independent of the blending function since it should be effective throughout the fluid domain.


#### 3.3.2. Turbulent Wall Function for Boiling Flow

A wall function for boiling flow is obtained by considering that the presence of the bubbles on the wall forces the liquid into a similar flow pattern as that observed in single-phase turbulent flow with wall roughness (e.g., Ramstorfer et al. [52] and Končar and Borut [53]). The latter is described by a modified law of the wall [54, 55]
(29)uuτ=1κln(yδ)+B=1κln(ys)+B*,
where uτ=(τW/ρ) is the friction velocity, δ=μ/(ρuτ) is the viscous length scale, s is the hydrodynamic roughness length scale, and κ≈0.41 is the von Karman constant. For flow over smooth walls, B≡const≈5.5, while for flow over rough walls, B* is a function of s/δ which in the limit of large s, that is, for the so-called fully rough walls, tends to a constant value of 8.5.The hydrodynamic roughnesss may be related directly to the bubble size and nucleation site density as [10]
(30)s=CSNdW3.The constant of proportionality is not known from theoretical considerations at present, but a value ofCS=1.0 has been found to yield good results.As shown in [10], the representation of radial velocity profiles is greatly enhanced by employing this two-phase boiling wall function over the often used single-phase wall function.


### 3.4. Modeling of the Momentum Transfer

For momentum exchange between the phases, the Ishii and Zuber [56] drag law was used. Furthermore, a lift force with sign reversal according to Tomiyama et al. [57] and a Favre averaged turbulence dispersion force [58] were included. As noted in Krepper and Rzehak [10], the effect of an additional wall force is small. Therefore, in the present study this force was neglected. In the following expressions, the forces for the dispersed gaseous phase are given.


#### 3.4.1. Drag Force

The volumetric source of gaseous momentum due to drag exerted by the liquid is given by(31)Fdrag=-34dBCDρLαG|uG-uL|(uG-uL).
The drag coefficient CD is calculated according to Ishii and Zuber [56] considering three different shape regimes as
(32)CD=max(CD,sphere,min(CD,ellipse,CD,cap)),
where
(33)CD,sphere=24Re(1+0.1Re0.75),CD,ellipse=23Eo,CD,cap=83.


#### 3.4.2. Lift Force

A lift force due to interaction of the bubble with the shear field of the liquid was first introduced to two-fluid simulations by Zun [59]. The corresponding momentum source is given by
(34)Flift=-CLρLαG(uG-uL)×rot(uL).
The classical lift force, which has a positive coefficient CL, acts in the direction of decreasing liquid velocity. In case of co-current upwards pipe flow this is the direction towards the pipe wall. Experimental [57] and numerical [60] investigations showed that the direction of the lift force changes its sign, if a substantial deformation of the bubble occurs. From the observation of the trajectories of single bubbles rising in simple shear flow of a glycerol water solution, the following correlation for the lift coefficient was derived:
(35)CL={min[0.288tanh(0.121Re),f(Eo⊥)]Eo⊥<4,f(Eo⊥)for4<Eo⊥for4<10,-0.2710<Eo⊥,withf(Eo⊥)=0.00105Eo⊥3-0.0159Eo⊥2-0.0204Eo⊥+0.474.
This coefficient depends on the modified Eötvös number given by
(36)Eo⊥=g(ρL-ρG)d⊥2σ.
Here d⊥ is the maximum horizontal dimension of the bubble. It is calculated using an empirical correlation for the aspect ratio by the following equation [61]:
(37)d⊥=dB1+0.163Eo0.7573.
For the water-air system at normal conditions, CL changes its sign at dB=5.8 mm which was confirmed by investigations of polydispersed upward vertical air/water bubbly flow [62, 63]. For R12 this value is decreased substantially to about 1.0 mm at 2.62 MPa.


#### 3.4.3. Turbulent Dispersion Force

The turbulent dispersion force is the result of the turbulent fluctuations of liquid velocity. Burns et al. [58] derived an expression by Favre averaging the drag force as
(38)Fdisp=-34CDαLdB(uG-uL)μeffσTD(1αL+1αG)∇αG.
Here, σTD=0.9 is an empirical parameter.As noted in Krepper and Rzehak [10], the effect of an additional wall force is small. Therefore, in the present study, this force was neglected.


## 3.1. Modeling of Boiling at a Heated Wall

In boiling, heat is transported from the hot wall to the fluid by several different mechanisms. On parts of the wall, where no bubbles reside, heat flows directly to the subcooled liquid in the same way as that in single-phase flow. On parts of the wall where bubbles grow, heat is consumed by the vapor generation which occurs at the so-called nucleation sites. Moreover, there is a liquid mixing mechanism due to the bubbles which leave the wall. As a consequence of the recirculation around the detaching bubble, a cold liquid from the bulk of the flow is brought into contact with the hot wall which leads to additional cooling. This mechanism, which is obviously not present in single-phase flows, is termed quenching.Accordingly, the given external heat fluxQtot, applied to the heated wall, is written as a sum of three parts:
(1)Qtot=QC+QQ+QE,
where QC, QQ, and QE denote the heat flux components due to single-phase turbulent convection, quenching, and evaporation, respectively. The individual components in this heat flux partitioning are then modeled as functions of the wall temperature and other local flow parameters. Once this is accomplished, (1) can be solved iteratively for the local wall temperature TW, which satisfies the wall heat flux balance. Denoting the fraction of area influenced by the bubbles as AW, the heat flux components are expressed as discussed in the following.The turbulent convection heat flux is calculated in the CFX model version (see Wintterle [19]) in much the same way as for a pure liquid flow without boiling but multiplied by the fraction of area unaffected by the bubbles; that is,
(2)QC=(1-AW)hC(TW-TL).
Here hC is the heat transfer coefficient which is written using the temperature wall function T+(y+) known from Kader [20] as
(3)hC=ρCPuτT+,
where nondimensional variables (indicated by superscript “+”) and the friction velocity uτ are defined as usual. Note that (2)-(3) may be evaluated at any location y, provided it is used consistently whenever a variable depends on position.Q
Q is represented in terms of the quenching heat transfer coefficient hQ:
(4)QQ=AWhQ(TW-TL).
A grid independent solution for QQ is obtained by evaluating the nondimensional temperature profile at a fixed value of y+.The evaporation heat fluxQE is obtained via the evaporation mass flux at the wall:
(5)QE=m˙WHLG,
where HLG is the latent heat of evaporation and the generated vapor mass m˙W is expressed in terms of the bubble diameter at detachment dW, bubble generation frequency f, and nucleation site density N as
(6)m˙W=ρGπ6dW3fN.
The bubble size at detachment is known to depend on a large number of variables, namely, the liquid subcooling, the flow rate, the heat flux, the wall superheat, the fluid properties at the system pressure, and the contact angle between the fluid and the wall. Although several differing proposals have been made (e.g., Ünal [21], Winterton [22], Zeng [23], Kandlikar et al. [24], Klausner et al. [25], Situ et al. [26], Basu et al. [27], and references therein), a generally valid empirical or theoretical correlation covering all of these dependencies is not available at present.Some more readily usable results are possible for restricted conditions, for example, for water at atmospheric pressure. However, for the present application to R12, unfortunately even such more specialized results are not available. Yet more complex is the dependence on the contact angle, which varies significantly not only with the material combination of working fluid and heating surface but also depending strongly on ill-characterized factors like surface roughness and contaminants (e.g., Griffith and Wallis [28] and Hong et al. [29]). For the present application, the value of the contact angle is not available.With respect to the other variables, it may be noted that for each individual test, the values of flow rate and heat flux are constant along the pipe. The available data on axial dependence of the wall superheat (cf., Section5) suggest that this is also the case for the wall superheat except probably close to the beginning of the heated section. Thus further simplification may be sought by fixing these conditions and focusing on the dependence with respect to the liquid subcooling.Common practice in engineering simulations of flow boiling largely relies on an experimental investigation of the bubble size at detachment for water at atmospheric pressure by Tolubinsky and Kostanchuk [30], where the observed dependence on the liquid subcooling can be fitted to a correlation
(7)dW=drefe-(Tsat-TL)/ΔTrefd.
Here again TL is obtained by evaluating the nondimensional temperature profile of Kader [20] at a fixed value of y+. The values of the fit parameters dref and ΔTrefd then must be expected to depend on all the other variables. By suitable calibration of these fit parameters, it has been shown for a number of systems that reasonable agreement with data can be obtained [6, 10].The approach taken in the present investigation is to allow a different calibration for different values of flow rate and heat flux. Eventually it will turn out that within certain limits the same calibration can be used even for different conditions. Detailed values will be given in Section4.2.1.Concerning the nucleation site density, most available correlations are expressed in the form of power laws depending on the wall superheat as(8)N=Nref(TW-TsatΔTrefN)p.
Notably, the only variable that appears is the wall superheat, although a recent compilation [31] shows that vastly different parameter values are required to match different data sets. This may be rationalized based on the microscopic theory of nonclassical heterogeneous nucleation [32]. Accordingly, bubble growth is initiated from preexisting nuclei trapped in surface corrugations of micrometer size as sketched schematically in Figure 2.Figure 2
Sketch of bubble nucleus.An important variable characterizing the thermodynamic stability of such nuclei is the local equilibrium temperature at the curved liquid vapor interface [33] which may be roughly identified with the wall superheat; namely,
(9)TW-Tsat=2RnσTsatρGHLG.
Due to the small spatial scales involved, this process may be expected to be independent of other parameters like liquid subcooling, flow rate, and heat flux. This provides some rationale for the appearance of the wall superheat as the only variable in correlations for the nucleation site density.The diversity of specific expressions for these correlations is likely related to the geometry of corrugations on a given heating surface which depends strongly on the processes that were used to finish the surface. These processes are very diverse and in most boiling experiments not specifically controlled. Moreover, a characterization of the resulting surface topography is rarely available.Parameter variations showed that the assumed nucleation site density has almost no influence on the calculated liquid temperature, a small influence on the calculated gas volume fraction, but a strong influence on the calculated wall superheatTW-Tsat. Due to the common case of missing information on the nucleation site density, the approach taken in the present investigation is to make use of the measured wall temperature to calibrate the correlation equation (8). Since this calibration depends only on properties of the heating surface, it can then be used for all different conditions, as will be corroborated by the results. Detailed values will be given in Section 4.2.1.In terms of bubble detachment diameter and nucleation site density given by (7) and (8), the wall area fraction influenced by vapor bubbles AW is given by
(10)AW=π(adW2)2N.
Here, a is the so-called bubble influence factor, for which a value of 2 is commonly used [4, 5]. Direct experimental evidence concerning this quantity is rather scarce. Probably the most relevant source is Han and Griffith [34] who in some “rough” experiments determined the hydrodynamic disturbance caused by lifting a spherical particle from a horizontal surface and found that it has a range of twice the size of the particle. A similar size has been claimed by Cieslinski et al. [35] from PIV measurement of the flow field around departing bubbles, although the quality of the images presented is rather poor.SinceAW=1 corresponds to the case where the whole surface is under the influence of bubbles, AW as calculated by (10) has to be limited to values smaller than this limit. Moreover, it should be kept in mind that already as AW approaches 1, the assumptions of the model are not really satisfied anymore.The bubble detachment frequencyf is given according to Cole [36] as a function of the detachment size dW as
(11)f=4g(ρL-ρG)3CDdWρL.
Relations of this type have been reviewed critically by Ivey [37] and Ceumern-Lindenstjerna [38].The quenching heat transfer coefficient is calculated using the analytical solution for one-dimensional transient conduction, as suggested by Mikic and Rohsenow [39]:
(12)hQ=2πftwaitkLρLCP,L,
where twait is the waiting time between the departure of a bubble and the appearance of the next one at the same nucleation site. A simple assumption by Kurul and Podowski [4, 5] is adopted here, too, where the waiting time takes a fixed fraction of the bubble departure period:
(13)twait=0.8f.
Data supporting this simple estimate come from the work of del Valle and Kenning [40] which, however, is limited to the regime where the heat flux is 75% of the critical heat flux and larger.


## 3.2. Population Balance Approach to Bubble Size Distribution

To describe polydispersed flows within a purely Eulerian approach, a number of different (multiple) bubble size groupsi=1⋯M is considered, each representing bubbles of typical size di. The fraction of gas volume contained in each bubble size group is denoted as αi so that the total gas volume fraction is given by
(14)αG=∑i=1Mαi.
It is useful to also define occupation numbers fi=αi/αG giving the contribution of each size group to the total gas volume fraction. Obviously we then have ∑ifi=1.For each size group, the equation of mass conservation assumes the form(15)∂∂t(αiρG)+∇(αiρGuG)=Γitopo+Γiphase,
where the right-hand side gives the net source of mass for group i which results from topological changes due to coalescence and breakup as well as phase change due to condensation and evaporation. Models for these two contributions will be discussed in the following subsections.For the homogeneous MUSIG model only one momentum and energy equation for the total amount of vapor is considered as well as the conservation equations of the liquid of course. In these equations, the total gas volume fractionαG is calculated according to (14) above. In addition, also the bubble size dB appears which is taken in the Sauter sense representing the interfacial area AI=6αG/dB. In order to preserve this interpretation, dB is calculated from the occupation number and bubble size for each group as
(16)dB=(∑i=1Mfidi)-1.
The advantage of the homogeneous MUSIG approach is that a large number of bubble size groups can be considered while keeping the computational effort within reasonable bounds. On the other hand, profound effects of bubble size are missed entirely like, for example, the change in the sign of the lift force as discussed in Section 3.4.2. To capture such phenomena, provision has to be made that bubbles of different size may move with different velocities. To overcome these limitations, the inhomogeneous MUSIG model [41, 42] was developed and applied to boiling phenomena in Krepper et al. [11]. For the present investigation, however, it will be seen that this more elaborate treatment is not needed.


### 3.2.1. Coalescence and Breakup

The net mass source for size groupi due to bubble coalescence and breakup can be expressed as the sum of bubble birth rates due to the breakup of larger bubbles from groups j>i to group i and coalescence of smaller bubbles from size groups j<i to group i as well as bubble death rates due to breakup of bubbles from size group i to smaller bubbles in groups j<i and the coalescence of bubbles from size group i with bubbles from any other group to even larger ones which belong to groups j>i. That is,
(17)Γitopo=Bibreak-Dibreak+Bicoal-Dicoal.
The birth and death rates in turn are commonly expressed in terms of the coalescence and break-up kernels such that
(18)Bibreak=ρGαG∑j>ib(mj,mi)fjmimjΔmi,Dibreak=12ρGαG∑j<ib(mi,mj)fiΔmj,Bicoal=12ρG2αG2∑j≤i∑k≤ic(mj,mk)Xjkifjfkmj+mkmjmk,Dicoal=ρG2αG2∑jc(mi,mj)fifj1mj.
Here Xjki is an approximation to the delta-function δ(mj+mk-mi). For the break-up and coalescence kernel functions b and c, the commonly used break-up models according to Luo and Svendsen [43] and the coalescence models of Prince and Blanch [44] are applied in the present work but were adjusted by factors to match the measured bubble sizes. In this way the applicability of the general framework is demonstrated, but of course further developments will be necessary to improve the physical models and overcome such tuning procedures.


### 3.2.2. Condensation and Evaporation Including Boiling at the Wall

When condensation or evaporation occurs, the volume fraction in size groupi changes for two reasons: (i) mass is transferred directly between the bubbles and the liquid and (ii) since due to this direct mass transfer, the bubbles are shrinking or growing, they may subsequently belong to a different size group.Written as a source term for size groupi the direct mass transfer to the liquid is given by
(19)Γ~i=-AI,iHLGhL,i(TL-Tsat),
where similar to previous work [10] the assumption has been made that the gas is at saturation temperature. The total source terms for size class i including also the ensuing change of bubble size, that is, Γiphase in (15), have been derived recently by Lucas et al. [45] as
(20)Γiphase={mimi-mi-1Γ~i-mimi+1-miΓ~i+1forΓ~i<0,thatis,condensation,mimi-mi-1Γ~i-1-mimi+1-miΓ~iforΓ~i>0,thatis,evaporation,
where mi=ρπd3/6 is the mass of each bubble in size group i. Basing the calculation on bubble mass rather than size for compressible flows has the advantage that since mass is conserved, no extra terms arise in the equations. Conversion to the corresponding bubble size which depends on the local density can be done straight forwardly as needed. For incompressible flows, no differences between mass- and size-based groups arise.In principleΓ~i should be evaluated with the group size di, but for practical reasons, an approximation is used where a direct transfer term Γ~j is calculated according to (21) only for each velocity group j using the Sauter mean diameter of (16) and the area based fraction thereof is used for Γ~i; that is,
(21)Γ~i=AI,iAI,jΓ~j.
In this way, the size dependency of the factor AI,i in (21) is treated exactly, but that of the factor hI,i is not. The liquid side heat transfer coefficient, finally, is calculated according to Ranz and Marshall [46] as
(22)hL=kLdBNu=kLdB(2+0.6Re1/2Pr1/3).
In addition to the source terms for the continuity equations for the bubble size groups, there is also a mass source for the liquid phase continuity equation which is given by
(23)ΓL=-∑iΓ~i.
Moreover, corresponding secondary sources appear in the momentum and energy equations.A validation of the above procedure against experimental data has been given by Krepper et al. [47].To include the generation of vapor bubbles at the wall, an additional source term,Srpi, is included in (19) following [48]. This source term applies only to the equation corresponding to the size group whose diameter is the closest to the bubble detachment diameter dW. It is given by the evaporation mass flux computed in the wall heat partitioning distributed evenly throughout the grid cells adjacent to the heated wall; that is,
(24)Srpi=m˙WSV,
where m˙W is given by (6) and S and V are wall surface area and volume of the corresponding grid cell, respectively.


## 3.2.1. Coalescence and Breakup

The net mass source for size groupi due to bubble coalescence and breakup can be expressed as the sum of bubble birth rates due to the breakup of larger bubbles from groups j>i to group i and coalescence of smaller bubbles from size groups j<i to group i as well as bubble death rates due to breakup of bubbles from size group i to smaller bubbles in groups j<i and the coalescence of bubbles from size group i with bubbles from any other group to even larger ones which belong to groups j>i. That is,
(17)Γitopo=Bibreak-Dibreak+Bicoal-Dicoal.
The birth and death rates in turn are commonly expressed in terms of the coalescence and break-up kernels such that
(18)Bibreak=ρGαG∑j>ib(mj,mi)fjmimjΔmi,Dibreak=12ρGαG∑j<ib(mi,mj)fiΔmj,Bicoal=12ρG2αG2∑j≤i∑k≤ic(mj,mk)Xjkifjfkmj+mkmjmk,Dicoal=ρG2αG2∑jc(mi,mj)fifj1mj.
Here Xjki is an approximation to the delta-function δ(mj+mk-mi). For the break-up and coalescence kernel functions b and c, the commonly used break-up models according to Luo and Svendsen [43] and the coalescence models of Prince and Blanch [44] are applied in the present work but were adjusted by factors to match the measured bubble sizes. In this way the applicability of the general framework is demonstrated, but of course further developments will be necessary to improve the physical models and overcome such tuning procedures.


## 3.2.2. Condensation and Evaporation Including Boiling at the Wall

When condensation or evaporation occurs, the volume fraction in size groupi changes for two reasons: (i) mass is transferred directly between the bubbles and the liquid and (ii) since due to this direct mass transfer, the bubbles are shrinking or growing, they may subsequently belong to a different size group.Written as a source term for size groupi the direct mass transfer to the liquid is given by
(19)Γ~i=-AI,iHLGhL,i(TL-Tsat),
where similar to previous work [10] the assumption has been made that the gas is at saturation temperature. The total source terms for size class i including also the ensuing change of bubble size, that is, Γiphase in (15), have been derived recently by Lucas et al. [45] as
(20)Γiphase={mimi-mi-1Γ~i-mimi+1-miΓ~i+1forΓ~i<0,thatis,condensation,mimi-mi-1Γ~i-1-mimi+1-miΓ~iforΓ~i>0,thatis,evaporation,
where mi=ρπd3/6 is the mass of each bubble in size group i. Basing the calculation on bubble mass rather than size for compressible flows has the advantage that since mass is conserved, no extra terms arise in the equations. Conversion to the corresponding bubble size which depends on the local density can be done straight forwardly as needed. For incompressible flows, no differences between mass- and size-based groups arise.In principleΓ~i should be evaluated with the group size di, but for practical reasons, an approximation is used where a direct transfer term Γ~j is calculated according to (21) only for each velocity group j using the Sauter mean diameter of (16) and the area based fraction thereof is used for Γ~i; that is,
(21)Γ~i=AI,iAI,jΓ~j.
In this way, the size dependency of the factor AI,i in (21) is treated exactly, but that of the factor hI,i is not. The liquid side heat transfer coefficient, finally, is calculated according to Ranz and Marshall [46] as
(22)hL=kLdBNu=kLdB(2+0.6Re1/2Pr1/3).
In addition to the source terms for the continuity equations for the bubble size groups, there is also a mass source for the liquid phase continuity equation which is given by
(23)ΓL=-∑iΓ~i.
Moreover, corresponding secondary sources appear in the momentum and energy equations.A validation of the above procedure against experimental data has been given by Krepper et al. [47].To include the generation of vapor bubbles at the wall, an additional source term,Srpi, is included in (19) following [48]. This source term applies only to the equation corresponding to the size group whose diameter is the closest to the bubble detachment diameter dW. It is given by the evaporation mass flux computed in the wall heat partitioning distributed evenly throughout the grid cells adjacent to the heated wall; that is,
(24)Srpi=m˙WSV,
where m˙W is given by (6) and S and V are wall surface area and volume of the corresponding grid cell, respectively.


## 3.3. Modeling of Turbulence in a Two-Phase Flow

Liquid-phase turbulence is modeled by a shear stress transport (SST) model [15, 16] with additional source terms describing the effects of the bubbles as detailed in Section 3.3.1. The reverse effect exerted by the turbulent fluctuations on the bubbles is taken into account by the turbulent dispersion force as will be described in Section 3.4.3. In addition, a turbulent wall function for boiling flows was employed [10] which is summarized in Section 3.3.2. Turbulence in the gas phase is neglected based on the small density as argued in Troshko and Hassan [49].


### 3.3.1. Bubble-Induced Turbulence

The effects of the dispersed gas bubbles on the turbulence in the continuous liquid phase are modeled by introducing appropriate source terms in thek- and ε- or ω-equations.The bubble-induced sourceSLk in the k-equation describes additional generation of turbulent kinetic energy due to the bubbles. Assuming that all energy lost by the bubble due to drag is converted into turbulent kinetic energy in the wake of the bubble, this source can be calculated as
(25)SLk=|Fdrag·(uG-uL)|,
where Fdrag is the drag force given by (31) and uG-uL is the relative velocity.The bubble-induced sourceSLε in the ε-equation describes additional dissipation of turbulent kinetic energy due to the bubbles. Following the same dimensional reasoning used in deriving the single phase ε-equation, this term is postulated proportional to the source in the k-equation divided by a suitable time scale τ; that is,
(26)SLε=Cε3SLkτ.
Unfortunately no theoretical derivation is available for the time scale and a number of different expressions have been proposed. For the present simulation, we follow Lee et al. [50], Pfleger and Becker [51], and Politano et al. [18] who assumed that the usual turbulence time scale may be used; that is,
(27)τ=kε.
The coefficient Cε3=1.92 was found to give satisfactory results.An equivalent source term for theω-equation is derived from SLk and SLε as
(28)SLω=1CμkSLε-ωkSLk,
where Cμ=0.09, a standard parameter in the k-ε turbulence model. In the SST model, this term is used independent of the blending function since it should be effective throughout the fluid domain.


### 3.3.2. Turbulent Wall Function for Boiling Flow

A wall function for boiling flow is obtained by considering that the presence of the bubbles on the wall forces the liquid into a similar flow pattern as that observed in single-phase turbulent flow with wall roughness (e.g., Ramstorfer et al. [52] and Končar and Borut [53]). The latter is described by a modified law of the wall [54, 55]
(29)uuτ=1κln(yδ)+B=1κln(ys)+B*,
where uτ=(τW/ρ) is the friction velocity, δ=μ/(ρuτ) is the viscous length scale, s is the hydrodynamic roughness length scale, and κ≈0.41 is the von Karman constant. For flow over smooth walls, B≡const≈5.5, while for flow over rough walls, B* is a function of s/δ which in the limit of large s, that is, for the so-called fully rough walls, tends to a constant value of 8.5.The hydrodynamic roughnesss may be related directly to the bubble size and nucleation site density as [10]
(30)s=CSNdW3.The constant of proportionality is not known from theoretical considerations at present, but a value ofCS=1.0 has been found to yield good results.As shown in [10], the representation of radial velocity profiles is greatly enhanced by employing this two-phase boiling wall function over the often used single-phase wall function.


## 3.3.1. Bubble-Induced Turbulence

The effects of the dispersed gas bubbles on the turbulence in the continuous liquid phase are modeled by introducing appropriate source terms in thek- and ε- or ω-equations.The bubble-induced sourceSLk in the k-equation describes additional generation of turbulent kinetic energy due to the bubbles. Assuming that all energy lost by the bubble due to drag is converted into turbulent kinetic energy in the wake of the bubble, this source can be calculated as
(25)SLk=|Fdrag·(uG-uL)|,
where Fdrag is the drag force given by (31) and uG-uL is the relative velocity.The bubble-induced sourceSLε in the ε-equation describes additional dissipation of turbulent kinetic energy due to the bubbles. Following the same dimensional reasoning used in deriving the single phase ε-equation, this term is postulated proportional to the source in the k-equation divided by a suitable time scale τ; that is,
(26)SLε=Cε3SLkτ.
Unfortunately no theoretical derivation is available for the time scale and a number of different expressions have been proposed. For the present simulation, we follow Lee et al. [50], Pfleger and Becker [51], and Politano et al. [18] who assumed that the usual turbulence time scale may be used; that is,
(27)τ=kε.
The coefficient Cε3=1.92 was found to give satisfactory results.An equivalent source term for theω-equation is derived from SLk and SLε as
(28)SLω=1CμkSLε-ωkSLk,
where Cμ=0.09, a standard parameter in the k-ε turbulence model. In the SST model, this term is used independent of the blending function since it should be effective throughout the fluid domain.


## 3.3.2. Turbulent Wall Function for Boiling Flow

A wall function for boiling flow is obtained by considering that the presence of the bubbles on the wall forces the liquid into a similar flow pattern as that observed in single-phase turbulent flow with wall roughness (e.g., Ramstorfer et al. [52] and Končar and Borut [53]). The latter is described by a modified law of the wall [54, 55]
(29)uuτ=1κln(yδ)+B=1κln(ys)+B*,
where uτ=(τW/ρ) is the friction velocity, δ=μ/(ρuτ) is the viscous length scale, s is the hydrodynamic roughness length scale, and κ≈0.41 is the von Karman constant. For flow over smooth walls, B≡const≈5.5, while for flow over rough walls, B* is a function of s/δ which in the limit of large s, that is, for the so-called fully rough walls, tends to a constant value of 8.5.The hydrodynamic roughnesss may be related directly to the bubble size and nucleation site density as [10]
(30)s=CSNdW3.The constant of proportionality is not known from theoretical considerations at present, but a value ofCS=1.0 has been found to yield good results.As shown in [10], the representation of radial velocity profiles is greatly enhanced by employing this two-phase boiling wall function over the often used single-phase wall function.


## 3.4. Modeling of the Momentum Transfer

For momentum exchange between the phases, the Ishii and Zuber [56] drag law was used. Furthermore, a lift force with sign reversal according to Tomiyama et al. [57] and a Favre averaged turbulence dispersion force [58] were included. As noted in Krepper and Rzehak [10], the effect of an additional wall force is small. Therefore, in the present study this force was neglected. In the following expressions, the forces for the dispersed gaseous phase are given.


### 3.4.1. Drag Force

The volumetric source of gaseous momentum due to drag exerted by the liquid is given by(31)Fdrag=-34dBCDρLαG|uG-uL|(uG-uL).
The drag coefficient CD is calculated according to Ishii and Zuber [56] considering three different shape regimes as
(32)CD=max(CD,sphere,min(CD,ellipse,CD,cap)),
where
(33)CD,sphere=24Re(1+0.1Re0.75),CD,ellipse=23Eo,CD,cap=83.


### 3.4.2. Lift Force

A lift force due to interaction of the bubble with the shear field of the liquid was first introduced to two-fluid simulations by Zun [59]. The corresponding momentum source is given by
(34)Flift=-CLρLαG(uG-uL)×rot(uL).
The classical lift force, which has a positive coefficient CL, acts in the direction of decreasing liquid velocity. In case of co-current upwards pipe flow this is the direction towards the pipe wall. Experimental [57] and numerical [60] investigations showed that the direction of the lift force changes its sign, if a substantial deformation of the bubble occurs. From the observation of the trajectories of single bubbles rising in simple shear flow of a glycerol water solution, the following correlation for the lift coefficient was derived:
(35)CL={min[0.288tanh(0.121Re),f(Eo⊥)]Eo⊥<4,f(Eo⊥)for4<Eo⊥for4<10,-0.2710<Eo⊥,withf(Eo⊥)=0.00105Eo⊥3-0.0159Eo⊥2-0.0204Eo⊥+0.474.
This coefficient depends on the modified Eötvös number given by
(36)Eo⊥=g(ρL-ρG)d⊥2σ.
Here d⊥ is the maximum horizontal dimension of the bubble. It is calculated using an empirical correlation for the aspect ratio by the following equation [61]:
(37)d⊥=dB1+0.163Eo0.7573.
For the water-air system at normal conditions, CL changes its sign at dB=5.8 mm which was confirmed by investigations of polydispersed upward vertical air/water bubbly flow [62, 63]. For R12 this value is decreased substantially to about 1.0 mm at 2.62 MPa.


### 3.4.3. Turbulent Dispersion Force

The turbulent dispersion force is the result of the turbulent fluctuations of liquid velocity. Burns et al. [58] derived an expression by Favre averaging the drag force as
(38)Fdisp=-34CDαLdB(uG-uL)μeffσTD(1αL+1αG)∇αG.
Here, σTD=0.9 is an empirical parameter.As noted in Krepper and Rzehak [10], the effect of an additional wall force is small. Therefore, in the present study, this force was neglected.


## 3.4.1. Drag Force

The volumetric source of gaseous momentum due to drag exerted by the liquid is given by(31)Fdrag=-34dBCDρLαG|uG-uL|(uG-uL).
The drag coefficient CD is calculated according to Ishii and Zuber [56] considering three different shape regimes as
(32)CD=max(CD,sphere,min(CD,ellipse,CD,cap)),
where
(33)CD,sphere=24Re(1+0.1Re0.75),CD,ellipse=23Eo,CD,cap=83.


## 3.4.2. Lift Force

A lift force due to interaction of the bubble with the shear field of the liquid was first introduced to two-fluid simulations by Zun [59]. The corresponding momentum source is given by
(34)Flift=-CLρLαG(uG-uL)×rot(uL).
The classical lift force, which has a positive coefficient CL, acts in the direction of decreasing liquid velocity. In case of co-current upwards pipe flow this is the direction towards the pipe wall. Experimental [57] and numerical [60] investigations showed that the direction of the lift force changes its sign, if a substantial deformation of the bubble occurs. From the observation of the trajectories of single bubbles rising in simple shear flow of a glycerol water solution, the following correlation for the lift coefficient was derived:
(35)CL={min[0.288tanh(0.121Re),f(Eo⊥)]Eo⊥<4,f(Eo⊥)for4<Eo⊥for4<10,-0.2710<Eo⊥,withf(Eo⊥)=0.00105Eo⊥3-0.0159Eo⊥2-0.0204Eo⊥+0.474.
This coefficient depends on the modified Eötvös number given by
(36)Eo⊥=g(ρL-ρG)d⊥2σ.
Here d⊥ is the maximum horizontal dimension of the bubble. It is calculated using an empirical correlation for the aspect ratio by the following equation [61]:
(37)d⊥=dB1+0.163Eo0.7573.
For the water-air system at normal conditions, CL changes its sign at dB=5.8 mm which was confirmed by investigations of polydispersed upward vertical air/water bubbly flow [62, 63]. For R12 this value is decreased substantially to about 1.0 mm at 2.62 MPa.


## 3.4.3. Turbulent Dispersion Force

The turbulent dispersion force is the result of the turbulent fluctuations of liquid velocity. Burns et al. [58] derived an expression by Favre averaging the drag force as
(38)Fdisp=-34CDαLdB(uG-uL)μeffσTD(1αL+1αG)∇αG.
Here, σTD=0.9 is an empirical parameter.As noted in Krepper and Rzehak [10], the effect of an additional wall force is small. Therefore, in the present study, this force was neglected.


## 4. Model Setup

The models described in the previous section present a rather general framework that can be specialized to the simulation of any subcooled flow boiling problem. The necessary specifications to simulate the DEBORA test cases will now be described. These comprise prescription of flow domain and boundary conditions and specification of grid and bubble classes as well as calibration of model parameters.The simulations are performed by ANSYS CFX 13. For details of the numerical procedures, we refer to the user guide [64].


### 4.1. Geometry and Boundary Conditions

The tests were simulated in a quasi-2D cylindrical geometry, that is, a narrow cylindrical sector with symmetry boundary conditions imposed on the side faces. The validity of this simplification has been verified by grid resolution studies and by comparison to a 3D simulation representing a 60° sector of the pipe.At the upstream end of the pipe, an unheated flow development zone has been added to obtain at the beginning of the heated section a fully developed turbulent flow that is independent of the detailed conditions imposed at the inlet. The required length of this flow development zone was determined by examining the development of velocity and temperature as well as turbulent kinetic energy and dissipation rate to ensure that these did not vary anymore in the axial direction before entering the heated section. At the outlet at the top, a pressure boundary condition was imposed.On the heated walls, boundary conditions for mass and energy equations are provided by the heat flux partitioning discussed in Section3.1. It remains to specify boundary conditions for the gas and liquid momentum equations. Since the gas volume fraction in our model represents bubbles that have detached from the wall, an appropriate boundary condition for the gas phase is the free slip condition. For the liquid phase, we argue that bubbles which have not left the wall are still attached to their respective nucleation site. Hence they restrain the liquid motion in the same way as the solid wall itself does. Therefore, we choose a no-slip condition for the liquid phase. While this issue does not appear to have received due attention in the literature, the results to be presented justify our choice as preliminary working solution until a more thorough investigation becomes available.All of these two-phase flow simulations have been carried out on a quite coarse grid for which the center of the grid cells adjacent to the wall has a nondimensional coordinate ofy+≈200. For the test P26-G2-Q74-T16, a grid refinement study was performed which showed no change of the results until this value of y+ has decreased to about 70. For still smaller values, no convergence could be achieved. This is a well-known problem of the Kurul and Podowski [4, 5] wall boiling model where all vapor generation occurs in the grid cell adjacent to the wall.


### 4.2.  Calibration of Parameters

#### 4.2.1. Wall Boiling Model

The necessity to recalibrate the boiling model parameters to the working fluid and heating surface of the experiment was discussed in a previous study [10].For the model of bubble detachment size in (7), we proceed as follows. The outermost measurement points of the experimental bubble size profiles are taken as representative of the detachment size. The liquid subcooling at the measurement location is determined from the averaged measured liquid temperature profile for the the test series P26-G2-Q74-Txx. For the two other series P26-G3-Q118-Txx and P26-G3-Qxx-T28, unfortunately no temperature measurements are available. Hence, the averaged liquid subcooling was determined from preliminary CFD calculations. The values of detachment size and liquid subcooling for all tests are listed in Table 3. For each of the three series of tests, values for the reference bubble size, dref, and the reference subcooling, ΔTrefd, were then found by fitting the expression (7) to the values of Table 3. The resulting correlations are shown in Figure 3 and the derived model parameters are collected in Table 4. It may be seen that the calibration parameters depend on the flow rate but not on the heat flux.Table 3
Parameters used for the bubble detachment model.
T
AV (K)
T
sub (K)
d
W
(
R
) (mm)
P26-G2-Q74-T16
359.5
0.1
0.247
P26-G2-Q74-T18
359.2
0.35
0.237
P26-G2-Q74-T20
359.0
0.62
0.228
P26-G2-Q74-T26
356.4
3.23
0.165
P26-G2-Q74-T28
351.6
8.03
0.154
P26-G3-Q118-T16
359.3
0.3
0.152
P26-G3-Q118-T20
358.5
1.1
0.152
P26-G3-Q118-T24
355.9
3.7
0.134
P26-G3-Q118-T28
352.5
7.1
0.124
P26-G3-Q129-T29
353.8
5.8
0.142
P26-G3-Q139-T28
355.7
3.9
0.163
P26-G3-Q148-T28
356.4
3.3
0.161Table 4
Adapted parameters for different test series.
Tests
T
sub
Equation (7) and Figure 3
Equation (8) and Figure 4
F
B
F
C
d
ref
Δ
T
ref
d
N
ref
Δ
T
ref
N
P26-G2-Q74-Txx
28, 26, 20, 18, 16
0.24
45
3.0
e
+
7
10
0.03125
0.5
28,
0.015
1.0
P26-G3-Q118-Txx
24,
0.155
31
3.0
e
+
7
10
0.015
1.0
20,
0.03125
1.0
16
0.03125
0.5
P26-G3-Q129-T28
0.0075
1.0
Q139,
28
0.155
31
3.0
e
+
7
10
0.0075
1.0
Q148
0.0
1.0Bubble size at detachment according to (7) fitted to measurements.
(a)
P26-G2-Q74-T16 to T28dref=0.24 mm, ΔT=45K
(b)
P26-G3-Q118-T16 to T28dref=0.155 mm, ΔT=31K
(c)
P26-G3-Q129-T29 to Q148dref=0.155 mm, ΔT=31KAs discussed in Section3.1 calibration of the model for nucleation site density in (8) is possible by matching the measured axial profiles of wall superheat. Such profiles are available only for the test series P26-G2-Q74-Txx, but as discussed in Section 3.1, the nucleation site density may be expected to depend on parameters subcooling, flow rate, and heat flux only indirectly through the wall superheat. Moreover, the sensitivity of the final results with respect to the exponent p=1.85 is rather low; so only the prefactor Nref will be adjusted for a value of the reference temperature ΔTrefN=10K. Wall temperature profiles for two selected tests, P26-G2-Q74-T18 and P26-G2-Q74-T26, are shown in Figure 4. A value of Nref=3.0e7 m−2 gives good agreement with the respective data.Figure 4
Wall superheat depending on the reference nucleation site densityNref.


#### 4.2.2. Coalescence and Break-Up Models

In the present work, bubble coalescence and breakup are described by the models proposed by Prince and Blanch [44] and by Luo and Svendsen [43]. To obtain agreement with the measurements, efficiency factors FC and FB were introduced and calibrated to match the measured radial bubble size profiles. Bubble size profiles for selected tests and several values of the calibration parameters FC and FB are shown in Figure 5. The procedure was successful so far as a good match could be obtained in all cases, but a different set of calibration parameters was necessary to this end (see Table 4). Therefore, this procedure has to be considered only as a first step and further model development needs to be done.Figure 5
Radial profiles of bubble size for different coefficientsFB and FC.As discussed in Section3.4.2, the lift force changes its sign at a certain bubble size. The value of this critical bubble size can be calculated as 1.0 mm for the material properties of R12 at a pressure of 2.62 MPa. Examination of the measured bubble sizes shows that these are lower than the critical value for all of the tests investigated here. Therefore, it is possible to use the simpler and computationally less expensive homogeneous MUSIG approach. For each test considered here, an equidistant bubble size discretization employing 15 size groups in the range of 0 to 1.5 mm was applied.


## 4.1. Geometry and Boundary Conditions

The tests were simulated in a quasi-2D cylindrical geometry, that is, a narrow cylindrical sector with symmetry boundary conditions imposed on the side faces. The validity of this simplification has been verified by grid resolution studies and by comparison to a 3D simulation representing a 60° sector of the pipe.At the upstream end of the pipe, an unheated flow development zone has been added to obtain at the beginning of the heated section a fully developed turbulent flow that is independent of the detailed conditions imposed at the inlet. The required length of this flow development zone was determined by examining the development of velocity and temperature as well as turbulent kinetic energy and dissipation rate to ensure that these did not vary anymore in the axial direction before entering the heated section. At the outlet at the top, a pressure boundary condition was imposed.On the heated walls, boundary conditions for mass and energy equations are provided by the heat flux partitioning discussed in Section3.1. It remains to specify boundary conditions for the gas and liquid momentum equations. Since the gas volume fraction in our model represents bubbles that have detached from the wall, an appropriate boundary condition for the gas phase is the free slip condition. For the liquid phase, we argue that bubbles which have not left the wall are still attached to their respective nucleation site. Hence they restrain the liquid motion in the same way as the solid wall itself does. Therefore, we choose a no-slip condition for the liquid phase. While this issue does not appear to have received due attention in the literature, the results to be presented justify our choice as preliminary working solution until a more thorough investigation becomes available.All of these two-phase flow simulations have been carried out on a quite coarse grid for which the center of the grid cells adjacent to the wall has a nondimensional coordinate ofy+≈200. For the test P26-G2-Q74-T16, a grid refinement study was performed which showed no change of the results until this value of y+ has decreased to about 70. For still smaller values, no convergence could be achieved. This is a well-known problem of the Kurul and Podowski [4, 5] wall boiling model where all vapor generation occurs in the grid cell adjacent to the wall.


## 4.2.  Calibration of Parameters

### 4.2.1. Wall Boiling Model

The necessity to recalibrate the boiling model parameters to the working fluid and heating surface of the experiment was discussed in a previous study [10].For the model of bubble detachment size in (7), we proceed as follows. The outermost measurement points of the experimental bubble size profiles are taken as representative of the detachment size. The liquid subcooling at the measurement location is determined from the averaged measured liquid temperature profile for the the test series P26-G2-Q74-Txx. For the two other series P26-G3-Q118-Txx and P26-G3-Qxx-T28, unfortunately no temperature measurements are available. Hence, the averaged liquid subcooling was determined from preliminary CFD calculations. The values of detachment size and liquid subcooling for all tests are listed in Table 3. For each of the three series of tests, values for the reference bubble size, dref, and the reference subcooling, ΔTrefd, were then found by fitting the expression (7) to the values of Table 3. The resulting correlations are shown in Figure 3 and the derived model parameters are collected in Table 4. It may be seen that the calibration parameters depend on the flow rate but not on the heat flux.Table 3
Parameters used for the bubble detachment model.
T
AV (K)
T
sub (K)
d
W
(
R
) (mm)
P26-G2-Q74-T16
359.5
0.1
0.247
P26-G2-Q74-T18
359.2
0.35
0.237
P26-G2-Q74-T20
359.0
0.62
0.228
P26-G2-Q74-T26
356.4
3.23
0.165
P26-G2-Q74-T28
351.6
8.03
0.154
P26-G3-Q118-T16
359.3
0.3
0.152
P26-G3-Q118-T20
358.5
1.1
0.152
P26-G3-Q118-T24
355.9
3.7
0.134
P26-G3-Q118-T28
352.5
7.1
0.124
P26-G3-Q129-T29
353.8
5.8
0.142
P26-G3-Q139-T28
355.7
3.9
0.163
P26-G3-Q148-T28
356.4
3.3
0.161Table 4
Adapted parameters for different test series.
Tests
T
sub
Equation (7) and Figure 3
Equation (8) and Figure 4
F
B
F
C
d
ref
Δ
T
ref
d
N
ref
Δ
T
ref
N
P26-G2-Q74-Txx
28, 26, 20, 18, 16
0.24
45
3.0
e
+
7
10
0.03125
0.5
28,
0.015
1.0
P26-G3-Q118-Txx
24,
0.155
31
3.0
e
+
7
10
0.015
1.0
20,
0.03125
1.0
16
0.03125
0.5
P26-G3-Q129-T28
0.0075
1.0
Q139,
28
0.155
31
3.0
e
+
7
10
0.0075
1.0
Q148
0.0
1.0Bubble size at detachment according to (7) fitted to measurements.
(a)
P26-G2-Q74-T16 to T28dref=0.24 mm, ΔT=45K
(b)
P26-G3-Q118-T16 to T28dref=0.155 mm, ΔT=31K
(c)
P26-G3-Q129-T29 to Q148dref=0.155 mm, ΔT=31KAs discussed in Section3.1 calibration of the model for nucleation site density in (8) is possible by matching the measured axial profiles of wall superheat. Such profiles are available only for the test series P26-G2-Q74-Txx, but as discussed in Section 3.1, the nucleation site density may be expected to depend on parameters subcooling, flow rate, and heat flux only indirectly through the wall superheat. Moreover, the sensitivity of the final results with respect to the exponent p=1.85 is rather low; so only the prefactor Nref will be adjusted for a value of the reference temperature ΔTrefN=10K. Wall temperature profiles for two selected tests, P26-G2-Q74-T18 and P26-G2-Q74-T26, are shown in Figure 4. A value of Nref=3.0e7 m−2 gives good agreement with the respective data.Figure 4
Wall superheat depending on the reference nucleation site densityNref.


### 4.2.2. Coalescence and Break-Up Models

In the present work, bubble coalescence and breakup are described by the models proposed by Prince and Blanch [44] and by Luo and Svendsen [43]. To obtain agreement with the measurements, efficiency factors FC and FB were introduced and calibrated to match the measured radial bubble size profiles. Bubble size profiles for selected tests and several values of the calibration parameters FC and FB are shown in Figure 5. The procedure was successful so far as a good match could be obtained in all cases, but a different set of calibration parameters was necessary to this end (see Table 4). Therefore, this procedure has to be considered only as a first step and further model development needs to be done.Figure 5
Radial profiles of bubble size for different coefficientsFB and FC.As discussed in Section3.4.2, the lift force changes its sign at a certain bubble size. The value of this critical bubble size can be calculated as 1.0 mm for the material properties of R12 at a pressure of 2.62 MPa. Examination of the measured bubble sizes shows that these are lower than the critical value for all of the tests investigated here. Therefore, it is possible to use the simpler and computationally less expensive homogeneous MUSIG approach. For each test considered here, an equidistant bubble size discretization employing 15 size groups in the range of 0 to 1.5 mm was applied.


## 4.2.1. Wall Boiling Model

The necessity to recalibrate the boiling model parameters to the working fluid and heating surface of the experiment was discussed in a previous study [10].For the model of bubble detachment size in (7), we proceed as follows. The outermost measurement points of the experimental bubble size profiles are taken as representative of the detachment size. The liquid subcooling at the measurement location is determined from the averaged measured liquid temperature profile for the the test series P26-G2-Q74-Txx. For the two other series P26-G3-Q118-Txx and P26-G3-Qxx-T28, unfortunately no temperature measurements are available. Hence, the averaged liquid subcooling was determined from preliminary CFD calculations. The values of detachment size and liquid subcooling for all tests are listed in Table 3. For each of the three series of tests, values for the reference bubble size, dref, and the reference subcooling, ΔTrefd, were then found by fitting the expression (7) to the values of Table 3. The resulting correlations are shown in Figure 3 and the derived model parameters are collected in Table 4. It may be seen that the calibration parameters depend on the flow rate but not on the heat flux.Table 3
Parameters used for the bubble detachment model.
T
AV (K)
T
sub (K)
d
W
(
R
) (mm)
P26-G2-Q74-T16
359.5
0.1
0.247
P26-G2-Q74-T18
359.2
0.35
0.237
P26-G2-Q74-T20
359.0
0.62
0.228
P26-G2-Q74-T26
356.4
3.23
0.165
P26-G2-Q74-T28
351.6
8.03
0.154
P26-G3-Q118-T16
359.3
0.3
0.152
P26-G3-Q118-T20
358.5
1.1
0.152
P26-G3-Q118-T24
355.9
3.7
0.134
P26-G3-Q118-T28
352.5
7.1
0.124
P26-G3-Q129-T29
353.8
5.8
0.142
P26-G3-Q139-T28
355.7
3.9
0.163
P26-G3-Q148-T28
356.4
3.3
0.161Table 4
Adapted parameters for different test series.
Tests
T
sub
Equation (7) and Figure 3
Equation (8) and Figure 4
F
B
F
C
d
ref
Δ
T
ref
d
N
ref
Δ
T
ref
N
P26-G2-Q74-Txx
28, 26, 20, 18, 16
0.24
45
3.0
e
+
7
10
0.03125
0.5
28,
0.015
1.0
P26-G3-Q118-Txx
24,
0.155
31
3.0
e
+
7
10
0.015
1.0
20,
0.03125
1.0
16
0.03125
0.5
P26-G3-Q129-T28
0.0075
1.0
Q139,
28
0.155
31
3.0
e
+
7
10
0.0075
1.0
Q148
0.0
1.0Bubble size at detachment according to (7) fitted to measurements.
(a)
P26-G2-Q74-T16 to T28dref=0.24 mm, ΔT=45K
(b)
P26-G3-Q118-T16 to T28dref=0.155 mm, ΔT=31K
(c)
P26-G3-Q129-T29 to Q148dref=0.155 mm, ΔT=31KAs discussed in Section3.1 calibration of the model for nucleation site density in (8) is possible by matching the measured axial profiles of wall superheat. Such profiles are available only for the test series P26-G2-Q74-Txx, but as discussed in Section 3.1, the nucleation site density may be expected to depend on parameters subcooling, flow rate, and heat flux only indirectly through the wall superheat. Moreover, the sensitivity of the final results with respect to the exponent p=1.85 is rather low; so only the prefactor Nref will be adjusted for a value of the reference temperature ΔTrefN=10K. Wall temperature profiles for two selected tests, P26-G2-Q74-T18 and P26-G2-Q74-T26, are shown in Figure 4. A value of Nref=3.0e7 m−2 gives good agreement with the respective data.Figure 4
Wall superheat depending on the reference nucleation site densityNref.


## 4.2.2. Coalescence and Break-Up Models

In the present work, bubble coalescence and breakup are described by the models proposed by Prince and Blanch [44] and by Luo and Svendsen [43]. To obtain agreement with the measurements, efficiency factors FC and FB were introduced and calibrated to match the measured radial bubble size profiles. Bubble size profiles for selected tests and several values of the calibration parameters FC and FB are shown in Figure 5. The procedure was successful so far as a good match could be obtained in all cases, but a different set of calibration parameters was necessary to this end (see Table 4). Therefore, this procedure has to be considered only as a first step and further model development needs to be done.Figure 5
Radial profiles of bubble size for different coefficientsFB and FC.As discussed in Section3.4.2, the lift force changes its sign at a certain bubble size. The value of this critical bubble size can be calculated as 1.0 mm for the material properties of R12 at a pressure of 2.62 MPa. Examination of the measured bubble sizes shows that these are lower than the critical value for all of the tests investigated here. Therefore, it is possible to use the simpler and computationally less expensive homogeneous MUSIG approach. For each test considered here, an equidistant bubble size discretization employing 15 size groups in the range of 0 to 1.5 mm was applied.


## 5. Reference Simulations

In this section, a more detailed discussion is given for a few selected tests, namely, P26-G2-Q74-T16, P26-G3-Q118-T16, and P26-G3-Q148-T28. As summarized in Table1, these tests have the same pressure but different subcooling, flow rate, and heat flux. A comparison of measured and calculated profiles for gas volume fraction, gas velocity, bubble size, and liquid temperature at the end of the heated length (x=3.5m) are shown in Figures 6, 7, and 8. Unfortunately, temperature measurements are available only for the first of these.Comparison of measured and calculated values for P26-G2-Q74-T16:P=2.62 MPa, G=2·103 kg m−2 s−1, Tin=70.7°C.
(a)
Gas volume fraction
(b)
Gas and liquid velocities
(c)
Radial bubble size profile
(d)
Bubble size distribution at four points
(e)
Temperature
(f)
Cross-sectionally averaged valuesComparison of measured and calculated values for P26-G3-Q118-T16:P=2.62 MPa, G=3·103 kg m−2 s−1, Tin=70.73°C.
(a)
Gas volume fraction
(b)
Gas and liquid velocities
(c)
Radial bubble size profile
(d)
Bubble size distribution at four points
(e)
Temperature
(f)
Cross-sectionally averaged valuesComparison of measured and calculated values for P26-G3-Q148-T28:P=2.62 MPa, G=3·103 kg m−2 s−1, Tin=58.2°C.
(a)
Gas volume fraction
(b)
Gas and liquid velocities
(c)
Radial bubble size profile
(d)
Bubble size distribution at four points
(e)
Temperature
(f)
Cross-sectionally averaged valuesFor the first two test cases, experimental gas volume fractions exhibit s-shaped profiles with an inflexion point which is absent for the third test case. In the simulations, the gas volume fraction profiles for all cases have a slope that decreases monotonously from the wall towards the center of the pipe. Cross-sectionally averaged values agree with the experimental ones for the first two cases but deviate for the last one. There does not appear to be a simple explanation for these deviations at hand. Gas velocities (part b of the figures) are slightly overpredicted for all cases. Average bubble size profiles (part c of the figures) agree with the data to a varying degree which is not surprising due to the rather crude modeling approach for coalescence and break-up rates. The calculated bubble size distributions at several radial positions (part d of the figures) show a rather narrow distribution near the pipe wall which broadens towards the center. This is due to the present modeling where all bubbles are generated with the same value ofdW at the wall. Where available, the temperature profiles show good agreement between simulation and experiment.For all test cases, overall reasonable agreement between experiment and simulation is found. The relative deviations in all other aspects of the comparison are in the 20 to 30% range which is comparable to other state-of-the-art investigations in the field.


## 6. Parametric Variations

In this section, tests in two series with varying subcooling for otherwise fixed parameters are compared, namely, P26-G2-Q74-T16 to P26-G2-Q74-T28, P26-G3-Q118-T16 to P26-G3-Q118-T28, and P26-G3-Q129-T28 to P26-G3-Q148-T28. The resulting changes in the profiles for gas volume fraction, bubble size, and liquid temperature at the end of the heated length (x=3.5 m) as the subcooling, considering the heat flux is varied, are shown in Figures 9, 10, and 11, respectively, where the left column gives the measured and the right columns the calculated values. For the latter two series, unfortunately no temperature measurements are available.Figure 9
Gas volume fraction, bubble size, and liquid temperature profiles for test series P26-G2-Q74-Txx.Figure 10
Gas volume fraction, bubble size, and liquid temperature profiles for test series P26-G3-Q118-Txx.Figure 11
Gas volume fraction, bubble size, and liquid temperature profiles for test series P26-G3-Qxx-T28.Thermal conditions in the liquid at the measurement location depend on all three parameters inlet subcooling, flow rate, and heat flux. The calculated temperature profiles compare quite well with the experimental ones for the case where the latter are available. As a reference the calculated profiles are shown also for the two cases where data are not provided.For all tests, the bubble size increases with increasing the distance from the wall despite the fact that the bulk liquid is subcooled. Clearly this must be due to coalescence of the bubbles. This phenomenon could not be captured by the monodisperse model approach of Krepper and Rzehak [10]. For both test series upon decreasing the subcooling, a strong increase of bubble size is observed in the center of the pipe while the detachment size changes only relatively little. The present polydisperse modelling approach describes this behaviour at least qualitatively correct.Comparing the radial gas volume fraction profiles with decreasing subcooling, a broadening of the wall peaking profile can be observed for both test series. This effect is due to a lower condensation rate in the bulk liquid as the subcooling decreases. Again this behaviour is described qualitatively correct by the present modeling approach where bubbles of all sizes move with the same velocity. A change from wall to core peaking profiles as for the test series considered by Krepper et al. [11] is not observed here. This matches with the explanation that such a transition is related to the sign change in the lift force as discussed in Section 3.4.2. For the presently investigated tests, the bubble size remains below this value while for the cases of Krepper et al. [11] it became larger for the lower values of subcooling.


## 7. Summary and Conclusions

Boiling at a heated wall has been simulated by an Euler/Euler description of two-phase flow combined with a heat flux partitioning model describing the microscopic phenomena at the wall by empirical correlations adapted to experimental data. Such an approach was previously used and adjusted to boiling experiments with water at a pressure of several MPa. We here have investigated the applicability and necessary readjustments for similar tests using R12 at the DEBORA facility. At the same time, the bubble size distribution in the bulk was described by a population balance approach by coupling the wall boiling model with the MUSIG model.A critical review of the detailed correlations used in previous work shows that some of the parameters used in these correlations have to be carefully recalibrated for the present applications. The DEBORA tests provide a large body of information that can be used to this end. Quantities with a strong influence on the amount of produced steam are the bubble size at detachment and the nucleation site density. The former can be taken straight forwardly from the measurements. On the latter, unfortunately no direct information is available; however, by matching the temperature of the heated wall, this gap can be closed. Previous work [10] has shown that the recalibration results in different values for different pressure levels. The present results suggest that for the bubble detachment size, a different calibration is needed for different flow rates but not for different heat fluxes while for the nucleation site density, the same calibration can be used independently of either liquid subcooling, flow rate, and heat flux.The measured gas bubble size profiles show an increase of the bubble size with increased distance from the heated wall. A monodispersed treatment is not able to capture this phenomenon, but including polydispersity by means of a MUSIG approach and suitable models especially for bubble coalescence this phenomenon can be described. In contrast to a previously investigated series of test cases [11], as the subcooling is decreased only a broadening of the wall peaking gas volume fraction profile occurs, but no transition to a core peaking profile is observed. Again this phenomenon could not be captured by a model with a monodisperse bubble size but can be described using a homogeneous MUSIG approach since the bubbles remain small enough so the sign change in the lift force does not occur.A complete polydispersed description requires that processes of coalescence/breakup and condensation/evaporation must be modeled explicitly. For the latter, a suitable model is readily obtained from first principles with the aid of a heat transfer correlation like that of Ranz and Marshall [46]. Unfortunately the former are the more important processes, and for these the situation is much less clear. In the present work, the commonly applied models for bubble coalescence according to Prince and Blanch [44] and for bubble breakup according to Luo and Svendsen [43] were used as a first step. To reach a fair agreement with the measurements, calibration factors had to be introduced, which need to be adjusted according to the experimental condition. In this way the suitability of the general model framework could be demonstrated in principle. For a trusted prediction, further development of the coalescence and break-up models is necessary.Bubble coalescence and breakup are heavily influenced by two-phase flow turbulence. Unfortunately in the literature, only few measurements of turbulent characteristics of two phase flow can be found and even less when boiling occurs. Furthermore, models of bubble-induced turbulence working well for air/water flow may fail for steam/water flow at higher pressure or for refrigerants. In the present work, the selection of a specific model for bubble effects on the turbulence is confirmed mainly by plausibility of the final results. Hence, a more systematic investigation of approaches to modeling bubbly turbulence would be desirable.Finally, looking carefully at the figures showing the gas volume fraction profiles in the near wall region, the calculated gas volume fraction is systematically too large. Reasons could be a missing force pushing the bubbles away from the wall or the neglect of swarm effects in the models of drag and lift forces even at gas volume fractions around 50%. Furthermore, the application of the simple heat transfer correlation of Ranz and Marshall [46] might be questionable. In transferring models used successfully for adiabatic air/water flows to the DEBORA tests, it should be noted that bubbles are much smaller here which may require changes beyond simple recalibration of parameters.Overall, our results confirm the great potential of the Euler/Euler two-phase flow and heat flux partitioning models for the simulation of subcooled flow boiling in industrial applications while at the same time highlighting the need for specific model improvements in order to achieve highly accurate quantitative predictions.



---

*Source: 687494-2013-01-22.xml*