Chemical Engineering Science, Vd. Printed in Great Britain.
45, No.
11. pp. 3367-3371.
NONLINEAR
1990.
oco!z-2so9/90 $3.00 + 0.00 0 1990 Pcrgamoa Press plc
STABILITY AND OF ANNULAR
DYNAMIC FLOW
SIMULATION
DVORA BARNEA and YEHUDA TAITEL Faculty of Engineering,Department of Fluid Mechanics and Heat Transfer, Tel-Aviv University, Ramat-Aviv (Received 15 Nouember
69978, Israel
1989; accepted 3 January 1990)
Abstract-Formulation of the two-fluid model may yield multiplesteady-statesolutions for cocurrent and countercurrent annular flow. The nonlinear stability of these solutions was studied by examining the system response to finite disturbances. The results show that steady-state solutions that are linearly stable may be unstable for finite disturbances and therefore are not expectedto exist. It was found that only the solutions that are associated with a thinner film are stable.
INTRODUCI’ION
The solution of steady-state cocurrent and countercurrent annular flow, using a two-fluid model in which an effective interfacial shear stress is used to account for the wavy structure of the liquid film, may yield multiple steady-state solutions, some of which are unstable, and therefore would not exist. A stability examination is needed to determine the stability characteristics of these steady-state solutions. A linear stability analysis of annular flow using the two-fluid model, neglecting surface tension (Stewart and Wendroff, 1984; Ardron, 1980, Lyczkowski et al., 1978; Ramshaw and Trapp, 1978), shows, as well known, that the steady-state solutions are always unstable owing to a Kelvin-Helmholtz (K-H) type instability. Recently, Bamea and Taitel(1989) claimed that the fully two-fluid model is not suited to examine the stability of steady annular flow where the solution is obtained for an average film thickness, using an effective interfacial shear to account for the wavy interface. They suggested that the stability of a steady-state structure be examined using a formulation that ignores the K-H instability by eliminating the Bernoulli amplification from the transient equations. Bamea and Taitel (1989) presented dynamic formulations that satisfy this requirement and performed a linear stability analysis. It seems that these dynamic formulations are the ones that are suited to answer the question whether annular flow with an unstable wavy interface is realized physically. In the present work the linear stability analysis is briefly reviewed. Next, nonlinear stability is examined to show that solutions that may be linearly stable are unstable for large disturbances and therefore are not expected to exist.
two separate fluids with average and constant crosssectional areas and velocity profiles. The information regarding local gradients and interfacial distribution is essentially lost in this approach. Thus, effective wall and interfacial shear stresses are required in order to replace this information. It has been shown that multiple steady-state solutions may be obtained, some of which are unstable. Various modes of dynamic formulation have been presented by Bamea and Taitel(1989) for the stability analysis of these steady-state solutions. The transient formulation which is consistent with the aforementioned steady-state solutions is provided by the “twofluid model” (see Fig. 1). In this model transient momentum and continuity equations are applied to each phase, leading to four simultaneous partial differential equations for the liquid holdup a=, the average axial liquid and gas velocities uL and z+, and the pressure P (Barnea and Taitel, 1989). The linear stability analysis of these equations yields the following criterion for stability (Bamea and Taitel, 1989): 2
p2aLaG
-
UL)2
(1)
where
--_ aTiG
LINEAR
PLPG
+ -(lJG
uG
au‘2uG >
STABILITY
The linear stability of steady-state cocurrent and countercurrent annular flow has been analyzed recently by Bamea and Taitel (1989). The steady-state solutions were calculated by the commonly used integrat approach, which treats the gas and the liquid as
p=%+po. &L
OIG
qL is the interfacial shear stress needed for a steadystate flow for a given liquid velocity IJ&,and film
3367
DVORABARNEAand YEHUDA
3368
TAITEL
tions are written for the liquid phase for a finite pipe length 1 (Barnea and Taitel, 1989):
(5)
da, -= dt
where uLs is the liquid superficial flow rate (positive for upward flow and negative for downward flow). For countercurrent flow the sign of the first term on the RHS of eqs (5) and (6) should be reversed. Note that the steady-state solutions of eqs (5) and (6) are the same as those obtained from the two-fluid model. A linear stability for eqs (5) and (6) leads to simple criteria for the stability of steady-state annular flow (Barnea and Taitel, 1989):
Fig. 1. Geometry of annular flow.
thickness 6: (Pi -
P&Q
sin/i + y (3)
riL = si($+d)”
whereas ziG is the shear stress provided by the gas which depends on the gas velocity uG and the film thickness S: ?iG
=
fr;:A&
f (U, - wk)
(4)
Note that in a steady state t= = riG = zr. Criterion (1) indicates that the how is always unstable. The second term on the LHS can be observed as the well-known K-H instability of the interface. The first term is the additional effect of the wall and the interfacial shear stresses which always enhance the K-H instability. Obviously the existence of K-H instability in steady annular flow is expected 4 priori and steady annular flow is dynamically unstable with regard to its interface, resulting in an unstable wavy interface. However, our question of interest is whether this kind of steady annular flow is a stable structure with respect to its average film thickness and the average phase velocities, and whether it is physically possible. The dynamic formulations of the two-fluid model are therefore unsuitable for studying the stability of steady annular flow, where the solutions are obtained for an average film thickness using the effective interfacial shear to account for the wavy interface. In order to ignore the K-H instability, the Bernoulli amplification should be eliminated from the transient formulation. Barnea and Taitel (1989) suggested few options that satisfy this requirement and showed that annular flow can be linearly stable. We will further consider the option in which the film thickness is a function of time only. Using this approach transient continuity and momentum equa-
2
llrJ > 2
$I”,.
-= zj”_
I_
for cocurrent flow
(7)
for countercurrent flow. (8)
As one can see the stability criterion is not restricted to a certain form of the correlation of either the interfacial shear stress ti or the wall shear stress rr. It is felt that this kind of stability criterion is indeed the one that is realized physically. Note that this criterion was also obtained for other formulations that ignore the K-H instability (Barnea and Taitel, 1989). The stability conditions (7) and (8) can be easily applied by using the steady-state diagrams, where tr is plotted vs 6 (Figs 2 and 3). The solid lines represent the steady-state relation (3) while the dashed lines represent the shear stress provided by the gas [eq. (4)]. The intersections of the solid and the dashed lines yield the steady-state solutions. The linear stability of these solutions can be determined by inspection of the
.02 0
0
.02
.04
Fig. 2. Steady-state solution 0.1 MPa, 25”C, air-water,
8/D
.c%
and
.08 stability
2.5 cm 0.005(1 + 3006/D).
.I0 diagram:
diameter, f, =
Nonlinear stability and dynamic simulation of annular flow
3369
physically unrealistic even if it is linearly stable. The nonlinear stability analysis will shed some new light on the stability of the steady-states solutions in this region. .06
NO-NEAR
0
0
.02
.04
06
06
0.1
0.12
8/D Fig. 3. Steady-state solution and stability diagram: air-water, 0.1 MPa, 25”C, 2.5 cm diameter, f, = 0.05.
of the slopes of the solid and the dashed lines at the steady-state solutions. Figures 2 and 3 are examples for steady annular flow for two different interfacial friction factorsf,. The wall shear stress rL is related to the liquid average axial velocity uL by
behavior
XL =
f PLIULIUL L
2
where, for a smooth pipe, fL = CL(D,uL/vL)-“, CL = 0.046 and n = 0.2 for turbulent flow, and CL = 16 and n = 1 for laminar flow. Note that for Si = 4A,/&D, steady annular flow S, = 4A/D, D, = 4A,/SL = a,D, uL = uLs/uL and D,uL = Du,. in the case of countercurrent flow, two steady-state solutions may exist for any given pair of uLs and uGs. In Fig. 2, for example, solution A which corresponds to the thinner film is stable, according to criterion (S), while solution B, which corresponds to the thicker film, is unstable. Thus for the case of countercurrent flow one may have either one stable solution or no solution. For cocurrent flow two situations are presented for different chosen values of A_ When the Wallis (1969) correlation, X = O.OOS(l + 3006/D), is used (Fig. 2) a single stable solution is obtained for any gas and liquid flow rates (point C or D). When a constant value forfi is used,fi = 0.05 (Fig. 3) multiple solutions may occur: two stable solutions (points A and C) and an unstable solution (point B). Note that in the case of upward cocurrent annular flow the solid lines on the zi vs 6 curves display a minimum. The steady-state solutions along the branch to the left of the minimum are always stable, independent of the correlation for X (the slope of the curve is negative). The solution along the branch to the right of the minimum may be either linearly stable (Fig. 2, point C) or unstable (Fig. 3, point B) depending on fP Barnes and Taitel(1985) showed that the minimum in the z, vs 6 curves is associated with the change in the direction of the velocity profile in the liquid lilm. Thus, it seems that steady-state solutions to the right of the minimum where the film is quite thick may be
STABILITY
AND
DYNAMIC
SIMULATION
The transient formulation presented by eqs (5) and (6) are used in order to examine numerically the nonlinear stability of the steady-state solutions. Figure 4 is a phase plane diagram that shows the steady-state relations between uL and 6/D in countercurrent flow with the same conditions as those in Fig. 2. The dome-shaped solid line is the relation between uL and S/D that satisfies the momentum equation du,/dc = 0 [eq. (5)]. The broken line satisfies the continuity equation daJdt = 0 Eeq. (6)]. The intersections of the two lines yield two points which are the steady-state solutions, points A and B. The linear stability of these solutions were determined by the aforementioned method, and it was found that point A is linearly stable whereas point B is unstable. The nonlinear stability was studied by examining the system response to finite disturbances. Numerical runs were carried out by starting the dynamic imulation at different points. The results are shown in Fig. 5. It is clearly seen that the trajectories are attracted to the stable steady-state solution (point A) and repelled
from the unstable solution (point B). There is a wide region of attraction towards A to the “left” of B whereas most of the region to the “right” of point B is a “runway” zone. Thus, countercurrent solutions that are linearly stable are also stable for finite disturbances. It is also interesting to observe that the transient solution in the phase space approaches very quickly the dome-shaped solid line and follows it to the stable steady-state solution A (Figs 4 and 5). This means that the dynamics of the momentum equation are very fast and that one can further simplify the problem by considering a quasi-steady state for the momentum equation.
%fD Fig. 4. Phase diagram for uLs= - 0.1 m/s, twos= 10 m/s (case shown in Fig 2): (0) linearly stable steady state, (0) unstable steady state.
3370
DVORA
BARNEA
and YEHUDATAITEL
simulation for U- = - 0.1 m/s, uoS Fig. 5. Dynamic = 10 m/s (case shown in Fig. 2): (Cl) linearly stable steady
Fig. 7. Dynamic simulation for uLS = 0.1 m/s, uGS = IO m/s (case shown in Fig. 2): (Cl) linearly stable steady state.
state, (0) unstable steady state.
For the case of cocurrent flow a single steady-state solution is obtained when the Wallis correlation is used for f, (solutions C and D in Fig. 2). The two solutions were found to be linearly stable. However, finite disturbances adjacent to these steady states indicate that the solution that corresponds to the thin tilm (along the branch to the left of the minimum in Fig. 2) is indeed a stable steady state that attracts the trajectories (Fig. 6). The solution which corresponds to the thicker film is of a peculiar character (Fig. 7). It resembles the behavior of a stable focus near the vicinity of the steady solution (even though the eigenvalues of the amplification matrix are real). However, one can observe that this solution is unstable for finite disturbances. As shown (Fig. 7), a finite decrease in the film thickness leads to negative film velocities, resulting in the destruction of cocurrent flow. Thus the solutions along the branch to the right of the minimum are practically unstable. An interesting case is the one of cocurrent flow with a constant interfacial shear stress coefficient,f, = 0.05. Figure 3 shows that in this case one can obtain three steady-state solutions, A, B and C. In Fig. 8 these three steady-state solutions are given as the intersection of
IS-
-
c
I , I 1
3 ;
MOMENTUM
Ea.
--a-__
.04
06
.08
.I2
.I
.I4
.I6
WD Fig. 8. Phase diagram for uLs = 0.05 m/s, uGs = 15 m/s (case shown in Fig. 3): (Cl) linearly stable steady state, (0) unstable steady state.
J
w I.5
F 5
1.0
& 0.5
0
!E’c~ .: A
E
0
0.02
0.04
Q
0.02
0.02
0.1
0.12
0.14
OM
Fig. 9. Dynamic simulation for uLs = 0.05 m/s, ncs = 15 m/s (case shown in Fig. 3): ( q) linearly stable steady state, (0) unstable steady state.
.06
.06
.I
Fig. 6. Dynamic simulation for uU = 0.1 m/s, uoS = 25 m/s (case shown in Fig. 2): (Cl) linearly stable steady state.
the solid line duJdt = 0 and the broken line daJdt = 0. According to the linear analysis [eqs (7) and (8)3 A and C are stable solutions whereas 5 is unstable. Figure 9 shows the dynamic simulations for this case. As shown, point A is a stable steady state that attracts the trajectories whereas point B is clearly a strong repellent point. The steady state C, however, attracts only those trajeotories that are adjacent to the steady state C. Its characteristic behavior is similar to point C shown in Fig. 7, that was discussed before. Namely, in this case also, any finite disturbance will
Nonlinear stability and dynal nit simulation of annular flow cause
the solution
that
D
this is an unstable steady state that is not likely to occur. Thus, out of the three steady-state solutions, we expect only the steady-state solution that corresponds to the thinnest film thickness to exist (point A).
to
run
away.
This
suggests
r” B I
practically
SUMMARY
AND
CONCLUSIONS
For co- and countercurrent annular flow one may obtain several steady-state solutions, some of which are linearly unstable. The nonlinear stability of the steady-state solutions was examined via dynamic simulations. It was indeed found that some linearly stable solutions are unstable for finite disturbances and therefore are not expected to exist. In countercurrent flow for any given gas and liquid flow rates we may have either two steady-state solutions or no solution. For the case where two solutions exist both the linear and the nonlinear stability analysis show that the solution that corresponds to the thinner film is stable whereas the solution that corresponds to the thicker film is unstable. The numerical simulation shows that the stable solution is a strong attraction point whereas the unstable solution is a strong repellent point. For cocurrent flow it is shown that the linear stability analysis is not sufficient to determine the status of the steady-state solutions_ It is shown that solutions that are linearly stable may not be stable to finite disturbances. The steady-state solutions along the left branch of the ziL vs S curve (the solid lines in Figs 2 and 3) are always stable. The solutions along the branch to the right of the minimum are always unstable for finite disturbances. The locus of the minima is the line of neutral stability and the solutions to the right of this locus indicate transition from cocurrent annular flow. The above-mentioned conclusions regarding the stability of annular flow can be applied to detect the inception of flooding and flow reversal (Bamea and Taitel, 1985) and to detect the transition from annular to intermittent flow (Barnea, 1986). NOTATION
A
C
pipe cross-sectional area constant in the friction factor correlation
: t u
3371
pipe diameter hydraulic diameter for the liquid friction factor gravitational acceleration pipe length constant in the friction factor correlation wetted periphery time velocity in the axial direction
Greek letters ; s V
P
T
void fraction angle of inclination, film thickness kinematic viscosity density shear stress
positive
for upward flow
Subscripts am superscripts G gas IL s
interface liquid superficial
REFERENCES
Ardron. K. H., 1980, One-dimensional two-fluid equations for horizontal stratified two-phase flow. Inr. .I. Multiphase Flow 6.295-305. Bamea. D., 1986, Transition from annular flow and from dis&rsed bubble flow-unified models for the whole range of viue inclinations. Int. J. Multiphase Flow 12, 733-744.
- -
Commun.
Hear
Barnes, D. and Taitel, Y., 1985,Stability of annular flow. 1~. Mass
Transfer 12, 611-622. Y., 1989, Transient-formulafion modes and stability of steady-state annular flow. Gem. Engng Sci. 44, 325-332. Stewart, H. B. and Wendroff B., 1984, Two-phase flow: models and methods. J. cumput. Whys. 56, 363409.
Bamea, D. and Taitel,
Lyczkowski, R. W., Gidaspow, D., Solbrig, C. W. and Hughes, E. D., 1978, Characteristics and stability analyses of transient one-dimensional two-phase Bow equations and their finite difference approximations. Nucl. Sci. Engng 66, 378396. Ramshaw, J. D. and Trapp, I. A., 1978, Characteristics, stability. and short-wavelength phenomena in two-phase flow equation systems. Nucc SC;. Engng 66,93-1021 Two-phase Flow. Wallis. G. B.. 1969, One-dimensional McGraw-Hill, New York.