CFD for Subcooled Flow Boiling : Parametric Variations

We investigate the present capabilities of CFD for wall boiling.e computational model used combines the Euler/Euler two-phase �ow description with heat �ux 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 dichlorodi�uoromethane (R12) as the working �uid. is facilitated measurements of radial pro�les 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, �ow rate, and heat �ux. It is investigated to which extent a set of calibrated model parameters suffices to cover at least a certain parameter range.


Introduction
Subcooled �ow 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 signi�cant area.is occurs at the critical heat �ux 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 �uid, the critical heat �ux depends on the �ow parameters as well as the geometry of the �ow domain.e veri�cation of design improvements and their in�uence on the critical heat �ux requires expensive experiments.erefore, 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 �ux were developed and �tted to data obtained from experimental tests.ese have been implemented mainly in purpose-speci�c 1D codes and applied for engineering design calculations.However, these correlations are valid only in the limited region of �uid properties, working conditions, and geometry corresponding to the tests to which they were �tted.�sing large look-up tables based on a great number of experiments, a signi�cant range of �uid properties and working conditions can be covered.But this method is still limited to only that speci�c 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 �ux reliably.A precondition would be the complete understanding and simulation of boiling as a preliminary state towards critical heat �ux.
For engineering calculations, currently the most widely used CFD approach to model two-phase �ows with signi�cant volume fractions of both phases is the Eulerian two-�uid framework of interpenetrating continua (see, e.g., [1][2][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.ese exchange terms consist of analytical or empirical correlations, expressing the interfacial forces, as well as heat and mass �uxes, as functions of the average �ow parameters.Since most of these correlations are highly problem-speci�c, their range of validity has to be carefully considered and the entire model has to be validated against experiments.
For the case of boiling �ows, where heat is transferred into the �uid 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 �ow boiling of water at high pressure �owing 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 �ow in a nuclear reactor.However, typical �ow 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 signi�cant challenges for measurements.e use of refrigerants can greatly relieve this burden.In the French DEBORA tests [8,9], dichlorodi�uoromethane (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.e 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.is enables a measurement of radial pro�les for gas volume fraction, temperature, liquid and gas velocities, and bubble sizes which allows a stringent validation of CFD models.
e applicability of CFD models to the DEBORA tests was recently studied in Krepper and Rzehak [10,11] to demonstrate their general validity and identify speci�c weak points.By coupling a population balance approach to the wall boiling model, it was shown that a reasonable prediction of the measured radial pro�les 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, �ow rate, and heat �ux.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.e paper is organized as follows.A brief summary of the DEBORA test facility and selected data is given in Section 2. In Section 3 the details of the models used for the simulations are described.e presentation is adapted from our previous work [10,11] to make the paper self-contained.e setup of the simulations together with necessary calibrations of the models is described in Section 4. e 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 pro�les and the gas volume fraction pro�les provided the models are suitably calibrated.A discussion of the achievements so far and needs for further research are �nally offered in Section 7.

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, dichlorodi�uoromethane (R12) is heated over a pipe length of 3.5 m as sketched in Figure 1.e facility is operated at various pressures, �ow rates, heat �uxes, and subcoolings.e radial pro�les for gas volume fraction and gas velocity at the end of the heated length are measured by means of an optical probe.Furthermore, pro�les of bubble size at this position are available.In addition, radial liquid temperature pro�les as well as axial pro�les of the wall temperature are measured by thermocouples.Unfortunately, temperature-and phase-related measurements are not both available for all experimental conditions.Table 1 lists the parameters of the DEBORA tests selected for the present investigation.e �rst 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].e next three cases extend the range of inlet subcooling at the same �ow rate and heat �ux.e second series of four cases covers the same range of inlet subcooling but at different values of �ow rate and heat �ux.In the last group, the heat �ux is varied at constant values of inlet subcooling and �ow rate.For all investigated tests, the system pressure is kept at a �xed value of    MPa.e question is to what extent the model is capable to reproduce the data under all of these conditions with as little recalibration as possible.

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.e 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 pro�les.
For R12 liquid and vapor, the relevant material properties were taken from the National Institute of Standards and Technology (NIST) Standard Reference Database on ermophysical Properties of Fluid Systems (http://webbook.nist.gov/chemistry/�uid). e saturation properties for the value of pressure considered, that is,   .62 MPa, are listed in Table 2. Based on the NIST data, �uid property tables were generated and used in the calculations.

The Models
e general equations for diabatic two-phase �ow in the Euler/Euler framework of interpenetrating continua have been reviewed in many places before (e.g., [1][2][3]).erefore, 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).
e �rst ma�or 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 �ux 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 speci�c conditions under investigation.e basis on which such a recalibration can be carried out for varying conditions is discussed.e �nal 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 �ow regime, it is convenient to use an equivalent Sauter diameter and work with the bubble size.Obviously, in boiling �ows the bubble size may change due to both condensation/evaporation and bubble coalescence/breakup. e importance of taking into account the latter processes has been shown in Krepper et al. [11].is 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, bubbleinduced 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 �ux.
Turbulent �uctuations are modeled by a shear stress turbulence (SST) model according to Menter [15,16] applied to the liquid phase.is corresponds to a k- model near the walls and a k- model far from the walls.e 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 �ows 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 pro�les.
For momentum exchange between the phases, �nally, li and turbulent dispersion forces are included in the model in addition to the ubiquitous drag force.For adiabatic �ows, 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 in�uence turned out to be small.erefore, in the present study, this force was neglected.In general the applicability of a wall force to �ow boiling should be investigated further.

Modeling of Boiling at a Heated Wall.
In boiling, heat is transported from the hot wall to the �uid by several different mechanisms.On parts of the wall, where no bubbles reside, heat �ows directly to the subcooled liquid in the same way as that in single-phase �ow.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 �ow is brought into contact with the hot wall which leads to additional cooling.is mechanism, which is obviously not present in singlephase �ows, is termed quenching.
Accordingly, the given external heat �ux  tot , applied to the heated wall, is written as a sum of three parts: where   , the heat �ux components are expressed as discussed in the following.e turbulent convection heat �ux is calculated in the CFX model version (see Wintterle [19]) in much the same way as for a pure liquid �ow without boiling but multiplied by the fraction of area unaffected by the bubbles; that is, Here ℎ  is the heat transfer coefficient which is written using the temperature wall function  + ( + ) known from Kader [20] as where nondimensional variables (indicated by superscript "+") and the friction velocity   are de�ned as usual.Note that ( 2)-( 3) may be evaluated at any location , provided it is used consistently whenever a variable depends on position.  is represented in terms of the quenching heat transfer coefficient ℎ  : A grid independent solution for   is obtained by evaluating the nondimensional temperature pro�le at a �xed value of  + .
e evaporation heat �ux   is obtained via the evaporation mass �ux at the wall: where   is the latent heat of evaporation and the generated vapor mass ṁ is expressed in terms of the bubble diameter at detachment   , bubble generation frequency , and nucleation site density  as e bubble size at detachment is known to depend on a large number of variables, namely, the liquid subcooling, the �ow rate, the heat �ux, the wall superheat, the �uid properties at the system pressure, and the contact angle between the �uid 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 signi�cantly not only with the material combination of working �uid 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 �ow rate and heat �ux are constant along the pipe.e available data on axial dependence of the wall superheat (cf., Section 5) suggest that this is also the case for the wall superheat except probably close to the beginning of the heated section.us further simpli�cation may be sought by �xing these conditions and focusing on the dependence with respect to the liquid subcooling.
Common practice in engineering simulations of �ow 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 �tted to a correlation Here again   is obtained by evaluating the nondimensional temperature pro�le of Kader [20] at a �xed value of  + .e values of the �t parameters  ref and Δ ref  then must be expected to depend on all the other variables.By suitable calibration of these �t parameters, it has been shown for a number of systems that reasonable agreement with data can be obtained [6,10].e approach taken in the present investigation is to allow a different calibration for different values of �ow rate and heat �ux.�ventually it will turn out that within certain limits the same calibration can be used even for different conditions.Detailed values will be given in Section 4.2.1.
Concerning the nucleation site density, most available correlations are expressed in the form of power laws depending on the wall superheat as 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.is 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.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 identi�ed with the wall superheat� namely, Due to the small spatial scales involved, this process may be expected to be independent of other parameters like liquid subcooling, �ow rate, and heat �ux.is provides some rationale for the appearance of the wall superheat as the only variable in correlations for the nucleation site density.e diversity of speci�c 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 �nish the surface.ese processes are very diverse and in most boiling experiments not speci�cally controlled.Moreover, a characterization of the resulting surface topography is rarely available.
Parameter variations showed that the assumed nucleation site density has almost no in�uence on the calculated liquid temperature, a small in�uence on the calculated gas volume fraction, but a strong in�uence on the calculated wall superheat   −  sat .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 in�uenced by vapor bubbles   is given by Here,  is the so-called bubble in�uence 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 liing 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 PI� measurement of the �ow �eld around departing bubbles, although the quality of the images presented is rather poor.
Since   = 1 corresponds to the case where the whole surface is under the in�uence of bubbles,   as calculated by (10) has to be limited to values smaller than this limit.Moreover, it should be kept in mind that already as   approaches 1, the assumptions of the model are not really satis�ed anymore.
e bubble detachment frequency  is given according to Cole [36] as a function of the detachment size   as Science and Technology of Nuclear Installations Relations of this type have been reviewed critically by Ivey [37] and Ceumern-Lindenstjerna [38].e quenching heat transfer coefficient is calculated using the analytical solution for one-dimensional transient conduction, as suggested by Mikic and Rohsenow [39]: where  wait 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 �xed fraction of the bubble departure period: 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 �ux is 75% of the critical heat �ux and larger.

Population Balance Approach to Bubble Size Distribution.
To describe polydispersed �ows within a purely �ulerian approach, a number of different (multiple) bubble size groups  =    is considered, each representing bubbles of typical size   .e fraction of gas volume contained in each bubble size group is denoted as   so that the total gas volume fraction is given by It is useful to also de�ne occupation numbers   =   /  giving the contribution of each size group to the total gas volume fraction.Obviously we then have ∑    = .
For each size group, the equation of mass conservation assumes the form where the right-hand side gives the net source of mass for group  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   is calculated according to (14) above.In addition, also the bubble size   appears which is taken in the Sauter sense representing the interfacial area   = 6  /  .In order to preserve this interpretation,   is calculated from the occupation number and bubble size for each group as e 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 li 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.
Here   is an approximation to the delta-function   +   −  ).For the break-up and coalescence kernel functions  and , 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.

Condensation and Evaporation Including Boiling at the
Wall.When condensation or evaporation occurs, the volume fraction in size group  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 group  the direct mass transfer to the liquid is given by where similar to previous work [10] the assumption has been made that the gas is at saturation temperature.e total source terms for size class  including also the ensuing change of bubble size, that is, Γ phase  in (15), have been derived recently by Lucas et al. [45] as where   =  3 /6 is the mass of each bubble in size group .Basing the calculation on bubble mass rather than size for compressible �ows 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 �ows, no differences between mass-and sizebased groups arise.
In principle  Γ  should be evaluated with the group size   , but for practical reasons, an approximation is used where a direct transfer term  Γ  is calculated according to (21) only for each velocity group  using the Sauter mean diameter of (16) and the area based fraction thereof is used for  Γ  ; that is, In this way, the size dependency of the factor   in ( 21) is treated exactly, but that of the factor ℎ  is not.e liquid side heat transfer coefficient, �nally, is calculated according to Ranz and Marshall [46] as 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 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,  rpi , is included in (19) following [48].is source term applies only to the equation corresponding to the size group whose diameter is the closest to the bubble detachment diameter   .It is given by the evaporation mass �ux computed in the wall heat partitioning distributed evenly throughout the grid cells adjacent to the heated wall; that is, where ṁ is given by ( 6) and  and  are wall surface area and volume of the corresponding grid cell, respectively.

Modeling of Turbulence in a Two-Phase Flow.
Liquidphase 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.e reverse effect exerted by the turbulent �uctuations 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 �ows 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].

Bubble-Induced
Turbulence. e effects of the dispersed gas bubbles on the turbulence in the continuous liquid phase are modeled by introducing appropriate source terms in the k-and -or -equations.e bubble-induced source    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 where  drag is the drag force given by (31) and   −   is the relative velocity.e bubble-induced source    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, 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], P�eger and Becker [51], and Politano et al. [18] who assumed that the usual turbulence time scale may be used; that is, e coefficient  3 = .92 was found to give satisfactory results.

Science and Technology of Nuclear Installations
An equivalent source term for the -equation is derived from    and    as where   = 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 �uid domain.

Turbulent Wall Function for Boiling Flow.
A wall function for boiling �ow is obtained by considering that the presence of the bubbles on the wall forces the liquid into a similar �ow pattern as that observed in single-phase turbulent �ow with wall roughness (e.g., Ramstorfer et al. [52] and Končar and Borut [53]).e latter is described by a modi�ed law of the wall [54,55] where   = (  / is the friction velocity,  = /(   is the viscous length scale, s is the hydrodynamic roughness length scale, and   0.1 is the von Karman constant.For �ow over smooth walls,   const  5.5, while for �ow over rough walls,  * is a function of / which in the limit of large s, that is, for the so-called fully rough walls, tends to a constant value of 8.5.e hydrodynamic roughness s may be related directly to the bubble size and nucleation site density as [10]  =    3   .
e constant of proportionality is not known from theoretical considerations at present, but a value of   = 1.0 has been found to yield good results.
As shown in [10], the representation of radial velocity pro�les is greatly enhanced by employing this two-phase boiling wall function over the oen used single-phase wall function.

Modeling of the Momentum Transfer.
For momentum exchange between the phases, the Ishii and Zuber [56] drag law was used.Furthermore, a li 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.erefore, in the present study this force was neglected.In the following expressions, the forces for the dispersed gaseous phase are given.

Drag Force. e volumetric source of gaseous momentum due to drag exerted by the liquid is given by
e drag coefficient   is calculated according to Ishii and Zuber [56] considering three different shape regimes as (33)

Li Force.
A li force due to interaction of the bubble with the shear �eld of the liquid was �rst introduced to two-�uid simulations by Zun [59].e corresponding momentum source is given by e classical li force, which has a positive coefficient   , acts in the direction of decreasing liquid velocity.In case of cocurrent upwards pipe �ow this is the direction towards the pipe wall.Experimental [57] and numerical [60] investigations showed that the direction of the li 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 �ow of a glycerol water solution, the following correlation for the li coefficient was derived: Here  ⟂ is the maximum horizontal dimension of the bubble.It is calculated using an empirical correlation for the aspect ratio by the following equation [61]: For the water-air system at normal conditions,   changes its sign at   = 5.8 mm which was con�rmed by investigations of polydispersed upward vertical air�water bubbly �ow [62,63].For R12 this value is decreased substantially to about 1.0 mm at 2.62 MPa.

Turbulent Dispersion Force. e turbulent dispersion
force is the result of the turbulent �uctuations of liquid velocity.Burns et al. [58] derived an expression by Favre averaging the drag force as Here,  TD = 0.9 is an empirical parameter.
As noted in Krepper and Rzehak [10], the effect of an additional wall force is small.erefore, in the present study, this force was neglected.

Model Setup
e models described in the previous section present a rather general framework that can be specialized to the simulation of any subcooled �ow boiling problem.e necessary spec-i�cations to simulate the D�BORA test cases will now be described.ese comprise prescription of �ow domain and boundary conditions and speci�cation of grid and bubble classes as well as calibration of model parameters.
e simulations are performed by ANSYS CFX 13.For details of the numerical procedures, we refer to the user guide [64].

Geometry and Boundary
Conditions. e tests were simulated in a quasi-2D cylindrical geometry, that is, a narrow cylindrical sector with symmetry boundary conditions imposed on the side faces.e validity of this simpli�cation has been veri�ed 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 �ow development zone has been added to obtain at the beginning of the heated section a fully developed turbulent �ow that is independent of the detailed conditions imposed at the inlet.e required length of this �ow 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 �ux partitioning discussed in Section 3.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 le 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.erefore, 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 �ow 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 of  + ≈ 200.For the test P26-G2-Q74-T16, a grid re�nement study was performed which showed no change of the results until this value of  + has decreased to about 70.For still smaller values, no convergence could be achieved.is 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.

Calibration of Parameters
4.2.1.Wall Boiling Model.e necessity to recalibrate the boiling model parameters to the working �uid 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.e outermost measurement points of the experimental bubble size pro�les are taken as representative of the detachment size.e liquid subcooling at the measurement location is determined from the averaged measured liquid temperature pro�le 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.e 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,  ref , and the reference subcooling, Δ ref , were then found by �tting the expression (7) to the values of Table 3. e 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 �ow rate but not on the heat �ux.
As discussed in Section 3.1 calibration of the model for nucleation site density in ( 8) is possible by matching the measured axial pro�les of wall superheat.Such pro�les 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, �ow rate, and heat �ux only indirectly through the wall superheat.Moreover, the sensitivity of the �nal results with respect to the exponent    is rather low; so only the prefactor  ref will be adjusted for a value of the reference temperature Δ ref  0 K. Wall temperature pro�les for two selected tests, P26-G2-Q74-T18 and P26-G2-Q74-T26, are shown in Figure 4.A value of  ref  30 m −2 gives good agreement with the respective data.

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   and   were introduced and calibrated to match the measured radial bubble size pro�les.Bubble size pro�les for selected tests and several values of the calibration parameters   and   are shown in Figure 5. e 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).erefore, this procedure has to be considered only as a �rst step and further model development needs to be done.As discussed in Section 3.4.2, the li force changes its sign at a certain bubble size.e 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.erefore, 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.

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 Table 1, these tests have the same pressure but different subcooling, �ow rate, and heat �ux.A comparison of measured and calculated pro�les for gas volume fraction, gas velocity, bubble size, and liquid temperature at the end of the heated length (  3 m) are shown in Figures 6, 7, and 8. Unfortunately, temperature measurements are available only for the �rst of these.
For the �rst two test cases, experimental gas volume fractions exhibit s-shaped pro�les with an in�exion point which is absent for the third test case.In the simulations, the gas volume fraction pro�les 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 �rst two cases but deviate for the last one.ere does not appear to be a simple explanation for these deviations at hand.Gas velocities (part b of the �gures) are slightly overpredicted for all cases.Average bubble size pro�les (part c of the �gures) agree with the data to a varying degree which is not surprising due to the rather crude modeling approach for coalescence and breakup rates.e calculated bubble size distributions at several radial positions (part d of the �gures) show a rather narrow distribution near the pipe wall which broadens towards the center.is is due to the present modeling where all bubbles are generated with the same value of   at the wall.Where available, the temperature pro�les show good agreement between simulation and experiment.
For all test cases, overall reasonable agreement between experiment and simulation is found.e 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 �eld.

Parametric Variations
In this section, tests in two series with varying subcooling for otherwise �xed 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.e resulting changes in the pro�les for gas volume fraction, bubble size, and liquid temperature at the end of the heated length (  3 m) as the subcooling, considering the heat �ux is varied, are shown in Figures 9, 10, and 11, respectively, where the le column gives the measured and the right columns the calculated values.For the latter two series, unfortunately no temperature measurements are available.
ermal conditions in the liquid at the measurement location depend on all three parameters inlet subcooling, �ow rate, and heat �ux.e calculated temperature pro�les compare quite well with the experimental ones for the case where the latter are available.As a reference the calculated pro�les 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.is 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.e present polydisperse modelling approach describes this behaviour at least qualitatively correct.
Comparing the radial gas volume fraction pro�les with decreasing subcooling, a broadening of the wall peaking pro�le can be observed for both test series.is 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 pro�les as for the test series considered by Krepper et al. [11] is not observed here.is matches with the explanation that such a transition is related to the sign change in the li 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.

Summary and Conclusions
Boiling at a heated wall has been simulated by an Euler/Euler description of two-phase �ow combined with a heat �ux 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.e DEBORA tests provide a large body of information that can be used to this end.Quantities with a strong in�uence on the amount of produced steam are the bubble size at detachment and the nucleation site density.e 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.e present results suggest that for the bubble detachment size, a different calibration is needed for different �ow rates but not for different heat �uxes while for the nucleation site density, the same calibration can be used independently of either liquid subcooling, �ow rate, and heat �ux.
e measured gas bubble size pro�les 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 pro�le occurs, but no transition to a core peaking pro�le 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 li 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 �rst 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 �rst step.To reach a fair agreement with the measurements,      only few measurements of turbulent characteristics of two phase �ow can be found and even less when boiling occurs.Furthermore, models of bubble-induced turbulence working well for air/water �ow may fail for steam/water �ow at higher pressure or for refrigerants.In the present work, the selection of a speci�c model for bubble effects on the turbulence is con�rmed mainly by plausibility of the �nal results.Hence, a more systematic investigation of approaches to modeling bubbly turbulence would be desirable.Finally, looking carefully at the �gures showing the gas volume fraction pro�les 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 li 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 �ows to the DEB�RA tests, it should be noted that bubbles are much smaller here which may require changes beyond simple recalibration of parameters.
�verall, our results con�rm the great potential of the Euler/Euler two-phase �ow and heat �ux partitioning models for the simulation of subcooled �ow boiling in industrial applications while at the same time highlighting the need for speci�c model improvements in order to achieve highly accurate quantitative predictions.tot: Total : Wall +: Nondimensional.

F 1 :
Sketch of the DEBORA test geometry.

F 2 :
Sketch of bubble nucleus.

F 4 :
Wall superheat depending on the reference nucleation site density  ref .

F 11 :
�as �olu�e fraction� �u��le si�e� and li�uid te��erature �ro�les for test series �����������T��� 3.2.1.Coalescence and Breakup.e net mass source for size group  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    to group  and coalescence of smaller bubbles from size groups    to group  as well as bubble death rates due to breakup of bubbles from size group  to smaller bubbles in groups    and the coalescence of bubbles from size group  with bubbles from any other group to even larger ones which belong to groups   .at is, ,sphere , min  ,ellipse ,  ,cap  ,

T 3 :
Parameters used for the bubble detachment model.
F 9: �as volume fraction, bubble si�e, and li�uid temperature pro�les for test series ���-��-�7�-Txx.calibrationfactors 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.�ubblecoalescence and breakup are heavily in�uenced by two-phase �ow turbulence.�nfortunately in the literature, sat F 10: �as �olu�e fraction� �u��le si�e� and li�uid te��erature �ro�les for test series ��������118�T���