Population balance modelling for bubbly flows with heat and mass transfer

Chemical Engineering Science 59 (2004) 3125 – 3139

www.elsevier.com/locate/ces

Population balance modelling for bubbly &ows with heat and mass transfer G.H. Yeoha;∗ , J.Y. Tub a Australian b School

Nuclear Science and Technology Organisation (ANSTO), PMB 1, Menai NSW 2234, Australia of Aerospace, Mechanical and Manufacturing Engineering, RMIT University, Vic. 3083, Australia Received 1 June 2003; received in revised form 1 February 2004; accepted 1 April 2004

Abstract Population balance equations combined with a three-dimensional two-&uid model are employed to predict bubbly &ows with the presence of heat and mass transfer processes. Subcooled boiling &ow belongs to this speci3c category of bubbly &ows is considered. The MUSIG (MUltiple-SIze-Group) model implemented in CFX4.4 is further developed to account for the wall nucleation and condensation in the subcooled boiling regime. Comparison of model predictions against local measurements near the test channel exit is made for the radial distribution of the bubble Sauter diameter, void fraction, interfacial area concentration and gas and liquid velocities covering a range of di:erent mass and heat &uxes and inlet subcooling temperatures. Additional comparison was also performed against existing boiling model in CFX4.4 and the modi3ed model developed in our previous work (Int. J. Heat Mass Transfer 45 (2002) 1197). Good agreement is better achieved with the local radial bubble Sauter diameter, void fraction, interfacial area concentration and liquid velocity pro3les against measurements using the newly formulated MUSIG boiling model over the simpler boiling models. However, signi3cant weakness of the model is still evidenced in the prediction of the vapour velocity. Work is in progress to circumvent the de3ciency of the model by the consideration of additional momentum equations or an algebraic slip model to account for bubble separation. ? 2004 Elsevier Ltd. All rights reserved. Keywords: Population balance; Two-&uid model; Bubble coalescence; Bubble break-up; Wall nucleation; Condensation; Subcooled boiling &ow

1. Introduction 1.1. Population balance approach Application of the population balance approach towards better describing and understanding complex industrial &ow systems has received an unprecedented attention and acceptance. A population balance of any system is concerned with maintaining a record for the number of entities, which for bubbly &ows are bubbles, or drops; whose presence or occurrence may dictate the behaviour of the system under consideration. In addition to the motion of these entities through the state space, it is usual to encounter “birth” processes that create new entities and “death” processes that destroy existing ones. The birth and death processes may depend on the states of the entities created or destroyed with an associated phenomenology; coalescence, breakage, etc. are examples of such processes. A population balance model is, ∗ Corresponding author. Tel.: +61-2-9717-3817; fax: +61-2-97179263. E-mail address: [email protected] (G.H. Yeoh).

0009-2509/$ - see front matter ? 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2004.04.023

therefore, formulated based on the collective phenomenology contained in the displacement of entities through their state space and the birth and death processes that terminate entities and produce new entities. Lately, mounting interest on population balances have resulted in a number of signi3cant developments especially towards modelling bubbly &ows. Bubble column reactors, widely used in many chemical, petroleum, mining, food and pharmaceutical industries, are known as excellent reactors for processes that require large interfacial area for gas–liquid mass transfer and eGcient mixing for competing gas–liquid reactions. The rate of transport of the gas to the liquid often limits productivity and is therefore a critical design criterion. The uncertainty in bubble column design arises from a lack of fundamental understanding of the local hydrodynamics and rate processes, which govern bubble size and thus the interfacial area between phases. With the advancement of computer technologies, the quest for improved designs has paved the trend biased towards the use of numerical simulations. Numerical methods are gradually gaining acceptance as a powerful tool

3126

G.H. Yeoh, J.Y. Tu / Chemical Engineering Science 59 (2004) 3125 – 3139

for design of chemical reactors. Several studies have been conducted using the computational &uid dynamics (CFD) methodology (Krishna et al., 1999; Shimizu et al., 2000; Pohorecki et al., 2001; Olmos et al., 2001). The use of CFD and population balance models has shown to expedite a more thorough understanding of di:erent &ow regimes and further enhance the description of the bubble characteristics in the column volume for design, especially with the consideration of bubble coalescence and break-up mechanisms in the model simulations. Recently, Ramkrishna and Mahoney (2002) have highlighted a promising future towards handling two-phase &ow systems using the population balance approach. Along similar developments, a transport equation for the interfacial area, analogous to the Boltzman equation has been considered by Kocamustafaogullari and Ishii (1995), Wu et al. (1998), Hibiki et al. (2001) and Hibiki and Ishii (2002) to handle two-phase turbulent bubbly &ows. This transport equation may be regarded as another form of a population balance equation addressing the distribution of bubbles. The development undertaken principally aims to address the many shortcomings found in thermal-hydraulic system codes for nuclear safety analysis caused by empirical correlations based on traditional two-phase &ow regimes and regime transition criteria. The one-group equation has been applied to predict the interfacial area in bubbly &ows. Hibiki and Ishii (2000) extended the one-group interfacial area transport equation to two-group transport equations thereby treating the bubbles in two groups, which are the spherical/distorted bubble group and the cap/slug bubble group. 1.2. Phenomenology consideration for bubbly :ows with heat and mass transfer Although considerable e:orts have been invested to develop more sophisticated models for bubble migration, attention of the transport processes is still very much focused on isothermal bubbly &ow problems. Such &ows greatly simplify the formulation of mathematical models where the heat and mass transfer processes can be safely neglected. However, Fleischer et al. (1996) investigated a process involving the chemisorption of CO2 into NaOH in a bubble column, where it was clearly demonstrated that mass transfer is an important process. Applying population balance equations with the consideration of mass transfer, their numerical predictions established that enhanced mass transfer resulted in large temporal and spatial variations of bubble sizes and of void fraction and interfacial area. There is therefore an increasing need for further development towards a more robust mathematical model capable of handling complex phenomena associated with hydrodynamics, heat and mass transfer, and bubbles coalescence and break-up. Subcooled boiling &ow belongs to a speci3c category of bubbly &ows, which embraces all the com-

plex dynamic interaction of the aforementioned phenomena. Here, heterogeneous bubble nucleation occurs within small pits and cavities on the heater surface, designated as nucleation sites. These nucleation sites are activated when the temperature of the surface exceeds the saturation temperature of the liquid at the local pressure. If the temperature of the bulk &uid remains below saturation at the same location, the boiling process is known as subcooled &ow boiling. Subcooled &ow boiling is thus characterised by a high-temperature two-phase region near the heated wall and a low-temperature single-phase liquid away from the heated surface as depicted in Fig. 1. Subcooled boiling starts at a point called the onset of nucleate boiling (ONB). It continues downstream from the ONB point until the void fraction begins to increase sharply at a location called the net vapour generation (NVG). The NVG point is the transition between two regions: low void fraction region followed by a second region, in which the void fraction increases signi3cantly. In our comprehensive investigation on axial void fraction distribution in channels, good agreement of the boiling &ow model has been achieved against a wide range of experimental data (Tu and Yeoh, 2002; Yeoh and Tu, 2002). Improvements made to the boiling &ow model include modi3cations to inter-phase heat transfer, mean bubble diameter in the bulk liquid and wall heat partition model. Nevertheless, further investigations in Yeoh et al. (2002) against local radial measurements of Yun et al. (1997) and Lee et al. (2002) for a low-pressure subcooled boiling annular channel &ow revealed signi3cant weakness of the model, predominantly in the radial prediction of the bubble Sauter diameter, liquid and vapour velocities. The empirical correlation applied for our axial comparison exercise that determines the bubble size in the bulk subcooled liquid was derived to only predict the macroscopic consideration of the boiling phenomenon. It is therefore not surprising that numerical models that employed this sort of relationship are unable to adequately represent the complex mechanistic behaviours of bubble coalescence and condensation (microscopic in nature) as observed through experiments. Experimental observations by Lee et al. (2002) using high-speed photography (see Fig. 2) clearly depicted presence of large bubble sizes away from the heated wall. The vapour bubbles, relatively small when detached from the heated surface, were seen to increase in size due to bubble coalescence. However, the bubbles gradually decreased in size due to the increased condensation as they migrated towards the opposite end of the unheated wall of the annular channel. This phenomenon predominant in subcooled boiling &ows was further con3rmed by experimental observations of Gopinath et al. (2002) (see Fig. 3), which illustrates a bubble gradually being condensed in a subcooled liquid away from the heated surface. The absence of the bubble mechanistic behaviour, such as bubble coalescence clearly observed during experiments in the vicinity of the heated wall and the condensation on

G.H. Yeoh, J.Y. Tu / Chemical Engineering Science 59 (2004) 3125 – 3139

Single

3127

Subcooled Boiling

Void Fraction

Phase

Highly

Slightly

Subcooled

Subcooled

ONB

OVG

Height

q Uin

Fig. 1. Subcooled &ow boiling regions.

19 mm

Heated Inner Tube Measuring Plane

1670 mm long Heated Section

1610 mm

9.25 mm

Test Channel Inlet

Inlet Liquid Flow

Fig. 2. Schematic drawing of the test channel.

the unheated &ow region, signi3cantly compromised the model predictions especially the bubble Sauter diameter distribution. In the two-&uid CFD model, which is the most

commonly used macroscopic formulation of the thermo-&uid dynamics of the two-phase systems, the phasic interaction term appears in the 3eld equations. These terms represent the important contribution of the mass, momentum and energy transfers through the interface between the phases. An accurate determination of the bubble Sauter diameter is crucial as the bubble size in&uences the inter-phase heat and mass transfer through the interfacial area concentrations and momentum drag terms. Because of the successes employing the population balance approach especially the feasibility to accurately account for the merging and destruction characteristics of bubbles in &ow systems, the potential to implement and extend the modelling to predict the non-uniform bubble size distribution in subcooled boiling &ows is of enormous signi3cance. Such a capability does not exist in the current state-of-the-art. Therefore, a successful and speci3c development of the population balance approach for boiling &ows can contribute to a signi3cant improvement of the two-&uid boiling model formulation. The objectives of this present study are twofold: to develop a population balance approach using on the MUSIG (MUltiple-SIze-Group) model coupled with a complete three-dimensional &ow numerical simulation for subcooled boiling &ows at low pressures including the associated phenomenology of nucleation and condensation in the generic computer code CFX4.4; and to evaluate the new MUSIG boiling model through validation against experimental measurements. Comparisons of local quantities for a range of di:erent mass and heat &uxes and inlet subcoolings are performed against recent radial measurements of Yun et al. (1997) and Lee et al. (2002).

3128

G.H. Yeoh, J.Y. Tu / Chemical Engineering Science 59 (2004) 3125 – 3139 6

Continuity equation of vapour phase:

Measurement

@(g g fi ) * + ∇ · (g g u g fi ) = Si − fi Plg : @t

Sauter Diameter (mm)

MUSIG Boiling Model 5

Modified Boiling Model Default Boiling Model

4

Momentum equation of liquid phase: *

3

* * @(l l u l ) + ∇ · (l l u l u l ) @t

Qwall = 152.9 kW/m2 T = 96.6oC

2

*

inlet

1

= − l ∇P + l l g

Pinlet = 0.142 MPa G = 474.0 kg/m2s

*

0

0.2

(a)

0.4

0.6

0.8

1

6 MUSIG Boiling Model

(3)

*

@(g g u g ) * * + ∇ · (g g u g u g ) @t

Modified Boiling Model Default Boiling Model

4

*

Momentum equation of vapour phase:

Measurement

Sauter Diameter (mm)

*

+(Plg u g − Pgl u l ) + Flg :

(r-Ri)/(Ro-Ri)

5

*

+∇ · [l le (∇ u l + (∇ u l )T )]

0

*

= − g ∇P + g g g

3

Qwall = 197.2 kW/m2 T = 95.0oC

2

*

Pinlet = 0.137 MPa

*

0 0

0.2

*

+(Pgl u l − Plg u g ) + Fgl :

G = 714.4 kg/m2s 0.4

0.6

0.8

1

(r-Ri)/(Ro-Ri)

(b)

*

+∇ · [g ge (∇ u g + (∇ u g )T )]

inlet

1

Qwall = 251.5 kW/m T = 92.1oC

5

2

Measurement

@(g g Hg ) * + ∇ · (g g u g Hg ) @t

MUSIG Boiling Model Modified Boiling Model

=∇ · [g ge ∇Tg ] + (Pgl Hl − Plg Hg ):

Default Boiling Model

2

(5)

Energy equation of vapour phase:

G = 1059.2 kg/m2s

3

Energy equation of liquid phase:

=∇ · [l le ∇Tl ] + (Plg Hg − Pgl Hl ):

inlet

Pinlet = 0.143 MPa 4

(4)

* @(l l Hl ) + ∇ · (l l u l Hl ) @t

6

Sauter Diameter (mm)

(2)

(6)

In Eq. (1), Plg represents the mass transfer rate due to condensation in the bulk subcooled liquid expressed by

1

0 0

0.2

(c)

0.4

0.6

0.8

1

(r-Ri)/(Ro-Ri)

Fig. 3. Local mean radial pro3les of bubble Sauter diameter: (a) C1, (b) C2 and (c) C3.

2. Mathematical formulation 2.1. Flow equations The two-&uid model treating both the vapour and liquid phases as continua solves two sets of conservation equations governing mass, momentum and energy, which are written for each phase as: Continuity equation of liquid phase: * @(l l ) + ∇ · (l l u l ) = Plg : @t

(1)

Plg =

haif (Tsat − Tl ) ; hfg

(7)

where h is the inter-phase heat transfer coeGcient determined from Ranz and Marshall (1952) Nusselt number correlation and aif is the interfacial area per unit volume. The wall vapour generation rate is modelled in a mechanistic way derived by considering the total mass of bubbles detaching from the heated surface as Pgl =

Qe ; hfg + Cpl Tsub

(8)

where Qe is the heat transfer due to evaporation. This wall nucleation rate is accounted in Eq. (2) as a speci3ed boundary condition apportioned to the discrete bubble class based on the size of the bubble departure criteria on the heated surface. On the right-hand side of Eq. (2), Si is the

G.H. Yeoh, J.Y. Tu / Chemical Engineering Science 59 (2004) 3125 – 3139

additional source terms due to coalescence and break-up based on the formulation, which are described in the next section. The term fi Plg represents the mass transfer due to condensation redistributed for each of the discrete bubble classes. The gas void fraction along with the scalar fraction of each size group fi are related to the number density of the discrete bubble ith class ni (similarly to the jth class nj ) as g fi = ni vi . The size distribution of the dispersed phase is therefore de3ned by the scalar size group fi . The population balance equation for each of the discrete bubble classes ni is provided in the next section. Inter-phase transfer terms in the momentum and energy equations—Pkj and Fkj —denote the transfer terms from phase j to phase k. The mass transfer Plg is already given in Eq. (7) while the total interfacial force Flg considered in the present study includes the e:ects of: drag dispersion lift lubrication Flg = Flg + Flg + Flg + Flg :

(9)

The terms on the right-hand side of Eq. (9) are the drag force, lift force, wall lubrication force and turbulent dispersion force respectively. Detail descriptions of these forces can be found in Anglart and Nylund (1996) and Lahey and Drew (2001). Brie&y, Interphase momentum transfer between gas and liquid due to drag force is given by drag Flg =

* * * * 1 CD aif l | u g − u l |( u g − u l ) 8

and (10)

Lift force in terms of the slip velocity and the curl of the liquid phase velocity is described by *

*

and

lift lift Flg = −Fgl :

(11)

Wall lubrication force, which is in the normal direction away from the heated wall and decays with distance from the wall, is expressed by 
* * g  l ( u g − u l ) Ds * lubrication n Flg =− max 0; Cw1 + Cw2 Ds yw and

lubrication lubrication Fgl = −Flg :

(12)

Turbulence induced dispersion taken as a function of turbulent kinetic energy and gradient of the void fraction of the liquid yields in the form of: dispersion Flg = −CTD l $∇l dispersion dispersion Fgl = −Flg :

Recommended value for CTD according to Kurul and Podowski (1990) of 0.1 is used for the turbulent dispersion force. A two-equation $-% turbulence model is employed for the continuous liquid and dispersed vapour phases. The e:ective viscosity in the momentum and energy equations is taken as the sum of the molecular viscosity and turbulent viscosity. The turbulent viscosity is considered as the total of the shear-induced turbulent viscosity and Sato’s bubble-induced turbulent viscosity (Sato et al., 1981). The local bubble Sauter diameter based on the calculated values of the scalar fraction fi and discrete bubble sizes di can be deduced from: Ds =  i

1 : fi =di

(14)

2.2. Bubble mechanistic model According to Fleischer et al. (1996), the bubble size distribution is calculated with the following population balance equation: *

* * * @n(V; x ; t) + ∇ · ( u g n(V; x ; t)) = G(V; x ; t); @t

(15)

*

drag drag = −Fgl : Flg

lift Flg = g l CL (˜u g − u l ) × (∇ × u l )

3129

and (13)

The drag coeGcient CD in Eq. (10) has been correlated for several distinct Reynolds number regions for individual bubbles according to Ishii and Zuber (1979). The constant CL takes a value of 0.01 (Wang et al., 1987). The wall lubrication constants Cw1 and Cw2 as suggested by Antal et al. (1991) are −0:01 and 0.05, respectively.

where n(V; x ; t) is the bubble number density distribution per unit mixture and bubble volume, which is a function * for the spatial range x for a given time t and volume V . On * the right hand side, the term G(V; x ; t) contains the bubble source/sink rates per unit mixture volume due to the bubble interactions such as coalescence, break-up and phase change. For the case of nucleate boiling and condensation in a subcooled boiling &ow, the phase change term includes the rate of change of bubble population with speci3c volumes. Phenomenological models developed by Prince and Blanch (1990) and by Luo and Svendsen (1996) allowed detailed description of the mechanisms for coalescence and break-up of intermittent bubbles. To e:ectively employed the number density transport equation given in Eq. (15) coupled with the aforementioned phenomenological models of coalescence and break-up, the technique proposed by Kumar and Ramkrishna (1996) that allows the usage of variable N bubble size groups to reduce the numerical e:ort is adopted, viz.,    * @ni + ∇ · ( u g ni ) = Rj + (Rph )i ; (16) @t j 

i

where ( j Rj )i represents the net change in the number density distribution due to coalescence and break-up processes. This entails the use of a 3xed non-uniform volume distribution along a grid, which allows a large size range to be covered with a small number of sections and yet still o:ers good resolution. The discrete bubble class between bubble volumes vi and vi+1 is represented by the centre point of the

3130

G.H. Yeoh, J.Y. Tu / Chemical Engineering Science 59 (2004) 3125 – 3139

grid interval. Such discretisation of the population balance equation has been found to allow accurate determination of desired characteristics of  the distribution (Ramkrishna, 2000). The interaction term ( j Rj )i (=PC +PB −DC −DB ) contains the source rates of PC , PB , DC and DB , which are, respectively, the production rates due to coalescence and break-up and the death rate to coalescence and break-up of bubbles formulated as: 2.2.1. Break-up PB =

N 

.(vj : vi )nj

j=i+1

DB = . i ni

with

.i =

i 

.ki

(17)

k=1

2.2.2. Coalescence PC =

1 2

i  i 

/i; kl ni nj

with

k=1 l=1

/i; kl = /kl if vk + vl = vi

else

/i; kl = 0 if vk + vl = vi ; DC =

N 

/ij ni nj :

(18)

j=1

The term (Rph )i in Eq. (16) comprises the essential formulation of the source/sink rates for the phase change processes associated with subcooled boiling &ow. At the heated surface, bubbles form at activated cavities known as active nucleation sites. The bubble nucleation rate from these sites can be expressed as 0WN =

N

f2H ; AC

(19)

where N

, f, 2H and AC are the active nucleation site density, the bubble departure frequency from the active sites, the heated perimeter and the cross-sectional area of the boiling channel, respectively. Since the bubble nucleation process only occurs at the heated surface, this heated wall nucleation rate is not included in (Rph )i but rather speci3ed as a boundary condition to Eq. (16) apportioned to the discrete bubble class ni based on the bubble departure criteria on the heated surface. The bubble sink rate due to condensation in a control volume for each bubble class can be determined from: 0COND = −

h(Tsat − Tl ) ni AB : VB g hfg

(20)

Given that the bubble surface area AB and volume VB based on the bubble Sauter diameter are, respectively, 4Ds2 and

4Ds3 =6, equation (20) can be rearranged as
  6g h(Tsat − Tl ) 1 (Rph )i = 0COND = − ni  g g D s hfg   haif (Tsat − Tl ) 1 ni : =− g  g hfg

(21)

The break-up of bubbles in turbulent dispersions employs the model developed by Luo and Svendsen (1996). Binary break-up of the bubbles is assumed and the model is based on the theories of isotropic turbulence. The break-up rate of bubbles of volume vj into volume sizes of vi can be obtained as 1=3 1 (1 + 2)2 % 2 211=3 dj 2min   12cf 6 d2; ×exp − 7l %2=3 dj5=3 211=3

.(vj : vi ) =C (1 − g )nj



(22)

where 2 = =dj is the size ratio between an eddy and a particle in the inertial sub-range and consequently 2min =min =dj ; and C and 7 are determined, respectively, from fundamental consideration of drops or bubbles breakage in turbulent dispersion systems to be 0.923 and 2.0. The variable cf de2=3 notes the increase coeGcient of surface area: cf = fBV + 2=3 (1−fBV ) −1, where fBV is the breakage volume fraction. The coalescence of two bubbles is assumed to occur in three steps. The 3rst step involves the bubbles colliding thereby trapping a small amount of liquid between them. This liquid 3lm then drains until it reaches a critical thickness and the last step features the rupturing of the liquid 3lm subsequently causing the bubbles to coalesce. The collisions between bubbles may be caused by turbulence, buoyancy and laminar shear. Only the 3rst cause of collision (turbulence) is considered in the present model. Indeed collisions caused by buoyancy cannot be taken into account here as all the bubbles from each class have been assumed to travel at the same speed. Moreover, calculations showed that laminar shear collisions are negligible because of the low super3cial gas velocities considered in this investigation. The coalescence rate considering turbulent collision taken from Prince and Blanch (1990) can be expressed as 
tij 4 /ij = [di + dj ]2 (uti2 + utj2 )0:5 exp − ; (23) 4 8ij where 8ij is the contact time for two bubbles given by (dij =2)2=3 =%1=3 and tij is the time required for two bubbles to coalesce having diameter di and dj estimated to be {(dij =2)3 l =166}0:5 ln(h0 =hf ). The equivalent diameter dij is calculated as suggested by Chesters and Ho:man (1982): dij = (2=di + 2=dj )−1 . According to Prince and Blanch (1990), for air–water systems, experiments have determined h0 , initial 3lm thickness and, hf , critical 3lm

G.H. Yeoh, J.Y. Tu / Chemical Engineering Science 59 (2004) 3125 – 3139 Table 1 Diameter of each discrete bubble class Class no.

Central class diameter di (mm)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

0.503 1.040 1.644 2.265 2.889 3.512 4.141 4.771 5.402 6.033 6.665 7.297 7.929 8.562 9.194

thickness at which rupture occurs to be 1 × 10−4 and 1 × 10−8 m, respectively. The turbulent velocity ut in the inertial subrange of isotropic turbulence (Rotta, 1972) is given by: ut = 1:4%1=3 d1=3 . The implementation of population balance equations originally developed by Lo (1996) for the generic computer code CFX4.4 is reformulated to account for the wall nucleation and condensation in the subcooled boiling &ow regime. To account for the non-uniform bubble size distribution in the bulk subcooled liquid, in this present study, bubbles ranging from 0 mm to 9:5 mm diameter are equally divided into 15 size groups (see Table 1). Instead of considering 16 different complete phases, it is assumed that each bubble class travels at the same mean algebraic velocity to reduce the computational time. This therefore results in 15 continuity equations for the gas phase coupled with a single continuity equation for the liquid phase. 2.3. Wall heat partition model Various experimental and theoretical investigations for low-pressure subcooled boiling &ow (Judd and Hwang, 1976) suggest that the wall heat &ux Qw can be divided into three components: heat transferred by conduction to the superheated layer next to the wall (nucleate boiling or surface quenching), Qq ; heat transferred by evaporation or vapour generation, Qe ; and heat transferred by turbulent convection, Qc . The surface quenching heat &ux is determined through the relationship:  

2

Qq = √ kl l Cpl f Aq (Tw − Tl ); (24) 4 where Tw is the wall temperature; Aq , the fraction of wall area subjected to cooling by quenching is calculated from Aq = N

(4d2bw =4); and f is the bubble departure frequency

given by Cole (1960):  4g(l − g ) : f= 3dbw l

3131

(25)

The bubble departure diameter, dbw , is formulated from considering the balance of surface tension and buoyancy forces at low pressures (Kocamustafaogullari and Ishii, 1995): 1=2
0:9
 l − g 6 dbw = 2:496 × 10−5 ; : (26) g gV The density of active nucleation sites, N

, is obtained from Lemmert and Chwala’s (1977) correlation of data, which is expressed by N

= [210(Tw − Tsat )]1:805 :

(27)

The bubble contact angle ; is taken to be at 55◦ for the present investigation as suggested in Hsu and Graham (1976) for most industrial metals and water. The heat &ux due to vapour generation at the wall in the nucleate boiling region can be simply calculated from Bowring (1962):

4 Qe = N

f (28) d3bw g hfg : 6 The heat &ux according to the de3nition of local Stanton number St for turbulent convection is Qc = St l Cpl ul (Tw − Tl )(1 − Aq ):

(29)

It is noted that ul is the local tangential liquid velocity adjacent to the heated surface. 3. Experimental details Experimental data of local subcooled boiling &ow measurements performed by our Korean collaborators that are used for the current validation exercise have been obtained from Yun et al. (1997) and Lee et al. (2002). The experimental set-up consists of a vertical concentric annulus with an inner heating rod of 19 mm outer diameter. The heated section is a 1:67 m long Inconel 625 tube with 1:5 mm wall thickness and is 3lled with magnesium oxide powder insulation. The rod is uniformly heated by a 54 kW DC power supply. The outer wall is comprised of two stainless-steel tubes with 37:5 mm inner diameter, which are connected by a transparent glass tube so that visual observation and photographic recording are made possible. The transparent glass tube is 50 mm long and is installed just below the measuring plane. The measuring plane is located at 1:61 m downstream of the beginning of the heated section. Demineralised water was used as the working &uid. The test channel inlet temperature was measured using the calibrated platinum resistance temperature detector with the estimated error of ±0:2◦ C. The absolute pressure at measuring plane was measured within the uncertainty of ±1 kPa. In this work, local gas phase

3132

G.H. Yeoh, J.Y. Tu / Chemical Engineering Science 59 (2004) 3125 – 3139

Table 2 Experimental conditions for C1, C2 and C3 Run

Pinlet (MPa)

Tinlet (◦ C)

Tsub (inlet) (◦ C)

Qw (kW=m2 )

G (kg=m2 s)

C1 C2 C3

0.142 0.137 0.143

96.6 94.9 92.1

13.4 13.8 17.9

152.3 197.2 251.5

474.0 714.4 1059.2

parameters such as local void fraction, bubble frequency and bubble velocity were measured by a two-conductivity probe method while the Pitot tube was used to measure the local measurement of liquid velocity with a mean relative error of 3:0%. However, the uncertainty of the bubble Sauter diameters (assuming spherical bubbles) determined through the Interfacial Area Concentration (IAC), calculated using the measured bubble velocity spectrum and bubble frequency, was diGcult to ascertain and will, at present, be estimated to be lower than 27%. More details regarding the experimental set-up can be found in Lee et al. (2002). Experimental conditions that have been used for comparison with the simulated results are presented in Table 2. Fig. 2 shows the schematic drawing of the test channel.

to generate the three-dimensional mesh within the annular channel resulting in a total of 13(radial) × 30(height) × 3(circumference) control volumes. Since wall function was used in the present study, the normal distance between the wall and the 3rst node in the bulk liquid should be such that the corresponding y+ was greater than 30. Grid independence was examined. In the mean parameters considered, further grid re3nement did not reveal signi3cant changes to the two-phase &ow parameters. Convergence was achieved within 1500 iterations when the mass residual dropped below 1×10−7 . Global execution time on the Silicon Graphics machine was about 30 min.

5. Results and discussion 4. Numerical details The discrete bubble sizes prescribed in the dispersed phase were tracked by solving an additional set of 15 transport equations, which these equations were progressively coupled with the &ow equations during the simulations. Sensitivity study on the number of size groups was performed through the consideration of equally dividing the bubble diameters into 10, 15 and 20 size groups. The analysis revealed that no appreciable di:erence was found for the predicted maximum bubble Sauter diameter between the 15 and 20 bubble size groups. For the subdivision into 10 size groups, the maximum bubble Sauter diameter was under predicted by a maximum di:erence of 2%. In view of the computational resources and times, it was therefore concluded that the subdivision of the bubbles sizes into 15 size groups were deemed suGcient and hereafter in the following the computational results are all based on the discretisation of 15 bubble size groups. Solution to the two sets of governing equations for the balance of mass, momentum and energy of each phase was sought. The conservation equations were discretised using the control volume technique. The velocity-pressure linkage was handled through the SIMPLE procedure. The discretised equations were solved using Stone’s Strongly Implicit Procedure (Stone, 1968). Since the wall heat &ux was applied uniformly throughout the inner wall of the annulus, advantage of the annular geometrical shape was utilised by modelling one quarter of the annulus as the domain for simulation. A body-3tted conformal system was employed

The radial pro3les of the bubble Sauter diameter, void fraction, interfacial area, vapour and liquid velocities located at 1:61 m downstream of the beginning of the heated section are predicted through the two-&uid and MUSIG boiling models. In all the 3gures presented, the dimensionless parameter (r − Ri )=(Ro − Ri ) = 1 indicates the inner surface of the unheated &ow channel wall while (r − Ri )=(Ro − Ri ) = 0 indicates the surface of the heating rod in the annulus channel. The MUSIG boiling model predictions against local measurements are further accompanied by computational results determined through empirical relationship of Anglart and Nylund (1996) to determine the local bubble diameter. They have proposed to estimate the interfacial transfer terms through a bubble diameter relationship assuming a linear dependence with local liquid subcoolings, which it can be expressed by d=

d1 (; − ;0 ) + d0 (;1 − ;) : ;1 − ;0

(30)

This relationship is still currently being used and applied in many boiling studies through the CFX4.4 code (numerical investigations in Lee et al., 2002; KonXcar et al., 2004). In this present study, the existing model will be known as the Default boiling model. The empirical bubble departure correlation in the wall heat partition model employs the relationship of Tolubinsky and Kostanchuk (1970) (in built into the Default boiling model): dbw = min[0:0006 exp(−;=45); 0:0014]:

(31)

G.H. Yeoh, J.Y. Tu / Chemical Engineering Science 59 (2004) 3125 – 3139

Reference diameters of d0 and d1 in equation (30) corresponding to the reference subcooling temperatures at ;0 and ;1 are usually not known a priori. We have assumed the local bubble diameters were evaluated between d0 = 1:5 × 10−4 and d1 = 7:0 × 10−3 m. We further assumed that both of the reference diameters corresponded to identical reference subcooling temperatures of ;0 = 13:0 K and ;1 = −5 K. In Tu and Yeoh (2002), improvements were made to the subcooled boiling &ow model to better predict the axial void fraction distribution. The mean bubble diameter relationship of Zeitoun and Shoukri (1997) was incorporated (here referred to as the Modi3ed boiling model): Ds



6=gV

=

0:0683(l =g )1:326

: 149:2( = )1:326 Re0:324 Ja + Bo0:487l Reg 1:6

(32)

For the bubble departure diameter, the correlation of Unal (1976) is adopted: dbw =

2:42 × 10−5 p0:709 a √ ; b?

where (qw − hl ;)1=3 kl

a= 2C 1=3 hfg 4kl =l cpl g

(33) 

kw w cpw ; kl l cpl

b = ;=2(1 − g =l ); C=

?=

hfg l [cpl =(0:013hfg Pr 1:7 ]3 ; [6=(l − g )g]0:5   u 0:47 l for ul ¿ 0:61 m=s; 0:61 1:0

for ul ¡ 0:61 m=s:

The stated range of this correlation is: Pressure Wall heat &ux Liquid velocity Liquid subcooling

0:1 ¡ p ¡ 17:7 MPa 0:47 ¡ qw ¡ 10:64 MW=m2 0:08 ¡ ul ¡ 9:15 m=s 3:0 ¡ ; ¡ 86◦ C

5.1. Local distributions of the bubble Sauter diameter, void fraction and interfacial area concentration Fig. 3 illustrates the local radial bubble Sauter diameter distribution at the measuring plane of the heated annular channel. In all the three cases, the empirical correlations from the Default and Modi3ed boiling models grossly misrepresented the local bubble sizes. The gradual increase of the bubble Sauter diameters towards the heated wall with the highest bubble sizes predicted at the heated wall by the empirical relationships contradicted the local radial measurements. In the experiments of Lee et al. (2002), high-speed photography clearly demonstrated large bubble sizes were present away from the heated wall not at the heated wall. However, this trend was correctly predicted by the MUSIG boiling model. Good agreement was achieved against the

3133

measured bubble sizes for all the three experimental conditions. The predicted bubble diameter behaviour determined through the empirical correlation was de3cient due to the absence of properly accommodating the bubble mechanistic behaviour coalescence and collapse due to condensation, which was succinctly observed in experiments. Evidently, the bubble size determination in the bulk liquid core was not strictly dependent on local subcoolings alone. This relationship was seen to signi3cantly compromise the model predictions. Extending the use of this bubble diameter correlation to predict local bubble sizes is invalid. It was also observed in Lee et al. (2002) that the vapour bubbles, relatively small when detached from the heated surface, increased in size due to bubble coalescence as they migrated towards the centre of the &ow channel. The bubble departure diameter evaluated from Eq. (26) resulted in a bubble size of approximately 1:9 mm. In all the three cases, a maximum predicted bubble size was respectively obtained for C1, C2 and C3 of about 4.5, 4.0 and 3:8 mm through the MUSIG boiling model con3rming the experimental observations. It was also interesting to note that coalescence of bubbles occurred axially along the heated surface. Experiments by Bonjour and Lallemand (2001) and Prodanovic et al. (2002) have clearly indicated the presence of bubbles sliding shortly after being detached from the heated cervices before lifting into the liquid core. These upstream bubbles travelling closely to the heated wall have the tendency of signi3cantly colliding with any detached bubbles downstream and subsequently forming bigger bubbles due to the bubbles merging together. Here, simulations have determined a bubble Sauter diameter of 3 mm corresponding to the adjacent points along the heated wall for all the three experimental conditions C1, C2 and C3. The bubble sizes being substantially larger than the bubble departure diameter have demonstrated to some degree the capability of the MUSIG boiling model to capture the coalescence behaviour of the bubbles sliding along the heated surface. As the bubbles migrated towards the opposite end of the adiabatic wall, they decreased due to the increased condensation. Here, only the low-temperature single-phase subcooled water existed. The bubble Sauter diameter pro3les of the MUSIG boiling model clearly showed the gradual collapse of the bubbles and the absence of bubbles near the adiabatic wall of the test channel. However, the bubble Sauter diameters were found to be under predicted when compared against the measured bubble sizes during the condensation process (see Figs. 3(a) and (b)). Gopinath et al. (2002) has demonstrated that the original form of the Ranz and Marshall Nusselt number correlation has a tendency to over predict the condensation Nusselt number, which increases the interfacial heat transfer coeGcient. It was observed through Gopinath et al. (2002) experiments that the thermal layer around the bubble thickened as the bubble shrunk in size and the net e:ect of the thickening of the boundary layer resulted in a decrease in the overall condensation rate. This phenomenon is evidently absent in the current model as

3134

G.H. Yeoh, J.Y. Tu / Chemical Engineering Science 59 (2004) 3125 – 3139

succinctly demonstrated through the predicted bubble sizes. Another important insight to the e:ect of condensation revealed that more bubbles were condensed with a higher inlet subcooling condition as shown in Fig. 3(c). With increasing mass &uxes, interfacial heat transfer was further enhanced thereby resulting in more bubbles being condensed in the subcooled liquid core. Fig. 4 presents the locally predicted void fraction pro3les against radial measured values. The peak local void fraction was always observed in the vicinity of the heated surface in a typical subcooled boiling &ow. This high local void fraction found here was explicitly due to the large number of bubbles generated from the active nucleation sites on the heated surface. Here, large amount of bubbles was generated from these nucleation sites when the temperature on the heated surface exceeded the saturation temperature. As these bubbles reached a critical size, they detached and migrated laterally toward the subcooled liquid core under the competing process of bubble coalescence and condensation as aforementioned. The Default and Modi3ed boiling models under predicted the void fraction pro3les for all the three experimental conditions; the most severe being case C3 where the channel boiling &ow was subjected to high inlet subcooling and mass &ux conditions (see Fig. 4(c)). The use of reference diameters in Eq. (30), albeit its simplicity formulation and application, failed to o:er any signi3cant bene3ts due to the ad hoc speci3cations. In some boiling problems, the setting of proper reference subcooling limits needs revision, which are also not known a priori. More importantly, extending the use of empirical relationships such as Eqs. (30) and (32) for other types of boiling &ow regimes may not be con3dently applied beyond the subcooled bubbly &ow regime where the two-phase &ow structures transit from bubbly to slug or churn turbulent boiling &ows, and other geometries. Nevertheless, the MUSIG boiling model (fundamentally derived from population balance principles) has the capacity of accommodating di:erent range of bubble sizes and mechanisms that may be present within the boiling liquid. It therefore presents enormous potential of possibly tracking the transition from one &ow regime to another and mechanistically predicting the bubble sizes associated for each of the boiling &ow regimes. This approach may well replace traditional &ow regime maps and regime transition criteria. For example, numerical studies of adiabatic bubbly &ows in bubble columns conducted through Olmos et al. (2001) have demonstrated the capability of the MUSIG model to predict the evolution of bubble sizes between two domains. In these two domains, the characteristics of the bubbles are typical of the homogeneous and transition regimes. Fig. 5 describes the local Interfacial Area Concentration (IAC) radial distribution. The IAC can be determined through the relationship: aif =

6 : Ds

(34)

The measured radial data followed the similar trend as the void fraction distribution in Fig. 5. Overall, better agreement between the measured and predicted IAC was achieved using the MUSIG boiling model. Here again, the Default and Modi3ed boiling models under predicted the IAC by an unacceptable margin. Based on a recent investigative study by Hibiki and Ishii (2003), they established the signi3cance of active wall nucleation site density linking to the prediction of the IAC. A generalised expression of the active nucleation site density on the heated surface covering a wide range of conditions is still not available and requires more analysis especially for large mass &uxes. One complexity arises in performing such experiments is in determining the population distribution of cavities, which they cannot be directly determined by measurement. They have to be indirectly inferred from experimental data. The local predictions of the void fraction and bubble Sauter diameter and consequently IAC immediately adjacent to the heated surface were determined through imposed boundary conditions at the heated surface. Therefore, they were greatly in&uenced by the density of the active wall nucleation sites and also by the bubble departure size at the wall, which explained the discrepancies seen in Fig. 5 for the Default and Modi3ed boiling models. Based on these results, it is clear that the applicability of Unal’s bubble departure relationship in the Modi3ed boiling model cannot be extended beyond the heat &ux range that it has been correlated. Work is in progress to adopt a more fundamental model through the formulation of appropriate force balances that are acting on the bubbles before departing from the heated surface under subcooled condition. This bubble departure model will accommodate a wider range of wall heat &uxes and &ow conditions in subcooled boiling &ows. 5.2. Local distributions of the vapour and liquid velocities The radial pro3les of the axial component of the local vapour velocity are shown in Fig. 6 while Fig. 7 presents the radial pro3les of the local liquid velocity for experimental conditions C1, C2 and C3. The vapour velocity was greater than the liquid velocity due to buoyancy force caused by density di:erence. As was observed in the experiment, the vapour velocity was higher at the centre than the velocities near the heating rod. This was probably due to the buoyancy e:ect being enhanced for the migration of the large bubbles there. However, the vapour velocity predicted by the MUSIG boiling model along with the Default and Modi3ed models showed higher velocities approaching the heated boundary. The MUSIG boiling model vapour velocities in the vicinity of the heated surface were rather similar to those of the simpler models for all the three cases because of the assumption that each bubble class travelled at the same mean algebraic velocity. Within the channel space, di:erent size bubbles are expected to travel with di:erent speeds. Additional momentum equations or an algebraic slip model could be

G.H. Yeoh, J.Y. Tu / Chemical Engineering Science 59 (2004) 3125 – 3139

3135

0.5

Qwall = 152.9 kW/m2 T = 96.6oC

0.4

inlet

Pinlet = 0.142 MPa

Void Fraction

G = 474.0 kg/m2s 0.3

Measurement MUSIG Boiling Model 0.2

Modified Boiling Model Default Boiling Model

0.1

0 0

0.2

0.4

0.6

0.8

1

(r-Ro)/(Ri-Ro)

(a) 0.5

Qwall = 197.2 kW/m2 T = 95.0oC

0.4

inlet

Void Fraction

Pinlet = 0.137 MPa G = 714.4 kg/m2s 0.3

Measurement 0.2

MUSIG Boiling Model Modified Boiling Model Default Boiling Model

0.1

0 0

0.2

0.4

0.6

0.8

1

(r-Ro)/(Ri-Ro)

(b) 0.5

Qwall = 251.5 kW/m2 T = 92.1oC

0.4

inlet

Void Fraction

Pinlet = 0143 MPa G = 1059.2 kg/m2s 0.3

Measurement MUSIG Boiling Model

0.2

Modified Boiling Model Default Boiling Model 0.1

0 0

(c)

0.2

0.4

0.6

0.8

(r-Ro)/(Ri-Ro) Fig. 4. Local mean radial pro3les of void fraction: (a) C1, (b) C2 and (c) C3.

1

3136

G.H. Yeoh, J.Y. Tu / Chemical Engineering Science 59 (2004) 3125 – 3139 1.4

Qwall = 152.9 kW/m2 Tinlet = 96.6 C 600

Pinlet = 0.142 MPa G = 474.0 kg/m2s Measurement

400

MUSIG Boiling Model Modified Boiling Model Default Boiling Model

200

0.2

0.4

0.6

0.8

Default Boiling Model

0.8 0.6

Qwall = 152.9 kW/m2 o Tinlet = 96.6 C Pinlet = 0.142 MPa G = 474.0 kg/m2s

0.4

0

1

0.2

0.4

800

0.6

0.8

1

(r-Ri)/(Ro-Ri)

(a)

(r-Ri)/(Ro-Ri)

(a)

1.8

Qwall = 197.2 kW/m2

Measurement 1.6

Tinlet = 95.0oC 600

Vapour Velocity (m/s)

Interfacial Area Concentration (1/m)

Modified Boiling Model 1

0 0

Pinlet = 0.137 MPa G = 714.4 kg/m2s

400

Measurement MUSIG Boiling Model Modified Boiling Model

200

Default Boiling Model

MUSIG Boiling Model Modified Boiling Model

1.4

Default Boiling Model 1.2 1 0.8

Qwall = 197.2 kW/m2

0.6

Tinlet = 95.5 C Pinlet = 0.137 MPa G = 714.4 kg/m2s

o

0.4 0.2

0 0

0.2

0.4

0.6

0.8

1

(r-Ri)/(Ro-Ri)

(b)

0

0.8

1

(r-Ri)/(Ro-Ri)

1.8

Vapour Velocity (m/s)

2

G = 1059.2 kg/m s Measurement MUSIG Boiling Model Modified Boiling Model Default Boiling Model

200

0.6

Measurement

o

Tinlet = 92.1 C Pinlet = 0.143 MPa

400

0.4

2

Qwall = 251.5 kW/m2 600

0.2

(b)

800

Interfacial Area Concentration (1/m)

MUSIG Boiling Model

0.2 0

MUSIG Boiling Model Modified Boiling Model

1.6

Default Boiling Model 1.4 1.2

Qwall = 251.5 kW/m2 o Tinlet = 92.1 C Pinlet = 0.143 MPa

1 0.8 0.6

0

G = 1059.2 kg/m2s

0.4 0

(c)

Measurement

1.2

o

Vapour Velocity (m/s)

Interfacial Area Concentration (1/m)

800

0.2

0.4

0.6

0.8

1

(r-Ri)/(Ro-Ri)

0

(c)

0.2

0.4

0.6

0.8

1

(r-Ri)/(Ro-Ri)

Fig. 5. Local mean radial pro3les of interfacial area concentration: (a) C1, (b) C2 and (c) C3.

Fig. 6. Local mean radial pro3les of vapour velocity: (a) C1, (b) C2 and (c) C3.

proposed to resolve the problem. Work is currently in progress to overcome this de3ciency of the two-&uid and MUSIG boiling models. The consideration of additional momentum equations to cater for each of the 15 bubble classes would increase the computational resources tremendously and deem impractical. Ongoing investigations are currently undertaken to test a pertinent choice of two or three dominant groups of bubbles transformed into the Eulerian phases to suGciently accommodate the hydrodynamics of wide bubble size distributed bubbly &ows. Also,

the process of developing of an algebraic slip model is in progress to account for the proper evaluation of slip velocities due to bubble separation. For the algebraic slip model, the terminal velocities for each of the bubbles can be determined possibly through an algebraic relationship suggested by Clift et al. (1978). Nevertheless, in Fig. 7, the MUSIG boiling model predictions of the liquid velocities were better represented when compared against the experimental measurements.

G.H. Yeoh, J.Y. Tu / Chemical Engineering Science 59 (2004) 3125 – 3139 1.4

Qwall = 152.9 kW/m2

Liquid Velocity (m/s)

1.2

Measurement MUSIG Boiling Model Modified Boiling Model Default Boiling Model

Tinlet = 96.6oC 1

Pinlet = 0.142 MPa G = 474.0 kg/m2s

0.8 0.6 0.4 0.2 0 0

0.2

0.4

0.6

0.8

1

(r-Ri)/(Ro-Ri)

(a) 1.8

Measurement MUSIG Boiling Model Modified Boiling Model Default Boiling Model

Liquid Velocity (m/s)

1.6 1.4 1.2 1 0.8

Qwall = 197.2 kW/m2 Tinlet = 95.0oC Pinlet = 0.137 MPa G = 714.4 kg/m2s

0.6 0.4 0.2 0

0.2

0.4

0.6

0.8

1

(r-Ri)/(Ro-Ri)

(b)

3137

sizes in the subcooled boiling &ow was distributed according to the division of 15 diameter groups through the formulation of a MUSIG model. Each of them experiencing coalescence and break-up phenomena has been considered. The MUSIG boiling model was developed to account for the wall nucleation or vapour generation on the heated surface and condensation process in the subcooled liquid core combined with the bubble coalescence of Prince and Blanch (1990) and bubble break-up of Luo and Svendsen (1997). Comparison of the predicted results was made against recent local measurements of Yun et al. (1997) and Lee et al. (2002). Additional comparison was also performed against the existing boiling model in CFX4.4 and the modi3ed boiling model developed in Tu and Yeoh (2002). Good agreement was achieved through the newly formulated MUSIG boiling model for the local bubble Sauter diameter, void fraction, IAC and liquid velocity pro3les. However, in the gas phase, since the assumption was invoked where each bubble class traveled at the same mean algebraic velocity in order to reduce the computational time and resources, signi3cant weakness of the model was evidenced in the prediction of the vapour velocity. Research is currently ongoing to consider additional momentum equations or develop an algebraic slip model to account bubble separation in order to yield a more realistic prediction of the vapour velocity.

2

Liquid Velocity (m/s)

1.8

Measurement MUSIG Boiling Model

1.6

Modified Boiling Model

Notation

Default Boiling Model

1.4 1.2 1

Qwall = 251.5 kW/m2 Tinlet = 92.1oC Pinlet = 0.143 MPa G = 1059.2 kg/m2s

0.8 0.6 0.4 0

(c)

0.2

0.4

0.6

0.8

1

(r-Ri)/(Ro-Ri)

Fig. 7. Local mean radial pro3les of liquid velocity: (a) C1, (b) C2 and (c) C3.

6. Conclusion A two-&uid model coupled with population balance approach is presented in this paper to handle bubbly &ows with the present of heat and mass transfer processes. The increase in complexity of modelling such &ows derives from the additional consideration of the gas or liquid undergoing a phase transformation. Subcooled boiling &ow belongs to a speci3c category of bubbly &ows with heat and mass transfer where it embraces all the complex dynamic interaction of the phenomena associated with hydrodynamics, heat and mass transfer, and bubbles coalescence and break-up. Modelling subcooled boiling &ows particularly at low pressures have been successfully demonstrated. The range of bubble

aif AB AC Aq cf C Cp CD CL CTD Cw1 ; Cw2 d dbw di ; dj dij d0 ; d1 DB DC Ds f fBV fi Flg drag Flg lift Flg

interfacial area concentration, m−1 bubble area, m2 cross-sectional area of boiling channel, m2 fraction of wall area subjected to quenching increase coeGcient of surface area constant in equation (22) speci3c heat, J=kg K drag coeGcient lift coeGcient turbulent dispersion coeGcient wall lubrication constants parent particle diameter, m bubble departure diameter, m daughter particle diameters, m equivalent diameter, m reference bubble diameters, m death rate due to break-up, m−3 s−1 death rate due to coalescence, m−3 s−1 bubble Sauter diameter, m bubble departure frequency, s−1 breakage volume fraction scalar fraction of each bubble size group total interfacial force, N drag force, N lift force, N

3138 lubrication Flg dispersion Flg g * g h G G(V; x; t) h0 hf hfg H k * n n(V; x; t) ni nj N N

P PB PC Qw Qc Qe Qq r RB Ri ; Ro Rph Si

St t tij T Tsat Tsub u * u ut v VB yw

G.H. Yeoh, J.Y. Tu / Chemical Engineering Science 59 (2004) 3125 – 3139

wall lubrication force, N turbulent dispersion force, N gravitational acceleration, m=s2 gravitational vector, m=s2 inter-phase heat transfer coeGcient mass &ux, kg=m2 s generation function, m−3 s−1 initial 3lm thickness, m critical 3lm thickness at rupture, m latent heat, J/kg enthalpy, J/kg thermal conductivity, W=m K normal to the wall surface bubble number density distribution, m−3 number density of the ith class, m−3 number density of the jth class, m−3 number of bubble classes active nucleation site density, m−2 pressure, Pa production rate due to break-up, m−3 s−1 production rate due to coalescence, m−3 s−1 wall heat &ux, W=m2 heat transferred by convection, W=m2 heat transferred by evaporation, W=m2 heat transferred by quenching, W=m2 radius, m bubble radius, m inner and outer radius of annular channel, m source/sink term due to phase change, m−3 s−1 source term due to coalescence and break-up, kg=m3 s stanton number thermo-&uid time scale, s coalescence time, s temperature, K saturation temperature, K subcooling temperature, K velocity, m/s velocity vector, m/s velocity due to turbulent collision, m/s volume corresponding to particle diameter d, m3 bubble volume, m3 adjacent point normal to the wall surface, m

Greek symbols  7 V % ; ;0 ; ;1 $  e

void fraction measured constant in Eq. (22) density di:erence = l− g , kg=m3 dissipation of kinetic energy, m2 =s2 bubble contact angle, rad reference subcooling temperatures = Tsat − Tl ; K turbulent kinetic energy, kg m2 =s2 size of an eddy, m e:ective thermal conductivity, W=m K

e 2  6 8ij P 0WN 0COND / .

e:ective viscosity, kg=m s size ratio between an eddy and a particle in the inertial sub-range density, kg=m3 surface tension, kg=s2 bubble contact time, s mass transfer rate, kg=m3 s bubble nucleation rate, m−3 s−1 bubble condensation rate, m−3 s−1 coalescence rate, m−3 s−1 break-up rate, m−3 s−1

Subscripts g gl l lg min w

vapour transfer of quantities from liquid phase to vapour phase liquid transfer of quantities from vapour phase to liquid phase minimum wall

References Anglart, H., Nylund, O., 1996. CFD application to prediction of void distribution in two-phase bubbly &ows in rod bundles. Nuclear Science and Engineering 163, 81–98. Antal, S.P., Lahey Jr., R.T., Flaherty, J.E., 1991. Analysis of phase distribution and turbulence in dispersed particle/liquid &ows. Chemical Engineering Communication 174, 85–113. Bonjour, J., Lallemand, M., 2001. Two-phase &ow structure near a heated vertical wall during nucleate pool boiling. International Journal of Multiphase Flow 27, 1789–1802. Bowring, R.W., 1962. Physical model based on bubble detachment and calculation of steam voidage in the subcooled region of a heated channel. Report HPR-10, Institute for Atomenergi, Halden, Norway. Chesters, A.K., Ho:man, G., 1982. Bubble coalescence in pure liquids. Applied Scienti3c Research 38, 353–361. Clift, R., Grace, J.R., Weber, M.E., 1978. Bubbles, Drops and Particles. Academic Press, New York. Cole, R., 1960. A photographic study of pool boiling in the region of the critical heat &ux. A.I.Ch.E. Journal 6, 533–542. Fleischer, C., Becker, S., Eigenberger, G., 1996. Detailed modelling of the chemisorption of CO2 into NaOH in a bubble column. Chemical Engineering Science 51, 1715–1724. Gopinath, R., Basu, N., Dhir, V.K., 2002. Interfacial heat transfer during subcooled &ow boiling. International Journal of Heat and Mass Transfer 45, 3947–3959. Hibiki, T., Ishii, M., 2000. Two-group area transport equations at bubbly-to-slug &ow transition. Nuclear Engineering and Design 202, 39–76. Hibiki, T., Ishii, M., 2002. Development of one-group interfacial area transport equation in bubbly &ow systems. International Journal of Heat and Mass Transfer 45, 2351–2372. Hibiki, T., Ishii, M., 2003. Active nucleation site density in boiling systems. International Journal of Heat and Mass Transfer 46, 2587–2601. Hsu, Y.Y., Graham, R.W., 1976. Transport Process in Boiling and Two-phase Systems. Hemisphere, Washington.

G.H. Yeoh, J.Y. Tu / Chemical Engineering Science 59 (2004) 3125 – 3139 Ishii, M., Zuber, N., 1979. Drag coeGcient and relative velocity in bubbly, droplet or particulate &ows. A.I.Ch.E. Journal 25, 843–855. Judd, R.L., Hwang, K.S., 1976. A comprehensive model for nucleate pool boiling heat transfer including microlayer evaporation. ASME Journal of Heat Transfer 98, 623–629. Kocamustafaogullari, G., Ishii, M., 1995. Foundation of the interfacial area transport equation and its closure relations. International Journal of Heat and Mass Transfer 38, 481–493. KonXcar, B., Kljenak, I., Mavko, B., 2004. Modelling of local two-phase &ow parameters in upward subcooled &ow boiling at low pressure. International Journal of Heat and Mass Transfer 47, 1499–1513. Krishna, R., Urseanu, M.I., van Baten, J.M., Ellenberger, J., 1999. In&uence of scale on the hydrodynamics of bubble columns operating in the churn-turbulent regime: experiments vs. Eulerian simulations. Chemical Engineering Science 54, 4903–4911. Kumar, S., Ramkrishna, D., 1996. On the solution of population balance equations by discretisation—I. A 3xed pivot technique. Chemical Engineering Science 51, 1311–1332. Kurul, N., Podowski, M.Z., 1990. Multi-dimensional e:ects in sub-cooled boiling. Proceedings of the 9th Heat Transfer Conference. Lahey Jr., R.T., Drew, D.A., 2001. The analysis of two-phase &ow and heat transfer using multidimensional, four 3eld, two-&uid model. Nuclear Engineering and Design 204, 29–44. Lee, T.H., Park, G.-C., Lee, D.J., 2002. Local &ow characteristics of subcooled boiling &ow of water in a vertical annulus. International Journal of Multiphase Flow 28, 1351–1368. Lemmert, M., Chwala, J.M., 1977. In&uence of Flow Velocity on Surface Boiling Heat Transfer CoeGcient. Academic Press and Hemisphere, New York, Washington. Lo, S., 1996. Application of population balance to CFD modelling of bubbly &ow via the MUSIG model. AEA Technology, AEAT-1096. Luo, H., Svendsen, H., 1996. Theoretical model for drop and bubble break-up in turbulent dispersions. A.I.Ch.E. Journal 42, 1225–1233. Olmos, E., Gentric, C., Vial, Ch., Wild, G., Midoux, N., 2001. Numerical simulation of multiphase &ow in bubble column reactors. In&uence of bubble coalescence and break-up. Chemical Engineering Science 56, 6359–6365. Pohorecki, R., Moniuk, W., Bielski, P., Zdrojkwoski, A., 2001. Modelling of the coalescence/redispersion processes in bubble columns. Chemical Engineering Science 56, 6157–6164. Prince, M.J., Blanch, H.W., 1990. Bubble coalescence and break-up in air-sparged bubble column. A.I.Ch.E. Journal 36, 1485–1499. Prodanovic, V., Fraser, D., Salcudean, M., 2002. Bubble behaviour in subcooled &ow boiling of water at low pressures and low &ow rates. International Journal of Multiphase Flow 28, 1–19. Ramkrishna, D., 2000. Population Balances. Theory and Applications to Particulate in Engineering. Academic Press, New York.

3139

Ramkrishna, D., Mahoney, A.W., 2002. Population balance modelling. Promise for the future. Chemical Engineering Science 57, 595–606. Ranz, W.E., Marshall, W.R., 1952. Chemical Engineering Progress 48, 141–148. Rotta, J.C., 1972. Turbulente Stromungen. B.G. Teubner, Stuttgart. Sato, Y., Sadatomi, M., Sekoguchi, K., 1981. Momentum and heat transfer in two-phase bubbly &ow—I. International Journal of Multiphase Flow 7, 167–178. Shimizu, K., Takada, S., Minekawa, K., Kawase, Y., 2000. Phenomenological model for bubble column reactors: prediction of gas hold-up and volumetric mass transfer coeGcients. Chemical Engineering Journal 78, 21–28. Stone, H.L., 1968. Iterative solution of implicit approximations of multidimensional partial di:erential equations. SIAM Journal of Numerical Analysis 5, 530–558. Tolubinsky, V.I., Kostanchuk, D.M., 1970. Vapour bubbles growth rate and heat transfer intensity at subcooled water boiling. Proceedings of the Fourth International Heat Transfer Conference, Vol. 5, Paper No. B-2.8. Tu, J.Y., Yeoh, G.H., 2002. On numerical modelling of low-pressure subcooled boiling &ows. International Journal of Heat and Mass Transfer 45, 1197–1209. Unal, H.C., 1976. Maximum bubble diameter, maximum bubble growth time and bubble growth rate. International Journal of Heat and Mass Transfer 19, 643–649. Wang, S.K., Lee, S.J., Lahey Jr., R.T., Jones, O.C., 1987. 3-D turbulence structure and phase distribution measurements in bubbly two-phase &ows. International Journal of Multiphase Flow 13, 327–343. Wu, Q., Kim, S., Ishii, M., Beus, S.G., 1998. One-group interfacial area transport in vertical bubbly &ow. International Journal of Heat and Mass Transfer 41, 1103–1112. Yeoh, G.H., Tu, J.Y., 2002. Implementation of a two-phase boiling model into the RELAP/MOD2 computer code to predict void distribution in low-pressure subcooled boiling &ows. Nuclear Science and Engineering 140, 182–188. Yeoh, G.H., Tu, J.Y., Lee, T.H., Park, G.-C., 2002. Prediction and measurement of local two-phase &ow parameters in a boiling &ow channel. Numerical Heat Transfer, Part A: Applications 42, 173–192. Yun, B.J., Park, G.-C., Song, C.H., Chung, M.K., 1997. Measurements of local two-phase &ow parameters in a boiling &ow channel. Proceedings of the OECD/CSNI Specialist Meeting on Advanced Instrumentation and Measurement Techniques. Zeitoun, O., Shoukri, M., 1997. Axial void fraction pro3le in low pressure subcooled &ow boiling. International Journal of Heat and Mass Transfer 40, 867–879.