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.
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 Section 2. 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.
Table 1 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 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.
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, MPa, are listed in Table 2. Based on the NIST data, fluid property tables were generated and used in the calculations.
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 flux , applied to the heated wall, is written as a sum of three parts: where , , and 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 , which satisfies the wall heat flux balance. Denoting the fraction of area influenced by the bubbles as , 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, 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 defined 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 profile at a fixed value of .
The evaporation heat flux is obtained via the evaporation mass flux 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 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., Section 5) 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 Here again is obtained by evaluating the nondimensional temperature profile of Kader [20] at a fixed value of . The values of the fit parameters and 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 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. 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.
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, 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 superheat . 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 is given by Here, 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.
Since corresponds to the case where the whole surface is under the influence 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 satisfied anymore.
The bubble detachment frequency is given according to Cole [36] as a function of the detachment size as 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]: where 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: 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 groups is considered, each representing bubbles of typical size . The 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 define 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 . In order to preserve this interpretation, is calculated from the occupation number and bubble size for each group as 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 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 . That is, The birth and death rates in turn are commonly expressed in terms of the coalescence and break-up kernels such that 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.
3.2.2. 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. The total source terms for size class including also the ensuing change of bubble size, that is, in (15), have been derived recently by Lucas et al. [45] as where is the mass of each bubble in size group . 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 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. The liquid side heat transfer coefficient, finally, 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, , 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 . 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, where is given by (6) and and 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 the - and - or -equations.
The bubble-induced source in the -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 is the drag force given by (31) and is the relative velocity.
The 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 -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], Pfleger and Becker [51], and Politano et al. [18] who assumed that the usual turbulence time scale may be used; that is, The coefficient was found to give satisfactory results.
An equivalent source term for the -equation is derived from and as where , 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] where is the friction velocity, is the viscous length scale, is the hydrodynamic roughness length scale, and is the von Karman constant. For flow over smooth walls, , while for flow over rough walls, is a function of which in the limit of large , that is, for the so-called fully rough walls, tends to a constant value of 8.5.
The hydrodynamic roughness may be related directly to the bubble size and nucleation site density as [10]
The constant of proportionality is not known from theoretical considerations at present, but a value of 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 The drag coefficient is calculated according to Ishii and Zuber [56] considering three different shape regimes as where
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 The classical lift force, which has a positive coefficient , 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: This coefficient depends on the modified Eötvös number given by 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 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 Here, 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 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 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 of . For the test P26-G2-Q74-T16, a grid refinement 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. 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, , and the reference subcooling, , 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.
| (a) P26-G2-Q74-T16 to T28 mm, |
| (b) P26-G3-Q118-T16 to T28 mm, |
| (c) P26-G3-Q129-T29 to Q148 mm, |
As discussed in Section 3.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 is rather low; so only the prefactor will be adjusted for a value of the reference temperature . Wall temperature profiles for two selected tests, P26-G2-Q74-T18 and P26-G2-Q74-T26, are shown in Figure 4. A value of m−2 gives good agreement with the respective data.
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 and 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 and 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.
As discussed in Section 3.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 Table 1, 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 are shown in Figures 6, 7, and 8. Unfortunately, temperature measurements are available only for the first of these.
(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 values
(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 values
(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 values
For 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 of 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 ( 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.
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.
Nomenclature
| : | Interfacial area density |
| : | Bubble-induced turbulence coefficient [17] model |
| : | Drag coefficient |
| : | Lift coefficient |
| : | Specific heat capacity at constant pressure |
| : | Turbulent dispersion coefficient |
| : | Virtual mass force coefficient |
| : | Wall force coefficient |
| : | Shear-induced turbulence coefficient ( model) |
| : | Bulk bubble diameter (m) |
| : | Bubble diameter perpendicular to main motion (m) |
| : | Bubble detachment diameter (m) |
| : | Pipe diameter (m) |
| : | Eötvös’number |
| : | Bubble detachment frequency (Hz) |
| : | Drag force |
| : | Lift force |
| : | Turbulent dispersion force |
| : | Virtual mass force |
| : | Wall force |
| : | Acceleration of gravity |
| : | Heat transfer coefficient for single-phase convection |
| : | Heat transfer coefficient for bulk evaporation/condensation |
| : | Heat transfer coefficient for quenching |
| : | Specific enthalpy |
| : | Jakob’s number |
| : | Thermal conductivity |
| : | Turbulent kinetic energy |
| : | Length scale (m) |
| : | Massflux |
| : | Morton’snumber |
| : | Nucleation site density () |
| : | Pressure (Pa) |
| Pr: | Prandtl number |
| : | Wall heat flux |
| : | Heat flux due to single-phase convection |
| : | Heat flux due to quenching |
| : | Heat flux due to evaporation |
| : | Radial coordinate (m) |
| Re: | Reynolds’ number |
| : | Hydrodynamic wall roughness (m) |
| : | Time (s) |
| : | Waiting time (s) |
| : | Temperature (K) |
| : | Saturation temperature (K) |
| : | Liquid subcooling (K) |
| : | Wall superheat (K) |
| : | Wall temperature (K) |
| : | Velocity (m ) |
| : | Friction velocity |
| : | Velocity scale (m) |
| : | Volume () |
| : | Axial coordinate (m) |
| : | Distance to the wall (m) |
| : | Volume fraction |
| : | Viscous length scale (m) |
| : | Turbulent dissipation rate |
| : | Temperature scale (K) |
| : | Dynamic viscosity |
| : | Kinematic viscosity |
| : | Density |
| : | Surface tension |
| : | Wall shear stress . |
| Bulk | |
| : | Due to single-phase convection |
| : | Due to evaporation |
| : | Gas |
| : | Interface |
| : | Liquid |
| : | Due to quenching |
| sat: | Saturation |
| tot: | Total |
| : | Wall |
| +: | Nondimensional. |
Acknowledgment
This work is funded by the Federal Ministry of Education and Research (Contract no. 02NUK010A).