Two-phase bubbly flow through vertical tube: Void-fraction distribution and velocity profiles

Physica A 192 (1993) 410-442 North-Holland

~

/

~

Two-phase bubbly flow through vertical tube" Void-fraction distribution and velocity profiles A.J.N.

Vreenegoor 1 and J.A.

Geurst 2

Department of Mathematics and lnformatics, Delft University of Technology, P.O. Box 356, 2600 AJ Delft, The Netherlands Received 10 August 1992

The macroscopic theory of two-phase bubbly flow developed recently by one of the authors is extended with a dispersive force that models the small irregular motions of the bubbles which, in the case of vertical two-phase flow, come into prominence near a solid boundary. It is demonstrated by means of analytical and numerical calculations how the dispersive force and the "lift" force already contained in the original theory combine to produce void-fraction distributions and velocity profiles for vertical pipe flow which are in good agreement with the experimental results obtained by Serizawa et al., Wang et al. and Kapteyn. In particular the "turbulent" character of the velocity profile for the liquid phase is reflected by the computations. In conformity with the experiments the void-fraction distribution displays a distinct maximum in the vicinity of the wall both for the cocurrent upward flow of a bubbly liquid/gas mixture and the collective rise of bubbles in a stagnant liquid.

I. Introduction It is a l r e a d y k n o w n for a long t i m e t h a t the b u b b l e s in t h e vertical u p w a r d flow o f a b u b b l y l i q u i d / g a s m i x t u r e t e n d to c o n c e n t r a t e n e a r the wall o f t h e t u b e (see refs. [1-4] a n d r e f e r e n c e s c o n t a i n e d in t h e s e p a p e r s ; similar p h e n o m e n a w e r e o b s e r v e d in a s y s t e m o f c o n v e r g i n g a n d d i v e r g i n g t u b e s e c t i o n s by d e J o n g [5]). S e v e r a l a t t e m p t s h a v e b e e n m a d e to e x p l a i n t h e o b s e r v e d v o i d f r a c t i o n d i s t r i b u t i o n by m e a n s o f l a t e r a l forces acting on the b u b b l e s . S o m e e x a m p l e s o f t h o s e forces a r e t h e s h e a r " l i f t " force c o n s i d e r e d b y Z u n , R i c h t e r a n d Wallis [6], t h e t u r b u l e n t p r e s s u r e in t h e liquid t a k e n into a c c o u n t by D r e w a n d L a h e y [7] a n d t h e n o n - d i s s i p a t i v e " l i f t " f o r c e u s e d r e c e n t l y by W a n g et al. [4] in a d d i t i o n to the R e y n o l d s stress. U p till n o w , h o w e v e r , d e s p i t e v a r i o u s a t t e m p t s to a n a l y s e vertical p i p e f l o w t h e o r e t i c a l l y (a concise s u r v e y with Present address: KSLA (Shell Laboratory Amsterdam), P.O. Box 3003, 1003 AA Amsterdam, The Netherlands. 2 Present address: Malvalaan 29, 5582 BC Waalre, The Netherlands. 0378-4371/93/$06.00 (~ 1993- Elsevier Science Publishers B.V. All rights reserved

A.J.N. Vreenegoor, J.A. Geurst / Two-phase flow through tube

411

references is given in ref. [8]) no satisfactory prediction of the distribution of the bubbles over a cross-section of the tube could be given by starting from first hydrodynamic principles (see the short review of bubbly flow in ref. [9]). It is one of the purposes of this paper to present such a derivation of the void-fraction distribution by using a properly extended form of the general theory for two-phase bubbly flow developed by Geurst [10]. In ref. [11] the two-phase flow equations derived in ref. [10] are examined in more detail. It is shown there, by going to the limit of vanishingly small values of the void fraction, that the theory implies the following equation for the non-dissipative motion of a single gas bubble through liquid: du 2 m 1 ~=(l+m,)

(O .V)u, ~-~+u, -ml(u 2-u,)×(V×u,).

(1)

In this equation u 1 represents the velocity of the liquid in the absence of the bubble, u 2 equals the velocity of the bubble when moving through the liquid and m I denotes the virtual-mass coefficient of the bubble (m 1 = 1/2 for a sphere). Note that, for mathematical convenience (see the subsequent sections), the subscript i is used in quantities like the velocity u i to refer to the liquid phase when i = 1 and to the gas phase when i = 2. The volume of the bubble is taken constant and the mass density of the gas is disregarded compared to the mass density of the liquid. External forces like the buoyancy force have been omitted. The last term at the right-hand side of (1) shows that an inertial "lift" force is obviously contained in the theory of Geurst [10]. That may be understood in the following way. According to the analysis in ref. [11] the virtual-mass terms appearing in (1) should have an objective, i.e., material frame indifferent character. The "lift" force apparently equals minus the Coriolis force experienced by a bubble when its motion is considered with respect to a frame attached to the locally rotating liquid. The "lift" force therefore contributes to making the acceleration of the bubbles in (1) material frame indifferent [12]. The equation of motion (1) may accordingly be written in the form mlare

I =

,9

~-~ + U 1

.V) U 1 ,

(2)

where are I represents the acceleration of the bubble relative to the rotating liquid. See section 2, where the theory of Geurst [10] is reviewed in a form suitable for the purposes of this paper. The "lift" force is discussed extensively in recent papers by Auton [13], Drew and Lahey [12] and Auton, Hunt and Prud'homme [14]. In addition to the inertial "lift" force and well-known forces like the viscous

A.J.N. Vreenegoor, J.A. Geurst / Two-phase flow through tube

412

stress of the bubbly liquid/gas mixture, the buoyancy force, the Levich-Moore drag force and the dissipative Fax6n force (see section 4), some other force appears to be required in order to explain the fact that the maximum of the void-fraction profile in cocurrent upward flow occurs at a certain distance from the wall of the tube. At the location of a solid boundary the probability to encounter a gas bubble vanishes. The void fraction accordingly assumes a zero-value at a solid wall. It is known from experiments that very large gradients of the void fraction occur in the immediate vicinity of the wall. The corresponding characteristic lengths appear to be of the order of magnitude of the diameter of the bubbles. The condition required for the application of a macroscopic theory, viz. local equilibrium of the distribution of the gas bubbles in a volume element, seems therefore violated in a relatively thin layer adjacent to the wall. It is possible to cope with that situation by admitting additional terms in the energy density of the two-phase mixture that depend on the gradient of the void fraction. We accordingly extend the theory developed in ref. [10] by introducing the following additional contributions AF and AK to, respectively, the free energy density F and the kinetic energy density K: AF=

1

2

~gllPe(Va)

2

,

A K = ½pe[(mll 3 - m ± ) w w + m ; w21]: V a V a .

(3) (4)

Here a denotes the fraction of the volume occupied by the bubbles (void fraction), w -- u: - u 1 represents the drift velocity of the bubbles with respect 2 3 3 to the liquid, while g~lPe, p t r n , and p~rnll are relatively small material coefficients which characterise the bubbly mixture. "Relatively small" here means that the expressions (3) and (4) contribute substantially only when Va has the order of magnitude of the reciprocal of a bubble radius. Since the expressions for AF and AK model in a macroscopic way deviations of the distribution of the gas bubbles from local equilibrium, the right-hand side of (3) may be attributed partly to the entropy portion - T S of the free energy density F. A distribution of the bubbles, which deviates appreciably from equilibrium as a result of a large gradient of the void fraction, has a relatively small probability, in particular at large values of the void fraction due to the excluded volume effect. The corresponding contribution to the entropy density is therefore negative. The major part of AF, however, comprising the energy associated with the small irregular motions executed by the bubbles during steady flow, has to be attributed to the internal energy portion U of the free energy density. The irregular motions will be accompanied by capillary waves on the surfaces of the bubbles. The coherence length associated with these fluctuations will have the order of magnitude of a bubble radius. The

A.J.N. Vreenegoor, J.A. Geurst / Two-phase flow through tube

413

small random motions, which come into prominence in the neighbourhood of the wall because of the locally modified bubble distribution, have a dispersive effect on the macroscopic distribution of the bubbles (see section 5). It will be plausible from the preceding remarks that gl~ ~>0. As a result of the large gradient of the bubble concentration in the vicinity of a solid boundary the local bubble distribution is no more symmetric with respect to an imaginary plane that is locally parallel to the wall. A bubble in the neighbourhood of the wall is consequently subject to a net non-vanishing Bernoulli force which acts in a direction away from the boundary. In a similar way the atoms and molecules at a free surface of a liquid or a solid body experience a net non-vanishing van der Waals force directed towards the interior of the body. It will be clear that the net Bernoulli force should follow from an additional contribution to that part of the kinetic energy density of a bubble dispersion that is associated with the local "backflow" of the liquid around the bubbles. When terms at most quadratic in the gradient of the void fraction are taken into account, that additional contribution is represented by the right-hand side of (4). The expression for AK implies that the virtual-mass coefficient of a bubble dispersion depends on the gradient of the volume density of the bubbles. In view of (4) it takes the form of an anisotropic tensor given by [m(c~) + p~m±(Vot)z]l + p~(mll - m±)Vo~Vc~ .

(5)

The virtual-mass coefficient for uniform two-phase flow has been denoted by Experiments indicate that the dispersive action of the small irregular motions of the bubbles in the neighbourhood of the wall is probably the most prominent effect in the vertical two-phase flow of a bubbly liquid/gas mixture (see e.g. ref. [4]). It will therefore be sufficient for the purpose of this paper to consider only the modification brought about in the two-phase flow equations by the addition of expression (3) to the free energy density (see section 3). The possible importance of an expression like (3) for describing the macroscopic behaviour of a two-fluid system in the immediate vicinity of a solid boundary was recognized more than three decades ago by Ginzburg and Pitaevskii [15] in the case of superfluid 4He (He II). In He II the superfluid mass density, which vanishes at a solid boundary, attains its bulk value within a few atomic distances from the wall. This so-called "healing" effect was also studied by Hills and Roberts [16] and Geurst [17]. Expression (3) may play a role in smoothing out the sharp transition in a discontinuous void-fraction wave (concentration wave) (Wallis, private communication, 1986). The effect of dispersive forces on the propagative of void-fraction waves is now under

414

A.J.N. Vreenegoor, J.A. Geurst / Two-phase flow through tube

investigation. It is important to realise that already at the turn of the century van der Waals [18] and Korteweg [19] introduced terms including gradients of the mass density in the expression of the free energy density in order to arrive at a general theory of capillarity in which the changes of the mass density proceed continuously. Reynolds stresses will not be considered in this paper. In fact, the tendency of bubbles to concentrate near the boundary in upward cocurrent flow has been observed under laminar as well as turbulent conditions as mentioned by Beyerlein et al. [3]. Even in the case where gas bubbles are allowed to rise in a stagnant liquid, a maximum of the void fraction is found close to the wall of the tube [20,21]. It will be shown, however, in sections 5 and 6 that although Reynolds stresses are not taken into account, velocity profiles for the liquid phase are obtained which look very similar to the velocity profiles known from one-phase turbulent flow.

2. Equations governing non-dissipative bubbly flow The non-dissipative equations of motion for a bubbly liquid-gas mixture were derived in ref. [10] from a generalised form of Hamilton's principle of least action. Special attention was given to the virtual mass of a bubble dispersion. Hamilton's principle was formulated in terms of Euler coordinates. A similar variational principle formulated in terms of Lagrange coordinates was discussed in ref. [11]. Surface tension was not considered in these two papers. It was later taken into account by Geurst and Vreenegoor [22] in their analysis of flow induced bubble deformation. We review in this section the macroscopic theory just mentioned in a form suitable for the purposes of this paper. The reduced densities p~ and P2 of the liquid and the gas phase of a bubbly liquid/gas mixture are defined by p, = p~(1 - a ) ,

p2 -- p , - .

(6)

It is assumed that the liquid is incompressible (constant mass density Pc) and the gas satisfies the ideal gas law (pg--(RT/M)pg). Note that the reduced density Pi refers to the liquid and the gas by means of, respectively, i = 1 and i = 2. The volume density of the bubbles (void fraction) is denoted by a in accordance with the notation in section 1. The conservation of mass for the two phases and the conservation of the number of bubbles (bubble number density n) are expressed, respectively, by

A.J.N. Vreenegoor, J.A. Geurst / Two-phase flow through tube

OP----2+ V . (p,u,) = 0

i = 1, 2

Ot

'

415

(7) '

On + V.(nu2) = 0 . 3t

(8)

- -

It is usual to assume that the flow proceeds isothermally. The free energy density F*(p~, n) of the bubbly liquid/gas mixture may be decomposed into the sum F(p~) of the free energy densities of the two phases and the interracial energy density of the bubbles in the following way: F * ( p i , n) = F ( p i ) + no'47ra 2 .

(9)

The surface tension coefficient is denoted by 0", while the equivalent bubble radius a is determined by ( 3a ~1/3 a = \ 4--~-nn/ "

(10)

It is shown in ref. [10] that the free energy density F ( p l ) and the thermodynamic potentials ~s (i = 1, 2) satisfy 2

(11)

d F = ~, ~'£idpi , i=1

1 d~l = ~ dpg,

1

d/z2 = ~gg dpg,

(12)

2 E

(13)

fli ~'l~i = F + pg .

i=l

It follows from (9) that dF* = / x ~ dPl + ~2 dP2 + /*n d n ,

(14)

where ~

20" = ~, - - - , pe a

~

0" a -

a n

(15)

Furthermore, P l l ~ + p2tZ2 + n~n = F* + p*

(16)

with P* = Pg

20" a

(17)

A.J.N. Vreenegoor, J.A. Geurst / Two-phase flow through tube

416

The kinetic energy density K(pi, ui) of the two-phase mixture is composed of the kinetic energy densities of the separate phases and the kinetic energy density associated with the local backflow of the liquid around the bubbles in the case of relative flow. We accordingly have 2

,

K(p,

~PiUi + ½pem(a) w 2 ,

(18)

i=1

where w = u 2 - u I denotes the relative velocity. The virtual-mass coefficient is represented by m(a). The validity of expression (18), in particular the additivity of the contributions at the right-hand side, is discussed at some length in ref. [11]. In ref. [22] the virtual-mass coefficient is allowed to be a function of the Weber number defined by 2

We-

PeW o'/2a "

(19)

We assume here that We ~ 1, so that in a first order of approximation the virtual-mass coefficient depends only on the volume density of the bubbles. The dynamic equations governing the non-dissipative flow of a bubbly liquid/gas mixture are obtained from Hamilton's principle of least action in the form t1

f dt f dV L(pi, n, ui) = O, to

(20)

.0

where the Lagrangian density L(Pi, n, u~) is given by

L ( p i, n, ni) = K(pi, ui) - F*(pi, n) .

(21)

The variations in (20) are subject to the side conditions (7), (8) and the two Lin constraints expressed by [23]

-

-

Ot

+ V.(~ui) = 0

'

i = 1,2

'

(22)

Introducing the quantities ~Oi (i = 1, 2) according to ~ = pi~O~we obtain from (22) O

~-~+Ui

.V)~Oi = 0

,

i=1,2

•

(23)

A.J.N. Vreenegoor, J.A. Geurst / Two-phase flow through tube

417

Variational principle (20) yields the following equations of motion for the liquid and the gas: 0 "/11"1

0t

1 i + ~m ( a ) w 2 - ~ u1l ] -2u

+ V[u 1 • 1rI + ~

_ ~ z

lx(Vx~I)=O

n

°~ffT2'1- V(U2" '/'$'2'[- ~J'2 2U2) -[- -- V~l'n -- U2 X (V X '/i"2)= O. Ot P2

(24)

(25)

The generalized momenta ni (i = 1, 2) taken per unit mass are given by Pem(a) Pl

q'/'l = U l

W,

7T2=

Pem(a) P2

U2 + -

W-

(26)

We remark that 1

n

d/x~ = -- d p * , Pe

1

d/x2 + -- d/x. = -- d p * . P2 Pg

(27)

The variational procedure yields the following Clebsch representations: ~rI =V~oI + ~blVx1,

~'z =Vq~2 + qJzVx2 + n Vq~..

(28)

/9 2

The quantities q~i (i = 1, 2), q~, and xi(i = 1, 2) represent Lagrange multipliers associated, respectively, with the constraints (7), (8) and (22) [23]. The Clebsch representations are valid locally [24]. The equations expressing the conservation of total momentum and energy read, respectively, --

OP + V.//=0, 0t

(29)

0H + V.Q = O, Ot

(30)

--

where 2 P =

~, p,~, i=1

2 = F, p,u,,

(31)

i=l

2

11 = Z PiUi ~[Ti + i=l

H=K+F*,

pl ,

(32) (33)

418

A.J.N. Vreenegoor, J.A. Geurst / Two-phase flow through tube

~J'J'l ~- ~Jb~ Jr- ½m'(a)

Q = p,u,(u,.

(

.

w 2 --

½u~)

,2) "

+ P2//2 u2 " ~r2 + tx2 + -~2 I~ - ~ u2

(34)

The average pressure p of the two-phase mixture is given by d ( m ( o 0 ]w2

"

(35)

Expression (35) evidently includes the Bernoulli effects that are associated with the local backflow of the liquid around the bubbles in the case of a relative motion of the two phases. A detailed stability analysis of the system of two-phase flow equations determined by (7), (8), (24) and (25) shows that a uniform two-phase bubbly flow is marginally stable if and only if the virtual-mass coefficient satisfies m ( a ) = ½c~(1 -

a)(1 -

300 .

(36)

For details see refs. [10,22]. Expression (36) is valid when W e < 1. The virtual-mass coefficient may be regarded as an order parameter that represents the local microscopic distribution of the gas bubbles on a macroscopic level. The bubble distribution underlying (36) is apparently anisotropic [10]. When the mass ratio Pg/Pe is disregarded and the bubble radius a is taken constant (ping-pong ball approximation) we obtain #1 from (7), (8), (24) and (25) in the limit of vanishingly small values of the void-fraction equation (1) for the motion of a single gas bubble in liquid. Indeed, in the limit of an infinitely dilute dispersion of bubbles the expression (O/at +//2 "•)U2 may be identified with the acceleration d u z / d t of a single bubble. Equation (1) may be put in the form (2), where the left-hand side represents m 1 times the relative acceleration of a bubble with respect to a frame that is rigidly attached to the locally rotating liquid. The left-hand side is accordingly objective [7,11].

3. Additional terms required near a solid boundary

It was discussed at some length in the introductory section of this paper that the free energy density F*(Pi, n) of a two-phase bubbly mixture has to be extended to the free energy density F(Pi, n, Vp~) determined by (see (3)) ~ We take the opportunity to make the followingcorrection in expression (6.3) of ref. [11]: the factor preceding ~r. v2 should read m(a)/a -m'(c~). The correction does not have consequences for the subsequent results of the paper.

A.J.N. Vreenegoor, J.A. Geurst / Two-phase flow through tube 1

F(pi, n, Vpl ) = F * ( p i, n) + ~ g11(Vpl)

2

419

(37)

in order to take into account the small irregular motions of the bubbles observed in two-phase bubbly flow near a solid boundary. An alternative way of introducing the dependence on Vpl would have been to include it as an additional variable in the virtual-mass coefficient. For reasons of simplicity, however, we have chosen for the formulation according to (37). If g~l is constant, it follows from (37) that (38)

dF =/21 dpl +/x z do2 q- /'£n dn + V. (g,,Vpl dp,), where

~P

/zl - ~p, - / x ~ - v . ( g l l v p , ) .

(39)

When expression (37) is used in the Lagrangian density (21), it is recognised by comparing (14) and (38) that the equations of motion resulting from Hamilton's principle of least action (20) may be obtained from (24) and (25) by substituting the thermodynamic potential /Jq for /z~'. In that case the pressure p* appearing in (35) has to be replaced by/~ determined by (cf. (16)) P = P l ~ I q- P2,0"2 +

n/x.

-

(40)

t~ .

By virtue of (16), (37), (39) and (40) we have (41)

/~ = P * - g l l [ P l Ap, + ½(Vp,)2],

where A denotes the Laplacian operator. It is derived in a straightforward way that the momentum f l u x / / ( s e e (32)), the energy density H (see (33)) and the energy flux Q (see (34)) are extended, respectively, to 2 I I = ~ PiUi'l'gi "b p l -'F g l l V P l V P l , i=1

(42)

H = K + F,

(43)

Q

=

P,U,(Ul " ~T1 31- ~ I

(

-~ IFn'(O/) W2 _ ~Ul)l 2

n

+ P2U2 U2" ~rg2+ bt2 + ~

12)

I't'n -- 2U2 -- g l l

~701

oo, Ot

(44)

In view of expression (35) and the remark preceding (40) the pressure p of the

420

A . J . N . Vreenegoor, J . A . Geurst / Two-phase flow through tube

bubbly liquid/gas mixture is given by d (\ 1m_ ( a o~/ )]w 2 • P = P + ½Pt( l -- oll2 -~a

(45)

The two terms that appear at the right-hand side of (42) and (44) with the coefficient g~l originate from the divergence term in (38). Similar expressions have been derived in the macroscopic analysis of "healing" phenomena in He II near a solid boundary [17].

4. Dissipative effects and external forces

The terms which represent dissipative effects have to be introduced in the equations of motion in such a way that the conservation of mass and the conservation of total linear momentum are not violated. In spite of these restrictions, however, there may be some arbitrariness in the selection of the dissipative forces. That is the reason why the presentation in this section looks somewhat different from the discussion given in the appendix of ref. [10]. The different approaches, however, prove to be equivalent. The two-phase flow equations for the liquid and the gas are extended, respectively, to C]"J"$'I qt- V [ U l o ,jT1 Jr ~"~1 + ½ m ' ( a ) W 2

Ot

1 2

- ~ u l ] - ul x ( V x ~J'gl)

=IF+_I V.~.w _ _ l V . 7 . + 1 Pl Pl P ~ F1

(46)

and (~~"2

at

-

1 2

n

+ V(u2. ~'2 + ~2 - ~u2) + ~ V~. - u2 × ( v × ~'2)

1 P2

F---

1 P2

1 1 V'a'w-- V'% + F2. P

P2

(47)

The external force densities are represented by F 1 and F 2. The viscous stress tensor of the bubbly liquid/gas mixture is denoted by % while F and V. ~'w are the densities of mutual friction forces which act on the two phases during relative flow. It is easily recognised that the balance of total momentum corresponding to (46) and (47) is expressed by 0P Ot

--

+ V"//=

F 1 + F 2 ,

(48)

A.J.N. Vreenegoor, J.A. Geurst / Two-phase flow through tube

421

where the m o m e n t u m density P is given by (31) and the momentum flux H is determined by 2

171 = ~

p , u , cr i

(49)

+ g,, 17p,17p, + p l + r u .

i=1

When dissipative effects are taken into account the energy conservation equation (30) has to be generalised to an equation expressing a balance of energy of the following form: OH

--

Ot

+ 17. Q = - R o ,

(50)

where R d denotes the rate of dissipation. The energy density H is given according to (43) by 2 H = ~

1

i=1

2

~fliUi -F ½Pem(a)

W2

O"

+ F(pi) + 3 - a + ½gll(VPl

)2

It is easily verified that the energy flux Q and the rate of dissipation take the form Q = plUl[Ul"

'~1"1 -[- ~L~1"[- ½m'(a)

(

+ P2U2 u2" ~2 +1~2 + --gllVO1~t

(51)

.

a

n 02

Rd

should

W2 -- ~U,]' 2

I~n -- ~u2

+%'w+%'u

(52)

and R d = F - w - ~'w:Vw - % : V u .

(53)

T h e average mass velocity of the two-phase mixture is denoted by u, i.e., 2

pu = ~

i=1

Piui,

(54)

where p = p~ + P2. Note that the fragmentation and coalescence of bubbles, which may be treated as an independent dissipative process, have been disregarded. A similar remark applies to the diffusive flux of the bubbles with respect to the average mass velocity u 2. A second-order tensor like ~r can be decomposed in the following way: ~" = dev 7 + 1 tr0. ) ! + ½(¢ _ ~.T),

(55)

A.J.N. Vreenegoor, J.A. Geurst / Two-phaseflow through tube

422

where the deviatoric part is given by dev r = l ( r + r T ) - ~ tr(7) I .

(56)

We apply the decomposition (55) to the second-order tensors that appear in the bilinear form (53) for the rate of dissipation R d and use the thermodynamics of irreversible processes (see ref. [25]). On account of Curie's theorem for an isotropic medium (note that the anisotropy associated with the direction of Vp~ is disregarded as in the case of the virtual-mass coefficient) we introduce the following constitutive relations:

-dev~r w

\A2~

A22

_ ½(.rw _ ~.T)= A½ [Vw -

dev~rw (vw)T]

'

(58)

,

F= Bw, (-~ __1

(59)

t r 0 " ) ] = ( C'1 C21 tr0"w) /

C'2](V'u) C22 / \ V . w

(60) •

The material coefficients, which are allowed to depend on p~ and P2, should satisfy the Onsager relations A12 = A21

,

C12 = C21 .

(61)

According to the second law of thermodynamics the matrices Aij and C a (i, j = 1, 2) are positive semi-definite, while the material coefficients A and B can take only non-negative values. Since the constitutive relations have to fulfil the requirement of objectivity (material frame indifference; see, e.g. ref. [26]) the viscous stress tensor ~-, is symmetric. The equation that expresses the vanishing of the antisymmetric part of ~', has been deleted from the constitutive relations (57) to (60). With a view to the analysis of the vertical flow of a bubbly liquid/gas mixture to be given in the next section and in order to indicate clearly the physical meaning of the various material coefficients we introduce the following notations and simplifications: All = 2flu , A = 2rtV,

Al2=A21 =0,

A~2 = 2r/w ,

(62) (63)

A.J.N. Vreenegoor, J.A. Geurst / Two-phase flow through tube

1 B = r/CM ~ , CII = ~ u ,

423

(64) C22 = ffw.

C12 = C21 = 0 ,

(65)

The constitutive relations (57) to (60) now take the form = n.[vu

+ ( v . ) w] +

- 7 w = Bw[VW +

1 F = ~LM--~W. a

-

(Vw) T] + (ffw - ~n~)V. w l + nF[VW -- (Vw)T],

(66) (67) (68)

Note that the cross-effects determined by the (small) off-diagonal elements of the matrices A i j and Cij (i, j = 1, 2) have been disregarded. It will be clear from (66) that r/, and ~', represent, respectively, the dynamic and second viscosity of the bubbly mixture. The last term at the right-hand side of (67) determines a couple of the Fax6n type. When the interactions of the bubbles are neglected we may write nF = ne ~ ,

(69)

where T/e denotes an effective viscosity of the bubble dispersion. In a local Stokes flow, where a no-slip boundary condition prevails at the surface of the bubbles, it is true that ne = ¼ " ,

(70)

where 7/ represents the viscosity of the liquid (see e.g. ref. [27]). Expression (68) for F corresponds to the Levich-Moore mutual friction force. In the case where the interactions of the bubbles are disregarded, we have [28] rtLM= 9rla .

(71)

In view of the simplifications which have been introduced it will be clear that the viscosity coefficients in the constitutive relations (66) to (68) are not completely independent of the physical conditions under which the two-phase flow of a bubbly liquid/gas mixture proceeds. Their values may in fact be adjusted somewhat in order to meet the requirements of a particular type of flow. A mild form of that procedure is adopted in section 4, where the theoretical results of section 5 concerning vertical two-phase flow are compared

A.J.N. Vreenegoor, J.A. Geurst / Two-phase flow through tube

424

with some known experimental data. Matching of the theoretical and experimental results proves to be successful in producing reasonable values for the material coefficients.

5. Bubbly flow through vertical tube We consider a two-phase bubbly flow through a vertical cylindrical tube of circular cross-section (radius R0). It is assumed that the flow is steady, vertically uniform and rotationally symmetric (no swirl). A system of cylindrical coordinates (r, 0, z) is used with symmetry axis along the centre line of the tube. The external forces experienced by the liquid and the gas of the bubbly mixture result from gravity, i.e., i = 1, 2.

F i = -pigiz,

(72)

It will be assumed in the analysis of this section that the masses of the gas bubbles are approximately equal. In addition, the mass density pg of the gas is taken constant. The spherical bubbles are accordingly characterised by a constant radius a (ping-pong ball approximation). That approximation seems to be justified locally, when the Mach number itself and the Mach number divided by the local Froude number are small, i.e. w

2

1,

gh

~ 1,

(73)

where h denotes the height of the tube section. The ping-pong ball approximation is evidently not valid g l o b a l l y because the bubble radius will vary (weakly) with vertical position. Since we are only interested in local properties of vertical two-phase flow, that variation is disregarded here. It is also in accordance with the assumption of uniformity of the flow in vertical direction. In view of the assumptions of steadiness, vertical uniformity and rotational symmetry of the flow, all quantities with the exception of the pressure depend only on the radial coordinate r. In the case of axial flow the pressure p* should in addition be a function of the coordinate z. The mass conservation eqs. (7) and eq. (8) expressing the conservation of bubble number are satisfied when U i = ( 0 , O, Ui) T where u i = u i ( r ). Accordingly P/'¥i = ( 0 , 0 , "N'I)T with 7ri = "l'ri(r). The equations of motion (46) and (47) now take the following form in cylindrical coordinates:

A.J.N. Vreenegoor, J.A. Geurst / Two-phaseflow through tube

d [

12

-~r u , ~ r , - ~ u , + l m ' ( a )

-

dlr 1 Ul dr

-

W2

(74)

11 ~[

Pla z ~TLMW Pl r dr

1 1r dr~ ( ~Tur ~r) dTr2 d(12) dr u27r2- 2 u2 u2 (Pl

da)]

+Pc ; -~r g,,r-~r

1 Op* Pe Or '

0~_,+ 1 Pe Oz +

1 d(

+ P2)

-

425

dr

~w]

(~w + r/F)r --d-;r

- g = 0, 1

(75)

Op*

(76)

pg Or

lop* 1 1 1 d [ dw] pg OZ p2a2 T/LMW+ P2-- -r --dr (~/~+ "t/F)r 1

(Pl

+

+

1 d( du) P2) r dr ~7.r -~r - g =0"

(77)

In deriving these equations the constitutive relations (66) to (68) have been used. It follows from (74) to (77) that (78)

p* = - C z + h(r),

where C is a constant. By combining (75), (77) and (78) we obtain the equations C+1

d (

r ~

and

(~g)

du)

(79)

7q"r-~r - - ( P l + P 2 ) g = O

(

~o )[1

11~((

) dw)] (80)

= 0.

Eq. (79) expresses the balance of total linear momentum in the axial direction. Elimination of the pressure gradient between (74) and (76) yields m(a) du, ( a(1 - a) w -~r + m ' ( a ) + Pe ~r

r dr

gllr~r

m_~))

=0.

dw

de~

w -~r + ½rn"(a) w 2 --dr

(81)

426

A.J.N. Vreenegoor, J.A. Geurst / Two-phase flow through tube

Since

Pg/Pt is very small, eqs. (79) and (80) may be reduced to

C+I

d (

du,]

r dr rl.r d r / -

(82)

pe(1-t~)g=O

and

1

eta 2 ~TLMw

-- --

(

a r dr ( % + r / v ) r

=0.

(83)

Eqs. (81)-(83) constitute a seventh order system of ordinary differential equations for the determination of a, u~ and w as a function of r. Note the role of the axial pressure gradient C as one of the external parameters characterising the physical conditions of the flow. In view of axial symmetry the solutions have to fulfil at the tube axis the requirements

(-~)r=0

(dul~ dw dr /Ir=o=(d--~T)r=0=0"

=\

(84)

Furthermore

(85)

(U~)r=0=(u,)0,

where (Ul) 0 is the given value of the liquid velocity at the tube axis.It seems natural to impose at the wall of the tube the boundary conditions

(a),=Ro = 0 ,

excluded v o l u m e ,

(Ul)r=Ro = 0 ,

no slip,

(d r) r=R0 =0,

vanishing Fax6n couple.

(86)

When the viscosity coefficient fILM is expressed according to (71), we derive from (83) on account of the boundary conditions (84) and (86): w = w0 ,

(87)

where w o denotes the value of the relative velocity at the tube axis. Eq. (83) accordingly reduces to 9rl C = ~ wo . a

(88)

It will be assumed that the value of the dynamic viscosity rl, of the bubbly mixture is constant and that the virtual-mass coefficient m(o 0 is determined by

A.J.N. Vreenegoor, J.A. Geurst / Two-phase flow through tube

427

its expression (36) corresponding to marginal stability of a uniform two-phase flow. When the velocities and the radial coordinate are made dimensionless by means of, respectively, the constant relative velocity w 0 and the tube radius R 0, eqs. (81) and (82) take the form a

r

- - o~ e d- 151 -

r dr

=0

dr~

du I (1-3a)~+(4-9a)dr

dR

(89)

'

d [1 d ( d R ) ] -3'11 drr r drr r-d-r-r = 0 ,

(90)

where 9rtw 0 Oee=l

Re-

pega

2 ,

151--

a 2 Fr R o Re ~

2 1420 Fr = -ga - ,

PeaWo

,

(91) a2 gripe Yll = 2 R~ wZa2"

Eq. (89) expresses the balance of forces acting on the bubbly mixture in vertical direction. The term a e appearing in (89) clearly vanishes when the bubbles of a dilute dispersion attain their terminal velocity, while the last term takes account of the viscous shear force. The three terms in (90) represent, respectively, the lift force experienced by bubbles in a shear flow, the gradient of the Bernoulli pressure associated with the local backflow of the liquid around a bubble (see the third term in the expression between square brackets appearing in (46)) and the dispersive force resulting from the small irregular motions of the bubbles. By eliminating u I between (89) and (90), we finally obtain 1 d (4-9a O/

--

R e

--

151

--

r ~

\1-3~

1d{ 1

+ 152 - -

da) r

r

d [1 d ( da)]} = 0 , r

(92)

where 152 = "Yl 1 151" When a is small, the equations may be linearised in order to arrive at an approximate solution. The linearised versions of (90) and (92) read

d ( da)]

du~

da d [1 dr + 4 ~ r - Y l l d r r r drr r-d-rr

a_ae_4E

1 d (da) r drr r ~ - r

=0,

1 d { d [1 d ( d a ) ] } + ez-r drr rdrr r drr r-d-Tr =0.

(93)

(94)

428

A.J.N. Vreenegoor, J.A. Geurst / Two-phase flow through tube

The boundary conditions to be associated with these equations should follow from (84), (85) and (86). We first infer from (93) that dot

{~r[a~r (r ~r)])r~0 0

(95)

By integrating (94) and (93) we then obtain

j ( a - a e ) r d r = a E l r d r ~o - E2rdrrd[l~(do)] r drr r ~

,

(96)

d (r ~ ) ] r = l }

(97)

0

j Ulrdr=-{2ae+Y-~[~

da dr (r ~r)]r=1}r2

0

+ ( Y l l - - 1 6 " l ) r - ~ r +4"2rdrr

r drr r~-r

"

(98)

The average (u~) of the liquid velocity is accordingly given by 1

(Ul)=2

-~r r=l J Ulrdr= - n a e + ( Yll - 32"l -8"2) (~°) -~11-8,~)(~)~1 +s,2~:~, 0

d2a

[ d3t~ \

(99)

while the void fraction averaged over a cross-section of the tube is determined by 1

(a)=2

f otr dr = ae + 2(4~1 + ~z)(d~) r=l

0

dr2]r=l + \ dr3 J r=l l • It will be clear that the solutions of the fourth order linear differential equation (94) have to satisfy the boundary conditions (95) and (99). Note that expression (99) may be used as an alternative for boundary condition (85). In view of (84) and (86) the set of boundary conditions is completed by requiring that (d°)

r=0 = (a)r=l =0"

(101)

A.J.N. Vreenegoor, J.A. Geurst / Two-phase flow through tube

429

The quantity a e and the average (u I ) of the liquid velocity play the role of external physical parameters. It should be remarked in this connection that for physical reasons the function a(r) and its derivatives have to take finite values at r = 0 . Eq. (94) may be solved in closed form in the terms of Bessel functions of order zero and complex argument (see the appendix). We prefer here an asymptotic analysis. Experimental results (see e.g. ref. [4]) indicate that a reasonable value for ")Ill is given by 711 ~

(102)

10-2"

The following orders of magnitude are commonly encountered in vertical bubbly flow (see the next section): -a- ~ Ro

10_ 2 ,

-Re -~ Fr

10.

(103)

It follows from (91) and (103) that e I ~ 10 -5 .

(104)

In view of (102) and (104) the outer solution satisfying (93) and (94) is, in a first order of approximation, given by O~ =

--life,

U 1 =

(U,)

(105)

0 .

We now apply the coordinate transformation y = -ln r. Eq.

(94)

(106)

takes in that case the form d 2 (e2y d2n~) \ - 4e I e dy 2 dy 2 /

e 2 e 2y - -

d2¢~ + dy 2

2y - -

a

-- a e =

0,

(107)

while eq. (93) is transformed into du, dot dy + 4 ~ y - - Y l l

d { d2a] ~yy \e 2y dy 2 / = 0 •

(108)

By introducing the strained coordinate )7 by means of 1 / 4 --

y = ~52 y

(109)

430

A.J.N. Vreenegoor, J.A. Geurst / Two-phase flow through tube

we derive from (107) and (108) that the inner solution is, in a first order of approximation, determined by d4ot d y 4 + ot - ote = O,

dff 1 dy

(110)

d3ot ---0,

(111)

dfi 3

Ul=('~ll/E1)l/21gl . According to (102) and (104) we have (61/ Y11)1/2 3 × 10 -2. The boundary conditions to be associated with eq. (110) are expressed by (see (99))

where

(d2ot

(ot)/=o = 0 ,

\ dfi2/y=o = - ( u l )

, (112)

dot

ot - %-->0,

Yd--~--->O' when 37--->~.

The conditions at infinity provide for the matching of the inner and outer solutions. Note that the average (K1) of the strained liquid velocity has been introduced according to ( u 1) = (Yll/el ) ~/2 ( ul ). The solution of (111) requires the additional boundary condition (ai);=(, = O.

(113)

The inner solution determined by the differential equations (110), (111) with the boundary conditions (112), (113) is now given by a-ae ote

-

( cos y , - -(u,)0 sin y* )

/~1- (/~I)0-

(a,),,

where y* = (114) that

2-1/2fi =

(/~1 ) ~-- (/~1)0"

ote

(cosy*+

e -y*

, (114)

OLe

sin y*) e y*

--2-l/2e2 ~/4 In r. It is easily verified by means of (112) and

(115)

The equality expressed by (115) is valid in a first order of approximation. It results from the fact that in that order of approximation the velocity profile of the liquid may be considered as flat (turbulent-like profile; see the next section). In the same order of approximation we have according to (100):

A.J.N. Vreenegoor, J.A. Geurst / Two-phaseflow through tube (a) = a e .

431 (116)

It is easily verified by means of (99), (100) and (115) that in a first order of approximation (dZa] (/~1)0 = -- \ d)72];=0 •

(117)

Eq. (117) implies that for upward (downward) flow of the liquid the voidfraction profile at the wall of the tube is concave (convex) when considered as a function of y = - I n r. That property of vertical two-phase flow is also found experimentally. We refer the reader for more details to the next section, where the solutions of eqs. (89) and (90) satisfying some appropriate dimensionless forms of the boundary conditions (84), (85) and (86) are studied numerically together with their linear and asymptotic approximations. The numerical results are matched there with some known experimental data.

6. Numerical results and discussion

In section 5 eqs. (89) and (90) were linearised and solved by deriving the first-order asymptotic approximation (114). We now wish to deal with the solutions of (89) and (90) in more detail. After integrating the equations numerically we are able to compare the asymptotic solution, the numerical solution of the linearised equations and the numerical solution of the nonlinear equations. The correlation of theoretical and experimental results is analysed next. In order to substantiate the theory a rather complete picture of vertical two-phase flow is given by making a comparison with experiments on upward flow [2, 4], downward flow [4] and the collective rise of bubbles in a stagnant liquid [21]. A summary of the contents completes the paper. The numerical integration of (89) and (90) is not straightforward due to the singular behaviour at r = 0. In addition, the conditions (84), (85) and (86) constitute a boundary-value problem which requires an iterative numerical procedure. The first difficulty may be removed by means of the transformation y* =2-1/2fi = -2-1/2~21/41nr, which was also used in section 5. As a consequence, the wall of the tube is positioned at y* = 0 and the axis of the tube at y* = ~. Since the case y* = ~ can not be treated numerically, some boundary values have to be given at a finite value y~ of y*. The flat behaviour of the asymptotic solution suggests to take (Of)y.=y~ = ae(1 + 8 ) ,

(Ul)y*=Y~ = (Ul)° '

(dUl'~ \dy*/y*=y~ = 0 '

(118)

432

A.J.N. Vreenegoor, J.A. Geurst / Two-phase flow through tube

where 6 is a small parameter. At the wall of the tube we require that

(oe)y*=O = (Ul)y*=O = O,

(119)

in acccordance with (86). The transformed eqs. (89) and (90) are subsequently written as a set of five first-order differential equations. Robust methods, based on shooting techniques, may be applied to obtain numerical solutions of the set of equations that satisfies the boundary conditions (118) and (119). We have chosen for the D02HBF routine from the N A G library, which is described in detail by Gladwell [29]. The routine requires, at each of the two boundaries, initial estimates of the quantities that are left unspecified there. In the case of the transformed version of the eqs. (89) and (90) the initial estimates at y * = y~ follow from the outer solution given by (105). They read

(doe) =(d2oe ~ d ~ v*=y~ \ dy*2 / y*=y~ = 0 "

(120)

At the wall of the tube y * = 0 the asymptotic solution (114) is the used to obtain the estimates

(dye,)

((/~1)o)

y*=o

=a e 1+

oee

( d2oe '

\ d Y *2]y*=°

= --2(/,71) 0 , (121)

dy,/y.=o=(Ul)o.

.

1

(~)o

"

A numerical exploration showed that it is favourable to impose the boundary conditions (118) at y* = y~ = 12 with 6 = 10 -3. It also demonstrated that the direction of integration (shooting direction) should be taken from y* = y~ to the wall of the tube at y* = 0. The physical parameters which may take various values are given by e I and 3'11, defined in (91). It follows from (91) that el is determined by el

--

2 WO~ 1.25 X lO-4Wo,

PegRo

(122)

when the value of the tube radius reported by Serizawa et al. [2] (2R 0 = 57.15 mm) is substituted and 7, is taken, in a preliminary way, equal to the dynamic viscosity "0 of water at 20°C (77 = 10 -3 kg/ms). A moderate value (in m/s) for the difference velocity corresponds to the estimate (104) for el. That value also applies to the experiments by Wang et al. [4], since their tube

A.J.N. Vreenegoor, J.A. Geurst / Two-phase flow through tube

433

d i a m e t e r is equal to 60 m m . The equivalent spherical diameters 2a reported by Serizawa et al. [2] vary from 3.5 m m to 4 mm. They lead to the estimates given by (103). Wang et al. [4] do not discuss the bubble size. In view of the preceding analysis we take txe = 0.05, (u~)0 = 5, ~1 = 4 × 10 -5 and YH = 1 0 - 2 , which represent typical values that may be encountered in practice for the case of upflow. It follows from (122) that a reasonable value of the difference velocity w 0, viz. 0.3 m / s , is associated with the value of El . Since the liquid velocity is made dimensionless by means of w 0, (u~) 0 = 5 implies that the liquid flows with a velocity of approximately 1.5 m / s at the centre of the tube. T h a t is certainly an acceptable value. Fig. 1 gives the void-fraction distribution o~ and the velocity profile u~ for three different cases: (1) the numerical solution of (89) and (90), (2) the numerical solution of the eqs. (89) and (90) when linearised and (3) the asymptotic solution (114). The numerical solutions are obtained by means of the transformed equations with the b o u n d a r y conditions (118) and (119) as discussed before. Clearly a distinct m a x i m u m in the void fraction occurs near the tube wall. The value of the m a x i m u m may exceed twice the value of the void fraction at the centre of the tube. The velocity profiles are flat and look similar to turbulent profiles, although turbulence has not been taken into account. The asymptotic solution and the numerical solution of the linearised equations overestimate the behaviour of the numerical solution of the nonlinear equations. That effect increases with increasing txe. Use of the linearised instead of the original non-linear equations leads to a m a x i m u m in tx slightly closer to the wall. It m a y therefore be concluded that when the gas fraction is not too large the calculated void-fraction distribution is not very sensitive to the precise form of the nonlinear terms in re(a). Varying ~/11 and ~ shows that these p a r a m e t e r s have a distinct effect on the value and the position of the m a x i m u m in the gas fraction a : increasing Yl~ decreases the m a x i m u m and moves it away from the

.2 Or.

Ul10t3

a

0

0

b

1

p

0

Fig. 1. (a) Void fraction a and (b) dimensionless liquid velocity u, as functions of r in vertical upward flow for a e = 0.05, (ul) 0 = 5, ~1 = 4 x 10-5, y. = 10-2; (1): numerical solution of eqs. (89) and (90), (2): numerical solution of eqs. (89) and (90) when linearised, (3): asymptotic solution (114).

434

A.J.N. Vreenegoor, J.A. Geurst / Two-phase flow through tube

wall of the tube, increasing E1 increases the maximum and also moves it away from the wall. Both effects will be used to obtain agreement with experiments reported in the literature. Note that in the present theory the void fraction can not exceed 1/3, since at that value the virtual-mass coefficient re(a) given by (36) changes sign. All theoretical curves to be presented below are obtained by numerically integrating the transformed version of the nonlinear eqs. (89) and (90) subject to the modified boundary conditions (118) and (119). Serizawa et al. [2] report on an experimental investigation of vertical upward flow. They measured the bubble velocity in addition to the liquid velocity. In the case of low void fractions their measurements indicate that a constant difference velocity w 0 is a good approximation. From their fig. 14(A) we obtain ae =0.035, w 0 = 0 . 1 4 m / s and u 1 = 0 . 9 4 m / s at the axis of the tube, which yields (ul) 0 -~7. We now vary e 1 and Yll in such a way that the position and magnitude of the maximum of a are close to their measured values. That is achieved by taking e 1 = 6.7 x 10 -5 and Yll = 2 x 10 -2. Fig. 2 shows that for those values of the physical parameters the theoretical and experimental curves for the void fraction nearly coincide. The calculated velocity profile for the liquid, however, is flatter than the measured one. The matching is clearly not as good as in the case of the void fraction. The value of the dynamic viscosity % of the bubbly mixture that corresponds to the value of e 1 is given by r/, = 3.8 × 10 -3 kg/ms. Obviously, the presence of the bubbles increases the viscosity ( r / = 10 -3 kg/ms for pure water at 20°C). This may be expected since the small irregular motions of the bubbles result in an additional contribution to the dynamic viscosity comparable with the eddy viscosity known from the theory of turbulence. Wang et al. [4] investigated both upflow and downflow. They report on measurements of the void fraction, the liquid velocity and turbulent fluctuations in the liquid. The gas velocity, however, is not measured. It is therefore not possible to obtain an estimate of the difference velocity w 0 directly from

Ol

0

0 I

r

0

r'

0

Fig. 2. (a) Void fraction a and (b) dimensionless liquid velocity u~ as functions of r in vertical u p w a r d flow for a e = 0.035, ( u l ) o = 7, E1 = 6.7 × 10 -5, %~ = 2 × 10-2; : numerical solution of eqs. (89) and (90), 0 : experiments (ref. [2], fig. 14(A)).

A.J.N. Vreenegoor, J.A. Geurst / Two-phase flow through tube

435

their figures. They do present, however, the values of the superficial liquid velocity j~ and the superficial gas velocity ]2 which are determined by j, = ((1 - a ) u , > ,

jz(aU2),

(123)

where 1

(') = 2 f - r d r .

(124)

0

Adding j~ and J2 and assuming a constant difference velocity we derive

(125)

J, +J2 = ( u , } + Wo(~) •

Relation (125) is used to obtain an approximate value for the difference velocity w 0. We selected figs. 6 (void fraction) and 11 (liquid velocity) from ref. [4] for a comparison with our theoretical results on upward flow. Wang et al. [4] report Jl = 0.71 m / s and J2 = 0.10m/s. By integrating the measured values for a and u~ we estimated the difference velocity by means of (125) to be approximately given by w 0 = 0.10 m/s. We also determined j~ from the measured values and arrived at j~ = 0.72 m/s, in accordance with the reported value. Fig. 6 in ref. [4] shows that a e =0.1 while the value 0.84m/s of the liquid velocity u~ at the tube axis, which is inferred from fig. 11 of ref. [4], yields (ul) 0 = 8.4 by using the estimated value of the difference velocity. Fig. 3 displays the good agreement that is obtained with regard to the void fraction and the liquid velocity when e~ = 1 . 2 x 10 -4 and VII 1.6× 10 -2. By choosing a slightly different value for the difference velocity according to w 0 = 0 . 1 0 2 m / s but keeping (Ul) 0 = 8.4 the computed values of the liquid velocity are increased =

9.8 U1

.3 OC 0

1

p

0

p

0

Fig. 3. (a) Void fraction a and (b) dimensionless liquid velocity u~ as functions of r in vertical u p w a r d flow for a e = 0.1, ( u l ) o = 8.4, ~ = 1.2 x 10 -4, %~ = 1.6 x 10-2; : numerical solution of eqs. (89) and (90), O: experiments (ref. [4], figs. 6 and 11).

436

A.J.N. Vreenegoor, J.A. Geurst / Two-phase flow through tube

only slightly, while the superficial velocities take the values j~ = 0.683 m/s and J2 = 0 . 1 0 4 m / s . These values deviate not more than four percent from the values reported in ref. [4]. The dynamic viscosity flu of the bubbly mixture that corresponds to the value of E1 is equal to 10 -2 kg/ms. The value of 7/, therefore increases as a result of the presence of the bubbles (see the remarks made above). A comparison of figs. 2 and 3 shows that the velocity profile measured by Wang et al. [4] is much flatter than the profile obtained by Serizawa et al. [2]. It is therefore not surprising that Wang et al. [4] observe turbulent fluctuations in the liquid. The magnitude of the fluctuations does not exceed 10 percent of the value of the liquid velocity at the tube axis. Since the fluctuations measured by Wang et al. [4] increase with the value of the void fraction, they may be related to the presence of the bubbles. The fluctuations are also larger near the wall of the tube. Wang et al. [4] therefore conclude that the turbulence level in the liquid consists of wall induced and bubble induced turbulence. It would be interesting to know more about the origin of the turbulent fluctuations. We recall that the energy associated with the irregular bubble motions which are particularly prominant in the neighbourhood of the wall is incorporated in the present theory by means of an additional term in the free energy density that is given by (3). That term might therefore be considered to represent some form of "wall induced turbulence". The "bubble induced turbulence" distinguished by Wang et al. [4] might be associated with the kinetic energy density of the local "backflow" of the liquid around the bubbles. Although the liquid velocity of that "backflow" does not show fluctuations in time in a frame of reference moving with the drift velocity of the bubbles, it does exhibit such fluctuations in a laboratory system of coordinates as a result of the irregular distribution of the bubbles in space. Apparently, turbulence does not necessarily have to be included explicitly in the form of a Reynolds stress in order to arrive at the trends which are observed in the experiments. The void-fraction profiles presented in figs. 1, 2 and 3 satisfy in the neighbourhood of the wall the inequality Aa < 0, where A denotes the Laplacian operator. The profiles when considered as a function of y = - l n r are therefore concave near the wall in the case of upflow. That behaviour was already concluded from the asymptotic analysis in section 5 (see expression (117)). Since as a consequence A& > 0 in the vicinity of the boundary it follows from expression (41) for/~ that wall effects decrease the pressure p determined by (45). Wang et al. [4] also investigated downflow. Their figs. 9 (void fraction) and 13 (liquid velocity) were selected for a comparison with the present theory. We discuss the case where the superficial velocities have both a negative sign and are given by Jl = - 0 . 7 1 m/s and J2 = - 0 . 1 0 m / s . Wang et al. [4] use positive

A.J.N. Vreenegoor, J.A. Geurst / Two-phase flow through tube

437

values for Jl and J2 in view of the fact that their test section was turned upside down to measure downflow. By using their m e a s u r e m e n t s and expression (125) in the way like was done earlier in the case of upflow, the difference velocity w 0 is estimated to be 0.38 m / s , a somewhat large value to our opinion. In addition, it follows from the integration of the measured values that Jl = - 0 . 6 4 m / s , a value which deviates considerably from the value reported in ref. [4]. We do not have an explanation for this discrepancy. We will assume here that Jl has the value which follows from the measured gas-fraction distribution and velocity profile, i.e., Jl -- - 0 . 6 4 m/s. From the graphs [4] it is inferred that a e = 0.17 and ul = - 0 . 7 8 m / s at the tube axis. By varying el and ")/11 and performing computations with several difference velocities we arrived at a satisfactory agreement with the experimental results when (ul) 0 = - 3 , e 1 = 1.25 x 10 -4 and y~ = 1.9 x 10 -2, as demonstrated in fig. 4. The corresponding difference velocity and dynamic viscosity are given by w 0 = 0.26 m / s and 71, = 4.2 x 10 -3 kg/ms. The superficial velocities computed from the numerical results are given by j~ = - 0 . 6 5 m / s and J2 = - 0 . 0 8 m / s , which again d e m o n s t r a t e s that jl cannot be equal to the reported value of - 0 . 7 1 m/s. The m e a s u r e d void-fraction distribution shows a weakly oscillating behaviour before going to zero at the wall. The computed distribution behaves similarly although it is hardly visible in the plot of fig. 4. Both the computed and the m e a s u r e d liquid velocity profiles display a minimum near the tube wall. The fact that the minimum of the computed curve is lower explains why the calculated value j ~ = - 0 . 6 5 m / s is somewhat smaller than the value jl = - 0 . 6 4 m / s determined from the measurements. Vertical two-phase flow in a stagnant liquid was investigated by Kapteyn [21]. That type of flow is characterised by j~ = 0 m / s and forms the transition between upflow and downflow. Kapteyn [21] measured a void-fraction distribution for a superficial gas velocity J2 = 0.015 m/s. The pipe radius R 0 was equal to 40 ram. From his measurements we obtained a e = 0.067. Starting with

0

U1 0(,

I

p

0

-3.9

p

0

Fig. 4. (a) Void fraction a and (b) dimensionless liquid u 1 as functions of r in vertical downward flow for a~ = 0.17, (ul) 0 = -3, et = 1.25 x 10-4, y, = 1.9 × 10 z; : numerical solution of eqs. (89) and (90), 0: experiments (ref. [4], figs. 9 and 13).

438

A.J.N. Vreenegoor, J.A. Geurst / Two-phase flow through tube

( u l ) 0 = 0.5 we varied E1 and YH until a satisfactory agreement was achieved. Subsequently we lowered (ul) 0 to the value 0.225 which hardly affected the void-fraction distribution but made the dimensionless quantity j l / w o approximately equal to zero (t~(10 4)). The corresponding values of EI and 711 are given by el = 2.5 x 10 5 and 711 = 7 x 10 -2 . Fig. 5 shows the void-fraction distribution and the velocity profile for those values. The dimensionless superficial gas velocity computed from the numerical results is given by j2/wo = 6.6 × 10 2. Taking into account that J2 = 0.015 m / s we therefore arrive at a difference velocity w 0 = 0 . 2 3 m / s , which is an acceptable value. The dynamic viscosity of the bubbly mixture is determined by ~/, = 1.7 x 10 - 3 kg/ ms. F r o m fig. 5 we see that the liquid velocity shows an interesting behaviour. At the tube axis the velocity is positive but near the wall there is a local backflow where u~ becomes negative. There is even a region near the wall where u 1 < - 1 , which means that the bubbles attain a negative velocity in this region and move downwards. Kapteyn (private communication, 1989) has observed that effect experimentally for average gas fractions ( a ) varying from 0.2 to 0.3. The comparison with the measurements given above yields a mixture viscosity ranging from r/,, = 1.7 × 10 - 3 k g / m s (a e = 0.067) to "0, = 10 2 k g / m s ( a e = 0.1). The spreading seems somewhat large and requires some comment. T h e dynamic viscosity 7/, was determined from the measurements and may therefore contain, according to our analysis, possible contributions arising from turbulence or the Fax6n force. In addition, the viscosity was calculated by means of (122). It is accordingly quite sensitive to the precise value of the difference velocity w,. That quantity, however, had to be estimated since it was not reported with the measurements (with the exception of ref. [2]). Since we took, furthermore, the L e v i c h - M o o r e mutual friction force to be linear in a according to (71), eq. (83) had the simple solution w = w 0. As a result (81)

.08 OC

Illlt~..i:iii i i i

°/ a

b 0

f,

0

Fig. 5. (a) Void fraction a and (b) dimensionless liquid velocity u 1 as functions of r in a stagnant liquid for ao = 0.067, (u l ) . = 0.225, E~ = 2.5 × 10 5, %~ = 7 × 10 2 : numerical solution of eqs. (89) and (90), I: experiments (ref. [21], fig. IV-12).

439

A.J.N. Vreenegoor, J.A. Geurst / Two-phase flow through tube

and (82) could be reduced to (89) and (90). At moderate values of the void fraction, however, the nonlinear behaviour of the mutual friction force as a function of a should be taken into account. A constant difference velocity is then no more a solution of (83) and the original set of eqs. (81) to (83) has to be used for a more detailed analysis. It will be clear that the pressure gradient C, the bubble radius a and the difference velocity w would have to be measured in that case. It may be expected that the last term in eq. (83) as well as the contribution in (81) that is associated with the Bernoulli pressure (the second term on the left-hand side which contains d w / d r ) might have a reducing effect on the range of values of ~/,. We finally summarise the contents of the paper. The macroscopic theory of two-phase bubbly flow developed in ref. [10] was extended with a dispersive force according to (3). That term models the energy of the small irregular motions of the bubbles which, in the case of vertical two-phase flow, come into prominence in the neighbourhood of a solid boundary. It is was demonstrated by means of analytical and numerical calculations that the dispersive force and the "lift" force already contained in the original theory (see (1)) combine to produce void-fraction distributions and velocity profiles for vertical flow in a cylindrical tube which are in good agreement with experimental results obtained by various investigators. In particular the "turbulent" character of the velocity profile for the liquid phase is reflected by the computations. The calculated void-fraction distributions display, in conformity with the experiments, a distinct maximum in the vicinity of the tube wall both for the cocurrent upward flow of a bubbly liquid/gas mixture and the collective rise of bubbles in a stagnant liquid. Other flow configurations may be dealt with in a similar way. In particular the effect of dispersive forces on the propagation of void-fraction waves in two-phase bubbly flow is now been investigated. In a slightly extended form the results of the present paper may be found in chapter V of ref. [30]. For other recent developments in bubbly two-phase flow one might consult refs. [31-35].

Appendix. Solution of eq. (94) satisfying appropriate boundary conditions By introducing the differential operator D and the function w ( r ) according to r -~r

r

w=o~-ot

e

the differential equation (94) takes the form

(A.1)

440

A.J.N. Vreenegoor, J.A. Geurst / Two-phase flow through tube

(1 - 4elD

+ e2D2)w

=

0

(A.2)

or 2

I~ (1 + / 3 , D ) w = 0,

(A.3)

k=l

where /3, = -2E, + ( - 1 ) k i(E2 - 4 ~ ) 1/2 ,

k = 1, 2.

(A.4)

Solutions w, of (1 + / 3 , D ) w , = 0, (k = 1, 2), clearly also satisfy (A.3). Since in general w k = AkJo(fl,

-1/2 r), + B k Y o ( f l ; l / 2 r )

,

k = 1,2,

(A.5)

[AkJo( p~-1/2 , r)x + B , Y o ( f l , -1,2 r l l ,

(A.6)

we have 2

2

W= E

Wk = E

k=l

k=l

provided fll # f12" Here J0 and 110 represent the well-known Bessel functions of order zero. Requiring that w is not singular at r = 0, we immediately infer from (A.6) that 2 O/ -

O/e = E A k J o ( ~ k - 1/2 r ) . k=l

The coefficients A k are determined by the boundary conditions at a = 0 a t r = l we have

(A.71 r =

1.

Since

2

--ae = E

Akgo(~;l/2) •

(A.8)

k=l

The other boundary condition is related to the value of the exterior parameter (u I). According to (98) f Ulrdr

da d = - [2a~ + ½Y,, (Da)r=l] r~ + (Y,, - 16e,)r d r + 4e2r drr ( D a ) .

o It follows from (A.9) and (A.7) that

(A.9)

A.J.N. Vreenegoor, J.A. Geurst / Two-phaseflow through tube

441

1

l (Ul)

= f ulrdr o 2 =-2ot~ + ~ Ak[ l"Yll[~kX Jo(flk 1/2) k=l h_ (4ERflk-1 -- "Y1, + 16E1)fl-kl/2Jl(fl;1/2)]"

(A.10)

T h e c o e f f i c i e n t s A i ( k = 1, 2) f o l l o w f r o m ( A . 8 ) a n d ( A . 1 0 ) .

References [1] [2] [3] [4] [5] [6] [7] [8] [9]

[10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]

[23] [24] [25]

K. Kobayashi, Y. Iida and N. Kanegae, Bull. JSME 13 (1970) 1005. A. Serizawa, I. Kataoka and I. Michiyoshi, Int. J. Multiphase Flow 2 (1975) 235. S.W. Beyerlein, R.K. Cossman and H.J. Richter, Int. J. Multiphase Flow 11 (1985) 629. S.K. Wang, S.J. Lee, O.C. Jones Jr. and R.T. Lahey Jr., Int. J. Multiphase Flow 13 (1987) 327. J.J.Th. de Jong, Master's Thesis, Laboratory for Aero- and Hydrodynamics, Delft University of Technology, Delft, 1987. I. Zun, H.J. Richter and G.B. Wallis, The Transverse Migration of Bubbles in Vertical Two-phase Flow (Dartmouth College, Thayer School of Engineering, Hanover, NH, 1975). D.A. Drew and R.T. Lahey, J. Fluid Mech. 117 (1982) 91. J.M. Delhaye, Int. Chem. Eng. 3, 23 (1983) 385. G.K. Batchelor, in: Theoretical and Applied Mechanics, Proceedings of the XVIIth International Congress of Theoretical and Applied Mechanics, P. Germain, M. Piau and D. Caillerie, eds. (North-Holland, Amsterdam, 1989) p. 27. J.A. Geurst, Physica A 129 (1985) 233. J.A. Geurst, Physica A 135 (1986) 455. D.A. Drew and R.T. Lahey Jr., Int. J. Multiphase Flow 13 (1987) 113. T.R. Auton, J. Fluid Mech. 183 (1987) 199. T.R. Auton, J.C.R. Hunt and M. Prud'homme, J. Fluid Mech. 197 (1988) 241. V.L. Ginzburg and L.P. Pitaevskii, Zh. Eksp. Teor. Fiz. 34 (1958) 1240 (Sov. Phys. JETP 7 (1958) 858). R.N. Hills and P.H. Roberts, Int. J. Eng. Sci. 15 (1977) 305. J.A. Geurst, Phys. Rev. B 22 (1980) 3207. J.D. van der Waals, Verhandel. Koninkl. Akad. Wetensch. Amsterdam 1, 8 (1893) 3 (J.S. Rowlinson, J. Stat. Phys. 20 (1979) 197). D.J. Korteweg, Arch. Need. Sci. Ex. Nat. 2 (6) (1901) 1. I. Zun, Int. J. Multiphase Flow 6 (1980) 583. C. Kapteyn, Ph.D. Thesis, Twente University, Enschede, 1989. J.A. Geurst and A.J.N. Vreenegoor, in: Proc. First Int. Conf. on Industrial and Applied Mathematics (ICIAM 87), Contributions from the Netherlands, A.H.P. van der Burgh and R.M.M. Mattheij, eds. (Dutch Mathematics Centre, Amsterdam, 1987) p. 335. J.A. Geurst, Physica A 152 (1988) 1. R. Salmon, in: Annual Review of Fluid Mechanics, vol. 20, J.L. Lumley, M. Van Dyke and H.L. Reed, eds. (Annual Reviews, Palo Alto, 1988) p. 225. S.R. de Groot and P. Mazur, Non-Equilibrium Thermodynamics (North-Holland, Amsterdam, 1962).

442

A.J.N. Vreenegoor, J.A. Geurst / Two-phase flow through tube

[26] C. Truesdell and W. Noll, in: Encyclopedia of Physics, vol. III/3, S. Fliigge, ed. (Springer, Berlin, 1965). [27] J. Happel and H. Brenner, Low Reynolds Number Hydrodynamics, second revised edition (Noordhoff, Leiden, 1973). [28] V.G. Levich, Physicochemical Hydrodynamics (Prentice-Hall, Englewood Cliffs, N J, 1962). [29] I. Gladwell, in: Codes for Boundary Value Problems in Ordinary Differential Equations, B. Childs, M. Scott, J.W. Daniel, E. Denman and P. Nelson, eds., Lecture Notes in Computer Science, vol. 76 (Springer, Berlin, 1979). [30] A.J.N. Vreenegoor, Ph.D. Thesis, Delft University of Technology, Delft (1990). [31] G.B. Wallis, Multiphase Sci. Technol. 5 (1989) 239. [32] G.B. Wallis, in: Two-phase Flows and Waves, D.D. Joseph and D.G. Schaefer, eds. (Springer, New York, 1990) p. 150. [33] L. van Wijngaarden, Int. J. Multiphase Flow 17 (1991) 809. [34] J.A. Geurst, Int. J. Multiphase Flow 17 (1991) 815. [35] C. Pauchon and P. Smereka, Int. J. Multiphase Flow 18 (1992) 65.