slug flow regime transition

Nuclear Engineed.ng Nuclear Engineering and Design 151 (1994) 15-39

ELSEVIER

andDesign

Analysis of the bubbly/slug flow regime transition* S. K a l k a c h - N a v a r r o ,

R . T . L a h e y , Jr., D . A . D r e w

Center for Multiphase Research, Rensselaer Polytechnic Institute, Troy, NY 12180-3590. USA

Received 25 September 1993

Abstract An extended two-fluid model that includes the dynamics of bubble cluster formation and breakup was developed and used successfully to predict the bubbly/slug flow regime transition. In particular, a linear void wave stability analysis of the extended two-fluid model was perfotTned and it was found that an instability of the void wave associated with the bubble clusters triggers a bubbly/slug flow regime transition. The predictions obtained for the change in sign of the void wave damping coefficient of this model, and tiros for the onset of flow regime transition, agreed very well with the available experimental data. Another important result of this work is that it provides a mechanistic transport model for the prediction of interfacial area density, bubble number density, and mean bubble radius, parameters which are required for closure of two-fluid models.

1. Introduction

When vapor and liquid flow in a vertical duct, the distribution of the phases can assume different configurations which are called flow regimes. If the particular flow regime present in a given twophase system is known, the most accurate formulation to analyze the thermal-hydraulic behavior of the system is called a two-fluid model. A twofluid model considers each phase independently. including a set of mass, momentum and energy equations for each phase, and it represents the phasic interactions using interracial transfer terms. The interracial transfer terms are modeled using constitutive laws which are flow regime depen-

* Mark Mills award paper. Elsevier Science S.A. SSDI 0029.5493( 94J 00846-Q

dent. This is why a knowledge of the flow regime is of great importance if one is to model mechanistically two-phase flow. However, the state-ofthe-art understanding of flow regime transition is quite limited, thus simplified flow regime maps and/or correlations hae been widely used to infer flow regime transition. Unfortunately, these maps are normally not mechanistically based and are known to be inadequate even for steady-state, let alone transient, applications. Nevertheless, if the physics that causes flow regime transition is included in the two-fluid model, the resultant extended two-fluid model should be able to predict row regime transition without the use of flow regime maps. Indeed, the objective of this paper is to present an extended two-fluid model which is capable of predicting the bubbly/slug flow regime transition. In bubbly flow, as the void fraction

16

S. Kalkach-Navarro et aL f Nuclear Engineering and Des(en 151 (19941 15 39

increases, the bubbles tend to group together and form bubble clusters. This phenomena has been previously observed in air-water mixtures by Saiz-Jabardo and Bour6 (1989) and by Park (1992) in air-oil mixtures. If the bubble clusters femain intact enough for the liquid film between the bubbles to drain away before the bubbles are knocked apart by liquid phase turbulence, they will coalesce and form larger bubbles, and this process will iead to the formation of Taylor bubbles (i.e. the transition to slug flow). Indeed, it appears that the dynamics of bubble cluster formation and breakup is responsible for the bubbly• slug flow regime transition. In order to predict the bubbly/slug flow regime transition, it is necessary to model the dynamics of the formation and breakup process of bubble clusters. A new mathematical modeling approach was developed to introduce the relevant physics, and tl;e resulting transport equations were coupled with a state-of-the-art two-fluid model to generate an extended two-fluid model which was found to be capable of predicting the bubbly/slug flow regime transition. It has been previously proposed by Bour6 (1988) that void wave instability may lead to a bubbly/ si~.g flow regime transition. Indeed, the results obtained from stability analysis of the model developed herein show that an instability of the void wave associated with the bubble clusters signals a bubbly/slug flow regime transition. Significantly, the predictions for the onset of the flow regime transition obtained from a linear stability analysis of the extended two-fluid model presented herein agree very well with the available experimental data. This is an important result since it provides a mechanistic approach for prediction of the bubbly/slug flow regime transition without the shortcomings of flow regime maps or correlations. A two-fluid model is based on a detailed treatment of the interactions at the interface. Such treatment is necessary since mass, momentum and energy are exchanged from one phase to the other through the interface. Moreover, this is why the interfacial area density is one of the most important parameters in a thermal-hydraulics analysis of two-phase flow systems based on two-fluid models. Thus, another important result of the

present research is that it allows one to predict the evolution of interfacial area density, mean bubble radius and bubble number density in two-phase flows.

Bubbly chlsters As noted previously, it has been observed by Park (1992) that bubbles in air-oil mixtures have a tendency to group together and form clusters. The same phenomenon was observed in air-water mixture~ by Saiz-Jabardo and Bout6 (1989). However, owing to the much shorter lifetime of bubble clusters in turbulent air-water flows, SaizJabardo and Bour6 were not able to quantify their effect on flow regime transition. Other authors have also observed the importance of bubble clusters (Smereka, 1993); (Sangani, 1991). In an air-oil mixture, the turbulent breakup of the bubble clusters is suppressed and the drainage time for the interstitial liquid is relative long, thus Park (1992) was able to measure the rise velocity of the clusters for different values of void fraction, obtaining a correlation for the dependence of the rise velocity of the bubble clusters on void fraction. Typical bubble clusters for air-oil bubbly flow in an annular conduit are shown in Fig. I. Coalescence between a bubble and a bubble cluster, or between two bubbles, in turbulent flows occurs in three stages. The first stage corresponds to a collision, where it is important to note that the velocity difference in the wake behind the bubble cluster enhances the rate of bubble collisions and, therefore, the formation and growth of bubble clusters. The second stage is the draining of the liquid film between the adjacent bubbles in a cluster until it reaches a critical thickness ho and the final stage is the rupture of the liquid film resulting in bubble coalescence. If a cluster remains intact long enough for the liquid film between the bubbles to drain to thickness hc before being knocked apart by liquid phase turbulence, the bubbles will coalesce and form larger bubbles. This process will eventually lead to the formation of Taylor bubbles and thus a transition to slug flow. This situation before and alter bubble coalescence can be seen in Fig. 1 and 2 respectively.

S. Kalkach-Navarro et al./ Nuclear Engineering and Design 151 (1994) 15-39

17

the void fraction ~, the time z required for drainage of the liquid film between adjacent bubbles and the time T that the clusters remain intact before being knocked apart by liquid phase turbulence.

. . .i!i

': z

•

2. Description of the mathematical model

~\

.+:

i

+2.

Fig. 1. Typical bubble clusters (Park, 1992).

r . •

,2+~!! i,

~

'? ',2

~'i!

,:'

ti ':

i:

+

/

? 'f /,

In order to model the behavior of the distribution of sizes present in bubbly two-phase flow, it is necessary to account for bubble cluster formation and breakup. Given the fact that the breakup of individual bubbles is not expected to be of importance for the case we are considering, only bubble cluster breakup will be considered. This implies that there should be a minimum volume above which a bubble cluster, instead of an individual bubble, exists, and therefore this should be the minimum volume at which breakup will occur. In particular, the minimum volume for breakup v~ is that of the largest individual bubble. This situation suggests that we can divide the total volume range into two groups: group 1 includes all the bubbles with a volume smaller than the critical volume v~ and thus comprises individual bubbles, whereas group 2 includes all the bubble clusters which have a volume larger than t~. This method of dividing the system in two regions is similar to the multigroup method of neutron transport theory (Lamarsh, 1972). Now, let fj(v, r, t) be the bubble ( i f j = 1), or bubble cluster ( i f j = 2), size distribution function, such that f.(t,,r, t)dv is the number of bubbles ( j = 1), or bubble clusters ( j = 2), per unit volume, at vectoral position r and time t, having a volume between v and v + dr. Since the size distribution function of the individual bubbles, f ~ ( v , r , t ) , is only valid for 0 < v < v~, and the size distribution function of the bubble clusters is only valid for v 1> t'c, the total size distribution f u n c t i o n f ( v , r, t) is given by f ( v , r, t) =.fl(v, r, t)H(v¢ - v) +f,(v, r, t)H(r - t~c)

Fig. 2. Formation of Taylc,~ bubbles (Park, 1992).

From this discussion, it is clear that in order to predict the bubble cluster breakup/formation rate ratio 7, the significant physical properties should be

(I) where H(v) is the Heaviside step operator. Thus, it can be observed that the total size distribution function is equal to the distribution of the

18

S. Kalkach-Nnrarro et al. / Nuclear Engineering and Design 151 (1994) 1 5 - 3 9

individual bubbles if v < vc or the bubble clusters if v ~ v~. Significantly, the total number density N ' , total mean radius/~[', interfaciai area density A~", and the local void fraction ~, are related to the zeroth, I 3, 2 and first moments of the distribution respec~, tively, and are given by

Substituting Eqs. (10) and (11) in Eq. (9) we obtain

(2)

In particular, the interfacial area density AI"--(36r01/3M2/3, which is given by the ~th moment, is given by

N"(r, t) =

f-

f ( v , r, t) dv

Rb(r, t) = R~'(r, t)/N"'(r, t)

(3)

M,2 =

v ~ ( v , r , t) dv

(11)

c

Mi = Mit + M,~

A~" - A~'~'+ A~

(12)

(13)

where -\4n}

Jo v'/3f(v'r't) dv

(4)

and A i'""(r, t) = (36rt) I/3 ~ " v2/~(v, r, t) dv

(5)

• (r, t) = / d

(6)

tf(v, r, t) dv O

We also note that the average volume of the bubbles is given by

fo

(7)

-~- ~ / N "

Ub = r : ~

J0 f ( v , r, t) dv Similar expressions can be obtained for the number density Nj", mean radius /~bj, interfacial area density, A ~', and local void fraction ~j for each of the groups, by noting that in general a moment of order i, M,., is given by

f

(8)

v f ( v , r, t) dv

substituting Eq. (I) into Eq. (8) we obtain M, =

fo

"~fl (v, r, t) dv +

f

vf2(v, r, t) dv

(9)

c

If we define Mq as the ith moment off(v, r, t) for the jth group ( j = I, 2), then Mil =

and

0f(v, r, t) ~ + V. If(v, r, t)uBl = B(v, r, t) + C(v, r, t)

(14)

" vf(v, r, t) dv

M, =

This implies, as expected, that the total interfacial area density is just the interfaciai area density of the individual bubbles plus the interfacial area density of the bubble clusters. The transport equation to describe the bubble size distribution function can be obtained by performing a balance of the loss and gain due to breakup and formation of the bubbles or bubble clusters, and is given by

fo

v ~ (v, r, t) dv

(10)

where B(v, r, t), C(v, r, t) are the total rates of bubble cluster breakup and formation respectively, and u~(v) is the velocity of the bubbles (if v < vc) or bubble clusters (if v >/v~), and is given by ug(v) = uel H(vc - v) + ug2H(v - vc)

(15)

where ugt is the average velocity of the individual bubbles and ug2 is the average velocity of the bubble clusters. It can be assumed that during bubble cluster formation or breakup, the total volume of gas does not change. For example, the grouping of a bubble of volume u with a bubble of volume (v - u ) restalts in a bubble cluster of volume v. The probability of the formation of a bubble Ouster involving two bubbles, or bubble clusters, is assumed to be proportional to the product of their distribution functions and a bubble duster formation kernel. The total rate of formation will be given by the integral of all possible combinations, and can be expressed as

S. Kalkach-Navarro et al. [ Nuclear Engineering and Design 151 (1994) 15-39

C(v, r, t) = ~ -

c(u, v - u)f(u, r, t)f(v, - u, r, t) du c(u, v)f(u, r, t)f(v, r, t) dv

(16)

where the bubble cluster formation kernel c(u, v) is the volume per unit time involved in bubble cluster formation. The first term in Eq. (16) represents the gain in the number of bubble clusters of volume v due to the grouping of smaller bubbles. The second term in Eq. (16) represents the loss in the number of bubbles or bubble clusters of volume v due to grouping of a bubble of volume v with any other size bubble. The factor of one-half in the first term avoids double counting. Next, the change of f(v, r, t) due to bubble cluster breakup is given by B(v, r, t) =

bffu, v)f(u, r, t) du -

b t(v, u)f(u, r, t) du

(17)

19

and break-up gain kernels over all the volumes that gain from break-up. The gain in the number density of bubbles of volume u, from the break-up of bubbles or bubble clusters of volume v, is given by Ng'ain(u, v) = bffu, v)f(v, r, t) du

Similarly, the gain in the void fraction of bubbles or bubble dusters of volume u, generated from the break-up of bubble clusters of volume v is given by ~g,i,(u, v) = ubg(u, v)f(v, r, t) du

(19)

The total void fraction gained by all bubbles or bubble clusters of volume smaller than v, owing to the breakup of bubble clusters of volume v, is ~g,~,(v) =

ubg(u, Of(v, r, t) du

(20)

The loss in number density of bubble clusters of volume v, due to their break-up to form bubbies of volume u, is Nios~(u, v) = bl(v. u)f(v, r, t) du

Where the bubble cluster breakup loss kernel b~(v, u) is the number of bubble clusters of volume v, per unit volume per unit time, that break-up to produce bubbles or bubble clusters of volume u. Similarly, the bubble cluster break-up gain kernel bg(u, v) is the number of bubbles or bubble clusters of volume u, per unit volume per unit time, that are gained from the break-up of bubble clusters of volume v. This equation implies that the gain o f f ( v , r, t) due to the break-up of bubble clusters of volume v depends on the distribution function of the bubble clusters with volumes larger than v being shattered, multiplied by a break-up gain rate kernal bg(u,v). It also implies that the loss of f ( v , , , t), due to break-up of bubble clusters of volume v to produce smaller bubbles or bubble clusters, depends on the distribution function of bubbles of volume v, multiplied by a break-up loss rate kernel b~(v, u). A relationship between these two kernels can be obatined by using the condition that the fraction gained from break-up has to be equal to that lost in break-up. This condition can be expressed by integrating the break-up loss

(18)

(21)

The loss in the void fraction of bubble clusters of volume v to form bubbles or bubble clusters of volume u is • los~(u, v) = vbJ(v, u)f(v, r, t) du

(22)

Hence, tile total void fraction of bubble clusters of volume v lost to form smaller bubbles is ~o.,,,(v) =

vbl(v, u)f(v, r, t) du

(23)

Thus the continuity condition that the total void fraction lost by bubble clusters of volume v to form smaller bubbles is the same as the total void fraction gained by those smaller bubbles, can be expressed as • ~,i,(v) = ~lo.,s(V)

(24)

Substituting Eq. (20) and Eq. (23) into Eq. (24) we obtain

f0

ubg(u, v)f(v, r, t) du =

vbl(v, u)f(v, r, t) du (25)

S. Kalkach-Navarro et al. / Nuch'ar Engineering and Design 151 (1994) 15-39

20

Rearranging Eq. (26) we obtain

+ ~

c(u, v -- u ) f ( u , r, O f ( v , - u, r, t) du

~ [ub~(u, v ) f ( v , r, t) -- vbl(v, u)f(v, r, t)] du = 0 -

(26) Eq. (26) is satisfied by the vanishing of the integrand, thus ub g(u. v) hi(v, u) -- - -

(27)

c(u, v)f(u, r, t ) f ( v , r, t) du

(29)

It can be easily shown that the right-hand side of Eq. (29) conserves void fraction by multiplying it through by v and integrating over all volumes (i.e. by obtaining the first moment):

v

This is the relationship between bg(u, v) and bin(v, u) that has to be satisfied for the break-up kernals to conserve void fraction. The modeling of the kernels for bubble cluster formation and break-up will be discussed in the next section. If we substitute the bubble cluster formation and break-up rate expressions, Eq. (16) and Eq. (17), into the balance equation for the bubble size distribution function, Eq. (14), we obtain the following form of the Boltzmann transport equation: Oftv, r, t)

~

,fo ;

v C ( v , r, t) dv -- -2

c(u, v - u)

x f ( u , r, t ) f ( v - u, r, t) du dv -

v

c(u, v ) f ( u , r, t)

x f ( v , r, t) du dv

(30)

Changing the order of integration in the first term on the right-hand side of Eq. (30) we obtain

i

+ V- [f(v, r, t)ug]

v

vC(v, r, t) dv =

v

c(u, v - u)

o

x f ( u , r, t ) f ( v -- u, r, t) dv du =

bg(u, v)f(u, r, t) du -

bl(v, u)f(v, r, t) du

+ ~

~t

=

bg(u, v)f(u, r, t) du

--

I ~'~ bg(u, v)f(v, r, t) du Jo

c(u, v)f(u, r, t)

(31)

Performing a change of variables in the first term on the right-hand side of Eq. (31) we obtain (28)

Next, substituting Eq. (27) into the second term on the right-hand side of Eq. (28), the transport equation for the bubble size distribution function, in terms of the break-up gain kernel, is ~f(v, r, t) - - + V " [ f ( v , r , t)Ug]

v

x f ( v , r, t) d u dv

c(u, v - u)f(u, r, t ) f ( v - u, r, t) du

- Jo c(u, v)f(u, r, O f ( v , r, t) du

;f:

;

,;f/

vC(v, r, t) dv = -~

( x + u)c(u, x)

x f ( u , r, O f ( x , r, t) du d x -fffoVC(u,v)f(u,r,t) x f ( v , r, t) du dv

(32)

The first and second terms of Eq. (32) cancel each other, thus we can conclude that bubble cluster formation conserves void fraction. The bubble cluster break-up terms conserve void fraction since, as mentioned before, such a condition, Eq. (24L was imposed on them.

21

S. Kalkach-Navarro et aL / Nuclear Engineering and Design 151 (1994) 15-39

Bubble cluster break-up and formation kernels

I

In order to use the transport equation for the bubble size distribution function, Eq. (29), it is necessary to define the bubble cluster break-up and formation kernels c(u, v) and bg(u, v) respectively. It has been proposed by Bilicki and Kestin (1987) that an important fator in the grouping of two bubbles or bubble clusters is the velocity defect in the wake of the leading bubble (or bubble cluster), which accelerates the trailing bubbles and causes them to overtake the leading bubble(s). Since the bubble cluster formation kernel represents the volume per unit time involved in bubble cluster formation, it can be modeled by considering the velocity difference in the wake behind a leading bubble (or bubble cluster) and can be defined as the volume V,.~the trailing bubble (or bubble cluster) has to occupy at time t-dt to group with the leading bubble (or bubble cluster) at time t, per unit time (i.e. Va), multiplied by the probability for the bubble to be in the wake of the leading bubble or bubble cluster P. This probability is proportional to the area of the cross-hatched region shown in Fig. 3. Thus, the bubble cluster formation kernel for leading bubble (or bubble cluster) i is given by: c, = P I;',.~

(33)

Again, referring to Fig. 3, the volume that bubble 1 has to occupy at time t-dt to group with bubble 2 per unit time is, in terms of the equivalent radii of the bubbles rhl and rb2

~'~l = 2~

rAut(r ~ _ r z) i/~,dr

(34)

where r~ = (rb~ + r~2). The integration domain corresponds to the shaded area in Fig. 3, and Aut is the velocity difference in the wake, which is proportional to (Schlichting, 1979)

f ga'~U"f rb, ~2/3Cg.~

"" :"'tT,) t-r/ where Co is the drag coefficient, ~ is the distance between the centers of the leading and trailing bubbles or bubble clusters, and k~ is a

I I

/

Au e

lz I

\

r

Fig. 3. A bubble i:~ the wake of a leading bubbJe.

non-dimensional constant of proportionality, which together with all the other constants that follow (k2,k3), will be included in a bubble cluster formation rate constant co determined from experimental data. The drag coefficient was modeled using Harmathy's model (1960) for distorted bubbles: c . = 4[g_(p,- P,)1j ''2

(36)

where a is the surface tension, at is the void fraction, p~ and pg are the densities of liquid and gas respectively. The dependence of Co on ~t, and the constants involved in the expression for Cn, can be included in the bubble cluster break-up rate/formation rate ratio ~, and determined experimentally. Substituting Eqs. (35) and (36) into Eq. (34) and integrating, we obtain

,, f3~2/3k rgS(p=-- flg)20-]l/12

rb'(rb' -t-rb2) 4/3 (37)

22

S. Kalkach-Navarro et al. / Nuclear Engineering and Design 151 (1994) 15-39

The probability that a bubble will be in the wake of the leading bubble(s) is assumed to be proportional to the area occupied by the wake divided by the area of the conduit: p = k2

/~(u,v)

(38)

2 where De is the diameter of the conduit and b is the width of the wake at location r~ and is given by (Schlichting, 1979) b = ka cU3(rbl + rb2 ) 1/3r~il3

(39)

Thus, the bubble cluster formation kernel when bubble I is the leading bubble is, in terms of rbl and rb2

/4

Cl|(rbl,rb2 ) = t ~ / l : )

\s/3

r

1

co[r3l(rbl-k-rb2) 2J

(40)

where

(3

.

co = zn -----h-2--

(41)

We note that co has units of rate per unit area (l/sin 2) and contains the non-dimensional constants k,, k2, k3 which are unknown. This parameter can be determined experimentally, The total bubble cluster formation kernel must consider the two possibilities that, bubble 1 is the leading bubble (q,) and that bubble 2 is the leading one (qz), hence

c(r:, I , rb2) = Cl, (rbl, rbZ) + Cl2(rbl, rb2)

(42)

Substituting Eq. (40) into Eq. (42) we obtain c(rb,, rb2) =

~)

Co[(r~, + r~2)(rb, + rb2)2l (43)

Finally, transforming Eq. (43) in terms of bubble volumes u and v, we obtain c(u, v) = Co(U + v)(u '/~ + v'/3) 2

Fig. 4. Bubble cluster break-up kernel. the distribution function obtained experimentally (Kalkach-Navarro, 1992). It can be observed in Fig. 4 that the pair of bubbles, or smaller bubble clusters, of volumes v~ and u - v) is less probable owing to the break-up of a bubble cluster of volume u, than the pair of bubbles, or bubble clusters, of volumes v2 and u - 02. As mentioned earlier, break-up is only considered for bubble clusters and the minimum volume at which it occurs is the critical volume, v~. Thus, the parabolic, normalized distribution for the break-up gain kernel is given by bffu, v) = 12bo

v(u - v)

u3

H(v - v ~ )

(45)

where v~ is the critical volume above which a bubble cluster will break into individual bubbles, H(v - v~) is the Heaviside step operator and bo is the rate of bubble cluster break-up, which must be determined experimentally. Substituting the bubble cluster formation and break-up gain kernels just derived into Eq. (29), we obtain the final form of the transport equation for the probability distribution function: ~f(v, r, t) + V" If(v, r, t)ug] 5t = 12

""

v(u -

v)

bo~H(v-vc)f(v,r,t)

du

(44)

For the bubble cluster break-up kernel, it can be assumed that a bubble cluster of volume u will break into a pair of bubbles, or bubble clusters, with volumes v and u-v, with a probability which is symmetric about ,t/2. The parabolic distribution shown in Fig. 4 was assumed because it was a simple symmetric function and the function that best represented the behaviour of the moments of

- 12

I

'

u(v -

bo ~

u)

H(u - re)f (v, r, t) du

! "'-

+~

Co[u +(V--U)][Ul/3+(V --U) 1/312

x f ( u , r, t ) f ( v - u, r, t) d u -

co(u + v)

x [u 1/3 + v 1/312f(u, r, Of(v, r, t) du

(46)

23

S. Kalkach-Navarro et al. [ Nuclem Engineering and Design 151 (1994) 15-39

Owing to the obvious complexity of Eq. (46), it was not possible to solve it analytically; however, since we are interested primarily in the moments of this equation, we may proceed to obtain the moments for each of the groups. If we determine the zeroth moment of Eq. (46) for the bubbles in group 1 we obtain the conservation equation for the number density N'/' of the individual bubbles as ON';' 0---7- + v . [ N ;'u~,]

had to determine a conservation equation for each of the moments involved in these equations, we would require an infinite number of equations to solve for the moments of each group, since the conservation equation for the moment M o. has introduced higher moments. Thus, it is clear that it is necessary to m'..ke some assumptions in order to close the system with a finite number of equations. To this end, we can perform a Taylor series expansion around the average volume for the corresponding group in the definition of the moment Mo.. For example, let us consider

3 -co[M5/31 N'" = bo(6h'l-22vc: - 4M-3:, v c) i ~ N " i' + 2M4/3, R " + 2M4/3RI" + A~"o~A[': +,,,5/3

- (M5/3, N~" + 2M4/3, R'(' + A ["oq )/2]

(47)

Following the same procedure as used to obtain Eq. (47), the conservation equation for the void fraction of the individual bubbles is given by the first moment of Eq. (46): &q 0--t + V • [~u~,] =bo( 4M_22v 2 - 3M_3: v~)

Mil =

(51)

v f ( v ) dv

A Taylor series expansion around volume vo yields

(52)

v ~ ~ v~ + ivJo- ~(v - Vo)

Hence, Eqs. (51) and (52) imply MiI ~

(53)

[rio + iv~o- I(v - vo)lf(v) dv

Thus - co[Ms/3, N " + Msl3~j + 2M7/3, R " M i, ~. v~o( ! -- i) I ''~f ( v ) dv + ivy- i

+ 2M4/3M4/3, + A["M21 + aMs/3, - - ( M s m N'(' + 2My31R'~'+ ,41~M2,)/2]

(54)

(48)

Similarly, the conservation equation for the number density N~' of bubble clusters will be

Using the definitions of N'[' and ~1, (55)

M:l ~ Vto( 1 -- i)N~" + ivy-'~l

~N~' O--t- + V " [N~'ugz] = bo(N '" 2 + 6 M _ ,-~ , v 2 - 4 M _ 3 , v 3) - co[Ms~s,- N ,,, + Ms/3N'[' + 2 M 4 / 3 2 R " - 2Ma/3R'~' + Ai"ot,

. . . . . NI ~ --gAil. +(Ms/s, - 2M

vf(v) dv

3o

w " + A ' "il~l )]2] 4/311),1

If we take Ooas the average volume of the bubbles in group 1, tS), Eq. (7) implies

VO~-/~I =N.I, ~

(56)

Hence, Eqs. (55) and (56) yield

(49) Finally, the void fraction ~z of the bubble clusters can be determined using the fact that the total void fraction ~ is not changed by the formation or break-up of bubble clusters, thus ~2 = ~ - ~t

(50)

In Eqs. (47), (48) and (49) we note that new moments are generated by the process, thus if we

M;, ~

(1 - i)N';' + / I T , , l

~,

(57)

or, rearranging, Mi, .~ O] N'~"

(58)

In general, we have M~/ "~ O~Nj"

(59)

24

S. Kalkach-Navarro et al. / Nuclear Engineer#~g and Design 151 (1994) 1 5 - 3 9

where ~s is the average volume of the bubbles group j and is given by a definition analogous to that presented in Eq. (7). That is,

.%

aN,'

- -

+ 7 . [N'~'us~]

0t

= OO(OOt 2 " P C 2 " D¢ - -

c t 6ats/3 N ''~/3 + ~t5/3 N,,,E/3N~,

(6O)

0~ = N j,,

4~t ~'3 N~'4v~)

~¢513~",2/31ff", .L. ~ 4 / 3

Substituting Eq. (59) into Eq. (12), the ith moment M~ is given by

' ~ 413/~d',,,- 113N 113

"~-~"1 ' ' 1 +

(61)

" + v"i, N : "' M, = v]N~

In particular, from Eqs. (5) and (59), the total interfacial area density will be given by AI" = (36rt) 1"3(13~/3N'~"+ f~/3 N~') = A,, + A,2 (62)

~2

&Tin213 "w2

~t~/3N~'ll3~I + ~t~/3N'~'ll3~x,_)

aNi" a--~- + V. [N]"~,] =

ON;."

a--S--+ v . [N_;',,,_]

- Ce(4~s,/3N~'113-- 2~13N'(113)

(65)

and, f,;-r th,: volume fraction of the individual bubbles ~ we obtain

0 t + V' [~1%~1 l>o~4a , "pc., v~. -- 3~ j 3 N'~'4v~ ) t. 1 ~ , . ~ 8 1 3 1 U ' - 2 / 3

,.¢8/3 M m - 5 / 3 M m

~ . 4 [ 3 '" ~ d 2' ' - 113,..,4/3 q- 015/3N'~'-213011 -k- •~'~2 ~1 '"~./-mI I/3 M ' 1 1 3 , ~~'1,, 2 M m+ 2~713N'I"-4/3~I3N~ '2/3 + ,~2/3 -'2 ,,,,

+ ~13Ni"-"13~)

I

(66)

-- ~-ot" I " [~513M'",,I "~r'"

-513 -5/3 #' ,,, "~-4/3 M i n t - 1 / 3 ~t'"' -~(v~ N -,I +v2 N , . ) N I +~-Vl , , ~ V l -I

+ v 2-~13~v~ ...... + 2~]I3N~'(f~I3N'I''+ f~:3N~'). .

+ 01Ni"(f~/3N'('+ O~/3N~ ') +

+ 4~3N~'4v~)

S t e a d y - s t a t e s o h a i o n s o f the s y s t e m o f c o n s e r v a t u m equations

boN,. (6v,_ "Vc -- 4V~-3 V~)-

i ~ 2 1 3'" 1 1 II" g~Vl - ,leT" ~1 , I

(64)

Following the same procedure with the other balance equations we obtain the following for N~':

= bo(N~' - 6~-2N~'3v~

which agrees with Eq. (13). This approximation allows us to obtain the ith moments as a function of the independent variables N]", ~q, N~' and ~t_,, and thus to close the system with one algebraic equation, and three additional transport equations, in additions to those in a two-fluid model for phasic mass, momentum and energy. An analysis of the interdependence between these equations is performed later on. If we substitute the approximation given in Eq. {61) into the balance equations for the number density of the individual bubbles, Eq. (47), we obtain

1~I"'- I/3 ,.,, 1131t[,,'213

+

v- 2 N ,'~' )

-2/3~513N~"2]

In order to obtain the steady-state solution of the system of conservation equations, the equations were non-dimensionalized using the following scaled, non-dimensional variables: 6. .~, = -/

(67)

a =-

(68)

trc

(63)

Now, introducing the definition of the average volume for the individual bubbles el and for the bubble clusters fz given in Eq. (60), the balance equation for the number density of the individual bubbles N'( can be expressed in terms of the variables ~ , ~z, N~" and N;_" as follows:

• )

=

A ,,,

-7

Nt =

N'l"V~

(69) (70)

S. Kalkach-Navarro et al.

~,

I

where

(71)

_ N~'v~

?

~2 = ~ - ~J

where the parameter ~: is the ratio of the bubble cluster break-up and formation rates: (72)

7 = bn[3o

in which 30 = Co/V~.0 has been used. The conditions satisfied for a steady, fully developed flow are c3m~" = 0 ~t

(73)

V. [N~",~,] = 0

(74)

-

-

Substituting Eqs. (73) and (74) into Eq. (64) we obtain b oLO~2"lv2v~ r.- -'~L,,,,3 "~ •

A~-31tT,,,4.,'t'. ~,2 "~- " ~ . ~ - C o ( 6 ~ / S N ' t

,,,5/3 &U"-"/3 M'"

"~'1

J'll

"-~'2

"l/s

5 1 3 M ' ~ - 2 / 3 M'"

"1-~2" Jr2

"~'t

+ 2~ 4/s N"'t/sx I~l.~Ar,,,2/s _,_ ,~,,4/.~ xr .... i/s,,~ 2t13 IAT,,,213 . aVl T ~ ' ~ I lVl v2 + ~/SN;"'/s~, + ct~/SN'['l/~) = 0

(75)

Eq. (75) can be expressed in terms of the non-dimen~ional variables just defined as follows: (AZ -2~T"3 ~2

~'¢2

A~ - 3 ~ , " 4 " 1

--"v~2

av2

i~,~513 AT'"II3

I--%°'~'1

,¢,Sl.~ 10 .... 213 10 ,,,

+~1

z 5/3 gr,,,- 213 ~,,,

z ' 2 "-t-:z2 ~v2

'Vl

*Vl

~'1

,%¢4/3 ,~,,,- I/3,¢ its ~,,,2t3 j 9 ^4/'~I~T .... n/3~

"F~2

~'1 ' V l

a'2

"

T~I

"tVl

"

~,,,'~

^

Since the system of Eqs. (76) through (79) is complicated, it was not possible to find an analytical expression for the fully developed, steadystate solution; however, the original set of partial differential equations can be evaluated numerically by assuming a set of initial conditions and integrating until a fully developed, steady-state is reached. The results obtained for the total, individual bubble and bubble cluster number densities 29" N t , and N z , respectively, are shown in Fig. 5. For very low ~ the total number density .~" is approximately equal to the number density of the individual bubbles N'~". This is, the number density of the cluster/V~' is almost zero and the system is basically composed of individual bubbles. This indicates that for a very low void fraction, very high breakup rate bo, or very low bubble cluster formation rate co, bubble clusters will essentially not be present in the steady-state. This behavior is physically reasonable. As the scaled void fraction ~ is increased, /V" and AT'(' increase up to a maximum and then N;' starts decreasing as more bubble clusters are formed, and ~ ' continues increasing. This be^

pit

~

rtt

•

0,12

(76)

In the same way, the non-dimensional steady, fully developed number density of the bubble clusters is given by Eq. (65) as . . . . .

(79)

i/=, , ~ ,,21"~

"~2 "~v2

+ ~2/Slfl~"/sa~ + ~/s~'('~/sa.,) = 0

(]~

25

Nuclear Engineering a n d Design 151 (1994) 1 5 - 3 9

0°~

*

N~"

o

Nm

~

3 ~,,,4 X

- (4a~:.~Ifl~'t/s - 2~/s~i "~/~) = 0

(77)

Finally, the steady, fully developed, non-dimensional void fraction of the individual bubbles is given by Eq. (66): ( 4 ~ ~-2/~,,s_ _ -- 3 ~ . ~ , 4 ) _ . -

_ ( 6 ~ / . ~ ~,,~.... 213 +~t~s/3~'"~

413)V~'+~'2Z5/3~''2..... 2:3.,¢~1 Tz ' ~")Z 2 ± 4/3 ~ ' ' 2. . . . 113(~ 413

X /V ~"- 1/3 -,-~l-I9~713 ,~vl~r"'4/3'~l ~2113 ~v2~'2/3

+~,,_,~:/.~¢~,,t/s,~2~O,,2~,~,,i.... t + ~t~/'~lfl~"-2/s~.,.) = 0

(78)

Fig. 5. Steady-statenon-dimensionalnumberdensity /V't",/V~' and /~/7

S. Kalkach-Navarro et al. / Nuch, ar Engineering and Design 151 (1994) 15- 39

26

haviour continues up to a point where steadystate could no longer be obtained numerically. After performing a stability analysis of the system of conservation equations, it was found that this was the point at which the system of transport equations, Eqs. (50), (64), (65) and (66), which were decoupled from the two-fluid model became unstable (i.e. the phasic momentum equations were not needed because fully developed flow was assumed). A more detailed description of the stability analysis of the coupled system is presented later by relaxing the assumption of fully developed flow, and using the phasic conservation equations in Park's (1992) two-fluid model. The fact that an instability has been obtained by introducing the evolution of the interfacial geometry of the system, and the physics of bubble cluster break-up and formation, implies that this model may be capable of predicting a bubbly/slug flow regime transtition: we prove that this is indeed the case shortly. The steady-state interfacial area density A~", mean bubble radius density R~' and mean bubble radius Rb can be non-dimensionalized as •47'

=

n~'=

A ,,,,t. ' ~

7

7

Rh = r It.~

(8o)

(81)

(82)

The parameter ,41" is plotted against ~ in Fig. 6. As shown in this plot, this moment of the probability distribution function has similar behavior to the number density, but it has a maximum at a different location. In particular, Fig. 6 indicates that the interracial area density reaches a maximum at about ~ = 0.07, after which it decreases with ~. This is exactly the type of behaviour one would expect during flow regime transition from bubbly to slug flow. The steady-state solution for /~b is ~hown in Fig. 7; it increases relatively slowly for small and at ~ = 0.076 the bubble clusters start increasing their size dramatically up to a point where |he system becomes unstable. Again. this is in agree-

0.40

O°:~'

M

A7

0.20

0.10

0 00

0 02

0 04

0 06

0 OG

0 10

& Fig. 6. Steady-stale interfacial area density ,)1%

1,00.

0,80 ¸

0,60

0.40

0,20

o.%~

•

............................ o.~ 0.04

'6.~g

6.~

....

~'.

a

Fig. 7. Steady-state mean bubble radius R~.

ment with the behaviour expected during bubbly/ slug flow regime transition. In order to complete this model, it is necessary to quantify the unknown parameters used herein, 7 and e~, in terms of properties of the flow. This has been done empirically using a special airwater experiment eKalkach-Navarro, 1992). The results obtained for ;, are shown in Fig. 8, in which 7 is ploited against the local void fraction

27

S. Kulkach-Navarro et aL/ Nuclear Engineering and Design 151 (1994) 15-39

~3/2U2

r*_,o = 31.09 ~

0.0 ~

11.3

~/3D~/3

(84)

where pl is the density of the liquid, a is the surface tension, el is the energy dissipation and Db is the average diameter of the bubbles. For other liquids, one must multiply the righthand side of Eq. (84) by the ratio of the critical thickness of water divided by the critical thickness x; of the other liquid that is, multiply by (hc.,_o/ hcx) ~/2. For example, for the mineral oil used by Park (1992), we have (hcH,o/hcoil) I/2= 2.77. It is convenient to define a non-dimensional critical volume as

48.0

=

Fig. 8. The variation of ~, with ~ and r*.

and the non-dimensional liquid film drainage time 3" which is defined as T T* = -

g(p,

(85)

pg)_l

The dependence of this non-dimensional critical volume on void fraction and the non-dimensional drainage time is shown in Fig. 9.

(83)

T

The length of time T that the bubble cluster remains intact before being knocked apart by turbulence was modeled using the expression developed by Levich (1962) and used by Thomas (1981), and the liquid film drainage time was modeled using the model proposed by Doubliez (1991). The resulting expression for the nondimensional liquid film drainage time for the liquid entrapped between adjacent air bubbles is given for water by

3. Stability analysis and bubbly/slug flow regime transition predictions The conservation equations for the moments of the distribution function obtained in previous sections are now coupled with the phasic mass and the momentum conservation equations and the stability of the resulting system is analyzed. The one-dimensional, adiabatic, incompressible form of the mass and momentum conservation equations for the individual bubbles for

3.0

2.0 0.0

/ A LYo" Fig. 9. The variation of t,, with :~ and r*.

28

S. Kalkach-Nacarro et al. / Nuclear Engineering arm Design 151 (1994) 15 39

monodispersed two-component bubbly flow can be expressed as (Park, 1992) (~lPl) + ~ (~tPlUl) = 0

(86)

& (~qp~) + ~ (~lp~ug~) = 0

(87)

c~p,

~,

R~ c~,

.~,.d, + Mi d' --4 ~ + ~qp~g cos 0 + ,,,~

( 88,

I r,,. = ~_f,.,pdt~u,

Tgw = ~.fgwPg//g I llg I

+~

~APli

C2

~

~l ]

Dt

~q / Dt

&trj , ~] = B 2 u r l - - ~ - B l u h t3--

g¢ ~ , -- T hi 0.~

4-~

(89)

where the gas bubbles (g) are the dispersed phase and liquid (!) is the continuous phase, 0 denotes the angle between tne axial direction z and the gravity vector, and the superscripts (d) and (nd) denote the drag and non-drag forces respectively. The constitutive (i.e. closure) relations for this one-dimensional system of equations are as follows (Park, 1992):

1 CDAI'~ U~I 8 ~l~l

2

+

[~,(r~, + ~ + ~)1

+ ~ , p g g c o s O -'~("d',..,, - M ~ d'

(98)

where a v , . = [ D v u v / D t - D ~ , q / D t ] is the virtual mass acceleration, u. = u,,, - u , is the relative velocity of the bubbles with respect to the liquid phase, and Po is the local liquid pressure if the bubbles are not present. The values of the closure constants for spherical bubbles are (Park, 1992) Cvm =0.5, Cm, =0.1, Cm2 =0.1, Cp = 0.25, C~0.2, C~ = 0.3. The two momentum equations, Eqs. (88) and (89), can be combined to eliminate the pressure gradient:

(:~1 Pgltgt ) + ~22 (~1 pgUgl)

~?Pl

(97)

1

(99)

where p * = Pg/Pl, and BI and B2 are given by (Park, 1992):

(

B, = 2 Cp~ + Cr :~ - Ci + -

Bj = ( 2Cp

~t

(100) 2~j J

(2 "~-(XI Ci ~l ) C r ~- oq

Cm2~ ~lcq /

(!01)

The system composed of Eqs. (86) (87) and (99) can be recast in matrix form as D -~~ - + E ~ - = f

(102)

p~a~,. -- Cml ~1 tOI/drl 0url

~i

where ., ~ t I

-- CmzplUr~ 8-M.~

(90)

3 CD

li ~ - -8- - /~t, p~g~u~

(91)

P~ = Po - Cp p~~t~u ~

(92)

Apli = - CpPl~lli~l

(93)

as = -- Cifllu~l

(94)

~,g~= _ C , p , -~'- u .2

(95)

r ~R~ = - C~p~u~

(96)

\ ui The characteristics r of this system can be found by solving det(Dr - E) = 0

(103)

For the case of monodispersed bubbly flow, Eq. (88) yields two void wave characteristics (Park, 1992), which are shown in Fig. 10. These two characteristics correspond to a slow highly damped

S. Kalkach-Navarro et at. [ Nuclear Engineering and Design 151 (1994) 15-39

number density of the individual bubbles is transferred to the void fraction and number density of bubble clusters. In a similar way, when bubble clusters are broken up by turbulence, the void fraction and number density are transferred from the bubble clusters to the individual bubbles group. All these effects are included in the righthand side term in Eq. (104). From the momentum equation for the individual bubbles developed by Park (1992), Eq. (99), and using the continuity equations to eliminate the liquid velocity u~ from this equation, we obtain

Well-posed Region -1

~3

y

I

,

t

I l 1 8

-2

-3 ,' o.o

, 0.5

1 1.0

0.5 0.6 0.5

0.25 0.25 1.0

-0.2 -0.2 -0.2

-0.3 -0.3 -0.3

-0.1 -0.1 -0.1

eu~,

~ ~,

c~u~,

+ ' ' 3 ~7. + A 4 t ~

=A 5

(105)

A~ =

1+

A~ =

1+

(106)

-0.3 -0.3 -0.3

void wave and to a faster less damped void wave. The system was found to be well posed up to the void fraction at which the two eigenvalues intersect, which if the constitutive values for spherical bubbles are used, is • ~ 0.13 for air-water flows at standard temperature and pressure (STP). Interestingly, the onset of ill posedness is strongly dependent on the closure parameter Cp. As seen in Fig. 10, an increase in Cp produced a considerable enlargement of the well posed region. As discussed previously, when the appearance of bubble clusters and the dynamics of bubble cluster break-up and formation is taken into account, the one-dimensional continuity equation for the individual bubbles is given by, Eq. (66) as • , + ~ (at Ugl ) = S,t

~, At-~-+A2-~

where

Fig. 10. Systemcharacteristics(Park, 1992).

I 2 3

29

(104)

where, for simplicity, the right-hand side of Eq. (66) has been expressed symbolically as S~. This conservation equation if similar to Eq. (87) but it includes a source]sink of void fraction of the individual bubbles due to break-up and bubble cluster formation. When individual bubbles group together to form a cluster, the void fraction and

- B,

(108)

+ Btu~t

(Xl / CXl I

-- B2url -- Ba url oq As-

CD

4

(109) 2

2

8oqoqA;'~u~I +-D-HH(fwtut --fw~ugO

(110)

The interracial area density of the individual bubbles is obtained as the 3th moment of tile probability distribution function, and approximated using Eqs. (59) and (60):

Ai't =(367t)l/s6~/SN';'=(367t)t/sazl/3(NT')t/s

(1~1)

Thus, the interfacial area density of the individual bubbles depends on their void fraction and number density, and to be able to close the system it is necessary to include the conservation equations for the number density of individual bubbles and clusters and for the void fraction of the clusters, Eqs. (50), (64) and (65), respectively: --= + ~t

S~t ~ (~,u:_) = -

(112)

S. Kalkach-Nm~rro et al. / Nuclear E.gineering and Design 151 (1994) 15-39

30 30

50 o

faster celerity

• slower celerity 40 !

20

~ 30

Uclusters0

(cm/s)

C a

-

ut (~n/s)

10

Bubbly I

O

Transition

Slug

I

I

0.1

p

0.2

l

Fig. I I. Rise velocity o f bubble clusters in air-oil mixture (Park, 1992).

5N['Ot + ~

(N'n"%,) = SNI'

~Ng' @t" + ~ (N;'%,)

u~ = 0.01 =

SN~

(114)

(115)

where, for air-oil flows, A = 0.04385 m s - ' , and B = 10.73. When the void fraction was increased further, Taylor bubbles appeared; for air-oil flows Park (1992) found that their rise velocity u, was well predicted by (Cliff, 1978)

o'.,

o12

&

o.;

Fig. 12. Void wave speed in a i r - w a t e r speed in a i r - w a t e r mixture vs. :~ (Park, 1992).

(113)

where, SNT and S~.~. are the right-hand sides of Eqs. (64) and (65) respectively. Park (1992) measured void wave celerities in air-oil bubbly two-phase flow. An important feature of his results is that he was able to measure the void wave related to the bubble clusters. The rise velocity of the bubble clusters measured by Park is presented in Fig. 11 as a function of the void fraction. As shown in this figure, the rise velocity increases as the void fraction increases up to a value of void fraction where the transition to slug flow occurs (~ ~ 16% in air-oil flows at STP). When the Taylor bubbles are formed, the rise speed of these large bubbles remains fairly constant with void fraction. The correlation proposed by Park (1992) for bubble clusters was %2 = AeB~ + ul

olo

I

0.3

gD~(Pnltl

pg)

(116)

The formation of bubble clusters and their velocity was not measured directly by Bout6 (1988) because, owing to the relatively low viscosity of water, the dynamics of bubble cluster formation and break-up is much more rapid than in .*,he case of air-oil flows. He did detect, however, the presence of a second void wave for :~ > 0.255, which he called a mode-4 wave. The speed of these two waves is shown in Fig. 12. For a k water flow Bour6 (1988) proposed that the vapor bubble velocity in the slug flow regime was well predicted by u~ =

0.369[g( 1 - p*)Od] ,/2

(117)

(1 - ~ )

Given the fact that very little is known about the drag laws for bubble clusters, the cluster velocity given in Eq. ( I ! 5) was used for the analysis that follows. Using the fact that, as expressed in Eq. (115), the velocity of the bubble clusters depends on the void fraction in Eqs. (i12) and (114) respectively, we obtain ~-T +

u~; ~ + ~, ~

0z

=

s,~

(118)

S. Kalkach-Navarro et al. / Nuclear Engineering and Design 151 (1994) 15-39

and

31

teristics (i.e. eigenvalues, 2) call be found by solving

ON~'ot+ ug2 -~zON~;'+ N2" ,,, Oug 2 ~ 0(Oq0z+=2)-S~;._ (119)

det(B - 2A) = 0

(121)

where the derivative of the velocity of the bubble cluster with respect to the void fraction can be obtained using. Eq. ( 1 ! 5). The system that includes Eqs. (104), (105), (118), (113) and (119) models the dynamics of the flow including the effects of bubble cluster breakup and formation and represents a distinct improvement over previous two-fluid models; we refer to this model as the extended two-fluid model. If we express this system of equations in matrix form, we obtain the following:

Two of the eigenvalues obtained this way (A,~ and Aft) correspond to the characteristics (r f and r~- respectively), since the new terms due to bubble cluster formation and break-up, and the presence of the distrubution of bubble sizes, are algebraic and are thus contained in vector c. This implies that the well posedness of the system will be the same as that of the two-fluid model given in Eq. (102). Indeed, the other three characteristics are always real and are given by 23 - 2n,i. = us1

(122)

A 0t

24 ----2N.~= ug2

(123)

~

=e

(120)

~l/g 2

where

2s = 2~2 = %2 + a2 0-a

I:l

The first eigenvalue corresponds to the number density N'(' of the individual bubbles, the second corresponds to the number density N~', of the bubble clusters and the last corresponds to the void fraction ~2 of the bubble clusters. The characteristics can be non-dimensionalized

gl

°= J; "i (N 'J C

=

A5

as

--S:d

2" = (2 - ul)/ur

SN,i. SNS: Ai

A=

,42 0 0 0 1 0 0 0 1 0 0 0

[i ° °°i] u~x /I 3

B=

(124)

Bug2

a2 ~

0 ., ~//g2

t~ N2 ~

at A4 0

N;' 0

0

0

0

0

0

0

0

0

u~l

0

0

ug2

u~2 + ~

~Ug2

o ~Ug2

N~' O~

It should be noted that the first two equations ~; the array, in the absence of bubble clusters and assuming only the bubble size, reduce to the model used by Park (1992). The systems charac-

(125)

Typical behaviour of the non-dimensional characteristics as a function of the void fraction a t of the individual bubbles is shown in Fig. 13. The continuous lines correspond to the characteristics obtained using the extended two-fluid model given in Eq. (120) and the dashed lines correspond to the characteristics of the two-fluid model given in Eq. (102). As expected, two of the characteristics obtained for the present model are identical with those obtained by Park (1992). There exists an important difference in behavior, however, since as shown in Section 2, the steady-state void fraction of the individual bubbles starts decreasing when bubble clusters are formed, in such a way that the extended two-fluid model never reaches the value at wh'ch it would become ill posed. Fig. 14 shows the same characteristics but as a function of the total void fraction a. In this case we can observe that for small void fraction the

32

S. Kalkach-Nat'arro et al. / Nm'tear Engineering and Des(gn 151 (1994) 15-39

The value of 2*.r remains constant with respect to ~ since the relative velocity of the individual bubbles is constant, while 2"~ increases exponentially with void fraction, since that is the behavior of the relative of the bubble clusters. We note that 2*2 increases with respect to void fraction faster than 2",~ since it additionally includes the derivative of the relative velocity of the bubble clusters with respect to void fraction.

2.~0 -

X'-r ~

a_-

A"

'x /

L#war analysiss of void wave propagation

O.~

0.e4

O.~IE

....

0'. 112' ' '

e.16

Ct1

Fig. 13. Characteristics of the extended two-fluid model vs. ~t.

When the two-phase system represented by Eq. (120) is perturbed about a fully developed steadystate situation, the behavior of the perturbed variables can be analyzed. Let us define Xo as the steady-state value of an independent variable and ~Jx as a small perturbation. Thus, the value of the perturbed variable is given by x = xo + Jx

)

(126)

If we apply this perturbation to each of the independent variables of Eqs. (104). (105), (118), (113) and (119). and subtract the steady-state, the linearized system of equations that will be obtained is as follows

06~t

06~,

0-~- + %1,, ~

06u~, + ~,,, O_

= fit 6~, + P2&tgl + fl~3x2 + fix2 + f146N'l" +

e.

~6..... e'.~ ..... i¢.'io ..... e'.'¢s ..... i ~ ..... ~'.~ Ct

Fig. 14. Characteristics of the extended two-fluid model vs. ,,.

f156N~'

06~, 06%, A 06~, 06u~1 4 , - ~ - + A2--~--[-- + .~--~-z + A4 0z = flfO~X I -'1-/JTftlgl + fl83~2

+ flofN';'+ fl,ohN"' total void fraction is approximately the void fraction of the individual bubbles, thus the characteristic 2*,- of the slow void wave increases, as shown in Fig. 13, up to a point at which the void fraction of the clusters starts becoming significant and the void fraction of the individual bubbles decreases; thereafter 2",- starts decreasing also, in such a way that the characteristics of the fast waves (2",-) and the slow waves ( 2 * - ) do not cross each other. Hence, as noted above, the system ahvays well posed.

(127)

(128)

06~, Of~, 06u~2 0t + u g , , , ~ + ~2,, O: = flllf~Xl + fl123Hgl "~ fll3f~X2

+ fl,4fN] " +

[3~sfN'~'

(129)

-N"' 06u~1 0fiN';' (rN]")+ i. Oz +ugl,, 0z Ot = fllff~l + flt:fu~l + fllxf~, + fli,jfN]" + [12ofN~'

(130)

S. Kalkach-Navarro et al. / Nuclear Engineering and Design 151 (1994) 15-39

,,, 06u~2

(6N~') + Nzo ~

O6N~' 0:

4. Analysis of results and predictions of the data

+ Ug2o

=/~2,/i~, +//:2'~u,~ + ~23~2

where the constants/I,, are given in Appendix A. Now, expressing the linearized system containing Eqs. (127)-(131) in matrix form:

Ao--~- + Bo--~z = C'o6eP

33

Typical behaviour of the five damping coefficients with respect to the total void fraction obtained from the extended two-fluid model, is shown in Fig. 15. One of the roots, that associated with the slow void were C~, is much more damped than the other four. A blow-up to the other four roots is shown in Fig. 16. It can be observed that for air-oil flows the damping co-

(132)

500"O01

where

kt. t

(~c)

b-$o

The dispersion relation can be obtained by assuming a modal solution of Eq. (132) of the form 6@ = @' exp[i(k: - tot)]

. k~ . kin ~ . ~Ju;

(133)

(134)

-~I (l~-t)

- 1000.~ -*t,

where k = n/2 is the wave number, and to = 2nfis the angular frequency. In order that the system have a non-trivial solution, the following dispersion relation must be satisfied: - ~00

det[Ao-( k )Bo + iC'o/to]=O

. 00 0.

(135)

ot

Fig. 15, Damping coefficients vs. :~.

To obtain the spatial damping of the waves compare with the available experimental data, we may assume that to is real but k is complex (k = ka + ik 0. Thus, Eq. (134) can be rewritten as

10.00

0.00

6 ~ = ~ ' e x p ( - k , z ) exp[ika(: - Ct)]

(136)

where C = to/ka are the various celerities. One of the system's characteristics and the kinematic void wave speed are the infinite and the zero frequency limits of the void celerity respectively. It is seen in Eq. (136) that the stability of the waves is determined by the sign of kt; a negative sign implies that the amplffudes of the corresponding waves increase as the waves travel, and thus that the wave is unstable. The roots of Eq. (135) yield the celerities of the system. These roots were found numerically since it was not possible to evaluate Eq. (135) analytically.

-le.~

- 2 0 . ;~o -kt(mq)

- 4;~. 0 0

- 50.00

Fig. 16. k m , kl,2, ktnt, kln~ vs. a (Cr=O.25, jl = O . O m s - t ) .

S. Kalkach-Navarro et al. [ Nuclear Engineering and Design 151 (1994) 15-39

34

efficient associated with the propagating void perturbation of the bubble clusters C,2 becomes negative (i.e. amplification occurs) at a global void fraction of about 16%, while the other waves remain damped. This indicates that the void fraction of the bubble clusters becomes amplified and, as is explained later, this leads to a flow regime transition. If the void fraction is increased further, we obtain a point at which the system no longer has a steady-state solution (i.e. :~ = 23% in Fig. 16). This point corresponds to the value of ~ of 0.076 at which the uncoupled system did not have a steady state. The void fraction at which the C~_, wave becomes amplified is dependent on the value of the ratio of the bubble cluster formation and break-up rate )' and the physical properties of the fluids. Previous authors such as Bour6 (1988) and Park (1992) have measured the spatial damping coefficient. Bour~ utilized air-water bubbly flow and Park obtained his data for air-oil bubbly flow. They also measured the speed of the various void waves. Fig. 17 shows that data for the damping coefficient associated with bubble clusters as a function of the global void fraction at a frequency of 2 Hz (the dominant frequency measured by Park (1992)) for zero liquid velocity (j. = 0.0) in air-oil flows (Park, 1992). It can be observed that the v, ave is damped for a void fraction smaller than

16%, the point at which the wave has neutral stability. This value coincides with the observed bubbly/slug flow regime transition (Park, 1992). This suggests that, as previously proposed by Bour6 (1988), a propagating perturbation in the void fraction may induce an instability which leads to a bubbly/slug flow regime transition. If the void fraction is increased further, the wave becomes amplified; this means that the amplitude of the wave will increase as it travels, triggering a cascade in the formation of bubble clusters, resulting in a bubbly/slug flow regime transition. The extended two-fluid model's predictions of the damping coefficient k~ for an air-oil system with a frequency f = 2 Hz are shown in Fig. 17 also. As can be seen, these predictions are in very good agreement with the experimental data. Moreover, the point at which the system becomes unstable (i.e. where k~ = 0), and the flow regime transition was observed, was predicted very well. When the liquid velocity was increased, Park's (1992) measurements showed that the magnitued of the damping for low void fractions was decreased and the onset of the bubbly/slug regime transition was thus shifted to a smaller value of global void fraction, and the amplification of the unstable C~2 waves becomes greater. This behaviour is clearly shown by the data of Fig. 18. 3.00.

tS"~ 1

Park'sdata[1992] (Ji ==0.0. f = 2.0Hz)

o

2.~

• Park's data [1~2] (Jr =05.71 cm/s;f : 2.0Hz) o Bubbly Flow

o

• Slug Flow

)

/ -kt (ra-l)

0. ~0, "I

-~'~

.

-I . ~ -2.00 0.

..... ~' '¢~.... ""o'.~ ..... ' ~'.~' ..... ~ ' ~

u (1 Fig. 17. Prediction o f d a m p i n g for a i r - o i l C~2 wave (Jr = 0.0).

Fig.

18. Prediction

(jr = 5.71 cm s-~).

of

damping

['or

air-oil

C~,_ wave

S. Kalkach-Navarro et al. i Nuclear Engineering and Design 151 (1994) 15-39 :3.e~-

K.eo-

35

o Bour~'sData [1988] Ut = o.a6 r~e;t" = 5.0Hz)

o Bour~'a Data [1988] (j! : O . 3 6 m / n ; / : 3.SHz)

I.eO

-kl, a(m-l)

-kla 2 (m -l)

1

--.....

a /

. . . .

-1 .~0

- I .£0"

-3. e'~

-3.00

........ 0.ia

0.0o

, . . . . . . . . .

,

0.~

.......

'e'.~e.... o ~ .... o'.~o

-5.00

......... , ................... 0.00 0.10 0.20

~

,

0.4-0

G

G

Fig. 19. Prediction of damping for air-water (Jl = 0.36 m s -I, f = 3.5 Hz).

, .........

0.30

C,2 wave

Shown also are the predictions obtained from the extended two-fluid model, which again agree very well with the data. Thus, we can observe that the model is able to predict the effect of increasing liquid velocity. It is interesting to note that the two (dark) data points for global void fractions greater than 0.2 are marginally stable, but not predicted since they belong to the slug flow regime, and the extended two-fluid model for bubbly flow does not have a steady-state solution for that case. Finally, the predictions and the experimental data of Bout6 (1988) for air-water flows are shown in Figs. 19 and 20. It is very encouraging to note that the onset of the bubbly/slug flow regime transition (i.e. where k~ = 0) obtained from this extended two-fluid model agrees with the well known "rule-of-thumb" for low pressure airwater flows of ct = 30%. Park (1992) measured the void celerities corresponding to the individual bubbles and bubble clusters. The behaviour of these celerities with void fraction for j~ = 0.0362 m s- ~ is presented in Fig. 21 for a frequency of 2 Hz (the dominant frequency measured by Park (1992) for void wave propagation). The celerity of the individual bubble void wave C~ was predicted by the dispersion model of the two-fluid model given in Eq. (102).

Fig. 20. Prediction of damping for air-water (Jl = 0.36 m s - I f = 5.0 Hz).

50 f

C,2 wave

o fastercelerity • slowerceled~

40

/,,.---- gq.(1151 /

'30 Speed (ends)

o _ / , o)

/~al.(l16).

~o

20 igenvalue,3.+ (Ug-U¢)+ut

0

0

I

~-~-

•

~

I

a

I

0.1

ct

Fig. 21. Void wave speed in (Jl = 3.62 cm s -~) (Park, 1992).

KinematicVoidWave.a+

0.2 air-oil

)

I

0.3 mixture

vs.

~t

The faster wave corresponds to the bubble clusters and is represented by Eq. (115). At a global void fraction of about 16%, in air-oil there is a transition to slug flow, and owing to the appearance of Taylor bubbles the celerity becomes fairly constant and its value is given by Eq. (It6). Typical behavior of the various celerities obtained from the dispersion analysis of the exo tended two-fluid model is shown in Fig. 22. The fastest two ~vaves are those corresponding to the

36

S. Kalkach-Navarro et al. / Nuclear Engineering and Design 151 (1994) 15 39 0.2~

.

Paxk's C. Cg

_ gxtended two-fluidmodel

e" "=+ai

e. 16

a':le 1 C+]eHty Ira/s]

e'1~ !

0.12

0)1112 -e. lO-

.eo

" . . . ,'c,r -0. ~0,

e.~4 . . . . . . . . . ~. . . . . . . . . +. . . . . . . . . u.., o.ee {a.05 o.10 0.1m . . . . 6'. ~1~ o

Fig. 22. Comparison of void wave celerities from two-fluid and extended two-fluid models.

void fraction of the bubble clusters C~2 and the number density of tile bubble clusters CNS. This is reasonable since it was observed by Park (1992) that the velocity of bubble clusters was significantly larger than that of individual bubbles. The celerities obtained using the model developed by Park (1992), Eq. (102), are also presented in Fig. 22, where it can be observed that the extended two-fluid model yields two celerities (C~ and C~i ) that are very close to those obtained by Park. In particular, the show void wave C;i is almost unaffected by the new transport equations and the fast void wave C,7 is also very close for low void fractions. The only difference is that the faster void wave increases somewhat for longer void fraction, since in this case the presence of bubble clusters becomes significant. t is interesting to note that the temporal amplification factor can be obtained by assuming that k is real and co is complex (,~ =tOR + i~o~). Substituting this into Eq. (134) yields 6(!) = ~ ' exp(ohz) exp[ik(: - Ct)]

(137)

The temporal amplification coefficient ~), of the void fraction of the bubble clusters is shown in Fig. 23. As Can be seen, at a given global void fraction, to~,z becomes positive, again indicating instability of the bubble clusters and the onset of the bubble/slug regime transition at ~t ~ 16% for bubbly air-oil flows.

(z

Fig. 23. Temporal damping of c,: wave vs. ~e. According to Taylor's hypothesis, for a nondispersive wave the temporal and spatial damping are related as follows: kl = -- ~ol/C~

(138)

It was shown numerically that the void waves for the individual bubbles and the bubble clusters obtained from the extended two-fluid model verified Taylor's hypothesis (i.e. Eq. (138) was satisfied), since the extended two-fluid model is fairly non-dispersive. From the results obtained using the extended two-fluid model it is possible to conclude that a properly formulated two-fluid model, which includes the physics of bubble cluster formation and break-up, is inherently capable of predicting the bubbly/slug flow regime transition. This is an important result, since it allows the prediction of the onset of a bubbly/slug flow regime transition without the need to rely on empirical flow regime maps.

5. Conclusions and recommendations

In order to be able to predict mechanistically a bubbly/slug flow regime transition, it was necessary to model the dynamics of the formation and break-up of bubble clusters. A mathematical

S. Kalkach-Navarro et ai. / Nuclear Engineering and Design 151 (1994) 15-39

modeling approach similar to the Boltzmann transport equation was developed to include the relevant physics. This model was then combined with a state-of-the-art two-fluid model to generate an extended two-fluid model which was found capable of predicting the bubbly/slug flow regime transition. The void wave dispersion relation was developed to quantify the behaviour of the measured void wave signals at different frequencies. The spatial damping coefficients that were predicted with the extended two-fluid model were in very good agreement with the experimental data obtained by Bour6 (1988) for air-water, and by Park (1992) .for air-oil bubbly flows. Significantly, the predicted void fraction at which the void wave associated with the bubble clusters became unstable agreed with the observed bubbly/ slug flow regime transition. Another important result of this work is that it allows one to model the interracial area density and mean bubble radius which are required for closure in two-fluid models. Acknowledgements

The authors wish to recognize the importance that the pioneering work of Dr. Novak Zuber has had on this and many other technical papers. Dr. Zuber was one of the first to attempt rigorous analysis of multiphase systems and his effects helped define the field and have inspired many students and researchers to push the frontiers of knowledge in muitiphase science. Novak, from all of us within the Center for Multiphase Research (CMR) at Renssclaer--Good health and long life! Last, but not least, we wish to acknowledge with gratitude the financial support given this research by Dr. Oscar Manley of the USDOEBAS.

37

J. A. Boure, Properties of kinematic waves in two-phase pipe flows consequences on the modeling strategy, Proc. European Two-Phase Group Meet., Brussels, 1988. R. Cliff, J. R. Grace and M. E. Weber, Bubbles, Drops and Particles, Academic Press, New York 1978. L. Doubliez, The drainage and rupture of a non-foaming liquid formed upon bubble impact with a free surface, Int. J. Multiphase Flow, 17(6) (1991) 783-803. T. Z. Harmathy, Velocity of large drops and bubbles in media of infinite or restricted extent, AIChe J. 6(2) (1960). S. Kalkach-Navarro, Mathematical modeling of flow regime transition in bubbly two-phase flow, PhD Thesis, Rensselaer Polytechnic Institute, Troy, NY, 1992. J. R. Lamarsh, Introduction to nuclear reactor theory, Addison-Wesley, Reading, MA, 1972, 2nd edn. V. G. Levich, Physicochemical Hydrodynamics, Prentice Hall, Englewood Cliffs, NJ, 1962. J.-W. Park, Void wave propagation in two-phase flow, PhD Thesis, Rensselaer Polytechnic Institute, Troy, NY, 1992. J. M. Saiz-Jabardo and J. A. Bour& Experiment on void fraction waves, Int. J. Multiphase Flow, 15, (1989) 260-493. K. Sangani and F. Didwania, Dynamic simulations of flows of bubbly liquids at large Reynolds number, J. Fluid Mech. 233, (1991) 65-69. H. Schlichting, Boundary Layer Theory, McGraw-Hill, New York, 1979, 7th edn. J. Smereka and K. Milton, Bubbly flow and its relation to conductio, in composites, J. Fluid Mech. 250, (1993) 307321. R. M. Thomas, Bubble cealescence in turbulent flows, Int. J. Muhiphase Flow, 7(6) (1981} 709-717.

Appendix A: linearization of the extended two-fluid model

The extended two-fluid model given by Eq. (105) can be linearized to obtain the system of linear equations comprising Eqs. (I 12)-(116). The parameters flj used there are given by: fl, = -OS~l -~

o

(AI)

f12 = OS~l%, o

(A2)

o f13 = OS~l ~

(A3)

References

0S~i 114= ~ o

(A4)

Z. Bilicki and J. Kestin, Transition criteria for two-phase flow patterns in vertical upward flow, ha. J. Muhiphase Flow, 13(3) (1987) 283-294.

115 = ON~--~ o

(A5)

S. Kalkach-Navarro et al. / Nuclear Engineering and Design 151 (1994) 15 39

38

~As

0As

Io fl25 = 5SN~' ~

(A7)

Now, substituting the expressions for S,I, SN'C, SN~ and As, we obtain the equations for //~ in terms of the steady-state variables:

(A8)

fll= _( 16~tS/aN,~,2/3+ ° /.... N,~,-213 ~ V 1 "51~ " 2tit"+ otsff3

(A25)

+ 8/3a~/3N,~,-I1.~i-i/3 + 14]3o~/3ot~/3N,~,2/3

~As I

(A9)

~A~ I0~

0

(AI0)

"{- ,~_2/3~r,-I/3,-; : ~2 ' v 2 ~1

o

fl,-, = ~

+ 5]3t~/30G)

(A26)

//2 = 0

(A27)

13

OS,d .o,,

(A6)

~ ) / ~ . 11/3~- 3

O~ 14/3~-4X

(All)

+ 8/3g~/3:t~/3N,~,-i/3+ 2/30t7/3N,:,-4/3~22/3

(AI2)

+ 2/36~u3~t~N': '-I + ~5/3N'~'-2/3

(A28)

~4 = (~/3N'1"-5/'3 + 5/3f~/3N~ ' + 2]3~'~/'3N~'-1/3~ 4/'3

0S, I

fl

13 ~"'="~2"0~2 O

8

0S, i

14~---~

0

(AI3)

_t_8/3~7/3~tl/3Na¢3 2 +~22/3 N 1'/'3-2 -5/3 2VI-I-2/3VI"0~2)

(A29)

(AI4) fls = 127(v~l/3f2~

0S, i

(AI5)

15~--'~

0SN;-

.a,,, = ~ 0 ~SN'I' fl19 = ON-~(' o OSu.,.

(AI9) (A20)

u2 !

--

~I(Ot~/'3N,1,,-2/'3

+ ~/..,w2 A~a.4/'3 ~v a! .... ± l/3f22/'3~t) I u3.l/a u2 T

16 ~t(1 _ ~ g ) 2

f17

4~1(1 - - ~ 1 ) 2

(A31)

COA[['frl

(A32)

/~8 = 0

(A33)

/~,=0

(A34)

fllo = 0

(A35)

fll,=(16~/3N'('-z/3+8/3~5/3N'~'+a2. OSN'~, =

(A21) o

~SN~'

(A22)

(A30)

7 coAi':-2~l

f12 =

(AI7) (A18)

vc

2/3 r,s/3 m,,, ")171a4/'3.I/'3M',,,2/'3 02 "tVl --'~l'~U2 ~"1 JVl

-

(A16)

0SN'r

,~1413~-3X

. + 8/3~/3N'~'-U3v]/3 +

+--~-2 ,,2

,v 2

5/3 ~ , 7 . , - 2 / 3

1/11"1 ~4/3N I/'3 ATm2/3 ,-'t-'~l

~'2

"2

-i +

(A36)

fit2 = 0 i~SN~ 0SN~ fl24 = ON--~" o

(A23)

ill3 = T(8UcI I/3/~23

-

(A37) --

90/4/'3~24) + 5/3~2/30g

+ 8/3g~/'3 ~4/'3N~"-I/3 + 2/3~ 7/3N'('-4/'3 t3~-2/'3

(A24)

5/'3 ,, M ' Is + 2/3621/'3~NS; '-I + ..,, ~m

2/3

(A38)

S. Kalkach-Navarro el al. / Nuclear Engineering and Design 151 (1994) 1 5 - 3 9 ill4

=

ff19

--( 4°t~/3 N'l '- 5/3 "t- 5]3v~/3 N'2 ' •~/'4,f,4/3Mm-I/3U4/3 ltl'li~7/3,vl/3Mm2/3 "Jt- :"l"~2 ~t2 v t di-Ul'-'u2 ~'2 ~v2 -]- ~213 N~'1/3~2 d- 2/3t3~5/3~ 2)

" ) / ~ 5 / 3 Jv Arm - - ~/'% 1 v5/3M ~v!" - 2 / 3 --~.,/..Jt,.i 2 + ot5/3N'2~-2/3 + 4 / 3 t X 4 / 3 N ~ ' - 1 / 3 ~ / 3

(A39) --

ill5 = -- 12~(Vc~/3022 -- I)14136 2 3) _

~6 m2/3

+ 2/3~t~/3N~ '- ~i317~-2/3 + 8[3fl/30t~/3 x N~ '2/3 q- ~2/3N~'~/3 -t- 2/3~'~/3ot2) fl~7 = 0 -

(A41)

"c

"JCu/'Ju2~/~'~I/3N A 1I/3 " 2~/ 3' 2~ 1

~2 , + 5[3~22/3N'( '

(A46)

fl22 = 0

(A47)

fl23 = ? 12(v~/3v2 - 3 _ ~"~~/3,~-4~2, -20]3ot~/3N'~ '1/3 (A48)

fl24 = 2/3~/3 N~ "-2/3

(A49)

fl25 = yv2/3( 1 -- 18vc2v22+ 16vc3~ -3) --4[3~/3N'~ '-2/3

-~- 213v~2/3~t~/3NI "-~/3_

+ 2[392 ~/3~ + ~2/3 N~,,I/3)

(A45)

[$21 = lO[3ot2/3 N'( '1/3

(A42)

.

(A44)

3) - (~:3N~"2/3

+4[3ot4/3N ..... 1/31/~ I/3 "3L 1/3622/3~Xl) (A40)

-2/3 ~r,,, --(lOot~/3N'~'n/3 + 5[31)! ~v 2

-

2 - 16v~u3~

1/3~2/30t2)

' ) / ~ -5/3 AT" "~11i;413N 113AT'213 - - /"IJU2 ~V l m ~l"U2 ~'1 ZVl

+ 1/3~2/3cg )

fl~6

'>/'1~4/3N M-,2/3 ~_ l ~'2I/3 1v2 T

~"l"u

/~20 = r(18v~/3~

~ (~/3 N .... 2/3 -- 2/31)2-S/3N~"

- - 2")/'1'74/3"~'11/3 ~'/'t' ~VlM"2/3 .3t_4/3ot4/3 N,(, -

39

(A43)

(A50)