Simulation input parameters for R12 tests
The predictions of the CFD Two Fluid Model for bubbly flows in ducts depend on constitutive relations for various interfacial and wall transfer terms. The objective is to use state-of-the-art data covering a wide range of experimental conditions to validate the constitutive relations for two cases of practical interest: adiabatic flows and subcooled and saturated boiling flows. The key to the validation is that all the data was obtained with either double conductivity or double optical needle probes to measure the bubble sizes.
This chapter is divided into two parts. The first part considers adiabatic flows where the interfacial and wall momentum transfer is represented by various forces acting between the two phases and between the wall and each phase. There is a considerable amount of work on the development of constitutive relations for these forces in the open literature. However, these relations have been derived under special flow conditions which may not hold in general. Hence, every force term needs to be calibrated and this introduces the closure coefficients which are dependent on the flow conditions. The air-water experimental validation and calibration of lift, wall and turbulent diffusion forces, which determine the distribution of void fraction in the direction normal to the flow, is carried out.
The second part describes the implementation and validation of heat and mass transfer models required for boiling flows. In particular the wall heat transfer model is based on splitting the wall heat flux into three components: single phase convection, evaporation at the bubble surface and quenching after bubble departure. The evaporation and quenching components depend on three parameters: bubble nucleation site density, bubble departure diameter, and bubble departure frequency. These parameters depend on complex physical phenomena that are not completely understood. However, a considerable amount of experimental data, scaling and various models exist for each one. Nevertheless, it is impossible to test these three separate phenomena simultaneously and in particular the bubble departure diameter remains an open problem (Colin et al., 2009). Therefore only the other two parameters were validated and the validation was performed with two data sets where the bubble diameters were measured: one with high pressure R12 and another with water at atmospheric pressure. The high pressure Freon data corresponds to a fluid-vapor density ratio which is equivalent to that of water-steam at 150 bar. Thus the liquid to vapour density ratio between Freon and water data varies over two orders of magnitude. The surface tension also varies over a wide range from 0.0017 to 0.057 N/m. The ratio of the flow channel hydraulic diameter to the bubble diameter encountered during these simulations varied from 4 to 40. As anticipated, the two-fluid model performs very well when bubbles are small compared to the channel width. However, appropriate modifications have to be made to the two-fluid model to account for cases where the bubbles are too large.
The constitutive relations remain the main limitation of any Two Fluid Model. In particular the present model lacks constitutive relations for the bubble size at the wall and in the bulk. Hopefully, current research such as the development of interfacial area transport models (e.g., Ishii & Hibiki, 2006) will help to resolve this shortcoming in the foreseeable future.
The commercial CFD code ANSYS CFX (version 11 for adiabatic and 12-beta for boiling flows) was used for all the simulations reported in this chapter.
3. Two fluid model mass and momentum equationsMass Conservation
The ensemble averaged two-fluid model continuity equations (Ishii & Hibiki, 2006) governing the motion of
For the bubbly flow analyzed during this study, the two-fluid model is comprised of two fields: liquid continuous (k = 1) and dispersed bubbles (k = 2) and the mass transfer across the interface is zero for adiabatic flows.
The interfacial momentum transfer force comprises of force terms due to drag, turbulent diffusion, lift and a wall force:
Other forces like the Virtual Mass force and Basset force are also present but are assumed to be negligible.
The drag force of the bubbles is given by,
The drag coefficient is given by the Ishii-Zuber correlation (1979). The relative velocity is given by
The lift force is given by (Auton, 1987):
The lift coefficient used for the calculations is
The wall force model accounts for the effect that keeps the centers of the bubbles no closer than approximately one bubble radius from the wall. This force is important when the lift force is present. It is given by (Antal et al., 1991):
The turbulent diffusion force is not yet well understood. The turbulent diffusion force model (Lopez de Bertodano, 1998) is based on a probability distribution function analysis of the two-fluid model (Reeks, 1991) which has the following form:
The time constant of the bubbles is derived from the equation of motion of a single bubble as:
The time constant of the turbulent eddies is proportional to the turnover time of the eddies and it is calculated according to the model by Mostafa and Mongia (1988):
The effect of eddy cross over is accounted for by a simple model when the particle is moving very fast with respect to the eddies, the time constant is given by:
This approximation has the right asymptotic behavior at both extremes when the bubbles are very small or very large. The model given by equations (7) to (11) may be simplified under several hypothetical asymptotic conditions. We consider two cases:
where the turbulent diffusion coefficient used for the present calculations is
The closure for the Reynolds stresses in equation (2) is based on the two-phase
The equivalent expression for the turbulent diffusivity of momentum of the liquid phase is given by Sato et al. (1981):
where the first term on the RHS corresponds to the
4. Adiabatic CFD model assessment with atmospheric air-water data
Serizawa et al. (1986) performed measurements of void fraction and interfacial area concentration in upward air water bubbly flow in a 30 mm diameter tube. A double conductivity probe was used to measure void fraction and interfacial area concentration. The radial profiles of void fraction and interfacial area concentration were obtained at 2.5 m axial distance from the inlet. The flow rates covered the range jL = 0.5 - 5 m/s and jG = 0 - 0.54 m/s.
The data set was used to compare the spatial distribution of void fraction. In particular the two-fluid model assessment was done with Experiment III data for jL ≤ 2 m/s, which was the range of interest for the present study.
The CFD calculations were performed with CFX-11. Geometric modeling and meshing was done using the meshing software GAMBIT. Following the CFX guidelines a one-node, 2-degree slice of the pipe was adopted as the solution domain as shown in Fig. 1. A hexahedral mesh with 30 uniform radial elements and 150 uniform axial elements was selected after performing a mesh sensitivity study. The mesh is shown in Fig. 1.
The two turbulent diffusion force models described in the previous section were compared. The effect of the turbulent diffusion models is shown in Fig. 2. The two models available in CFX correspond to the models discussed in the section describing the turbulent diffusion force of the two-fluid model. The model of Lopez de Bertodano et. al. (2006) (i.e., Case 1) is similar to the
The results with these two models did not differ much at jL = 1 m/s. A typical result at 1 m/s is shown in Fig. 2. However, Fig. 2 shows that at 2 m/s the
Calculations were performed for two sets of wall force parameters, viz. Antal et al. (1991) (
Figures 3 to 5 show the prediction obtained using the two wall force models. A general trend noticed here is that the averaged experimental void fraction is higher than the computed results for all cases and the deviation is more for jL = 1.5 m/s and 2 m/s cases. This has a significant impact and it is the main reason why data is consistently higher than the calculated results in the center of the pipe. One explanation for this systematic discrepancy could be that the volumetric gas flow rates were measured at the inlet pressure instead of atmospheric pressure, and in such case correcting the computations for the effect of pressure would give better results. However, this has not been done here.
It can be seen from Figures 3-5 that the effect of the wall force coefficients is significant both in the prediction of the wall peak and in the void distribution in the far wall region. In general, the wall peak data is bounded by Antal et al.’s coefficients from below and Krepper and Prasser’s coefficients from above.
Fig. 3 shows comparisons for jL = 1.0 m/s. At this lower liquid flux, the void distribution is primarily controlled by the lift and the wall forces. The higher void fraction case (Fig. 3) shows that there is a wall peak in the near wall region and relatively uniform void distribution in the far wall region which implies that the turbulent diffusion force has a
negligible effect. It is clear from the comparisons that Krepper and Prasser’s coefficients underestimate the void fraction in the region away from wall.
For jL = 1.5 m/s, Krepper and Prasser’s coefficients predict the wall peaks accurately (Fig. 4) but the void fraction in the far wall region is under predicted. For the jG = 0.083 m/s, the far wall region void fraction is well predicted with Antal et al’s coefficients.
For jL = 2 m/s, the data shows that the void peak moves further away from the wall than that in jL = 1.5 m/s. This effect cannot be predicted by the present model without making some significant changes to the lift coefficient and it is an unresolved problem of two-phase flow CFD as of now. The void peak magnitude and bulk void fraction is better predicted by Antal et al.’s coefficients (Fig. 5).
Thus, Antal et al.’s coefficients predict the wall peak magnitude and the trend of the void distribution in the far-wall region quite well. Furthermore, the far-wall results indicate that the turbulent diffusion force model is satisfactory. As a result of this study, following models were selected for the further analysis using One-Group IATE.
Turbulence Diffusion Force:
Wall Force Model Coefficients: Antal et al’s (1991) coefficients. (
5. Two fluid model energy equation
The interface to liquid heat transfer (
In the above equation,
Reynolds number Re
The interfacial mass transfer is given as,
Wall heat transfer is the most important aspect of boiling as it provides the sources for energy and phase mass balance equations. Most wall heat transfer models are based on partitioning the wall heat flux into three components. In particular, the multidimensional CFD model by Kurul and Podowski (1990) has been implemented in the CFX-12 beta version (Egorov & Menter, 2004) and some of the details of this implementation are reported in the next section.
The three components of the wall heat flux are:
Evaporation heat flux (
Forced convection to the liquid (
Wall quenching by liquid phase (transient conduction) (
The wall surface is assumed to be split into two parts (
The evaporation heat flux is obtained from the rate of vapor generation at the wall which is given as a product of mass of a bubble at detachment, the detachment (departure) frequency and the nucleation site density,
The specification of evaporation heat flux using above equation requires closure for bubble departure diameter, bubble departure frequency and the site density. The default model in CFX uses following closure relations:
The bubble departure diameter is given by the empirical correlation of Tolubinsky and Kostanchuk (1970), the Nucleation Site Density is obtained using the correlations of Lemmart and Chawla (1977) and the departure frequency is closed using the simplest available expression (Cole, 1960), which uses a characteristic bubble velocity (terminal velocity of bubble rise) and a characteristic bubble size (departure diameter).
The frequency of bubble departure obtained by Cole (1960) was derived under the assumption of bubble moving with terminal velocity once detached from the surface. This is a fair assumption for pool boiling scenario at low heat fluxes where bubble detachment is hydrodynamically governed, but may not be accurate for flow boiling at high heat fluxes where the thermodynamics governs the bubble growth and detachment. The closure relations are given below
An alternative to the Lemmart and Chawla (1977) correlation for the nucleation site density is the correlation of Kocamustaffaogullari and Ishii (1983) validated against water data for 0-50 bar. The Kocamustaffaogullari-Ishii (K-I) correlation relates the site density to the wall superheat as follows:
Under the following conditions:
The density dependent parameter in equation (29) is given as,
The effective wall superheat in equation (32) is usually less than the actual wall superheat. The reason for this is as follows. A bubble nucleated at the wall grows through a liquid film adjacent to the wall where a considerably high temperature gradient exists. So, in reality, it experiences a lower mean superheat than the wall superheat. In case of pool boiling, the difference is not significant and the superheat based on the wall temperature can be taken as the effective superheat. In case of forced convective boiling, there exists a steeper temperature gradient in the near wall field due to thermal boundary layer. Hence, the effective superheat experienced by the bubble is smaller than the actual wall superheat and it is given as,
The two-phase Reynolds number is given as,
The Martinelli parameter can be approximated as,
Thus, equations (29)-(40) represent the complete model of Kocamustafaogullari-Ishii (1983), which has been implemented in CFX-12 for the present study.
As mentioned previously, the quenching heat flux is the component of wall heat flux utilized to heat the cold liquid replacing the detached bubble adjacent to the heated wall. In order to evaluate this component, Mikic and Rohsenow (1969) used an analytical approach starting with transient conduction in semi-infinite medium with heated wall being the only boundary where temperature is specified.
The above equation has a solution of the form,
The heat flux at the wall boundary can be obtained as,
The above heat flux is averaged over a “characteristic” time period which is approximated by the inverse of bubble departure frequency. This leads to the original closure for quenching heat flux [Mikic & Rohsenow, 1969].
Kurul and Podowski (1990) made minor modification to the above expression with an assumption that quenching occurs between detachment of one bubble and appearance of next bubble (nucleation), and this time was assumed to be 80% of the detachment period. This is the final form of quenching heat transfer closure of the CFX-12 model.
This is evaluated under the standard assumption of a logarithmic temperature profile across the turbulent boundary layer (Kader (1981), ANSYS CFX Solver Theory Guide (2006)).
where, Δy is the distance of the wall adjacent node from the wall,
The convective heat flux component through the liquid area fraction is thus given as,
It should be noted here that all the wall heat flux model correlations above use a liquid temperature (
6. Boiling CFD Model assessment with high pressure R12 data
Morel et al. (2003) performed boiling experiments in a vertical pipe of 19.2 mm internal diameter and a length of 5 m. The pipe had a heated section of 3.5 m preceded by and followed by unheated lengths of 1 m and 0.5 m respectively. The input parameters for the four tests reported were as summarized in Table 1. The liquid-vapor density ratio in these experiments corresponds to that for water-steam at 95-150 bar.
Radial Profiles of Void Fraction, Liquid Temperature and Interfacial Area Concentration were measured. Dual sensor fiber optics probes were used to measure void fraction and interfacial area concentration. Miniature thermocouples were used to measure the fluid temperature.
The CFD calculations were performed with CFX-12. Geometric modeling and meshing was done using the meshing software GAMBIT. A 2 wedge of the pipe is used as the domain. The mesh had 20 radial and 250 axial uniformly spaced elements. Fig. 6 shows the domain and cross-section of the mesh used.
|Parameter||Test1 (Deb5)||Test2 (Deb6)||Test3 (Deb13)||Test4 (Deb10)|
|Mass Flux (kg/m 2 )||1986||1984.9||2980.9||2027.0|
|Inlet Subcooling ( o C)||18.12||16.11||18.12||23.24|
|Heat flux (W/m 2 )||73.89||73.89||109.42||76.24|
All the data sets include radial profiles of bubble diameter. Note that with Freon at these test pressures, the surface tension is 0.0017 N/m for 26 bar which results in very small bubble sizes (i.e., 0.5 mm). Instead of using bubble size correlations for the CFD calculations the mean of the measured radial profile value is used for the bubble diameter in the bulk, whereas, the value measured at the wall is used as the departure bubble size (Fig. 7).
It was found during all the cases that a better match with experimental data is obtained with the coefficient of the turbulent diffusion force
Fig. 8 shows the prediction of radial variation of void fraction obtained using two different site density models. The Lemmart-Chawla (1977) model which is used in CFX was found to under-estimate the void fraction in all the cases. Kocamustafaogullari-Ishii model (termed as
K-I Model in the figures below) predicted the data well in all the cases. The prediction however deteriorated at low pressure data (P = 14.59 bar) and this needs further investigation.
Fig. 9 shows the prediction of radial profiles of interfacial area concentration (IAC). The IAC is obtained using the expression for the Sauter mean diameter:
As the bubble diameter is assumed to be uniform, the deviation of the CFD values from the experimental data near the wall, where the bubble diameter is smaller, is justified.
Fig. 10 shows prediction of liquid temperature for the first two cases (for which data is available). The simulation results show a good match with the experimental data in the bulk. However, towards the wall, the results with Lemmart-Chawla correlation showed steeper rise in temperature which caused over-prediction in this region.
Fig. 11 shows the bubble diameter prediction using the CFX default bubble departure diameter model with the Kurul and Podowski (1990) correlation for the bubble size in the bulk of the flow. The bubble departure diameter and the bulk diameter predicted by the models are approximately three times higher than the actual value, leading to a very high rate of vapor generation. This results in over-prediction of void fraction as shown in Fig. 11.
7. Subcooled boiling CFD model assessment with atmospheric water dataExperimental data and CFD model
Lee et al. (2002) performed subcooled boiling experiments in a vertical annulus having a 19 mm diameter heating rod placed inside a tube with internal diameter of 37.5 mm. The total length of the test section was 2.37 m. The central heated length of 1.67 m was preceded and followed by unheated portions of 0.188 m and 0.518 m respectively as shown in Fig. 12.
Yeoh and Tu (2005) have reported two cases for this experiment where the bubble diameter measurements have been done. These two cases were selected for the analysis here. The input parameters for the two tests reported are summarized in Table 2 below.
|Mass Flux (kg/m 2 )||474||1059.2|
|Inlet Subcooling ( o C)||13.11||17.82|
|Heat flux (kW/m 2 )||152.9||251.5|
Radial profiles of void fraction, vapor velocity and liquid velocity were measured at 1.61 m downstream of the start of the heated section. Two-point conductivity probe was used to measure void fraction and bubble velocity. Pitot tube was used for liquid velocity measurements with a correction factor for two phase flow.
The CFD calculations were performed with CFX-12. Geometric modeling and meshing was done using the meshing software GAMBIT. A 2 wedge of the annulus is used as the domain. The mesh had 20 radial and 142 axial uniformly spaced elements. Fig. 12 shows the domain and cross-section of the mesh used.
The bubble diameter measurements for the two cases are shown below in Fig. 13. It must be noted here that although a variation in the bubble size has been measured along the radial direction, however, the bubble layer is one bubble thick and the bubble size is of the order of 4-5 mm which occupies almost half of the annulus gap. This means the departure size can be assumed the same as the bubble size, which is equal to the maximum size of bubble measured across the cross-section. Hence, a bubble size (bulk and departure) of 4.7 mm is used in the low heat flux case and that of 4.1 mm is used in high heat flux case.
The site density model comparison was done for these two cases also and the void fraction comparison is shown below in Fig. 14. It can be seen that the Lemmart-Chawla correlation over-predicted the void fraction at the wall by a significant margin. Although results using K-I model did not accurately match the void fraction distribution, the wall void fraction was consistently better predicted.
It is to be noted here that the K-I model did a better job even at the high pressure data and hence appears to be better than the Lemmart-Chawla model.
As stated earlier, in this problem, the bubble size is comparable to the annulus gap width and hence further modifications were made to the two-fluid model to improve the results:
Addition of a wall drag force: The bubble moves touching the wall and hence experiences the shear due to the liquid velocity gradient near the wall. In such cases the bubble is further slowed by an additional wall drag which is essentially an additional shear stress on the bubble which is moving in the liquid shear layer. A model for this wall drag force has been derived by Lahey et al. (1993) and is given as follows:
The value suggested for
Exclusion of Turbulent Dispersion: The turbulent dispersion force accounts for transverse dispersion of bubbles due to interaction with turbulent eddies. However, in the present case where there exists only one bubble across the cross-section, such a diffusion model is invalid and hence this force is turned off.
Fig. 15 shows the void fraction prediction for the two cases where the wall drag is included. It can be noticed that the wall drag effect is more in the high mass flux case. This is due to the fact that wall drag is comparable to the interfacial drag in this case.
Fig. 16 shows the void fraction predictions after turning off turbulent dispersion (
Fig. 17 shows the vapor velocity comparison with two simulations, viz., the original two-fluid model and with wall drag and no turbulent dispersion. The low flux case showed no significant difference in the three cases. A deviation near the outer wall in the no TD case can be neglected since this region has no vapor. The vapor velocity did not change much even in the high flux case but showed a slight reduction when wall drag was added.
With this benchmark the two-fluid model with the wall-boiling module was validated over a wide range of pressures (~1-150 water-steam equivalent bar). It is worth noting again that the low pressure case presented above has a large bubble size (i.e., half the annulus size) which makes the two-fluid model constitutive relations invalid in a strict sense. Despite this shortcoming, the results are acceptable.
Several CFD Two Fluid Model constitutive relations were assessed with experimental data over a wide range of conditions for adiabatic and boiling flows to provide a guideline to choose the constitutive models.
For the air-water adiabatic flow data Auton’s (1987) formulation of the lift force was found appropriate for the prediction of the void fraction profiles. The turbulent diffusion force given by the Lopez de Bertodano et al. (1994a) produced satisfactory results, especially at higher liquid fluxes. Two sets of coefficients for Antal et. al.’s (1991) wall force were compared and it was found that Antal et al.’s coefficients produced more consistent predictions and with the combination of selected lift and turbulent diffusion models, it produced overall better predictions throughout the pipe cross-section. The interfacial momentum models selected using adiabatic data were then used for the boiling flow comparisons with satisfactory results.
For the high pressure R12 boiling flow data the nucleation site density model by Kocamustaffaogullari and Ishii (1983) and the bubble departure frequency model of Cole (1960) were found to predict the data well. For the atmospheric pressure subcooled boiling water data set, where the bubble size is about half the flow channel width and where the thickness of the bubble layer is one bubble diameter, use of a wall drag force and elimination of the turbulent diffusion force improved the predictions.
The present model does not include constitutive relations for the bubble size at the wall nor in the bulk. These remain as open problems to be resolved in the future.
Antal S. Lahey R. Flaherty J. 1991Analysis of phase distribution in fully-developed laminar bubbly two phase flow.
ANSYS CFX Solver Theory Guide 2006ANSYS CFX Release 11.0.
Auton T. R. 1987The Lift Force on a Spherical Body in a Rotational Flow.
Chen J. C. 1966Correlation for boiling heat transfer to saturated fluids in convective flow.
Cole R. 1960A Photographic Study of Pool Boiling In the Region of the Critical Heat Flux.
Colin C. Legendre D. Yoshikawa H. Montout M. 2009Hydrodynamics of bubble detachment in convective boiling,
Egorov Y. Menter F. 2004Experimental Implementation of the RPI Wall Boiling Model in CFX-11. Ansys Cfx, Staudenfeldweg 12, 83624 Otterfing, Germany.
Ishii M. Hibiki T. 2006
Ishii M. Zuber N. 1979Drag coefficient and relative velocity in bubbly, droplet or particulate flows.
Kader B. A. 1981Temperature and concentration profiles in fully turbulent boundary layers.
Kocamustafaogullari G. Ishii M. 1983Interfacial area and Nucleation Site Density in boiling systems.
Krepper E. Prasser-M H. 2000Measurements and CFD Simulations of bubbly flow in a vertical pipe,
Kurul N. Podowski M. Z. 1990Multidimensional effects in forced convection subcooled boiling,
Lahey R. T. Jr Bertodano M. Jones O. C. Jr 1993Phase distribution in complex geometry conduits.
Lee T. H. Park G. C. Lee D. J. 2002Local flow characteristics of subcooled boiling flow of water in a vertical concentric annulus.
Lemmert M. Chawla J. M. 1977Influence of flow velocity on surface boiling heat transfer coefficient, In: Hahne, E., Grigull, U. (Eds.),
Lopez de Bertodano, M.A.,
Lopez de Bertodano. M. A. Lahey R. T. Jr Jones O. C. 1994aPhase distribution in bubbly two-phase flow in vertical ducts.
Lopez de Bertodano. M. A. Lahey R. T. Jr Jones O. C. 1994bDevelopment of a k-epsilon model for bubbly two-phase flow.
Lopez de Bertodano. M. 1998Two fluid model for two-phase turbulent jets.
Lopez de Bertodano. M. Ishii M. Sun X. Ulke A. 2006Phase distribution in the cap bubble regime in a duct.
& Mikic B.B. Rohsenow W.M. 1969A new correlation of pool-boiling data including the fact of heating surface characteristics.
Moraga F. J. Larreteguy A. E. Drew D. A. Lahey R. T. Jr 2006A center-averaged two-fluid model for wall-bounded bubbly flows.
Morel C. Yao W. Bestion D. 2003Three dimensional modeling of boiling flow for the NEPTUNE code.
Mostafa A. A. Mongia H. C. 1988On the interactions of particles and turbulent fluid flow.
Prabhudharwadkar D.M.; Bailey C.A.; & Lopez de Bertodano M.A. Buchanan J.R. 2009Two-Fluid CFD model of adiabatic air-water upward bubbly flow through a vertical pipe with a One-Group Interfacial Area Transport Equation,
Ranz W. E. Marshall W. R. 1952Evaporation from drops.
Reeks M. W. 1991On a kinetic equation for the transport of particles in turbulent flows.
Sato Y. Sadatomi M. Sekoguchi K. 1981Momentum and heat Transfer in Two-phase Bubble Flow-II.
Serizawa A. Kataoka I. Michiyoshi I. 1986Phase Distribution in Bubbly Flow, Data Set 24
Tolubinsky V. I. Kostanchuk D. M. 1970Vapor bubbles growth rate and heat transfer intensity at subcooled water boiling,
Tomiyama A. Tamai H. Zun I. Hosokawa S. 2002Transverse migration of single bubbles in simple shear flows.
Yeoh G. H. Tu J. Y. 2005A unified model considering force balances for departing vapour bubbles and population balance in subcooled boiling flow.