Open access peer-reviewed chapter

CFD Two Fluid Model for Adiabatic and Boiling Bubbly Flows in Ducts

Written By

Martin Lopez de Bertodano and Deoras Prabhudharwadkar

Published: January 1st, 2010

DOI: 10.5772/7097

Chapter metrics overview

4,983 Chapter Downloads

View Full Metrics

1. Introduction

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.


2. Nomenclature


3. Two fluid model mass and momentum equations

Mass Conservation

The ensemble averaged two-fluid model continuity equations (Ishii & Hibiki, 2006) governing the motion of nfields, k= 1 to n, has the following form:


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.

Momentum conservation


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.

Drag force

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υ_¯R=υ_¯2υ_¯1

Lift force

The lift force is given by (Auton, 1987):


The lift coefficient used for the calculations is CL= 0.1 (Lopez de Bertodano, 1992). This value is of the same order as the experimental value (Tomiyama et al., 2002) for a bubble in Couette flow (CL= 0.288).

Wall force

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):


wheren_andywallare the normal inward vector and the distance from the wall respectively. The values for the coefficients that are used in the CFD code CFX are:cw1= -0.01 andcw2= 0.05. Note that the ratio|cw2/cw1|determines the thickness of a layer next to the wall (in terms of the bubble diameter) in which the wall force is applicable.

Turbulent Diffusion force

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:


where λeis the Eulerian length scale of the eddies. The effective time constant is obtained combining equations (9) and (10) as:


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:

Case 1 -Forτcτd,τc1τc2it can be shown that equation (7) is equivalent to the turbulent diffusion force model obtained by Lopez de Bertodano et. al. (2006):


whereν2is the turbulent diffusivity of the bubbles. This equation was used by Lopez de Bertodano (2006) to model the diffusion of small bubbles by the very large vortices generated by cap bubbles in cap-bubbly flow. This model scales likeυ'(υ)and is adequate when the turbulence is bubble induced. However for shear induced turbulence a diffusion model with stronger flow dependence is needed.

Case 2 -Forτc1τc2it can be shown that equation (7) is equivalent to the model of Lopez de Bertodano et al. (1994a) for the diffusion of bubbles in turbulent duct flow, which is implemented in CFX as:


where the turbulent diffusion coefficient used for the present calculations is CTD= 0.25. This model scales likekυ'2(υ2)which provides more diffusion at higher flow velocities in the case of shear induced turbulence.

Turbulence Transport

The closure for the Reynolds stresses in equation (2) is based on the two-phase k-εmodel developed by Lopez de Bertodano et al. (1994b). This model assumes that the shear induced (SI) and bubble induced (BI) turbulent stresses are added together:


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 k-ε model for the shear induced diffusivity and the second term corresponds to Sato et al.’s (1981) model for the bubble induced diffusivity. The coefficient Cμ= 0.09 is the standard value for the k-ε model. However, the coefficient CDBcan be adjusted to fit the velocity profile data. A value of 0.6 is used for CDBin the present simulations. Standard k-ε model transport equations are used for the continuous phase. For the dispersed phase, the kinematic eddy viscosity is assumed to be equal to that of continuous phase, and hence no turbulence transport equations are needed.


4. Adiabatic CFD model assessment with atmospheric air-water data

Experimental data and CFD model

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.

Figure 1.

The Computational Domain and Typical CFX Mesh

Turbulent Diffusion Force Model Comparison

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 Favre Averaged Drag Modelof CFX and the Lopez de Bertodano Modelof CFX is the model discussed in Case 2.

Figure 2.

Turbulent Diffusion Model Comparison for jL = 1 m/s and 2 m/s

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 kdependent model (i.e., Case 2) with CTD= 0.25 is a significant improvement over the Favre type model (i.e., Case 1). The previously discussed inadequacy of the Case 1 model to produce sufficient diffusion for shear induced turbulence is evident in the figure and therefore the Lopez de Bertodano Modelof CFX was selected for all further simulations.

Wall Force Model Comparison

Calculations were performed for two sets of wall force parameters, viz. Antal et al. (1991) (Cw1= -0.01, Cw2= 0.05) and Krepper and Prasser (2000) (Cw1= -0.0064, Cw2= 0.016). The uncertainties in the wall force are very significant and these two sets of parameters are used as the two bounds. The bubble size used for the calculations is Db = 5 mm., which is bigger than the reported value of 3 mm. This is done so because the bubbles are oblate and the measurements correspond to the smaller diameter, whereas the distance from the wall is determined by the larger diameter.

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

Figure 3.

Wall Force Model Comparison (jL = 1 m/s)

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.

Figure 4.

Wall Force Model Comparison (jL = 1.5 m/s)

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).

Figure 5.

Wall Force Model Comparison (jL = 2 m/s)

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: Lopez de Bertodano Modelof CFX (CTD= 0.25)

Wall Force Model Coefficients: Antal et al’s (1991) coefficients. (cw1= -0.01, cw2= 0.05)

Lift Force: CL= 0.1.


5. Two fluid model energy equation

Energy conservation


where, λkis the thermal conductivity and νtkis the turbulent viscosity for phase k. Pr tis the turbulent Prandtl number, Hkiis the mean enthalpy of phase knear interfaces, and qkiis the interfacial heat flux whose closure is described in detail later. In the present calculations of subcooled boiling, the vapor generated was assumed to be at saturated temperature and hence was treated as isothermal phase.

Interfacial heat transfer

The interface to liquid heat transfer (k=l) is expressed as,


where, hliis the heat transfer coefficient between the liquid and the interface, closed as follows,


In the above equation, λlis the liquid thermal conductivity and dsis a length scale which is assumed equal to the bubble diameter. The Nusselt number (Nu) in the above closure equation is given by following expression (Ranz & Marshall, 1952):


Reynolds number Re bused in the above closures is defined as,


The interfacial mass transfer is given as,


where, Hlgis the latent heat of phase change.

Wall heat transfer

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 (qe) – The part of heat flux utilized in formation of vapor at the wall,

Forced convection to the liquid (qc),

Wall quenching by liquid phase (transient conduction) (qq) – This accounts for the heat transfer to the subcooled liquid that replaces the detached bubble at the wall.

The CFD Wall Boiling Model

The wall surface is assumed to be split into two parts (A1, A2) each under the influence of one phase. Fraction A2 is influenced by the vapour bubbles formed on the wall and participates in the evaporation and quenching heat transfer. Fraction A1 is the remaining part of the wall surface, (A1=1A2) and participates in the convective heat transfer to the liquid. A1 and A2 are related to the nucleation site density per unit wall area (nsite) and to the influence area of a single bubble forming at the wall nucleation site. The CFX-12 model assumes, that the diameter of the bubble influence zone is twice as big as the bubble departure diameter (a= 2).


Evaporation Heat Flux:

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

Ddep=Drefexp((TsatTliq)/Tref),Dref=0.6 mm,  Tref= 45 KE28
nsite=nref ((TwallTsat)/ΔTref)1.805,nref=7.9384×105 m2,  ΔTref= 10 KE29

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:




Rcis the critical radius of the surface cavity which represents a minimum cavity size which can be activated at a given wall superheat. Rcis given as,


Under the following conditions:


Rccan be simplified as,


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 multiplier Sin the above equation is the superheat suppression factor and is given as (Chen, 1966),


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.

Quenching Heat Flux:

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.


Single Phase Convection Heat Flux:

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, klis the liquid turbulence kinetic energy. In the above equation βis a function of liquid Prandtl number and Γ is a function of y*(ANSYS CFX Solver Theory Guide (2006)).

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 (Tliq). The single phase heat transfer correlation above uses the wall adjacent node temperature. However, in the bubble departure diameter and the quenching heat transfer correlations which were based on one-dimensional models, the liquid temperature refers to the bulk mean temperature. As a good estimate, CFX approximates this temperature with the temperature at a fixed y* (250). This is done in order to have a mesh size independent evaluation of wall heat flux partitions.


6. Boiling CFD Model assessment with high pressure R12 data

Experimental data and CFD model

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.

ParameterTest1 (Deb5)Test2 (Deb6)Test3 (Deb13)Test4 (Deb10)
Mass Flux (kg/m 2 )19861984.92980.92027.0
Pressure (bar)26.1526.1526.1714.59
Inlet Subcooling ( o C)18.1216.1118.1223.24
Heat flux (W/m 2 )73.8973.89109.4276.24

Table 1.

Simulation input parameters for R12 tests

Figure 6.

Domain and Mesh for R12 boiling simulation case


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 CTD= 0.5 instead of the value obtained for adiabatic air water flows, CTD= 0.25. This change may be attributed to bubble diameters in these cases which are one order of magnitude smaller than that in atmospheric air-water systems (Prabhudharwadkar et al., 2009). Therefore there are more interactions of the bubbles with smaller turbulent eddies. Two comparisons are presented below for this particular benchmark problem: the nucleation site density model comparison and the bubble size model comparison.

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

Figure 7.

Bubble diameter approximation for the simulations

Figure 8.

Void Fraction prediction using the two Nucleation Site Density models

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.

Figure 9.

Interfacial Area Concentration prediction using the two Nucleation Site Density models

Figure 10.

Liquid Temperature prediction using the two Nucleation Site Density models

Figure 11.

Bubble diameter prediction using Kurul-Podowski model and void fraction prediction.


7. Subcooled boiling CFD model assessment with atmospheric water data

Experimental 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.

Figure 12.

Domain and Mesh for water boiling simulation case

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 )4741059.2
Pressure (bar)1.421.43
Inlet Subcooling ( o C)13.1117.82
Heat flux (kW/m 2 )152.9251.5

Table 2.

Simulation input parameters for water boiling tests

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.

Figure 13.

Bubble Diameter data for the water boiling tests

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.

Figure 14.

Void Fraction prediction using different site density models (Water Data)

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.

Modifications of the Two-Fluid 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 CDWis 1.

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 (MTD) force. The peak void fraction increased substantially in both cases but it also resulted in slight narrowing of the vapor boundary layer near the boiling wall. Overall, the two modifications to the two-fluid model improved the results. However, these computations based on a model where the void distribution is the result of the balance between the lift and wall forces are not correct in a strict sense because the thickness of the bubble layer is determined instead by the diameter of one bubble. This fundamental problem has been addressed by Moraga et al. (2006).

Figure 15.

Void Fraction predictions after including the wall drag term

Figure 16.

Void Fraction predictions after excluding 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.

Figure 17.

Vapor velocity predictions without and with two-fluid model modifications


8. Conclusions

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.


  1. 1. AntalS.LaheyR.FlahertyJ.1991Analysis of phase distribution in fully-developed laminar bubbly two phase flow.International Journal of Multiphase Flow,175635652,0022-1120
  2. 2. ANSYS CFX Solver Theory Guide2006ANSYS CFX Release 11.0.
  3. 3. AutonT. R.1987The Lift Force on a Spherical Body in a Rotational Flow.Journal of Fluid Mechanics,183199213,1469-7645
  4. 4. ChenJ. C.1966Correlation for boiling heat transfer to saturated fluids in convective flow.Industrial and Engineering Chemistry Process Design and Development,53322329.0196-4305
  5. 5. ColeR.1960A Photographic Study of Pool Boiling In the Region of the Critical Heat Flux.AIChE Journal,64533534,0001-1541
  6. 6. ColinC.LegendreD.YoshikawaH.MontoutM.2009Hydrodynamics of bubble detachment in convective boiling,Proceedings of the 7th ECI International Conference on Boiling Heat Transfer (BOILING 2009), Florianopolis-SC-Brazil.
  7. 7. EgorovY.MenterF.2004Experimental Implementation of the RPI Wall Boiling Model in CFX-11. Ansys Cfx, Staudenfeldweg 12, 83624 Otterfing, Germany.
  8. 8. IshiiM.HibikiT.2006Thermo-fluid Dynamics of Two-phase Flow, Springer,0-38728-321-8States.
  9. 9. IshiiM.ZuberN.1979Drag coefficient and relative velocity in bubbly, droplet or particulate flows.AIChE Journal,255843855,0001-1541
  10. 10. KaderB. A.1981Temperature and concentration profiles in fully turbulent boundary layers.International Journal of Heat and Mass Transfer,24915411544.
  11. 11. KocamustafaogullariG.IshiiM.1983Interfacial area and Nucleation Site Density in boiling systems.International Journal of Heat and Mass Transfer,26913771397,0017-9310
  12. 12. KrepperE.Prasser-MH.2000Measurements and CFD Simulations of bubbly flow in a vertical pipe,AMIFESF Workshop, Computing Methods in Two-Phase Flow,18.
  13. 13. KurulN.PodowskiM. Z.1990Multidimensional effects in forced convection subcooled boiling,Proceedings of the Ninth International Heat Transfer Conference,2126, Jerusalem, Israel, Hemisphere Publishing Company.
  14. 14. LaheyR. T.JrBertodanoM.JonesO. C.Jr1993Phase distribution in complex geometry conduits.Nuclear Engineering and Design,1412177201,0029-5493
  15. 15. LeeT. H.ParkG. C.LeeD. J.2002Local flow characteristics of subcooled boiling flow of water in a vertical concentric annulus.International Journal of Multiphase Flow,28813511368,0301-9322
  16. 16. LemmertM.ChawlaJ. M.1977Influence of flow velocity on surface boiling heat transfer coefficient, In: Hahne, E., Grigull, U. (Eds.),Heat Transfer in Boiling,237247, Hemisphere Publishing Corporation,0-12314-450-7
  17. 17. Lopez de Bertodano,M.A.,Turbulent bubbly two-phase flow in a triangular duct, PhD Thesis, Rensselaer Polytechnic Institute, Troy, NY, U.S.A, 1992.
  18. 18. Lopez deBertodano. M. A.LaheyR. T.JrJonesO. C.1994aPhase distribution in bubbly two-phase flow in vertical ducts.International Journal of Multiphase Flow,205805818,0301-9322
  19. 19. Lopez deBertodano. M. A.LaheyR. T.JrJonesO. C.1994bDevelopment of a k-epsilon model for bubbly two-phase flow.Journal of Fluids Engineering,1161128134,0098-2202
  20. 20. Lopez deBertodano. M.1998Two fluid model for two-phase turbulent jets.Nuclear Engineering and Design,17916574,0029-5493
  21. 21. Lopez deBertodano. M.IshiiM.SunX.UlkeA.2006Phase distribution in the cap bubble regime in a duct.Journal of Fluids Engineering,1284811818.
  22. 22. MikicB.B.&RohsenowW.M.1969A new correlation of pool-boiling data including the fact of heating surface characteristics.ASME Journal of Heat Transfer,91245250,0022-1481
  23. 23. MoragaF. J.LarreteguyA. E.DrewD. A.LaheyR. T.Jr2006A center-averaged two-fluid model for wall-bounded bubbly flows.Computers & Fluids,354429461,0045-7930
  24. 24. MorelC.YaoW.BestionD.2003Three dimensional modeling of boiling flow for the NEPTUNE code.Proceedings of the 10th International Topical Meeting on Nuclear Reactor Thermal Hydraulics, Seoul, Korea.
  25. 25. MostafaA. A.MongiaH. C.1988On the interactions of particles and turbulent fluid flow.Int. J. Heat Mass Transfer,311020632075,0017-9310
  26. 26. PrabhudharwadkarD.M.;BaileyC.A.;Lopez de BertodanoM.A.&BuchananJ.R.2009Two-Fluid CFD model of adiabatic air-water upward bubbly flow through a vertical pipe with a One-Group Interfacial Area Transport Equation,Proceedings of the ASME 2009 Fluids Engineering Division Summer Meeting, Vail, Colorado, USA.
  27. 27. RanzW. E.MarshallW. R.1952Evaporation from drops.Chemical Engineering Progress,483141146,0360-7275
  28. 28. ReeksM. W.1991On a kinetic equation for the transport of particles in turbulent flows.Phys. Fluids A,33446456,1070-6631
  29. 29. SatoY.SadatomiM.SekoguchiK.1981Momentum and heat Transfer in Two-phase Bubble Flow-II.International Journal of Multiphase Flow,72179190,0301-9322
  30. 30. SerizawaA.KataokaI.MichiyoshiI.1986Phase Distribution in Bubbly Flow, Data Set24Proceedings of the Second International Workshop on Two-Phase FlowFundamentals,Troy, NY, USA.
  31. 31. TolubinskyV. I.KostanchukD. M.1970Vapor bubbles growth rate and heat transfer intensity at subcooled water boiling,Proceedings of the 4th International Heat Transfer Conference,5Paris, France. (PaperB-2.8).
  32. 32. TomiyamaA.TamaiH.ZunI.HosokawaS.2002Transverse migration of single bubbles in simple shear flows.Chemical Engineering Science,571118491858,0009-2509
  33. 33. YeohG. H.TuJ. Y.2005A unified model considering force balances for departing vapour bubbles and population balance in subcooled boiling flow.Nuclear Engineering and Design,2351012511265,0029-5493

Written By

Martin Lopez de Bertodano and Deoras Prabhudharwadkar

Published: January 1st, 2010