The importance of correct modeling of

bubble size and condensation in prediction of sub-cooled boiling flows

S. Lo*1, A. Splawski1 and B.J. Yun2

1CD-adapco, Didcot, Oxfordshire, UK 2Korea Atomic Energy Research Institute, Daejeon, Korea

Received: 9 January 2012, Accepted: 28 June 2012

Abstract

This paper describes the updating of the sub-cooled boiling model used in CFD codes with the more recent and better sub-models. The improved sub-models include: Hibiki and Ishii [1] correlation for nucleation site density, Kocamustafaogullari [2] correlation for bubble departure diameter and the S-gamma model of Lo and Rao [3] for bubble size distribution in the flow. The new model has been tested against measured data from Debora [4] and Bartolomei [5]. The results show that improvement in the bubble size prediction has the most significant impact on the accuracy of the model.

1. INTRODUCTION

Modeling of sub-cooled boiling flows using three dimensional CFD methods can be considered fairly well established as most commercial CFD codes offer such modeling capability. This paper describes the updating of the sub-cooled boiling model used by the author with the more recent and better sub-models. The improved sub-models include: Hibiki and Ishii [1] correlation for nucleation site density, Kocamustafaogullari [2] correlation for bubble departure diameter and the S-gamma model of Lo and Rao [3] for bubble size distribution in the flow. It is important to note that both the Hibiki-Ishii and Kocamustafaogullari correlations include experimental data from atmospheric pressure to high pressure flows. Hibiki-Ishii model ranges up to 198 bar and Kocamustafaogullai model up to 142 bar. These models are therefore suitable for analyzing flows at PWR operating pressure.

The new model has been tested against measured data from Debora experiments [4] and five sets of measured data from Bartolomei [5] for sub-cooled steam-water flows at 147 bar, close to PWR operating pressure. The results show that improvement in the bubble size prediction has the most significant impact on the accuracy of the model. The S-gamma model solves a transport equation for the bubble interfacial area density. It takes the bubble departure diameter and generation rate as boundary conditions and models the bubble size change due to evaporation, condensation, breakup and coalescence.

It has been observed that bubble size is a function of pressure. At higher pressures the vapor bubbles generated at the heated wall and moving in the bulk flow are very small. It is therefore important to take the bubble departure diameter as the boundary value in the calculations. The bubble size predicted by S-gamma is significantly smaller than the simple linear model of Kurul and Podowski [6] which is a function of liquid temperature only and has no dependency on pressure. With smaller bubbles in the S-gamma results, the interfacial area available for heat transfer is increased. The steam condensation rate in the sub-cooled area is greatly enhanced and that has led to much better match of the axial void profiles against the experimental data for all five Bartolomei cases studied.

2. MATHEMATICAL MODEL

The standard "two-fluid model" was used in modeling the boiling two-phase flows considered in this paper. In this model the conservation equations for mass, momentum and energy are solved for both phases.

'Corresponding Author: E-mail address: simon.lo@uk.cd-adapco.com

The conservation of mass for phase is:

where ak is volume fraction of phase k, pk is phase density, uk is phase velocity, m ki and m ik are mass transfer rates to and from the phase, and N is the total number of phases. The sum of the volume fractions is clearly equal to unity.

Xa = i (2)

The conservation of momentum for phase k is:

dt (pu )+V-K PUU ) -v-(a (t+Tk )) = -ak vp+a PkS+M (3)

where Tk and T'k are laminar and turbulence shear stresses, p is pressure and M is the sum of the interfacial forces that include drag, turbulent dispersion, virtual mass and lift forces. A full description of these forces can be found in Chapter 5 of the book by Azzopardi et al. [13]. These interfacial forces are described by a rather standard set of formulae and usually available as standard models in commercial CFD software, for example, see description in Chapter 13 of the STAR-CD software Methodology Manual [14].

The conservation of energy for phase k is:

^KpA )+V-{akPkuA )-V

at VT + V\ )

where hk is phase enthalpy, Xk is thermal conductivity, T is temperature, ¡it is turbulent viscosity, ah is turbulent Prandtl number and Q is the interfacial heat transfer and other heat sources. Since in sub-cooled boiling flows the vapour is always at saturation temperature, the vapour phase can therefore be treated as isothermal at saturation temperature without solving the energy equation. Further details of the "two-fluid model" used can be found in Lo [7].

2.1 Wall Boiling Model

At the heated wall, boiling occurs when the wall temperature exceeds the saturation temperature. The steam generation rate is determined by the wall heat partitioning model as follows,

qw = q+%+(5)

where, qw is the total heat flux from the wall, ql is the single phase convection heat flux that takes place outside the influence area of nucleation bubbles, qq is the quenching heat flux within the bubble influence area and qe is the evaporation heat flux. The convective heat flux is calculated using the standard wall function in the CFD code. The quenching heat flux is modeled as transient heat conduction in a semi-infinite slab according to Del Valle and Kenning [12].

The bubble influence area Ae is defined by,

nd2 / f\

A = F—d-N (6)

where, FA ,dd, N" are model constant, bubble departure size and active nucleation site density, respectively. FA = 2 is used in this study.

The evaporation heat flux is clearly an important parameter for an accurate prediction of the sub-cooled boiling flows and it can be expressed as,

q p h fN" (7)

e g g fg

where, pg is the steam density, hfg is the latent heat and f is the bubble departure frequency.

In most commercial CFD codes, Tolubinsky [8] bubble departure size model, Kurul & Podowski [6] active nucleation site density model and Cole's [9] bubble departure frequency model are adopted in their wall boiling models. A summary of these models can be found in Lo [7] and Tentner el al. [11]. These models have very simple forms and they do not reflect properly their dependency on the flow, pressure and fluid properties. In this paper, we have replaced these models with better ones.

2.2 Active nucleation site density

For active nucleation site density, Hibiki et al.'s [1] model was adopted. Characteristics of the Hibiki et al. model is that (a) it considers a boundary condition for a wall superheating, and (b) it is validated against a great number of experimental data. It has a wide range of applicability in terms of the mass flow, pressure and contact angle. Hibiki et al.'s correlation is expressed as follows,

exp \f(p+)^-\-l

where, 6 is the contact angle, /! = 0.722, Nn = 4.72 x 105 sites/m2, X = 2.50 x 10-6 m. Rc is a critical cavity radius given as,

R _ 2a{l + (pg / pf)}/ P (9)

c exp{hfg(T -T )/RTTa}-1

where, a is the surface tension, Pf is the liquid pressure and R is the gas constant based on the molecular weight of fluid.

The f (p+) in equation (8) is a function to consider the pressure effect on the active nucleation site density and given as,

f (p+) _ -0.01064 + 0.4824p+ - 0.22712p+2 + 0.05468p+3 (10)

where, p+ = log (p*), p* = Ap / pg

Hibiki et al.'s model is applicable in the range of 0.0 to 886 kg/m2s for mass flux, 0.101 to 19.8 MPa for pressure, 5 to 90o for contact angle and 1x104 to 1.51x1010 sites/m2 for active nucleation site density. For the calculation of Rc in equation (9), superheated liquid temperature Tg near the heated wall is required. However, this temperature is not available in the CFD calculation, Tg is therefore assumed to be the surface temperature at the heated wall, Tw in the present work.

2.3 Bubble departure size

For bubble departure diameter the empirical correlation of Kocamustafaogullari [2] is used.

d = 2.64 x\059

f Y-s f * ^ - wAp

An important point about this correlation is that it is based on steam-water data and a wide pressure range, from 0.067 to 141.87 bar.

2.4 Bubble size distribution

For the prediction of bubble size in the flow, the S Y (S-gamma) model was used [15]. This model represents the bubble size distribution using the Method of Moment approach. The bubble size distribution is assumed to have a pre-defined shape and this shape is retained during the process under investigation. A log-normal distribution is used in this study. The parameter S Y is a conserved quantity on the volumetric basis and is related to the moment M Y by:

Sy = nMY = n J dY P( d )d (d) (12)

From equation (12) we can obtain the bubble number density n (=S0) from the zeroth-moment; the interfacial area density at (=nS2) from the second-moment; and the void fraction a (=nS3/6) from the third-moment. The Sauter mean diameter, d32, can be obtained from:

d = s=a. (13)

32 S n S

A transport equation is solved for S Y as follow:

^^—L + V.(Pl3uS) = prn (5 +5 +5 +5 ) (14)

dt V' g g Y' ' g \ br cl mass boil)

On the right hand side of equation (14) are the sources for breakup, coalescence, mass transfer due to evaporation/condensation and wall boiling. Full description of the breakup and coalescence sources can be found in Lo [3], so will not be repeated here. The sources related to boiling mass transfers in the bulk fluid (mmass) and at the wall (mbail) can be written as follow:

s = ■

Y mm Ï aPe

6 mboi

•d Y

In equation (16) dw is the bubble departure diameter obtained from the bubble departure size model given in equation (11), hence linking the wall boiling model to the bubble size distribution model.

The Sauter mean diameter computed by equation (13) is used to compute quantities such as forces and the second-moment S2 is used for the interfacial area density ai in the heat and mass transfer terms.

Equation (14) for the second-moment is in fact equivalent to the Interfacial Area Transport (IAT) equation reported in the literature. The differences between the two models are in the source terms for breakup and coalescence. Note that several S Y equations can be solved together for different moments such that the model can be extended to more complex distribution functions. Also note that density variation is fully incorporated in the equation.

3. RESULTS

The model described above was tested against a large number of cases with experimental data in a systematic study of the model performance. The study includes grid sensitivity study, turbulence

model, wall treatment, and various aspects of the wall boiling model. These details are reported in Yun et al. [10] so will not be repeated here.

3.1 Debora test case

One test case from Yun et al. [10] is described here for discussion. The test case chosen, DEB10, is one of the DEBORA experiment reported by Yao and Morel [4]. The computed results using the model described above are labeled "S-gamma". For comparison the same calculations were repeated using the correlation from Kurul and Podowski [6] giving bubble size as a function of liquid sub-cooling. The results are labeled "Linear".

Figure 1 shows the results from both methods are in reasonable agreement with the experimental data for the radial distribution of void fraction at the measuring plane near the pipe exit. The advantage of the SY method is shown in Figure 2a. Since it includes the bubble departure diameter at the heated wall into account, it can capture the near wall profile better. The curves show there is substantial bubble growth due to coalescence and evaporation in the bulk fluid. The interfacial area density is better predicted by the SY method as shown in Figure 2b.

Figure 1 shows an over-prediction of void in the near wall region and under-prediction in the pipe centre region. With a small amount of void in the centre region the bubble size is correspondingly smaller as shown in Figure 2a. These results suggest that the radial migration of the bubbles as represented by the lift force is not strong enough. With adjustment of various model parameters described in [10], better results can be obtained as shown in Figure 3.

The significance of this work is that the Symodel provides a mechanistic approach to model the bubble size distribution in boiling two-phase flows accounting for breakup, coalescence, evaporation, condensation and vapor generation at the boiling wall. This model has shown to give better results for boiling two-phase flows.

Figure 1. Radial distribution of void fraction.

Figure 2a. Radial distribution of Sauter mean diameter.

Figure 2b. Radial distribution of interfacial area density.

Radial Position (m) Radial Position (m) Radial Position (m)

Figure 3. Radial distributions for DEB10 with improved model. Extracted from Yun et al [10]

3.2 Bartolomei test cases

Bartolemei et al. [5] carried out a systematic experimental study of steam-water flows in vertical pipes investigating the effect of pipe diameter (D), pressure (P), water flow rate (G), wall heat flux (Q) and inlet sub-cooling (AT). Five test cases conducted at 147 bar were selected for analyses in this paper, since this pressure is close to the PWR operating pressure that we are interested in. The operating conditions of these five cases are tabulated in Table 1. In these experiments the axial void profiles were measured and the results were plotted again thermodynamic quality as shown in Figure 4. Note that Figure 4 was extracted from the paper by Bartolemei [5] published in 1980. The paper was written in Russian, it was difficult to determine the accuracy of the measured data (i.e. no error bar). The vertical axis in Figure 4 is void fraction and the horizontal axis is thermodynamic quality as a measure of axial distance. The computed axial void profiles for the five test cases are shown in Figures 5, 6 and 7. The figures show that the "Linear" model consistently over-predicted the void in the pipe whereas the "S-gamma" results agree well with the measured data.

Table 1. Experimental conditions of Bartolemei cases

Cases P (bar) G (kg/m2s) Q (MW/m2) AT (K)

Bart22 □ 147.9 1878 0.42 16.43

Bart23 A 147.4 1874 0.77 27.47

Bart24 • 147.5 2123 1.13 48.59

Bart25 x 147.0 2014 1.72 63.38

Bart26 0 149.9 2012 2.21 144.51

A. Mfla piqxr (M'-O « Sr. K ?

□ ff¡7S /878 0,fZ 603

A /847 0,77 598

• m,7S 2/23 /,/3 f63

x WO W/t 7,72 S4S

o /4,39 2/7/2 ¿.2/ ÍS3 aa £

• • L X • A A

O noo O- x 0 X • • A ! A a

axoxcP x x x » ® * S- X • • ; , »V tft> u

— » •

-o,z -o,/ f} o ¿car„ q?

Figure 4. Experimental results of Bartolemei cases.

-0.15 -0.1 -0.05 0 0.05 -0.2 -0.1 0 0.1

Thermodynamic Quality Thermodynamic Quality

Figure 5. Comparison of results for Cases Bart22 and Bart23.

-0.2 -0.1 0 -0.3 -0.2 -0.1 0

Thermodynamic Quality Thermodynamic Quality

Figure 6. Comparison of results for Cases Bart24 and Bart25.

-0.2 -0.1 Thermodynamic Quality Figure 7. Comparison of results for Case Bart26.

Contour plots of the computed results are provided in Figures 8, 9 and 10 to illustrate the differences between the "Linear" and "S-gamma" results. Figure 8 shows that the bubble diameter computed by the "Linear" model is much larger than the "S-gamma" model. The "Linear" model calculates the bubble diameter based on the local water temperature and takes no consideration of the bubble departure diameter at the wall nor the time and conditions required for the bubbles to grow. Since the water temperature near the wall is high the model returns a larger bubble diameter there. The "S-gamma" model takes the bubble departure diameter as the boundary value and

computes the bubble size change due to coalescence, breakup, evaporation and condensation. Figure 8 shows that the bubble size computed by the "S-gamma" model is much smaller. With smaller bubbles, the interfacial area for heat transfer increases and the condensation rate is higher as shown in Figure 9. The net effect of higher condensation rate is the smaller void in the flow shown in Figure 10.

Figure 8. Bubble diameter: Left (Linear), Right (S-gamma).

Figure 9. Condensation rate: Left (Linear), Right (S-gamma).

CD-adapco

pro-STAR 4.14

5-Nov-l 0 Volume Fraction

PHASE г ITER = 1D41 LOCAL MX- 0,5071 LOCAL MN= 0.1Z15E-11

0.5071

0.4709

0.4347

0.3934

0.3622

0.3260

0.2808

0.2536

0.2173

0.1811

0.1449

0.1087

0.7244E-01

0.3622E-01

0.1215E-11

CD-adapco

pro-STAR 4.14

5-NOV-10 Volume Fraction

PHASE г ITER = 1010 LOCAL MX= 0.4402 LOCAL MN= 0.100DE-11

0.4402

0.4087

0.3773

0.3453

0.3144

0.2330

0.2515

0.2201

0.1386

0.1572

0.1 258

0.9432E-01

0.6288E-01

0.3144E-01

0.1 000E-11

Figure 10. Void distribution: Left (Linear), Right (S-gamma).

4. CONCLUSIONS

In this paper we have described the sub-cooled boiling model with the improved sub-models of Hibiki and Ishii [1] correlation for nucleation site density, Kocamustafaogullari [2] correlation for bubble departure diameter and the S-gamma model of Lo and Rao [3] for bubble size distribution in the flow. We have found that the bubble size model has the strongest effect on the accuracy of the computed results. The bubble diameter computed by the "Linear" model is too large which gives rise to over-prediction of the void in the flow. The bubble diameter computed by the "S-gamma" model is smaller leading to higher interfacial area and stronger condensation. The computed void is smaller and much closer to the measured data. This paper therefore demonstrated the importance of correct modeling of bubble size and condensation in prediction of sub-cooled boiling flows.

REFERENCES

[1] T. Hibiki, M. Ishii, Active Nucleation Site Density in Boiling Systems, Int. J. Heat Mass Transfer, Vol. 46, pp. 2587-2601 (2003).

[2] G. Kocamustafaogullari, Pressure dependence of bubble departure diameter for water, Int. Comm. Heat Mass Transfer, Vol. 10, pp. 501-509, 1983.

[3] S. Lo and P. Rao, Modelling of droplet breakup and coalescence in an oil-water pipeline, 6th International Conference on Multiphase Flows, ICMF 2007, Leipzig, Germany, July 9-13, 2007.

[4] W. Yao and C. Morel, Volumetric interfacial area prediction in upward bubbly two-phase flow, International Journal of Heat and Mass Transfer, Vol. 47, pp. 307-328 (2004).

[5] G.G. Bartolemei, G.N. Batashova, V.G. Brantov, et.al. In: Teplomassoobmen-IV. Vol.5. Minsk. ITMO AN BSSR Press, 1980, p.38.

[6] N. Kurul and M.Z. Podowski, Multidimensional effects in sub-cooled boiling, Proc. of 9th Heat Transfer Conference, Jerusalem, 1990.

[7] S. Lo, Modelling Multiphase Flow with an Eulerian Approach, VKI Lecture Series - Industrial Two-Phase Flow CFD von Karman Institute, May 23-27 (2005).

[8] V.I. Tolubinsky, D.M. Kostanchuk, Vapour Bubbles Growth Rate and Heat Transfer Intensity at Subcooled Water Boiling Heat Transfer, 4th International Heat Transfer Conference, Paris, vol. 5, Paper No. B-2.8 (1970).

[9] R. Cole, A photographic study of pool boiling in the region of the critical heat flux, AIChE J. Vol.6, 533-542 (1960).

[10] B.J. Yun, A. Splawski, S. Lo and C.H. Song, Prediction of a subcooled boiling flow with mechanistic wall boiling and bubble size models, CFD4NRS-3 Workshop, September 14-16, 2010, Washington DC, USA.

[11] A. Tentner, S. Lo, A. Ioilev, M. Samigulin, V. Ustinenko, Computational fluid dynamics modelling of two-phase flow in a boiling water reactor fuel assembly, Mathematics and Computation, Supercomputing, Reactor Physics and Nuclear and Biological Applications, Avignon, France, 12-15 September 2005.

[12] M.V.H. Del Valle and D.B.R. Kenning (1985), Subcooled flow boiling at high heat flux, Int. J. Heat Mass Transfer, Vol. 28, pp. 1907-1920.

[13] B.J. Azzopardi, R.F. Muddle, S. Lo, H. Morvan, Y.Y. Yan and D. Zhao, Hydrodynamics of Gas-Liquid Reactors -Normal Operation and Upset Conditions, Wiley, 2011.

[14] STAR-CD Methodology Manual, CD-adapco.

[15] S. Lo and D. S. Zhang (2009), Modeling of Break-up and Coalescence in Bubbly Two-Phase Flows, J. Comp. Multiphase Flows, Vol. 1, pp. 23-38.