1. No category

Chemical Engineering Science 58 (2003) 3765 – 3776
www.elsevier.com/locate/ces
Analytical and numerical evaluation of two-%uid model solutions for
laminar fully developed bubbly two-phase %ows
Osvaldo E. Azpitarte, Gustavo C. Buscaglia∗
Centro Atomico Bariloche and Instituto Balseiro, 8400 Bariloche, Argentina
Received 15 July 2002; received in revised form 22 February 2003; accepted 4 March 2003
Abstract
An analysis of the two-%uid model in the case of vertical fully developed laminar bubbly %ows is conducted. Firstly the phase distribution
in the central region of the pipe (where wall e4ects vanish) is considered. From the model equations an intrinsic length scale L is deduced
such that the scaled system reduces to a single equation without parameters. With the aid of this equation some generic properties of the
solutions of the model for pipes with diameter greater than about 20L (the usual case, since L is of the order of the bubble radius) are
found. We prove that in all physically meaningful solutions an (almost) exact compensation of the applied pressure gradient with
the hydrostatic force e4 g occurs (with e4 the e4ective density and g the gravity). This compensation implies %at void fraction and
velocity pro9les in the central region not a4ected by the wall, even when no turbulence e4ects are accounted for.
We then turn to consider the complete problem with a numerical approach, with the e4ect of the wall dealt via wall forces. The previous
mathematical results are con9rmed and the near-wall phase distributions and velocity pro9les are found. With the numerical code it is
also possible to investigate the regime in which the pressure gradient is greater than the weight of the pure liquid, in which case a region
of strictly zero void fraction develops surrounding the axis of the pipe (in upward %ow of bubbles).
Finally, the same code is used to study the e4ect of reducing the gravity. As g decreases, so does the relative velocity between the
phases, making the lift force increasingly dominant. This produces, in upward bubbly %ows, narrower and sharper void fraction peaks
that also appear closer to the wall.
? 2003 Elsevier Ltd. All rights reserved.
Keywords: Fluid mechanics; Multiphase %ow; Laminar %ow; Bubbly %ow; Modelling; Two-%uid model
1. Introduction
The two-%uid model (Ishii, 1975) is, at present, widely
used in the simulation of two-phase %ows. It emerges from
the exact conservation equations after ensemble averaging,
introducing unknowns such as the phase fraction k (x; t),
which is the probability of 9nding phase k at the position x
at time t. The ensemble averaging process needs to be complemented with several closure relations which are far from
established, either conceptually or quantitatively. For example, the lift term modeling transversal forces that act upon
bubbles or particles in the presence of a velocity gradient,
is not completely understood. The e4ective lift coe>cient
for bubbles that 9ts pipe-%ow data (CL 0:05; Lahey Jr.,
1990) is much smaller than that deduced from potential %ow
∗
Corresponding author.
E-mail address: [email protected] (G. C. Buscaglia).
0009-2509/$ - see front matter ? 2003 Elsevier Ltd. All rights reserved.
doi:10.1016/S0009-2509(03)00239-2
calculations (CL = 12 ; Drew & Lahey Jr., 1987), while CL
for particles depends on the rotational motion of the particle
around its center. Some authors, because of this uncertainty,
simply omit the lift term (Uchiyama, 1999). The e4ect that
solid walls exert on its surrounding bubbles is included
by means of a so-called wall force. Antal, Lahey Jr., and
Flaherty (1991) modeled this force considering spherical
rigid bubbles, but a recent model (Larreteguy, Drew, &
Lahey Jr., 2002) assumes deforming and elastic bubbles.
In this state of a4airs it is interesting to go back to very
simple %ows and try to extract as much information as
possible about the behavior of the model, starting with a
model that satis9es physical constitutive principles (Drew &
Lahey Jr., 2001). Perhaps the simplest %ow is laminar fully
developed %ow in a circular pipe. A decade ago Antal et al.
(1991) reported on fully developed solutions of the two-%uid
model equations. They solved the equations numerically for
some speci9c data so as to reproduce experimental data from
Nakoryakov et al. (1986). Our objective here is to go beyond
3766
O. E. Azpitarte, G. C. Buscaglia / Chemical Engineering Science 58 (2003) 3765 – 3776
that analysis and explore generic properties of the two-%uid
model equations (and solutions) under fully developed
conditions.
The two-%uid model is brie%y recalled in Section 2.
Section 3 is devoted to algebraic manipulations of the equations. We arrive at a single ordinary di4erential equation
(ODE), which is rendered non-dimensional by introducing
an intrinsic length scale L, of the order of the bubble’s radius. The solutions of this ODE are easily investigated using
Mathematica (Wolfram, 1991) since far from the wall no
non-dimensional number is involved. It turns out that for
pipes with diameter greater than about 20L no physically
meaningful solution exists other than %at void fraction and
velocity pro9les. This comes from an exact compensation
of the applied pressure gradient with the hydrostatic force
e4 g, where e4 is the e4ective density of the mixture and
g the acceleration of gravity.
Section 4 contains the results of a numerical approach
similar to that of Antal et al. (1991). In our case, the
presentation focuses on the con9rmation of the results of
Section 3 once wall e4ects are introduced by means of a
wall force. The numerical approach also allows us to extend
the study to the regimes where regions of strictly zero void
fraction are predicted.
Finally, Section 5 discusses the behavior of the two-%uid
model solutions under reduced gravity, which is a matter
of growing interest and conceptual importance because it
changes the relative weight of the interfacial forces, making
the lift term to become increasingly dominant. Some conclusions are drawn in Section 6.
We consider the mass and momentum conservation equations neglecting temperature e4ects (Drew & Passman,
1998),
@(k k )
+ ∇ · (k k˜vk ) = 0;
(1)
@t
@(k k˜vk )
+ ∇ · (k k˜vk˜vk )
@t
= ∇ · [k (Tk +
vG − ˜vL |2 I + B(˜vG − ˜vL )(˜vG − ˜vL )];
Re
L = −G L [A|˜
+ k k˜g + Mk ;
(2)
Mk = pki ∇k − ki · ∇k + Mk ;
ML = pLi ∇L − [L (∇vL + ∇vLT )
− L (A|˜vr |2 I + B(˜vr ⊗ ˜vr ))]∇L + ML
(Antal et al:; 1991);
Tk = −pk I + k [∇˜vk + (∇˜vk )T ];
(5)
ML = pLi ∇L − [L (ab|˜vr |2 I − a2 (˜vr · ˜vr ))]∇L
− L L ∇ · (˜vr ⊗ ˜vr )) + ML
(Drew & Passman; 1998);
(6)
where a = −9=20 and b = 3=20.
The pressure at the interphase (pLi ) can be derived from
the expression (Stuhmiller, 1977):
C = 14 :
(7)
We will consider drag, lift and wall forces, Mk = Mk D +
Mk L + Mk W , with
MGD = −MLD = −
3 G
CD L˜vr |˜vr |;
8 Rb
CD =
24
(1 + 0:1Re0:75 );
Re
Re =
2Rb L |˜vr |
;
m

G L |u |2
Rb


CW 1 + CW 2
nMw


Rb
y0
=



0
in which k refers to the phase (L for liquid, G for gas), k
is the volumetric (in fact, probabilistic) fraction of phase k,
˜vk the corresponding velocity, k the density, Tk the stress
tensor,
(4)
where subscript i refers to values at the interphase.
Di4erent models exist for Mk , in particular for the liquid
phase (k = L),
m =
(8)
(9)
L
;
(1 − G )
MGL = −MLL = −CL G L˜vr ∧ (∇ ∧ ˜vL );
MGW = −MLW
(3)
where customary values for A and B are 3=20 and 1=20,
respectively. We denote hereafter ˜vG −˜vL by ˜vr , and de9ne
vr = |˜vr |.
In Eq. (2), Mk stands for the interfacial momentum exchange, which relates to the actual force between phases
(Mk ) according to
pL − pLi = CL L |˜vr |2 ;
2. The two-uid model
Re
k )]
pk the pressure, and k the viscosity. Re
k is the Reynoldslike stress due to statistical %uctuations, modeled by
(Nigmatulin, 1979):
G L |u |2
Rb
CW 1 + CW 2
¿ 0;
if
Rb
y0
G L |u |2
Rb
CW 1 + CW 2
6 0;
if
Rb
y0
(10)
(11)
(12)
u = ˜vr − [nMw · ˜vr ]nMw ;
(13)
CW 1 = −0:1;
(14)
CW 2 = 0:147;
where y0 is the distance to the wall, nMw the exterior unit
normal and Rb the bubble radius.
O. E. Azpitarte, G. C. Buscaglia / Chemical Engineering Science 58 (2003) 3765 – 3776
Considering fully developed %ow in a vertical cylindrical
duct, the complete set of equations presented can be summarized in a system of 9ve di4erential equations corresponding to 9ve unknown functions, i.e. vL (r), vG (r), pL (r; z),
pG (r; z) and G (r), as follows:
• Momentum conservation equation—gas phase—radial
direction (r)
M
G
@pG
@vL
+ C L L vr
@r
@r
Rb
L vr2
= 0;
CW 1 + C W 2
+
Rb
(R − r)
(15)
where R is the internal radius of the duct.
• Momentum conservation equation—gas phase—axial direction (z)
M
@pG
3 CD
L vr |vr |:
= −G g −
@z
8 Rb
(16)
This equation is used for the calculation of the relative
velocity (vr ). This velocity will hereafter be considered
uniform. This is an excellent approximation because all
coe>cients in Eq. (16) are constants with the exception
of CD , which depends only weakly on G for diluted %ows
(see Eqs. (9)–(10)).
• Momentum conservation equation—liquid phase—radial
direction (r)
M
L
@
@pL
@vL
= − [AL L (1 − L )vr2 ] + CL G L vr
@r
@r
@r
+ CL (1 − G )vr2
+
G L vr2
Rb
@G
@r
CW 1 + CW 2
Rb
(R − r)
(17)
3 G
@pL 1 @
@vL
L
rL L
+
CD L vr |vr |
=
@z
r @r
@r
8 Rb
@vL @G
:
@r @r
(18)
• Jump condition at the interphase
pG − pL =
2&
− CL (1 − G )vr2 :
Rb
Let us de9ne
≡ e4 = L L + G G ;
O ≡
1
S
S
dS;
@P
@p
≡
+ g;
O
@z
@z
(20)
where S is the cross-sectional area of the pipe.
Eliminating the drag force from Eqs. (16) and (18), and
using L = ( − G )=(L − G ), one arrives at
@vL
@P
( − G ) L @
r
+ ( − )g
O =
@z
(L − G ) r @r
@r
+ (C2 − 1)L
@vL @G
:
@r @r
(21)
On the other hand, replacing Eqs. (17) and (19) into
Eq. (15) one obtains, after multiplication by L =(L vr ),
1
@
( − G )
(L − G ) (L − G ) @r
=
CL
@vL
(2C − A)vr @r
1
Rb
C W 1 + CW 2
+
Rb (2C − A)
R−r
−
(C1 + AG )
1
@
:
(2C − A) (L − G ) @r
(22)
@P ∗
@P
1
≡
;
@z
(L − G )g @z
• Momentum conservation equation—liquid phase—axial
direction (z)
M
− L L g + C2 L
3. System reduction, scaling and generic properties of the
solutions
The system is thus reduced to only two equations, which
motivates the introduction of the following non-dimensional
(scaled) variables:
@G
− C1 L vr2
:
@r
3767
(19)
L2 ≡
2L vr (C − A)
;
g(L − G )CL
v∗ =
CL
vL
;
vr 2(C − A)
r∗ ≡
r
;
L
R∗b ≡
Rb
;
L
D≡
R∗ ≡
E≡
CW 2
:
2(C − A)
(23)
R
;
L
CW 1
;
2(C − A)R∗b
(24)
Table 1 shows values of the intrinsic length scale L calculated for di4erent dispersed phases in stagnant water, assuming a constant bubble (or particle) radius of 0:5 mm and
a lift coe>cient CL = 0:1. To leading order, L depends only
3768
O. E. Azpitarte, G. C. Buscaglia / Chemical Engineering Science 58 (2003) 3765 – 3776
Table 1
Some typical values of the length scale L
Mixture
c (kg= m3 )
c (kg=(m s))
d (kg= m3 )
vr (m= s)
L (m)
Water–air
Water–CO2
Water–polystyrene
Water–wood
Water–oil
Water–glass
Water–steel
Water–vapor
998
998
998
998
998
998
998
958.12
0.001
0.001
0.001
0.001
0.001
0.001
0.001
0.000282
1.19
1.9022
55
545
874
2500
7900
0.598
0.1184
0.1183
0.1143
0.0714
0.0297
−0:153
−0:387
0.1545
1:557 × 10−4
1:557 × 10−4
1:573 × 10−4
1:793 × 10−4
2:211 × 10−4
1:442 × 10−4
1:07 × 10−4
0:964 × 10−4
Subscript c refers to the continuous phase and subscript d to the dispersed phase. The water-vapor system is considered at 100◦ C and atmospheric
pressure. The bubble (or particle) radius was assumed as 0:5 mm.
on the bubble radius and on the model constants C − A and
CL , and for reasonable values of these coe>cients is of the
order of the bubble radius.
From Eqs. (21)–(22), using the scaled variables in
Eqs. (23)–(24) and considering the wall-force terms as null,
the following ordinary di4erential equation is obtained:
@2 L
(L )2 ∗2 + (2 − C2 )L
@r
@L
@r ∗
2
(L )2 @L
(C1 + A)
− L +
r ∗ @r ∗
2(C − A)
2 1 @L
@2 L
@L
L
+ ∗2 − (C2 − 1)
r ∗ @r ∗
@r
@r ∗
@P ∗
− L ;
@z
(25)
O G )=(L −G ). The equation holds wherever
where L =(−
the wall force is not active (R − r ¿ − DL=E). Notice that
the right-hand side is a constant, independent of r ∗ .
We have 9xed the values of the constants A and B at their
usual values A = 3=20, B = 1=20, and considered two sets of
values for C1 and C2 , namely C1 = −3=20 and C2 = 1 taken
from Antal et al. (1991) (Model I), or C1 = −27=400 and
C2 = 0, taken from Drew and Passman (1998) (Model II).
Adopting, for example, the model of Antal et al., Eq. (25)
is reduced to
2
(L )2
@ L
+ L
@r ∗2
=
@P ∗
− L ;
@z
@L
@r ∗
2
+
@P ∗
− L = −L (0)(1 + ,);
@z
2
(L ) @L
− L
r ∗ @r ∗
(26)
which in some sense is the “universal” equation of laminar
fully developed bubbly %ow far from the wall, in the sense
that it contains no parameters except for the constant on the
right-hand side. Di4erent solutions of Eq. (26) are obtained
by changing two parameters: the pressure gradient and the
amount of gas in the %ow. For the analysis that follows we
(27)
9nally resulting in
2
2
@L
(L )2 @L
2 @ L
+
+
− L
(L )
L
@r ∗2
@r ∗
r ∗ @r ∗
= − L (0)(1 + ,):
+
=
will change the second parameter for the liquid fraction
at r ∗ = 0, L (0), and will rewrite the right-hand side of
Eq. (26) as
(28)
If, on the other hand, the model of Drew and Passman is
adopted, an analogous equation is derived:
2
2
@L
(L )2 @L
2 @ L
(L )
+
2
+
− L
L
@r ∗2
@r ∗
r ∗ @r ∗
2 1 @L
@L
33
@2 L
L
+
+ ∗2 +
80
r ∗ @r ∗
@r
@r ∗
= − L (0)(1 + ,):
(29)
We thus solve ODE (28)–(29) with initial conditions L (0)
given and L (0) = 0, for di4erent constant right-hand sides
which we parameterize with ,. Notice that
@pL
= −(0)g − ,[(1 − G (0))(L − G )g];
@z
(30)
so that , is a measure of the di4erence between the applied
pressure gradient and the e4ective weight of the mixture at
the center of the pipe.
Eqs. (28)–(29) are then solved using Mathematica
(Wolfram, 1991) for a wide set of values of the parameters
L (0) (equal to 1 − G (0)) and ,. Some solutions can be
found in Figs. 1–3. A generic fact about the biparametric
family of solutions of Eqs. (28)–(29) is that there exists a
limit value of r ∗ , which we will denote critical radius rc∗ ,
at which L becomes unphysical (either negative or greater
than one). In Fig. 4 we plot rc∗ as a function of |,| for different values of G (0). Consider for example the plot with
G (0) = 0:016. Remember that the length scale is of the
O. E. Azpitarte, G. C. Buscaglia / Chemical Engineering Science 58 (2003) 3765 – 3776
3769
1.0
0.9
0.8
0.7
0.6
Mode l I
0.5
0.0001
0.4
0.001
0.01
Mode l II
0.3
0.2
0.0001
0.1
0.001
0.0
0.01
0
2
4
6
8
10
12
14
r*
Fig. 1. G (r ∗ ) − , ¿ 0—upward %ow.
0.018
45
0.016
40
Mode l I
Upward flow
0.014
M ode l I
0.010
- 0.01
0.008
- 0.001
- 0.0001
0.006
M ode l II
0.004
- 0.01
0.002
- 0.001
Critical radius (r c*)
35
0.012
- 0.0001
0.000
0
1
2
3
4
5
6
7
8
9
30
0 .016
0 .1
25
0 .2
20
15
10
10
5
r*
0
Fig. 2. G (r ∗ ) − , ¡ 0—downward %ow.
1 E6
1 E5
1 E4
1 E3
0 .01
|λ|
Fig. 4. rc∗ as a function of |,| for di4erent G (0).
1. 0
Mode l I
0. 9
50
0.016
0. 8
45
0.1
0. 7
0.2
0. 6
Critical radius (r c*)
Mode l II
0.016
0. 5
0.1
0. 4
0.2
0. 3
0. 2
0. 1
40
Mode l I -upward flow
35
Mode l I -downward flow
Mode l II -upward flow
30
Mode l II -downward flow
25
20
15
10
0. 0
5
0
1
2
3
4
5
6
7
8
9
10
11
12
13
0
r*
1E 6
Fig. 3. G (r ∗ ) for di4erent G (0).
1E 5
1E 4
1E 3
0. 01
|λ|
Fig. 5. rc∗ as a function of |,|.
order of Rb , and to 9x ideas assume that the radius of the
pipe you are dealing with is 40 times Rb . This is indeed
a small-diameter pipe. However, from Fig. 5 it is evident
that unless |,| ¡ 10−6 there exists no solution to the governing equations (28)–(29) that remains between 0 and 1
throughout the section of the pipe. In all reasonable cases
(that is, pipes that are signi9cantly bigger than the bubble
radius) the band of “allowed” values for ,, which is centered at zero, is extremely narrow. What does this mean?
3770
O. E. Azpitarte, G. C. Buscaglia / Chemical Engineering Science 58 (2003) 3765 – 3776
It means that in fully developed laminar %ow the solutions
of the two-%uid model (the only physically meaningful solutions that exist) are such that the e4ective weight at the
center of the pipe (i.e., e4 ˜g) practically balances the applied pressure gradient. By “practically” we mean that the
di4erence, as in the example above, is one part per million
or less.
This striking result allows us to prove, as a corollary, that:
0.0
-0.1
Consequence 2. The %at void fraction pro9le away from the
wall of the pipe has a value of L (and thus of G ) that only
depends on the applied pressure gradient, the density of the
phases, and the gravity. In particular, the value of G does
not depend on (a) the gas %ow rate, (b) the viscosity of the
phases, or (c) the slip velocity. The proof is immediate from
Eq. (30).
Once L (r ∗ ) is obtained, the equation
1 @
@v∗
− L (0)(1 + ,) = L ∗ ∗ r ∗ ∗ − L ;
r @r
@r
(31)
with initial condition @vL∗ =@r ∗ (r ∗ = 0) = 0, yields the liquid
velocity up to a constant (the value at the center). It is
again solved using Mathematica, with some sample results
shown in Figs. 6 and 7. With a proof analogous to that of
Consequence 1 we get
Consequence 3. The %at velocity pro9les typical of
two-phase pipe %ow are the only possible solutions of the
two-%uid model (unless the diameter of the pipe is too
small). In particular, the %atness of two-phase laminar %ow
is not linked to turbulent-like e4ective stresses. It simply
arises from a balance of the pressure gradient with the
e4ective weight of the %uid column.
0 .0001
v*
0 .001
-0.3
0 .01
Mode l II
0 .0001
-0.4
Consequence 1. The %at void fraction pro9les typical of
two-phase pipe %ow are the only possible solutions of the
two-%uid model (unless the diameter of the pipe is too
small).
This can be proved simply by inspecting the solutions of
Eqs. (28)–(29) which have |,| ¡ 10−6 , and con9rming that
all of them are extremely %at (as functions of r ∗ ) until very
close to rc∗ , at which point they exhibit a rapid growth (if ,
is positive) or decrease (if , is negative) that take the value
of G above one or below zero. Another, more heuristic way,
of proving that void fraction pro9les are necessarily %at,
once we know that |,| ¡ 10−6 , comes from a term-by-term
analysis of Eqs. (28)–(29). The fact that |,| is small implies
the “practical” cancellation, at r ∗ = 0, of the fourth term
on the left-hand side of Eqs. (28)–(29) with its right-hand
side. This leaves us with the remaining terms, all of them
containing derivatives of L , equated to zero (“practically”)
and with the boundary condition L (0) = 0. The solution
to this problem is (again “practically”, in the sense already
explained) a constant, and thus a %at void fraction pro9le.
This is not the only consequence of , being necessarily
very small.
Mode l I
-0.2
0 .001
0 .01
-0.5
0
2
4
6
8
10
12
14
r*
Fig. 6. v∗ (r ∗ ) − , ¿ 0—upward %ow.
0 .024
.
0 .022
0 .020
0 .018
0 .016
Mode l I
0 .014
v *
0 .012
−0.01
0 .010
−0.001
−0.0001
0 .008
Mode I lI
0 .006
−0.01
0 .004
−0.001
0 .002
−0.0001
0 .000
0
1
2
3
4
5
6
7
8
9
r*
Fig. 7. v∗ (r ∗ ) − , ¡ 0—downward %ow.
Remarkably, though our analysis concerns the equation
without the wall-force terms, we can conclude about
near-wall void-fraction pro9les. The consequence is thus
independent of the wall–bubble interaction model:
Consequence 4. In upward bubbly %ows there must exist a
so-called “peak” in the void fraction in the close vicinity of
the wall.
The proof is not di>cult. It is evident that, for the gas-liquid
mixture to %ow upward, the applied pressure gradient
must overcome the horizontally-averaged weight ˜
Og, or, in
non-dimensional variables, we must have
−
@P ∗
¿ L :
@z
But, since , is very small, −@P ∗ =@z L (0), leading to
L (0) ¿ L and thus to G (0) ¡ G . On the other hand, we
know that G is (practically) constant except very near the
pipe wall. Thus, in most of the %ow G takes the value G (0)
which is less than the mean void fraction. It is obvious that
a signi9cant accumulation of bubbles must take place at the
wall, leading to a “peak” in the void fraction.
O. E. Azpitarte, G. C. Buscaglia / Chemical Engineering Science 58 (2003) 3765 – 3776
In the case C ¡ A, the a4ected terms in the scaled equation
change sign, yielding (Model I):
2
2
@L
(L )2 @L
2 @ L
−(L )
−
−
− L
L
@r ∗2
@r ∗
r ∗ @r ∗
= − L (0)(1 + ,);
(32)
which accumulates bubbles at the center of the pipe in upward %ow, just as if the lift were negative. It is remarkable
that coe>cients A (related to the e4ective Reynolds-like
stress) and C (related to the interphase pressure correction)
are mathematically linked in such a strong way. One cannot
switch on and o4 the coe>cients A and C arbitrarily, since
their di4erence is as important as the lift coe>cient itself.
We add to both equations a time derivative term. Although we are considering a stationary case, treating it as
transient allows the variables to adapt to their 9nal value in
a progressive and numerically more stable way.
We also add an arti9cial di4usion term to avoid undershoots that would lead to negative values of G .
The resulting system is, thus
L
@pL
+ (1 − G )L g + G G g
@z
@vL
@vL @G
1 @
r(1 − G )L
+ C2 L
=
:
r @r
@r
@r @r
(33)
Combining Eqs. (15), (17) and (19), and after mathematical
manipulation, the di4erential equation for G (r) results
1
1
@G 2
vr
− A + C1 − G
G
@r
2
5
2
vr
@vL
Rb
= − G CL vr
:
(34)
− G Cw1 + Cw2
@r
R − r Rb
The following boundary conditions are imposed to the
system:
@vL
(0) = 0;
@r
vL (R) = 0;
@G
(0) = 0;
@r
with the additional requirement:
1
G dS = G ; where G is given:
S S
(35)
(36)
4.1. Numerical method
Eqs. (33) and (34) have been slightly modi9ed in order
to simplify their numerical treatment.
@vL
@pL
+
+ (1 − G )L g + G G g
@t
@z
@vL
@vL @G
1 @
r(1 − G )L
+ C 2 L
=
;
r @r
@r
@r @r
(37)
1
1
@2 G
− A + C1 − G − K
2
5
@r 2
2
vr
@vL
Rb
− G Cw1 + Cw2
= − G C L v r
;
(38)
@r
R − r Rb
@G
@G 2
+ G
v
@t
@r r
4. Numerical simulation and parametric study
The analysis of the previous section refers to the region
where the wall forces are not active. Moreover, not all applied pressure gradients have been considered, since when
−@pL =@z is greater than L g the pressure gradient cannot be
balanced by adjusting the void fraction at the center, even
if the void fraction there is strictly zero. In the present section we report on some numerical studies made using 9nite
elements to understand further the behavior of the system.
Eliminating the drag term from Eqs. (16) and (18), the
di4erential equation for vL (r) results:
3771
where K is the arti9cial di4usion coe>cient.
The time-marching scheme is semi-implicit:
L
@pL
vLn − vLn−1
+
+ (1 − Gn−1 )L g + Gn−1 G g
Tt
@z
@vLn
@vn @n−1
1 @
n−1
r(1 − G )L
+ C2 L L G ;
=
r @r
@r
@r @r
(39)
n∗
Gn∗ − Gn−1
1
1 n−1
n−1 @G 2
v
− A + C1 − G
+ G
Tt
@r r
2
5
2
@2 Gn∗
vr
Rb
+ Gn∗ Cw1 + Cw2
−K
R − r Rb
@r 2
= − Gn−1 CL vr
@vLn−1
;
@r
(40)
where Gn∗ is an intermediate variable that, at every time step,
is corrected to satisfy condition (36):
Gn = Gn∗ + (G − Gn∗ ):
De9ning
G=
G(n−1) vr2
1
− A + C1
2
(41)
2
vr
Rb
H = Cw1 + Cw2
;
R − r Rb
1 (n−1)
− G
;
5
(42)
(43)
in Eq. (40), G plays the role of a “velocity” and H the role of
a “reaction”. The di4usion coe>cient K is thus selected as
usual for reaction-convection problems (e.g. Codina, 1998):
Gh
if Hh2 6 Gh;
K=
(44)
Hh2 if Hh2 ¿ Gh:
3772
O. E. Azpitarte, G. C. Buscaglia / Chemical Engineering Science 58 (2003) 3765 – 3776
0 .08
0 .10
0 .07
0 .09
0 .08
Liquid velocity (m/s)
0 .06
Mo de l
0 .05
Experimental data
0 .04
0 .03
0 .02
0 .07
0 .06
0 .05
Mode l
0 .04
Experimental data
0 .03
0 .02
0 .01
0 .01
0 .00
0 .0
0. 1
0 .2
0. 3
0 .4
0. 5
0 .6
0. 7
0 .8
0. 9
1 .0
0 .00
0 .0
0. 1
0 .2
0. 3
0 .4
r/ R
0. 5
0 .6
0. 7
0 .8
0. 9
1 .0
0 .8
0. 9
1 .0
r/R
Fig. 8. Void fraction pro9le.
Fig. 9. Liquid velocity pro9le.
∆
0 .06
z
In this article, from now on, we will use the Drew and Passman model (Model II). Although the Antal et al. model renders very similar results, it does not conserve the momentum
in the axial direction.
p / ρL g
− 0 .992
0 .05
− 0 .994
0 .04
− 0 .996
− 0 .998
0 .03
− 1 .000
0 .02
4.2. Comparison with experimental data
0 .01
In order to verify the numerical model, its results are
compared to the experimental data reported by Nakoryakov
et al. (1986), as previously done by Antal et al. (1991).
Important parameters of the experiment are: Flow:
bubbly, laminar, vertical, upward; Mixture: air–water;
Pressure: atmospheric; Pipe diameter: 1:5 cm; Bubble diameter: 0:87 mm; Reynolds number (liquid phase): 1267; Area
average void fraction (G ): 0.019.
The physical parameters are the following: L =
1000 kg=m3 ; G = 1:19 kg=m3 ; L = 0:001 kg=(m s);
@p=@z = −9649 Pa=m. For the simulation we used CL = 0:1,
vr = 0:1028 m=s, G = 0:019.
Fig. 8 shows the comparison of the void fraction pro9les.
As it was already pointed out, the pro9le in the central region of the pipe is %at (Consequence 1). The pro9le presents
a peak close to the wall (Consequence 4) before going to
zero under the e4ect of the wall force. The numerical void
fraction at the center of the pipe is G (0) = 0:0154332, so
that
−
Pa
@p
= 9649
≈ [(1 − G (0))L + G (0)G ]g
@z
m
= 9648:935
Pa
;
m
(45)
and Consequence 2 is con9rmed.
From Eq. (30) one gets , = 6:77 × 10−6 .
Fig. 9 shows the predicted liquid velocity pro9le, which
is indeed %at though no turbulence was considered.
0 .00
0 .0
0. 1
0 .2
0. 3
0 .4
0. 5
0 .6
0. 7
r/R
Fig. 10. Void fraction pro9les—upward %ow.
Finally, the average liquid velocity was calculated:
1
m
vL =
vL dA ≈ 0:08215 ;
(46)
A A
s
and the resulting liquid phase Reynolds number was Re =
1232:2, close to the experimental value (Re = 1267).
4.3. Parametric study
A parametric study was carried out using a general case,
in order to analyze the behaviour of the main variables as a
function of the applied pressure gradient and of the average
void fraction (G ).
The parameters used are the following: L = 1000 kg=m3 ;
G = 1:19 kg=m3 ; L = 0:01 kg=(m s); CL = 0:1.
Fig. 10 shows the void fraction pro9les as a function of
the applied pressure gradient, for upward %ow and constant
average void fraction G . As the absolute value of the applied
pressure gradient increases, the void fraction at the central
region of the pipe diminishes, the peak increases its height
and the radial position of the maximum moves towards the
wall. This behaviour is due to an increasing lift force that
obeys to an increasing liquid velocity gradient at the region
O. E. Azpitarte, G. C. Buscaglia / Chemical Engineering Science 58 (2003) 3765 – 3776
0 .05
∆
Liquid velocity (m/s)
0 .04
0 .03
z
p / ρL g
−0 .992
0 .02
−0 .994
−0 .996
−0 .998
0 .01
3773
Fig. 16 shows the total mass %ow rate as a function of
the applied pressure gradient, for di4erent values of G . As
the absolute value of the applied pressure gradient increases,
a proportional increase in the mass %ow rate is observed.
When the pressure gradient balances the speci9c weight of
the liquid phase, a change of regime occurs: the %at velocity
pro9le changes to a parabolic-type pro9le (see Fig. 11), and
the slope of the mass %ow rate increases.
−1 .000
−1 .003
5. E%ect of reducing the acceleration of gravity
0 .00
0 .0
0. 2
0 .4
0. 6
0 .8
1. 0
r/R
Fig. 11. Liquid velocity pro9les—upward %ow.
0 .020
0 .018
0 .016
0 .014
∆
z
0 .012
p / ρL g
0 .010
−0 .985
0 .008
−0 .984
0 .006
−0 .983
−0 .982
0 .004
−0 .981
0 .002
0 .000
0 .0
0. 1
0 .2
0. 3
0 .4
0. 5
0 .6
0. 7
0 .8
0. 9
1 .0
r/R
Fig. 12. Void fraction pro9les—downward %ow.
close to the wall (see Fig. 11). From the results shown in
Fig. 10, another important conclusion can be drawn: when
−@p=@z is equal or greater than the speci7c weight of the
liquid phase (L g), the void fraction at the central region
becomes null.
Figs. 12 and 13 show the void fraction and velocity pro9les for several applied pressure gradients, for downward
%ows. The lift force accumulates the bubbles at the central
region, leaving a depleted zone close to the wall.
It has been previously demonstrated that , is necessarily
very small (Section 3). The parametric study allows us
to verify and extend this conclusion: , is very small
(|,| 6 10−6 ) when G (0) = 0 (−@p=@z 6 L g). When
G (0) = 0 (−@p=@z ¿ L g), , depends linearly on the pressure gradient, as can be inferred from Eq. (30). This is
shown clearly in Fig. 14 for di4erent values of G .
Considering the fact that , is very small, Eq. (30) reduces
to
@p
= g(L − G )G (0) − gL :
@z
(47)
Notice that, in this case, G (0) depends linearly on the applied pressure gradient (see Fig. 15).
The implemented 9nite element code is also used to analyze the solutions rendered by the presented two-%uid model
under reduced gravity conditions. Therefore, numerical solutions of a general case are obtained decreasing g, for constant G and wall shear stress (w ).
The acceleration of gravity is explicitly present in the
momentum conservation equation of the gas phase in axial
direction (Eq. (16)), used for the algebraic calculation of the
relative velocity. As shown in Fig. 17, the relative velocity
decreases when g decreases.
The in%uence of a decreasing g on the void fraction pro9le
can be deduced analyzing the equation corresponding to
G (r) (Eq. (34)). Dividing the whole equation by vr2 , only
the lift term remains divided by vr . Therefore, when g and vr
decrease, the lift term becomes increasingly dominant over
the other terms of the equation.
Figs. 18 and 19 show the void fraction pro9les predicted
for decreasing g: the relative increasing lift force causes the
bubbles to accumulate towards the wall for upward %ow,
and to accumulate at the central region of the pipe for downward %ow.
In the equation corresponding to vL (r) (Eq. (33)), g is
also explicitly present. Expressing the pressure gradient as
the sum of a frictional and a gravitational term:
@p
@p @p +
=
@z
@z F
@z G
@p =
− [(1 − G )L + G G ]g;
(48)
@z F
where
−
@p S
= w ;
@z F A
S: perimeter;
Eq. (33) can be rewritten as
@pL + (G − G )(L − G )g
@z F
1 @
@vL
r(1 − G )L
:
=
r @r
@r
(49)
(50)
Fig. 20 shows the velocity pro9les for di4erent values of g
in the case of upward %ow. The void fraction in the central
3774
O. E. Azpitarte, G. C. Buscaglia / Chemical Engineering Science 58 (2003) 3765 – 3776
0 .000
− 0 .002
− 0 .004
Liquid velocity (m/s)
− 0 .006
− 0 .008
∆
− 0 .010
z
p / ρL g
− 0 .012
−0 .985
− 0 .014
−0 .984
−0 .983
− 0 .016
−0 .982
− 0 .018
−0 .981
− 0 .020
0 .0
0. 2
0 .4
0. 6
0 .8
1. 0
r/ R
Fig. 13. Liquid velocity pro9les—downward %ow.
0 .0035
0 .0040
0 .0030
0 .005
0 .0035
0 .01
0 .0025
0 .0030
0 .02
0 .005
Mass flow rate (kg/s)
0 .01
0 .0020
0 .02
λ
0 .0015
0 .0010
0 .0005
0 .0000
0 .0025
0 .0020
0 .0015
0 .0010
0 .0005
0 .0000
0 .0005
0 .996
0 .997
∆
∆
Fig. 14. , for di4erent values of pressure gradient.
0 .998
z
0 .999
1 .000
1 .001
1 .002
1 .003
1 .004
− 0 .968
− 0 .970
− 0 .972
− 0 .974
− 0 .976
− 0 .978
− 0 .980
− 0 .982
− 0 .984
− 0 .986
− 0 .988
− 0 .990
− 0 .992
− 0 .994
− 0 .996
− 0 .998
− 1 .000
− 1 .002
− 1.006
− 1 .004
p / ρL g
z
p/ ρ L g
Fig. 16. Mass %ow rate for di4erent values of pressure gradient.
0 .025
0 .030
0 .025
0 .005
0 .020
0 .01
0 .020
Relative velocity (m/s)
0 .02
0 .015
0 .010
0 .005
0 .015
0 .010
0 .005
0 .000
− 0 . 968
− 0 . 970
− 0 . 972
− 0 . 974
− 0 . 976
− 0 . 978
− 0 . 980
− 0 . 982
− 0 . 984
− 0 . 986
− 0 . 988
− 0 . 990
− 0 . 992
− 0 . 994
− 0 . 996
− 0 . 998
− 1 . 000
− 1 . 002
− 1. 006
− 1 . 004
∆
z
p / ρL g
Fig. 15. G (0) for di4erent values of pressure gradient.
region of the pipe is smaller than the average void fraction
(G ¡ G ), and the second term on the left-hand side of Eq.
(50) is positive. As g decreases, so does the void fraction
in the central region, and the second term on the left-hand
0 .000
0
1
2
3
4
5
6
7
8
9
10
2
g (m /s )
Fig. 17. Relative velocity for di4erent values of g.
side increases. The result is a reducing velocity (g = 9:8, 7,
5 m=s2 ). Nevertheless, when the void fraction in the central
region becomes null, the second term on the left-hand side
decreases as g is reduced, and the result is an increasing
velocity (g=3, 1 m=s2 ). Notice that the slope of the di4erent
O. E. Azpitarte, G. C. Buscaglia / Chemical Engineering Science 58 (2003) 3765 – 3776
0 .05
3775
0 .000
0 .04
− 0 .005
2
0 .03
L iquid velocity (m/s)
g ( m/s )
1 .0
5 .0
9 .8
0 .02
0 .01
− 0 .010
2
g (m/s )
1 .0
− 0 .015
3 .0
5 .0
− 0 .020
7 .0
9 .8
0 .00
− 0 .025
0 .60
0 .65
0 .70
0 .75
0 .80
0 .85
0 .90
0 .95
0 .0
1 .00
0. 1
0 .2
0. 3
0 .4
r/R
0. 5
0 .6
0. 7
0 .8
0. 9
1 .0
r/R
Fig. 18. Void fraction pro9les for di4erent values of g—upward %ow.
Fig. 21. Liquid velocity pro9les for di4erent values of g—downward %ow.
0 .00012
0 .09
0 .08
0 .002
0 .00011
0 .07
0 .004
2
Mass flow rate (kg/s)
g (m/s )
0 .06
1 .0
0 .05
3 .0
0 .04
5 .0
0 .03
7 .0
9 .8
0 .02
0 .006
0 .00010
0 .00009
0 .00008
0 .00007
0 .01
0 .00
0 .00006
0 .0
0. 1
0 .2
0. 3
0 .4
0. 5
0 .6
0. 7
0 .8
0. 9
1 .0
0
1
2
3
4
r /R
5
6
2
7
8
9
10
g (m /s )
Fig. 19. Void fraction pro9les for di4erent values of g—downward %ow.
Fig. 22. Mass %ow rate for di4erent values of g—upward %ow.
−0.00050
0 .0030
−0.00055
0 .0025
0 .002
−0.00060
2
g (m/s )
0 .0015
1 .0
0 .0010
3 .0
5 .0
0 .0005
7 .0
Mass flow rate(kg/s)
Liquid velocity (m/s)
0 .004
0 .0020
−0.00065
0 .006
−0.00070
−0.00075
−0.00080
−0.00085
9 .8
−0.00090
0 .0000
0 .0
0. 1
0 .2
0. 3
0 .4
0. 5
0 .6
0. 7
0 .8
0. 9
1 .0
r /R
Fig. 20. Liquid velocity pro9les for di4erent values of g—upward %ow.
−0.00095
0
1
2
3
4
5
6
7
8
9
10
2
g ( m /s )
Fig. 23. Mass %ow rate for di4erent values of g—downward %ow.
pro9les are equal at the wall, corresponding to a constant
wall shear stress. In this case, @pL =@z|F = −10 Pa=m and
w = 0:025 Pa were used.
In the case of downward %ow, the void fraction in the
central region of the pipe is greater than the average void
fraction (G ¿ G ), and the second term on the left-hand side
of Eq. (50) is negative. As a consequence, |vL | at the central
region increases as g decreases (see Fig. 21). Notice, also,
that the slope of the di4erent pro9les are equal at the wall,
corresponding to a constant wall shear stress. For downward
%ow, @pL =@z|F = 78 Pa=m and w = −0:195 Pa were used.
Figs. 22 and 23 show the mass %ow rate as a function of
g for upward and downward %ow respectively.
3776
O. E. Azpitarte, G. C. Buscaglia / Chemical Engineering Science 58 (2003) 3765 – 3776
It is important to stress that the results obtained in this section are nothing but the numerical outcome of the two-%uid
model as presented in Section 2. We did not make g tend
to zero since the trend is already evident from the case analyzed. Comparison with appropriate experimental data will
reveal the improvements needed, if any, to accurately apply
the model in reduced-gravity conditions.
6. Conclusions
An analysis of the two-%uid model in the very simple case
of a laminar fully developed bubbly %ow in a circular pipe
has allowed us to identify several physical and mathematical
properties of the model.
All experimental results in bubbly %ows show %at velocity and void-fraction pro9les (Serizawa, Kataoka, &
Michiyoshi, 1986; Nakoryakov et al., 1986; Wang, Lee,
Jones, & Lahey Jr., 1987). In particular, the velocity pro9les are %atter than in single-phase %ows (Wang et al.,
1987). The mechanism that %attens the velocity pro9le in
single-phase turbulent %ow is, as is well known, the increase of the eddy viscosity with the distance from the wall.
Our results have shown that in vertical two-phase %ows an
additional mechanism is predicted by the two-%uid model,
namely the cancellation of the pressure gradient with the
e4ective speci9c weight of the %uid column. This mechanism is not related to the Reynolds-like stresses of Eq. (3),
and is a consequence of the model which is independent
of turbulence or %uctuations. In fact, notice that in Section
3 we prove that the %at velocity and void fraction pro9les
are the only valid solutions of the two %uid model. This
result remains valid in the case in which the Reynolds-like
stresses are neglected (A = 0).
The numerical results of Section 4 provide further insights
into the physics of the model. One observes how the lift force
tends to distribute the bubbles so as to balance the applied
pressure gradient. For high enough pressure gradients in
upward %ows this is no longer possible and a new %ow
regime appears with a pure-liquid core. In downward %ows,
on the other hand, the pure-liquid region occurs close to the
walls and grows thicker with increasing %ow rate.
Finally, the model predictions under reduced gravity have
been presented. The smaller the gravity, the greater the e4ect
of the lift force. As the lift becomes dominant, the bubbles
accumulate towards the wall (in upward %ows) so that very
high void fractions are to be expected.
Future work will be aimed at extending these results
to the turbulent case, coupling the two-%uid model with a
low-Reynolds-number k–j turbulent model.
Acknowledgements
GCB is also a fellow of CONICET, Argentina. Partial
support was received from ANPCYT, Argentina, through
grant PICT’99 No. 6337. Fruitful discussions with R.T. Lahey Jr., A. Larreteguy, F. Moraga and F. Bonetto are gratefully acknowledged. This paper is dedicated to the late Prof.
Ben C. Yen. He was an inspiration to all of us who had the
good fortune to be associated with him.
References
Antal, S., Lahey Jr., R., & Flaherty, J. (1991). Analysis of phase
distribution in fully developed laminar bubbly two-phase %ow.
International Journal of Multiphase Flow, 17(5), 635–652.
Codina, R. (1998). A comparison of some 9nite element methods for
solving the di4usion-convection-reaction equation. Computer Methods
in Applied Mechanics and Engineering, 156, 185–210.
Drew, D., & Lahey Jr., R. (1987). The virtual mass and lift force on a
sphere in rotating and straining inviscid %ow. International Journal
of Multiphase Flow, 13(1), 113–121.
Drew, D., & Lahey Jr., R. (2001). The analysis of two-phase %ow and
heat transfer using a multidimensional, four 9eld, two-%uid model.
Nuclear Engineering and Design, 204, 29–44.
Drew, D., & Passman, S. (1998). Theory of multicomponent %uids. In
Applied Mathematical Sciences, vol. 135. Berlin: Springer.
Ishii, M. (1975). Thermo-<uid dynamic theory of two-phase <ow:
Direction des Etudes et Recherches d’ElectricitXe de France: Eyrolles.
Lahey Jr., R. (1990). The analysis of phase separation and phase
distribution phenomena using two-%uid models. Nuclear Engineering
and Design, 122, 17–40.
Larreteguy, A., Drew, D., & Lahey Jr., R. (2002). A
particle-centre-averaged two-%uid model for wall-bounded bubbly
%ows. In Fifth international symposium on numerical methods for
multiphase <ow at the 2002 ASME <uids engineering division summer
meeting. New York: The American Society of Mechanical Engineers.
Nakoryakov, V., Kashinski, O., Kozemenko, B., & Gorelik, R. (1986).
Study of upward bubbly %ow at low liquid velocities. Izvestiya
sibirskogo Otdeleniya. Akademik Nauk SSSR, 16, 15–20.
Nigmatulin, R. (1979). Spatial averaging in the mechanics of
heterogeneous and dispersed systems. International Journal of
Multiphase Flow, 4, 353–385.
Serizawa, A., Kataoka, I., & Michiyoshi, I. (1986). Phase distribution in
bubbly %ow. Proceedings of the second international workshop on
two-phase <ow fundamentals, Data Set No. 24.
Stuhmiller, J. (1977). The in%uence of interfacial pressure on the
character of two-phase %ow model equations. ASME symposium on
computational techniques for non-equilibrium two-phase phenomena,
p. 118.
Uchiyama, T. (1999). Petrov–Galerkin 9nite element method for
gas-liquid two-phase %ow based on an incompressible two-%uid model.
Nuclear Engineering and Design, 193, 145–157.
Wang, S., Lee, S., Jones Jr., O., & Lahey Jr., R. (1987). 3-D
turbulence structure and phase distribution measurements in bubbly
two-phase %ows. International Journal of Multiphase Flow, 13(3),
327–343.
Wolfram, S. (1991). Mathematica. Reading, MA: Addison-Wesley.