Inclined two-layered stratified channel flows: Long wave stability analysis of multiple solution regions

International Journal of Multiphase Flow 62 (2014) 17–29

Contents lists available at ScienceDirect

International Journal of Multiphase Flow j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / i j m u l fl o w

Inclined two-layered stratified channel flows: Long wave stability analysis of multiple solution regions R. Kushnir, V. Segal, A. Ullmann, N. Brauner ⇑ School of Mechanical Engineering, Tel Aviv University, Tel Aviv 69978, Israel

a r t i c l e

i n f o

Article history: Received 22 October 2013 Received in revised form 20 January 2014 Accepted 21 January 2014 Available online 19 February 2014 Keywords: Stratified flow Concurrent Countercurrent Stability Multiple solutions Inclined channel

a b s t r a c t In the present work, the linear stability of two-layered stratified channel flows to long wave disturbances is studied. In particular, the study addresses the stability of laminar inclined counter-current and concurrent flows in the regions of multiple solutions for the holdup and pressure drop. The analysis is carried out by solving the Orr-Sommerfeld equations for two-plate geometry, through a formal power series in the wave number. The results are summarized in the form of stability boundaries on flow rate maps, which enable a systematic study of the effect of the system physical parameters on the stratified-smooth/ wavy transition in gas–liquid and liquid–liquid systems. It is demonstrated that for counter-current flow there is a region of low flow rates where the two solutions for the holdup are stable. Likewise, the results of concurrent gas–liquid upward flows reveal a region where all three solutions are stable. Moreover, it was found that the middle solution is always stable within the entire 3-s domain. Additionally, the analysis of the wave induced stresses in the axial direction reveals that the terms in phase with the wave slope should be considered in long wave stability analyses of stratified flows. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Stratified flow is a basic flow pattern in horizontal and inclined gas–liquid and liquid–liquid systems in a gravity field, where a continuous layer of a light phase flow on top of a heavier phase. Stratified two-phase flow regime is frequently encountered in various important chemical and industrial processes. For certain operating conditions however, interfacial instabilities can arise and may produce undesired effects. Therefore, the exploration of the interface stability dependence on the operational parameters is of practical importance. Exact solutions for steady laminar stratified flow in inclined pipes are available in the literature (e.g. Ullmann et al., 2004). However, an exact formulation of transient flow in pipe geometry is too complicated to conduct a rigorous stability analysis, and the stability boundaries are commonly predicted based on the simplified 1D two-fluid models. On the other hand, exact analysis of the flow in the simpler geometry of two plates can be conveniently carried out to provide insights into the mechanisms involved in the destabilization of stratified flows. Indeed, since the classical work of Yih (1967), the stability of stratified flow in the two-plate geometry has been extensively studied in the literature (e.g. Hoo-

⇑ Corresponding author. Tel.: +972 3 640 8127; fax: +972 3 640 7334. E-mail address: [email protected] (N. Brauner). http://dx.doi.org/10.1016/j.ijmultiphaseflow.2014.01.013 0301-9322/Ó 2014 Elsevier Ltd. All rights reserved.

per and Boyd, 1983; Yiantsios and Higgins, 1988; Charru and Fabre, 1994; Kuru et al., 1995). Although, a considerable amount of studies addressed horizontal flow (where the gravity driven multiple solutions are absent), only few referred also to inclined flows (Tilley et al., 1994; Boomkamp and Miesen, 1996; Amaouche et al., 2007; Vempati et al., 2010). In those few studies however, the issue of multiple solutions was not considered. Exact solutions for steady flow with a smooth interface indicates that in the counter-current region there exist always two possible solutions, whereas in the concurrent up-flow and downflow a triple solution is obtained in a limited range of the flow parameters (Ullmann et al., 2003a,b, 2004). The introduction of multiple solution regions on flow pattern maps of various twophase systems shows the practical significance of multiple solutions, and that their boundaries may be associated with flow pattern transition. The feasibility of multiple holdups was also verified experimentally (Ullmann et al., 2003a,b). However, stability analysis is required to determine the range of parameters where the assumption of a smooth interface is valid. Such an analysis may also rule out the feasibility of part of the solutions in the multiple solution regions. The stability of the flow with respect to long wave disturbances is of particular interest since long wave is an inherent approximation of two-fluid models. In the literature of gas–liquid and liquid– liquid in channels, the instability of the stratified flow pattern is

18

R. Kushnir et al. / International Journal of Multiphase Flow 62 (2014) 17–29

usually associated with to the Kelvin–Helmholtz (K–H) mechanism, and the stability analysis is carried out based on the transient two-fluid model equations. The K–H mechanism attributes the growth of interfacial disturbances to the streams inertia forces, which give rise to wave induced pressure fluctuation in phase with the wave height. However, it has been long recognized that the K– H mechanism is not the one responsible for wind-wave generation. The literature on wind waves attributes the energy transfer from the wind to the waves to the Miles-Phillips theory (e.g. Miles, 1962) and to the sheltering mechanism (Jeffreys, 1925), where wave induced pressure fluctuation in phase with the wave slope are responsible for the wind drag and momentum transfer to the waves. According to Miles-Phillips theory, pressure fluctuation in phase with the wave slope results from wind-wave interactions within the critical layer above the water surface, where the air velocity is lower than the wave speed. Deformation of the air flow within the critical layer results in a low pressure at the leeward side and high pressure on the windward side of the wave. According to the Jeffery’s sheltering mechanism, the low pressure at the leeward side is attributed to separation of the air flow over the crest, which may occur if the air flow is faster than the wave. However, since during the initial stage of the wave growth, the wind at the wave surface is slower than the wave, the sheltering mechanism is considered to be relevant only at some mature stage of the wave growth, when the local slope becomes larger than a critical value (e.g., Banner and Melville, 1976; Kawai, 1982; Kharif et al., 2008). In the general case of two-phase flow in channels, it is not apparent whether the wave velocity is faster or slower than the interfacial velocity. Moreover, it is not obvious which of the phases is the faster one and thus dominates the interfacial interactions that lead to instability and wave growth. Nevertheless, it is possible that the interaction between the wave and the flow fields in the two layers will give rise to wave induced stresses in phase with the wave slope. Exploring this possibility in the limit of the long wave approximation is of significance for further analysis and understanding to what extent the above wind-wave generating mechanisms are of importance also in destabilizing the interface between two laminar layers and in the framework of two-fluid models. This study was undertaken with the purpose of exploring the linear stability of two-layered stratified channel flows to long wave disturbances at operational condition associate with multiple solutions. The analysis was carried out for the general case of countercurrent and concurrent flows, by solving the well-known Orr-Sommerfeld equations for the two-plate geometry, through a formal power series in the wave number. It provided a closed form solution for the eigenvalues and eigenfunctions. The results are summarized in the form of stability maps showing the stable and

unstable ranges of the various system parameters. The solution is also utilized for deriving the expressions for calculations of the wave induced tangential and normal stresses. 2. Formulation of the problem Consider the flow of two immiscible, incompressible fluids, labeled j = 1, 2, flowing in an inclined channel (0 6 b 6 p/2) as shown in Fig. 1. The flow, assumed isothermal and two dimensional, is driven by an imposed pressure gradient and a component of gravity in the ^ x direction. For the indicated coordinate system, the dimensionless continuity and momentum equations governing the flow in each phase are:

@uj @ v j þ ¼0 @x @y

ð1aÞ

@uj @uj @uj q @pj 1 þ uj þ vj ¼ 1 þ @t @x @y r qj @x Re2 þ

mj @ 2 uj @ 2 uj þ m2 @x2 @y2

!

sin b Fr2

ð1bÞ

@v j @v j @v j q @pj 1 þ uj þ vj ¼ 1 þ @t @x @y r qj @y Re2 

mj @ 2 v j @ 2 v j þ m2 @x2 @y2

!

cos b Fr2

ð1cÞ

These equations are subject to the no-slip impermeable wall boundary conditions

u2 ¼ v 2 ¼ 0

y ¼ 1; y ¼ n;

ð2aÞ

u1 ¼ v 1 ¼ 0

ð2bÞ

and to the following boundary conditions at the interface, y = g(x, t)

u1 ¼ u2 ;

vj ¼

v1 ¼ v2

ð2cÞ

@g @g þ uj @t @x

ð2dÞ

" ( )#   2 ! ml @u @ v @g @u @ g 1 þ 4 ¼0 @y @x @x @x l1 @x "

ð2eÞ

!(  2 !   )# @u @g @u @ v @ g pþ 1 þ þ @y @x @x l1 1  ð@ g=@xÞ2 @x @x  2  We1 @ g=@x2 2 ¼ 3=2 ð1 þ ð@ g=@xÞ2 Þ ml

Re1 2

ð2fÞ

q2

ρ 2, μ 2 q1

ρ 1, μ 1

ηh2

g yˆ

Eqs. (2c)–(2f) represent kinematic and dynamic conditions (see Joseph and Renardy, 1993 for more details), where [f] denotes the jump f2–f1 across the interface in any quantity f. The dimensionless variables and parameters are defined as follows:

^ j ; v^ j Þ=ui ; ðuj ; v j Þ ¼ ðu

ˆx

h2

m ¼ ll1 ; 2

r ¼ qq1 ; 2

^Þ=h2 ; ðx; yÞ ¼ ð^x; y

n ¼ hh12 ;

Re2 ¼ umi h22 ;

t ¼ ^tui =h2 ; u2

Fr2 ¼ ghi2 ;

^j =q2 u2i pj ¼ p We2 ¼

q2 h2 u2i r

ð3Þ

h1

β Fig. 1. Schematic description of two-layer flow configuration in an inclined channel.

^ j ; v^ j , and p ^j are the velocities components and pressure of where u fluid j, lj, mj and qj are the corresponding dynamic viscosity, kinematic viscosity and density, and ^t, g and r denote the time, gravitational acceleration, and interfacial surface tension, respectively. As seen, the length, velocity, time and pressure scales are h2, ui, h2/ui

19

R. Kushnir et al. / International Journal of Multiphase Flow 62 (2014) 17–29

and q2 u2i , respectively, where ui represents the basic flow interface velocity, and hj is the basic flow layer height of fluid j (see Fig. 1).

1

0.8

Y=-500

3. The basic flow

∼ h

With reference to Fig. 1, the basic flow, assumed laminar and ^ j varying only steady, is parallel to the ^ x-axis, with the velocity u ^. Under these conditions, the exact solution of the differenwith y tial system (1)-(2) is as follows:

g¼0

ð4aÞ

v j ¼ 0;

700

0.6

3.1. Velocity profiles and pressure gradient

uj  U j ¼ 1 þ aj y þ bj y2 ;

-7 100

j ¼ 1; 2

e  ðm þ nÞð1  rÞ dpi dP 2 ð1 þ nrÞ Y  ¼ nð1 þ nÞð1  rÞ dx Re2 dx

!

0.2

0.4

0.15

0.2

0.05

100 -500

0 -0.2

0 -150

-50

ð4bÞ ð4cÞ

-7

0.1

50

-0.1

150

700 0

0.1

0.2

250

X2 Fig. 2. Single, double and triple holdup solution as calculated for the base flow (m = 2).

where

e a2 m  n2 þ n Y ; a2 ¼ 2 n þn m e e mþnY m þ n þ nY ; b2 ¼  b1 ¼  2 2 ðn þ nÞm n þn nð1  rÞRe2 sin b e Y ¼ 2Fr2

a1 ¼

ð5aÞ ð5bÞ ð5cÞ

e are based on the holdup and interfacial velocity, Note that n and Y which are implicitly dependent on the specified system parameters and operational conditions. The latter are used to define the Martinelli and inclination parameters:

X2 

b ^xÞ ðd P=d 1s ¼ mq; b ^xÞ ðd P=d 2s

q¼

q1 ; q2

Y

q2 ðr  1Þg sin b b ^xÞ ðd P=d 2s

ð6Þ

b ^xÞ ¼ 12l q =H3 is Here qj is the feed flow rate of phase j and ðd P=d j j js the corresponding superficial pressure drop for single phase flow in the channel, where H = h1 + h2. Given the parameters m, Y, and X2 ~ ¼ h1 =H ¼ n=ð1 þ nÞ, is obtained (or q), the lower phase holdup, h by solving the following implicit algebraic equation (e.g., Ullmann et al., 2003a):

~ 2 ½ð1þ2hÞmþð1mÞ ~ ~ ~ ~ 2 ½ð32hÞmþð1mÞ ~ ~2  hð4 hÞ h h mqð1 hÞ 3 ~3 ð1 hÞ ~ ½hþmð1 ~ ~ hÞ 4h

~2 þ 4mð1  mÞh ~ þ m2 1  ½2ð1  2h ~þh ~2 Þ þð4  12m þ 6m2 Þh ~ h ~ 1 ð11Þ ð1 þ hÞ

¼0 ð7Þ

e is obtained by Then, Y

~  1Þ2 ð4h ~ þ m  2mh ~ þ mh ~2  h ~2 Þ Yðh e Y ¼ ~4  3Y h ~3 þ 3Y h ~2  Y h ~  1Þ ðY h

Multiple solutions for the holdup is an inherent characteristic of inclined stratified flow (e.g., Ullmann et al., 2003a): Double holdups are obtained the whole range of countercurrent flow, while in concurrent flow a multiple (triple) holdup is predicted only in a limited range of operational conditions. For example, Fig. 2 shows the holdup vs. X2 for different values of the inclination parameter Y. Two solutions for the holdup are obtained in countercurrent flow, where the light phase is flowing upward (X2 < 0, Y < 0), triple solution for the upward co-current flow (X2 > 0, Y < 0) are typically obtained for small positive values of X2(q2 > > q1), and triple solution for concurrent down-flow (X2 > 0,Y > 0) are typically obtained at relatively high X2(q2 < < q1). Evidently, there is no solution for the countercurrent cases where the light phase is flowing downward (X2 < 0, Y > 0). For a specified value of Y, the range of X2 corresponding to multiple solutions is bounded by the values of X2 corresponding 2 ~ ¼ 0. Using Eq. (7) for calculation of this derivative to dX =dh yields:

~ 4  2h ~3 Þ þ 2mðm  1Þh ~  m2  Y ¼ ½ð1  2m þ m2 Þðh ~ 4 þ ð6 þ 10m  4m2 Þh ~3  ½ð1  2m þ m2 Þh

Y 

3.2. Multi-holdup regions

ð8Þ

The corresponding pressure gradient and interfacial velocity can be calculated by:

b e  d P=d^x  q2 g sin b P b ^xÞ ðd P=d

Eqs. (7) and (11) can be solved simultaneously to yield the X2 vs. Y relationship that forms the boundaries of multiple solution regions ~ Eq. (11) yields Y, and these are for a given system. Given m and h, then substituted into Eq. (7) to obtain the corresponding X2. The so-obtained regions of multiple solutions in both co-current and countercurrent flow configurations are depicted in Fig. 3. 4. Linear stability 4.1. The differential system governing stability

2s

¼

~ 2  4mhð1 ~  hÞ ~ h ~2 3mqð1  hÞ 2 ~  hÞ ~ ½ð1 þ 2hÞm ~ þ ð1  mÞhð4 ~  hÞ ~  3h ~ 4hð1

~  hÞðY ~ h ~  PÞ e ui 6hð1 ¼ ~ ~ U 2s mð1  hÞ þ h where Ujs = qj/H is the superficial velocity of fluid j.

ð9Þ

ð10Þ

A linear stability of the basic flow solution to infinitesimal, twodimensional disturbances was carried out in the routine manner. The basic velocities and pressure fields are perturbed by infinitesimal disturbances, whereby uj ¼ U j þ u0j ; v j ¼ v 0j ; pj ¼ Pj þ p0j , and g = g0 represents the dimensionless disturbance of the interface. The disturbance velocities are represented by the corresponding stream function, u0j ¼ @ Uj =@y; v 0j ¼ @ Uj =@x. All perturbation quantities are assumed to have exponential growth, hence

20

R. Kushnir et al. / International Journal of Multiphase Flow 62 (2014) 17–29

1000

m =2 countercurrent (light phase down) 500

concurrent down 1-S

Y

0-S

3-S

0

countercurrent 0-S

concurrent up

0-S

-7.5

1-S

2-S

-5

2-S

3-S

-10 -0.02 -0.01

-500 -100

-50

0

50

X

100

0

1-S

0.01 0.02

150

200

2

Fig. 3. Multiple-solution regions for concurrent and counter-current flows (m = 2).

0

1 0 1 Uj /j ðyÞ B 0C B C @ pj A ¼ @ Pj ðyÞ AeikðxctÞ H g0

ð12Þ

where /j, Pj, and H are the amplitudes of the perturbations. Upon substitution of these disturbances in the linearized momentum Eqs. (1b) and (1c) and elimination of Pj, the well-known Orr-Sommerfeld equations are obtained:

/ijv  2k /00j þ k /j ¼ ikRe2   i m2 h 2  ðU j  cÞ /00j  k /j  /j U 00j j ¼ 1; 2 2

4

mj

ð13Þ

where the superscripts in /j and Uj denote differentiation with respect to y. In the framework of linear analysis, a Taylor expansion in g around y = 0 is used to formulate the corresponding boundary conditions at y = 0 (see, Segal, 2008, for more details). The resulting set of boundary conditions read:

/2 ¼ /02 ¼ 0; at y ¼ 1

ð14aÞ

/1 ¼ /01 ¼ 0; at y ¼ n

ð14bÞ

/1 ¼ /2 ; at y ¼ 0

ð14cÞ

/01 þ

/1 / U 0 ¼ /02 þ 2 U 02 ; at y ¼ 0 c1 1 c1 2

/002 þ k /2 þ U 002

  /2 / 2 ¼ m /001 þ k /1 þ U 001 1 ; at y ¼ 0 c1 c1

    2 0 2 0 000 m /000 1  3k /1  /2  3k /2      þ ikRe2 r ðc  1Þ/01 þ U 01 /1  ðc  1Þ/02 þ U 02 /2 h i / 2 1 2 ; at y ¼ 0 ¼ ikRe2 ðr  1ÞFr1 2 cos b þ k We2 c1

ð14dÞ

ð14eÞ

then the phase speed of the perturbation, while kcI is the growth rate. Given m, n, r, b, Fr2, Re2, and We2, the wave speed c has to take on certain values in order that /j will not be identically zero. These values are the eigenvalues of the problem, and the corresponding /j solutions are the eigenfunctions. The stability of the basic flow is determined by the sign of cI. Neutral stability corresponds to cI = 0. Note however, that positive (negative) value of cI does not necessarily indicate instability (stability), since the velocity is scaled by ui, which may attain negative values (in co-current up-flows and in countercurrent flows). Hence, when ui > 0, instability is indicated by cI > 0, and vice versa (see also Eq. 27 below). Note that in this study a two-dimensional disturbance is considered. Yet, it was shown that Squire’s theorem holds at least for horizontal stratified flows (e.g., Hesla et al., 1986). Accordingly, the flow is considered to be more stable with respect to 3-three-dimensional perturbations, and therefore it is sufficient to consider two-dimensional perturbations in order to determine the stability boundary. 4.2. Long wave analysis In order to find the eigenvalues of this differential system for an arbitrary wave number, a numerical method must be adopted. However, as previously shown (e.g. Yih, 1967) interfacial instability due to viscosity and density stratification may be conveniently studied via long-wave analysis. 4.2.1. Asymptotic solution In the long-wave asymptotic limit, the eigenfunctions and the eigenvalue are represented as a regular perturbation series in k: 2

3

/j ¼ /j;0 þ k/j;1 þ k /j;2 þ k /j;3 þ    2

j ¼ 1; 2

ð15Þ

3

c ¼ c0 þ kc1 þ k c2 þ k c3 þ   

ð16Þ

When these expansions are substituted into Eqs. (13) and (14) and terms of like order in k are grouped together, a simplified set of fourth-order equations and boundary conditions for /j,M are obtained (M is the order of the solution). These equations may be solved sequentially to determine /j,M, and an eigenvalue relation for cM at each order of the approximation. For long-wavelength (k  1), it suffices to consider only the first two terms in Eqs. (15) and (16). In the zero order approximation all terms containing k in the differential system are ignored, whereby /j = /j,0 (j = 1, 2), and c = c0. For this case (assuming k Re2  1) Eq. (13) becomes: iv /iv 1;0 ¼ 0 and /2;0 ¼ 0

ð17Þ

The boundary conditions (14) reduce to

/2;0 ¼ /02;0 ¼ 0; at y ¼ 1 /1;0 ¼ ð14fÞ

Note that the same set of equations and boundary conditions are obtained by using a variable transformation, which can be repre_ _ _ _ _ sented by: x ¼ x ; y ¼ f ðy ; gðx ; t ÞÞ; t ¼ t , where the old variables are at the left hand side, and the new variables on the right. In this _ new coordinates the interface location is invariant (namely y ¼ 0 _ _ for y ¼ gðx ; t Þ), and there is no need to carry a Taylor expansion (Segal, 2008). The differential system governing the stability problem consists of (13) and (14). The temporal growth of the disturbances is studied by taking k as a dimensionless real wave number (k = 2ph2/k, with k being the wave length), and allowing the dimensionless wave speed c to be complex, thereby c = cR + icI. The quantity cR is

/01;0

ð18aÞ

¼ 0; at y ¼ n

ð18bÞ

/1;0 ¼ /2;0 ; at y ¼ 0 /1;0 0 /2;0 0 /01;0 þ U ¼ /02;0 þ U ; at y ¼ 0 c0  1 1 c0  1 2   /2;0 /1;0 /002;0 þ U 002 ; at y ¼ 0 ¼ m /001;0 þ U 001 c0  1 c0  1

ð18cÞ ð18dÞ ð18eÞ

000 m/000 1;0  /2;0 ¼ 0; at y ¼ 0

ð18fÞ

The solution of Eq. (17) are third-order polynomials:

/j;0 ¼ c0  1 þ B1j;0 y þ B2j;0 y2 þ B3j;0 y3 ;

j ¼ 1; 2

ð19Þ 

B01;0 ; B02;0



were Without loss of generality the arbitrary constants set to c0  1, thus (18c) is satisfied. It worth emphasizing that in previous published analyses (e.g. Yih, 1967; Yiantsios and Higgins, 1988) B01;0 and B02;0 were set to 1. However, that choice resulted in

21

R. Kushnir et al. / International Journal of Multiphase Flow 62 (2014) 17–29

an apparent singularity of the solution under conditions where c0 = 1 (i.e., when the wave velocity equals the basic flow interfacial velocity, ui), which arises for vanishing c0  1 in (18d) and (18e). Such conditions were referred to in the literature as critical points pffiffiffi (Yih, 1967). In horizontal flows, c0 = 1 when m = 1 or m ¼ n, in which case c1 = 0 (Eq. (25) below), and therefore these points correspond to neutral stability. In inclined flows, however, the value of c1 at the critical points is not necessarily zero, and c0 = 1 introduces difficulties in the computation of c1. This is resolved by taking B01;0 ¼ B02;0 ¼ c0  1 (instead of 1). The solution obtained for the  other six constants Bl1;0 ; Bl2;0 ; l ¼ 1; 2; 3 that satisfy the boundary conditions is given in Appendix A, and the solution obtained for c0 is given by:

e þ 2nðm  n2 Þðm  1Þ ðn  2n þ 2mn  mÞn Y c0 ¼ 1 þ n4 þ 4n3 m þ 6n2 m þ 4nm þ m2 2

ð20Þ

Hence, for the zero order approximation, c0 is real. In fact, it represents the velocity of a kinematic wave, which can also be obtained directly from the basic flow solution Eqs. (7)–(10):

  1 dU 1s c0 ¼ ~ U ui dh

1s þU 2s ¼const

~ @F=@ h U 2s ¼ ðq þ 1Þ@F=@q þ @F=@Y ui

ð21Þ

where the function F represents the left hand side of Eq. (7). Note e ¼ 0), the result reduces to that obtained by that, for b = 0 (i.e. Y Yiantsios and Higgins (1988) for horizontal Poiseuille flows. The first order approximation is obtained when all terms containing k2 and higher order are ignored whereby /j = /j,0 + k/j,1 (j = 1, 2), and c = c0 + kc1. The substitution of these expansions into (13) and (14) yields the following differential system: v /i1;1

1

h

c0 Þ/001;0

/1;0 U 001

¼ iRe2 m r ðU 1   h i v ¼ iRe2 ðU 2  c0 Þ/002;0  /2;0 U 002 /i2;1

i

Thus, in the first-order solution the wave speed is complex, with cR = c0 and cI =  ikc1. Therefore, the stability of the flow with respect to long wave disturbances can be determined from (25). Finally, the kinematic boundary condition is used to obtain the corresponding solution for the wave amplitude, H:

H¼

/1 ð0Þ /2 ð0Þ ¼ c1 c1

The validity of this asymptotic solution has been confirmed through the comparison with the numerical results of Tilley et al. (1994), as seen in Table 1. 4.2.2. Wave-induced stresses The obtained solution can further be used to evaluate the wave induced velocity components, shear stresses at the interface and at the channel walls, and pressure difference across the interface. The wave induced velocity components and the wave profile are represented by the real part of the solution:

  u0j ¼ /0jR cos h  /0jI sin h eðkcI tÞ

j ¼ 1; 2

v

j ¼ 1; 2

0 j

¼ kð/jI cos h þ /jR sin hÞe

g ¼ ðHR cos h  HI sin hÞe 0

/2;1 ¼

/jR ¼ /j;0 ; HR ¼

s ¼

/1;1 ¼ /01;1 ¼ 0; at y ¼ n

ð23bÞ

/1;1 ¼ /2;1 ; at y ¼ 0

ð23cÞ

c1 /01;0 þ ðc0  1Þ/01;1 þ /1;1 U 01 ¼ c1 /02;0 þ ðc0  1Þ/02;1 þ /2;1 U 02 ; at y ¼ 0 ð23dÞ c1 /002;0 þ ðc0  1Þ/002;1

  ¼ m c1 /001;0 þ ðc0  1Þ/001;1 ; at y ¼ 0

  0 000 0 m/000 1;1  /2;1 þ iRe2 ðr  1Þ ðc 0  1Þ/2;0 þ U 2 /2;0 ! 2 /2;0 ðr  1Þ cos b k þ ; at y ¼ 0 ¼ iRe2 Fr2 We2 ðc0  1Þ

ð23eÞ

/j;1

  ¼ iRe2 B1j;1 y þ B2j;1 y2 þ B3j;1 y3 þ wj ðyÞ ;

j ¼ 1; 2

The coefficients Blj;1 ðl ¼ 1; 2; 3Þ and the particular solutions wj are given in Appendix A. The value of c1 is then calculated from

c1 ¼

  iRe2 B12;1  B11;1 ðc0  1Þmðn þ 1Þn eÞ ðm  1Þðm  n2 þ n Y

¼

ð25Þ

ð27bÞ

0

h2

lj ui  h2

0

!

0

ð/00jR

HR U 00j Þ cos h

þ

y¼0



   /00jI þ HI U 00j sin h

y¼0

eðkcI tÞ

ð28Þ

The first (cos h) term on the right-hand side (RHS) of Eq. (28) is the interfacial shear stress in phase with the wave elevation, g  ðHR cos hÞekcI t , while the second term (sin h), is in phase with the wave slope, @ g=@x  ðkHR sin hÞekcI t , and can be represented b i @ g=@x, where by K



bi ¼  K

lj ui /00jI þ HI U 00j  h2

 

kHR

ð29Þ y¼0

b i can be calculated based on the solution obtained for The value of K either of the phases (j = 1, 2):

ð23fÞ

ð24Þ

HI ¼

;

@ 2 U @u @ v g 2j þ j þ j @y @y @x

lj ui

bi ¼  K

Notice that, the term containing the surface tension was retained (in spite that it contains a high order wave number) to account for surface tension effects in case We2 = O(k2). The solutions of (22) are the following polynomials (the arbitrary constants B0j;1 , j = 1, 2 were set to zero):

h ¼ kðx  c0 tÞ

ðc0 1ÞcI ðc0 1Þ2 þc2I

The dimensional tangential wave induced shear stress at the interface (only linear terms in k are retained) is:

^0i

ð23aÞ

ð27aÞ

ðkcI tÞ

/jI ¼ ik/j;1 ;

ðc0 1Þ2 ðc0 1Þ2 þc2I

ð22Þ

¼ 0; at y ¼ 1

ðkcI tÞ

where

with the boundary conditions

/02;1

ð26Þ

  2lj ui ic1 bj Re2 B2j;1 þ h2 c0  1

ð30Þ

Similarly, the wave induced dimensional wall shear stresses can be calculated by

s^0w1 ¼ 



l1 ui @u01 h2

@y

þ

@ v 01 @x

 ; y¼n

s^0w2 ¼



l2 ui @u02 h2

@y

þ

@ v 02 @x

 ð31Þ y¼1

Using Eq. (27) the corresponding terms of the wave induced (lower and upper) wall shear stresses in phase with the wave slope, b 1 @ g=@x; K b 2 @ g=@x are identified. The resulting coefficients K b 1 and K b K 2 are: 00  b 1 ¼  l1 ui /1I  K ; h2 kHR y¼n

where

00  b 2 ¼ l2 ui /2I  K h2 kHR y¼1

ð32Þ

22

R. Kushnir et al. / International Journal of Multiphase Flow 62 (2014) 17–29

Table 1 Calculated eigenvalues: comparison between the m ¼ 100; r ¼ 2; We2 ¼ ð1 þ nÞRe22 =25000; k ¼ 0:001=ð1 þ nÞ.

a

present

asymptotic

solution

and

Tilley

et

al.

(1994)

numerical

results,

for

n

b

Re2

Fr2

Asymptotic (cR, cI)a

Numeric (cR, cI)a

1/2 1/2 1/2 1/2 3 3 3 3

0 0 1 1 0 0 1 1

4.3269 4.3269 5.0742 11.8 26.2 26.2 26.69 31.098

0.025275 0.002527 0.034759 0.018797 17.5727 1.75727 18.236 2.47573

(0.25306345, 0.00012602292) (0.25306345, 0.00007233734) (0.29693843, 0.00012364063) (0.69181325, 0.00005630773) (6.218727, 0.001810625) (6.218727, 0.001573444) (6.118395, 0.0017395574) (5.215411, 0.000966384)

(0.25306349, 0.00012602286) (0.25306349, 0.00007233730) (0.29693846, 0.00012364057) (0.69181321, 0.00005630771) (6.218724, 0.001810619) (6.218724, 0.001573439) (6.118392, 0.0017395518) (5.215409, 0.000966382)

(cR, cI) are scaled by the characteristic velocity l1/(q1H) (instead of ui).

 00  ! !1 1 w1 y¼n ic1 b1 1  B21;1 2 B21;1 ðc0  1ÞRe2 B21;1 ! !1 b2 B32;1 1 ðw002 Þy¼1 ic1 b2 K 1 ¼ 1þ3 2 þ bi B2;1 2 B22;1 ðc0  1ÞRe2 B22;1 K b1 K ¼ bi K

1  3n

B31;1

þ

ð33aÞ

ð33bÞ

b i; K b 1; K b 2 are independent on k indicates that the The fact that K components of the wave induced wall and interfacial shear stresses which are in phase with the wave slope should not be ignored in the framework of long wave stability analysis, where terms of order k are retained. Wave induced axial drag in phase with the wave slope can also result from the axial gradient of the wave induced pressure differ  ence across the interface, q2 u2i =h2 @ðp1  p2 Þ=@xjy¼g . The term in b p @ g=@x, where phase with the wave slope, can be represented by K 2

b p ¼ q2 ui K h2

2

ðr  1Þ cos b k þ HR Fr2 HR We2

! ð34Þ

Inspection of Eq. (34) shows that in the limit of long wave approximation, the wave induced term in phase with the wave slope includes the gravity and surface tension stabilizing terms (in case of r P 1 and r – 0). These terms are commonly considered also in the framework of TF models. The above findings regarding the significance of the tangential and normal wave induced shear stresses in the framework of long wave stability analysis of laminar stratified flows are of importance for the assessment of the accuracy of results obtained via stability analysis of simplified two-fluid models.

(a) 1000

m

10 1

where the function f(m, n, r) is linear with respect to r(see Appendix A). Eq. (35) is valid also in the case of r = 1, except that f(m, n, r) is replaced with f(m, n,1). Thus, in addition to m = 1, the neutral stabil~cr ¼ ncr =ð1 þ ncr Þ ity boundary corresponds also to a critical holdup, h (Yiantsios and Higgins, 1988) and a critical flow rate ratio, qcr, which are determined solely by the viscosity ratio and are given by:

~cr ¼ h

pffiffiffiffiffi m pffiffiffiffiffi ; 1þ m

qcr ¼

pffiffiffiffiffi mþm pffiffiffiffiffi 1þ m

ð36Þ

The neutral stability curves in the (q, m) plane under zero gravity conditions are shown in Fig. 4a, where stable regions are denoted by ’’S’’, and unstable by ‘‘U’’. As seen, the stable region for m > 1 (higher viscosity of the lower phase) corresponds to q > qcr ~cr < h ~ < 1Þ, while for m < 1 the stable region is for ðh ~ mmin = 1.8593 corresponds to r ? 0, and the stable region diminishes as r increases. The maximal value of r (<1) for which density stratification can stabilize the interface (the Sr regions in

Sr

u

0.5

S 0.4

r=0.1 r→0

u

r=10

S

0.3 0.2

Sr

0.01 0.001 0.01

ð35Þ

(b) 0.6

r →∞ 0.1

cI c1 ¼ ¼ ðm  1Þðm  n2 Þf ðm; n; rÞ kRe2 iRe2

rmax

100

4.2.3. Stability under zero gravity conditions Zero gravity condition is a limiting case of particular interest. It corresponds to either g = 0, or to r = 1. It follows from Eq. (25) that in the limit of g ? 0(Fr2 ? 1):

0.1 0

0.1

1

q

10

100

1

10

100

1000

m

Fig. 4. Stability under zero gravity conditions. (a) long-wave neutral stability curves in the (q, m) plane; (b) maximal r for which a stable region (Sr) can be obtained for a specified m > 1.

23

R. Kushnir et al. / International Journal of Multiphase Flow 62 (2014) 17–29

Fig. 4a) is dependent on the viscosity ratio, and is shown in Fig. 4b. Likewise, the dashed curve of m < mmax = 1/1.8593 in Fig. 4a corresponds to r ? 1, and the stable region diminishes as r decreases. The minimal value of r(>1) in this case can be deduced from Fig. 4b, whereby rmin(1/m) = 1/rmax(m). 5. Results and discussion There are several mechanisms by which a flat interface solution can become unstable. In addition to the shear flow instability, which is encountered also in single-phase Poiseuille flow when the velocity exceeds a critical value (and leads to transition to turbulent flow), in the case of two-phase flow, instability of the interface can evolve due to the viscosity and/or density stratification. Viscosity difference between the fluids (m – 1) results in a jump in the velocity gradient U 0j across the interface, while density difference in inclined flows (r > 1 and b – 0) results in a jump in the curvature U 00j of the primary flow velocity profiles across the interface, such that mU 001 – U 002 (see Eq. (5b)). It has been shown that both jumps lead to energy transfer from the main flow to the disturbed flow, and are referred to in the literature as the ‘viscosity induced’ and ‘gravity induced’ instabilities, respectively (e.g., Hooper and Boyd, 1983; Boomkamp and Miesen, 1996). These are the dominant mechanisms for the growth of long wave disturbances (Segal, 2008). Note that the shear mode instability leading to transition to turbulence in either of the phases is associated with short waves and is not captured by the long wave stability analysis. However, the results of the long wave analysis provide insight to the stability of smooth stratified flow in a variety of two-phase systems, which include gas–liquid and liquid–liquid horizontal and inclined flows under normal gravity and in zero gravity conditions. 5.1. General considerations Only three dimensionless parameters (Y, m, q) are needed to determine the basic flow characteristics. However, neglecting the stabilizing term of surface tension, six dimensionless parameters are needed to determine the stability (e.g. Y, m, X2, r, b, Re2s or Fr2s, where Rejs = qjUjsH/lj and Frjs ¼ U 2js =ðgHÞ are the superficial Reynolds and Froude numbers). Conveniently, when the flow rate of the lighter phase is constant, the parameters used are Y and Re2s, while for constant flow rate of the heavier phase these two are replaced by Y/X2 and Re1s. The effect of varying light phase flow rate is demonstrated in Fig. 5. The conditions are of a constant and positive heavy phase

feed flow rate (i.e. downward flow, U1s > 0) and a variable light phase flow rate. For a given set of m, r, b, Re1S and Y/X2, the holdup can be plotted as a function of 1/X2 (/q2). The figure discloses the entire holdup curve, calculated from Eq. (7) for m = 2 and Y/ X2 = 104.4. As seen in the figure, the plot unveils the number of possible holdup solutions. The corresponding holdup vs. (cI/k)(ui/ U1s) is also plotted to determine the stability of the solutions (values indicated on the upper abscissa). Note that (cI/k)(ui/U1s) > 0 corresponds to growing disturbances. ~ ¼ 1) at It is seen that there are two points (in addition to h which the net flow rate of the light phase is zero. The correspond~ 01 and h ~ 02 . The holdup ing holdups of which will be referred to as h ~ 01 < h ~ 0), corresponds to ~
(c I /k )(u i /U2s ) -5000 1

(c I /k )(u i /U1s ) -300 1.2

-200

-100

0

100

200

1

1000

3000

5000

0.8

Y /X =104.4 β =1.8° r =1.02 m =2 Re1S=14.4

0.8

-1000

300

2

∼ h 02

-3000

r =2 β =0.5°

0.6

r =1.01 r =1.8 r =1.5 r =1.2 β =66.5° β =0.6° β =1.0° β =2.6°

∼

h

0.4

∼ h 0.6 0.4

h∼ f

0.2

∼ h 01

0.2

0 -150

0 -30

Y =500 m =2 ReS2=15

-20

-10

0

1/X

10

20

30

2

Fig. 5. The effect of varying light phase flow rate on both the holdup and cI.

-100

-50

0

X

50

100

150

2

Fig. 6. The holdup vs. the heavy phase flow rate (lower scale) for constant Y and m, and the effect of varying r and b on the corresponding normalized perturbations growth rate (upper scale).

24

R. Kushnir et al. / International Journal of Multiphase Flow 62 (2014) 17–29 1.2 (a)

(b)1.2 ∼ h02

1

β=1o r=1.1 Re1s=14.4

0.8

∼

h

β=1o

Y/X =268.9

s

u

∼

h ∼ hf

0.4

0.2

∼ h01

0.2

u

Fr1s=7×10-6

0.6

0.4

0 0.01

m=1 Re1s=14.4

0.8

2

0.6

∼ h02

1

s

u

∼ hf ∼ h 01

0 0.1

1

m

10

100

1

3

5

r

7

9

~ plane; (b) in the ðr; hÞ ~ plane. Fig. 7. Long-wave stability diagrams for slightly inclined channel. (a) in the ðm; hÞ

tion may cause instability mainly when the heavy phase is the viscous one (m > 1). Note however, that instability can also be caused by density stratification. To isolate the density effects, neutral stability curves were calculated and plotted in the (r, h) plane for m = 1, as seen in Fig. 7b. Evidently, density stratification alone can cause instability, even at small inclinations. It is worth noting that when the viscosity of the heavy phase increases (m > 1), the unstable regions of Fig. 7b increase. However, when the heavy phase is less viscous (m < 1), the unstable regions of Fig. 7b shrink. The above stability diagrams enable the demonstration of the effect of the various dimensionless system parameters on the stability of inclined stratified flows. However, from the practical point of view, given a particular two-phase system, it is of interest to use the results of the stability analysis in order to identify the region where stable stratified flow can be expected on a conventional flow pattern map for the operational window U2svs. U1s. To elucidate the effect of the inclination, the stability boundaries are first depicted on flow pattern maps of horizontal flows. 5.2. Stability maps for horizontal flows Stability of gas–liquid horizontal flow has been extensively studied in the literature. Here, it is of interest to demonstrate the stabilizing effect of gravity, by showing the stability boundary on flow pattern map, in comparison to the zero gravity stability boundary corresponding to q = qcr. Fig. 8(a)–(c) shows the stability boundary for gas–liquid systems (r 1) at different density ratios, viscosity ratios and channel heights, respectively. In the limit of r = 1, the stability boundary would correspond to the line of q=qcr. As the viscosity ratio considered is larger than 1, the entire region below the qcr boundary would be unstable for r = 1. It is worth emphasizing that since the basic flow solution for horizontal flow is independent of the phases’ densities, the flow rates corresponding to the conditions of q = qcr, are independent of the value of r. As q = qcr also correspond to equal phases’ average velocities, for any r, the region below q = qcr corresponds to average gas velocity higher than the average liquid velocity, whereas in the region above q = qcr the average liquid velocity is higher than the average gas velocity. The stabilizing effect of gravity for r > 1 is clearly shown in Fig. 8a by the extension of the stable region to q’s that are lower than qcr, hence to higher gas velocities. The larger is r, the wider is the stable region for q < qcr. In this region (low liquid flow rates), the critical superficial gas velocity is scaled approximately with [(r  1)Hg]1/2 (i.e., corresponds to a constant gas Froude number, U2s[(r  1)Hg]1/2). It is worth noting, however, that the critical superficial gas velocity at low liquid flow rates is not constant,

but decreases with further decreasing the liquid flow rate, below the values shown in Fig. 8(a)-(c). On the other hand, in the region of high liquid and gas flow rates, where the effect of gravity diminishes as compared to viscous effects, the stability boundary for r > 1 approaches the zero gravity stability boundary. Hence the stability boundary in the region of high liquid flow rates is not scaled with the liquid Froude number. The prediction of stable stratified flow in the region of low gas flow rates, independently of the liquid flow rate, implies that the transition from stratified flow to elongated-bubble (or slug) flow is not associated with the long-wave stability boundary. The viscosity ratio effect on the stability is demonstrated in Fig. 8b. It is shown that the critical gas flow rate at low liquid flow rates is mildly affected by the viscosity ratio. The viscosity effect is significant only at high liquid flow rates, where the stability boundary approaches the zero-gravity limit. Obviously, in this region, an increased m is associated with an increase of qcr, hence the stable region diminishes and is restricted to lower critical gas flow rates. Likewise, the stable region decreases with smaller channel heights, as can be seen from Fig. 8c. This is consistent with the observation made with respect to Fig. 8a, that at low liquid flow rates the critical superficial gas velocity is scaled approximately with [(r  1)Hg]1/2. Hence, smaller channel heights are associated with lower critical gas velocities. Similar behavior is observed in Fig. 9a, showing the stability boundary for horizontal liquid–liquid flow in case the heavier phase is the more viscous one (m = 2). As the density ratio is typically much lower than in gas–liquid systems, the stabilizing effect of gravity is milder, and the expansion of the stable region to q < qcr is limited to lower light phase flow rates. However, in many practical cases of oil–water systems, the lighter phase is the more viscous (m < 1). In such systems, the stable region for zero gravity (or r = 1) conditions is opposite to that obtained for m > 1, and is in the region of q < qcr (see Fig. 9b). With increasing r, the stable region extends to higher flow rates (and average velocities) of the heavy phase. In all horizontal systems, the stability boundary approaches the zero-gravity limit (equal average velocities) at high flow rates. The significance of the wave induced interfacial shear compob i @ g=@x, on the stability boundnent in phase with the wave slope, K ary, can be evaluated by referring to the study of Brauner and Moalem Maron, 1993. In that study, such a term was included in the stability analysis of the two-fluid (TF) model, and was represented by Chq(U1  U2)2@ g/@x (q is the density of the faster phase, Uj is the basic flow average velocity of fluid j, and Ch is an empirical b i; Ch ¼ K b i =ðqðU 1  U 2 Þ2 Þ. The parameter), hence in terms of K inclusion of that term in the TF model for turbulent gas–liquid flows was shown in that study to give rise to a destabilizing term,

25

R. Kushnir et al. / International Journal of Multiphase Flow 62 (2014) 17–29

(a) 10

(a)

β =0

10

m=55 H=0.02 m 1

10 -1

=7 q cr

.4

2 16

10-2

10-1

10

r=100 10-3 -3 10

r=1.25 m=2 H= 0.02 m

u

S

U1S [m/s]

U1S [m/s]

1

β =0

10

-2

10

-1

1

-2

=1 r

.41

42

qc

1000 10-3 10-3

10

10

-2

U2S [m/s]

(b) 10

(b)

10

r=1000 H=0.02 m 1

1

-1

1

10

1

10

β =0 r=1.25 m=0.5 H= 0.02 m

u

u

S

55 110 m=550

62 4.1 .488 2 7 = q cr r=10 .416 q c =7 q cr

U1S [m/s]

U1S [m/s]

10

U2S [m/s]

β =0

10 -1

u

S

10-2

10-1

10

S

-2

=0

.70

71

q cr 10-3 -3 10

10

-2

10

-1

1

10

10-3 10-3

U2S [m/s]

(c) 10

U1S [m/s]

1

Jh ¼ C h

u

0.05 m

S

0.02 m

10 -1

16 7.4

H=0.005 m

2

10-2

10

-3

10

-1

Fig. 9. Long-wave neutral stability curves for horizontal liquid–liquid system: (a) the heavier phase is the more viscous one (m = 2); (b) the light phase is the more viscous one (m = 0.5).

r=1000 m=55

10-3

-2

U2S [m/s]

β =0

= q cr

10

10

-2

10

-1

1

10

U2S [m/s] Fig. 8. Long-wave neutral stability curves for horizontal gas–liquid system at different: (a) density ratios; (b) viscosity ratios; (c) channel heights. (The dashed line represents the critical flow rate ratio of zero gravity conditions, Eq. (36)).

which was related to Miles-Phillips and Jeffreys mechanisms that are responsible for wind generated waves. Its impact on the stability boundary was represented by Jh (its ratio with respect to the gravity stabilizing term), where in the two-plate geometry is given by:

bi q ðU 1  U 2 Þ2 1 K 1 ¼ ~ ~ ~ ~  hÞ Dq Hg cos b hð1  hÞ DqHg cos b hð1

ð37Þ

It was shown in Brauner and Moalem Maron (1993) that a value of Jh of the order of 1 indicates that the interface stability is to a large extent affected by the gas–liquid interactions that give rise to an interfacial shear component in phase with the wave slope, and not only by the inertia of the two phases that are responsible for the K–H mechanism. Fig. 10 shows the value of Jh obtained by Eq. (37) along the neutral stability boundary for two typical cases corresponding to gas–liquid and liquid–liquid stratified flows (studied in Figs. 8 and 9). The results indicate that the destabilizing mechab i @ g=@x must be accounted for when analyzing the stanism due to K bility of laminar stratified flows even in the long wave approximation, in particular in the case of thin viscous layers. Obviously, the impact of the wave-induced wall-shear components in b 1 and K b 2 ) should also be phase with the wave slope (due to K considered. 5.3. Stability maps for inclined flows Long-wave stability analysis has been conducted for inclined gas–liquid and liquid–liquid flows, in particular for operational

26

R. Kushnir et al. / International Journal of Multiphase Flow 62 (2014) 17–29

1

(a)

β =0

1

H=0.02 m 10

0.8

H=0.0144 m

β =10°

-1

r =1.0335 m =1.52

10-2

U1S [m/s]

0.6

Jh 0.4

HPD

u

SHPD

f looding

10-3 10-4

0.2

r=1.25, m=2

LPD

S

10-5

S LPD

r=1000, m=55 0 10

-3

10

-2

10

-1

1

10-6

10

10-6

U1S [m/s] Fig. 10. Values of Jh along the neutral stability boundary as a function of the heavy phase superficial velocity for horizontal gas–liquid (r = 1000, m = 55), and liquid– liquid (r = 1.25, m = 2) flows.

(b)

5.3.2. Concurrent flows The interpretation of the results of the stability analysis in the region of triple solution (3-s) in concurrent down-flow and up-flow to obtain the stable regions on flow maps is rather complex. For example, Fig. 12 shows the long-wave stability regions for slightly upward inclined air–water flow. Compared to horizontal flow (Fig. 8a), the stable region for low gas and liquid flow rates is drastically reduced, whereby the flow becomes unstable for much lower gas flow rates. In this region of low gas flow rates, there is only a single solution of high holdup. At high liquid flow rates, where gravity effect becomes negligible, the stability boundary approaches the critical flow rate ratio line of zero gravity, q = qcr. With increasing the air flow rate (at relatively low liquid flow rates), a stable region re-appear, in the triple solution region, or in its vicinity. In the 3-s region (confined by the red boundaries), an additional two solutions for the holdup are obtained, of a low holdup and a medium holdup, and the stability boundaries of all

10

-4

10

-3

10

-2

10

-1

1

1

H=0.0144 m

β =26°

-1

r =1.0335 m =1.52

U1S [m/s]

10-2 HPD

10

u

-3

f looding

SHPD 10-4

10-5

S

SLPD

LPD

10-6 10-6

10

-5

10

-4

10

-3

10

-2

10

-1

1

-U2S [m/s] Fig. 11. Long-wave stability maps for countercurrent liquid–liquid flow (q2 = 916.6 kg/m3, l2 = 2.4  103 Pa s). HPD- high holdup of the heavy phase mode, LPD- high holdup of the light phase mode. (a) b = 10°; (b) b = 26°.

10

β =0.1° r=1000 m=55 H= 0.02m

1

-U1S [m/s]

5.3.1. Countercurrent flows Long-wave stability regions for countercurrent liquid–liquid flow in a b = 10° and 26° inclined channel, are shown in Fig. 11. The liquid properties are the same as those used in the experiments described in Ullmann et al. (2003a) in 1.44 cm pipe. The neutral stability boundary of the lower heavy phase holdup solution is denoted by LPD, and that of the high holdup by HPD (see Section 5.1). The upper bound for countercurrent flow is also shown in the figure (the flooding curve). As observed, there is a region of low flow rates where both solutions are stable, and another region where both solutions are unstable. This is in accordance with the experimental results, where for 10° inclination both, the LPD and HPD modes for this fluid system were observed to be smooth stratified flow in the operational window of U1s = 0– 0.005 m/s, U2s = 0–0.008 m/s (Zamir, 2003). The region where the LPD mode is stable extends to higher light phase flow rates compared to the HPD solution, whereas the region where the HPD is stable extend to higher heavy phase flow rates. Similar behavior is observed for 26° inclination, though for both LPD and HPD branches, the flow becomes unstable for lower flow rates compared to the case of 10° inclination.

-5

-U2S [m/s]

10

conditions associated with multiple holdup solutions. It includes countercurrent flows for which two different solutions exist, and triple solutions region of concurrent flows.

10

10

-1

10

-2

10

-3

S

upper solution

= q cr 10

-4

10

-5

16 7.4

2

S

U 3

lower solution 1

10

-5

10

-4

10

-3

10

-2

10

-1

1

10

-U2S [m/s] Fig. 12. Long-wave stability map in the single and triple solution regions for slightly upward inclined air–water flow.

the 3 solutions should be carefully examined. To elucidate the results obtained, the triple solution region stability map is plotted along with holdup curves for constant water flow rates, in

27

R. Kushnir et al. / International Journal of Multiphase Flow 62 (2014) 17–29

(a)

(a) 100 1

U1S [m/s]

3

10

upper & middle solutions

β =5°

10

r=1000 m=55 H=0.02 m

U

1 0.1

1

10

(b) 10

-2

U 10

-3

1

-4

β =5°

lower solution

r=1000 m=55 H=0.02 m

S

r=1000 m=55 H= 0.02m

10-5 0.1

-1

100

U2S [m/s]

β =0.1° 1

U

1

10

10

S lower solution

U1S [m/s]

Fig. 13a and b, respectively. As observed, in the 3-s region, a stable loop-shaped region appears. The middle holdup solution is stable within the entire 3-s region. Between line 1 and the l.h.s boundary of the 3-s region, the flow is stable for the lowest holdup solution; it becomes stable also for the upper holdup solution, between line 3 and the RHS of the 3-s boundary. Hence, even in gas–liquid systems there are regions where all three solutions are stable. The stable region extends also outside the 3-s region, where there is again a single solution. In this region, the holdup decreases sharply with small changes of the gas flow rate (e.g., U1S =  2.5  103 m/s in Fig. 13). It is worth emphasizing, that the region where stratified smooth or wavy flow is experimentally observed in upward inclined systems is in the vicinity of the 3-s region (see discussion in Ullmann et al., 2003b; Ullmann et al., 2004). As observed four regions emerge within the 3-s domain. In region I all the three solutions are stable, in region II the middle and upper solutions are stable, in region III only the middle solution is stable, and in region IV the middle and lower solutions are stable. The results become more comprehendible if the stability boundaries are examined while following the holdup variations at constant water flow rates (Fig. 13b). For example, for the holdup curve of U1S =  2.5  105 m/s, it is clearly seen that for gas velocities within region IV both the middle and lower solutions are stable. With increase of the gas flow rate (region III), the lower holdup becomes unstable and only the middle solution is stable. With a further increase of the gas flow rate (region II), the upper solution

1

10

100

-U1S [m /s]

-2

-U1S=2.5×10 m /s

10

U2S [m/s]

-2

Fig. 14. Long-wave stability maps in the single and triple solution regions for slightly downward inclined air–water flow. (a) high water flow rates; (b) low water flow rates.

-3

2.5×10 m /s

10

U

-3

S

-4

2.5×10 m /s

10

3

-4 -5

2.5×10 m /s

IV 1

10

I

III II

-5

1

(b) 0.5

10

-U2S [m /s] β =0.1° r=1000 m=55 H= 0.02m

0.4

-2

U1S=-2.5×10 m /s

0.3

~ h 0.2

-3

-2.5×10 m /s

0.1

-4

-2.5×10 m /s -5

-2.5×10 m /s

becomes stable as well. Eventually, for flow rates higher than the right boundary of the 3-s region, these two (upper and middle) stable solutions no longer exist, and the single (lower) solution is unstable. Long-wave stability regions for downward inclined air–water flow are shown in Fig. 14a and b for high and low liquid flow rates, respectively. Compared to horizontal flow, at low gas flow rates, the downward inclination reduces drastically the critical liquid flow rate for which the interface is stable. The critical liquid flow rate is constant for relatively wide range of gas flow rates (Fig. 14b), which is in accordance with experimental data for slightly inclined downward flow (e.g., Barnea et al., 1982). However, with further increase of the liquid flow rate, the region of 3-s is reached (Fig. 14a), where additionally to the lower holdup solution, the middle and upper solutions for the holdup are obtained. For these two additional modes the flow is stable (except for a very limited region near the tip of the 3-s region). Moreover, the solution corresponding to low holdup (which is unstable at lower liquid flow rates), becomes stable in a narrow range of liquid flow rates within the 3-s region. For even higher liquid flow rates (above the 3-s region), the single high holdup solution is stable.

0 1

-U2S [m /s]

10

Fig. 13. Triple solution region of slightly inclined air–water flow: (a) long-wave stability map; (b) holdup curves for constant water flow rate. The bold (dashed) curves denote stable (unstable) conditions, respectively.

6. Conclusions The linear stability with respect to long wave disturbances of laminar countercurrent and concurrent stratified flows in two-

28

R. Kushnir et al. / International Journal of Multiphase Flow 62 (2014) 17–29

plate geometry was investigated. The analysis was carried out by solving the well-known Orr-Sommerfeld equations. A closed-form solution for the eigenvalues and eigenfunctions was obtained through a formal power series in the wave number. The results are summarized in the form of stability boundaries on flow rates maps. Special attention is given to the stability of inclined flows in the regions of multiple solutions in countercurrent and concurrent flows. It is demonstrated that for countercurrent flow there is a region of low flow rates where both solutions are stable. The region where the lower holdup solution is stable extends to higher light phase flow rates, whereas the region where the upper holdup solution is stable extends to higher heavy phase flow rates. Likewise, the results of concurrent gas–liquid upward flows reveal a region where all three solutions are stable. Moreover, it was found that the middle solution is always stable within the entire 3-s domain. Another important aspect of the present study is deriving the expressions for calculations of the wave induced tangential and normal stresses. In two-fluid models, it is commonly assumed that the local shear stresses are adequately represented by quasi-steady models. Consequently, in the stability analysis, only the components of the stresses in phase with the wave height are considered. Indeed, the values of the shear stresses are insignificantly affected by the wave slope in the limit of long wave disturbances. Nevertheless, the expressions obtained in this study for the wave induced wall and interfacial shear stresses enable examining the components of the wave induced stresses, which are proportional to the wave slope. The proportionality coefficients were shown to be independent of the wave number. Therefore, the components of the stresses in phase with the wave slope are of the same order of the gravity force stabilizing term, and should not be discarded even in the framework of long wave stability analysis. The results from the exact long wave stability analysis obtained in this study can be used to identify the necessary modifications that should be introduced in two-fluid models in order to improve the prediction of the stability boundary of smooth stratified flow via the simplified one dimensional approach.

B22;0 ¼

þ

B32;0 ¼

    2 2 4 a1 B21;0 3ðc0 1ÞB31;0 b1 B11;0 y5 r 4 ðc0 1Þ B1;0 þb1 y w1 ðyÞ ¼  þ m 12 60 # a1 B31;0 y6 b1 B31;0 y7 þ þ 60 210     2 ðc0 1Þ B2;0 þb2 y4 a2 B22;0 3ðc0 1ÞB32;0 b2 B12;0 y5 þ w2 ðyÞ ¼ 12 60 a2 B32;0 y6 b2 B32;0 y7 þ þ 60 210 ðA:7Þ

B21;0 ¼

4

2

ð1  mÞ½mnð4n þ 2n2  3mÞ  4m2 þ n5  mðn þ 1Þð4mn þ 6mn2 þ n4 þ m2 þ 4n3 mÞ

ðA:1Þ

ðA:8Þ

B21;1

ðA:9Þ

B31;1

  2 0  n w2  n2 w2 þ mnw01 þ mw1  2n2 G B12;0  B11;0    i ¼ h n2 2mðn þ 1Þ B12;0  B11;0 þ ðm  n2 Þ B22;0  mB21;0   2 0  n w2  2n2 w2  n2 w01  2nw1  n2 G B22;0  mB21;0    i þ h n2 2mðn þ 1Þ B12;0  B11;0 þ ðm  n2 Þ B22;0  mB21;0

2ð1  mÞ½mn2 ð1  nÞ  m2 þ n5  mnðn þ 1Þð4mn þ 6mn2 þ n4 þ m2 þ 4n3 mÞ

B31;0

B32;0 ¼ m

B12;0

e ½n2 ð3 þ 2n  2n2 Þ  mð1 þ 2n þ 4n3 þ m þ 3n2 Þn Y ¼ nðn þ 1Þð4mn þ 6mn2 þ n4 þ m2 þ 4n3 mÞ ð1  mÞðm2 þ 4n3 m þ 2mn2  4n5  3n4 Þ þ nðn þ 1Þð4mn þ 6mn2 þ n4 þ m2 þ 4n3 mÞ

2mB31;1

w02

ðA:11Þ

þ w2 þ 2G

ðA:12Þ

B32;1 ¼ mB31;1  G ( ) 2   1 cos b k 1 ðr  1Þ  ðc0  1Þ B2;0 þ a2 þ G¼ 6 Fr2 We2

ðA:13Þ

¼



ðA:14Þ

In Eqs. (A.8)–(A.14), w1 and w01 are evaluated at y =  n, w2 and w02 at y = 1, and a1, a2, b1, b2 are the coefficients of the dimensionless velocity profile of the basic flow (see Section 3.1).

e ½mnð2  3n þ 7n2  2n3 þ mÞ  2n4  m2  Y mnðn þ 1Þð4mn þ 6mn2 þ n4 þ m2 þ 4n3 mÞ þ

2w1 n w01 w1 3 ¼ 2nB1;1 þ þ 2 n n

B11;1 ¼ n2 B31;1 þ w01 þ

B12;1 ¼ mB31;1 þ w02  2w2  G

e ½mnð4  3n  2n  n þ 2m þ 3mnÞ  n  2m  Y ¼ mðn þ 1Þð4mn þ 6mn2 þ n4 þ m2 þ 4n3 mÞ þ

e þ ð1  mÞðn4 þ m2  2mn2 Þ ½4mn þ ðn2 þ mÞð1 þ mÞn Y nðn þ 1Þð4mn þ 6mn2 þ n4 þ m2 þ 4n3 mÞ

A.2. First order solution coefficients

B22;1 B11;0

ðA:5Þ

ðA:10Þ

A.1. Zero order solution coefficients

3

2ð1  mÞðmn2  n3 m  m2 þ n5 Þ nðn þ 1Þð4mn þ 6mn2 þ n4 þ m2 þ 4n3 mÞ

ðA:6Þ

Appendix A. Asymptotic solution coefficients

2

e ½n3 ðn  1Þ þ mð2 þ 7n þ 3n2 þ 2n3 þ 2mÞn Y nðn þ 1Þð4mn þ 6mn2 þ n4 þ m2 þ 4n3 mÞ

ðA:2Þ

ðA:3Þ

A.3. Zero gravity solution The function f(m, n, r) appearing in Eq. (35), is linear with respect to r, and is given by the expression

f ðm; n; rÞ ¼

f1 ðm; nÞr þ f2 ðm; nÞ 420m2 ðn

2

þ 1Þ ðm2 þ 4mn þ 6mn2 þ 4mn3 þ n4 Þ

3

ðA:15Þ ðA:4Þ

where

R. Kushnir et al. / International Journal of Multiphase Flow 62 (2014) 17–29

f1 ðm;nÞ ¼ n14 þð62mÞn13 þð32m9Þmn12 þ120ðm1Þmn11 þð16m2 222m88Þmn10 ð444mþ60Þm2 n9 þð424272m2 þ646mÞm2 n8 þð1080m2 þ2344mþ224Þm2 n7 þð672m2 þ3021mþ1512Þm3 n6 þð1486m2 þ2086mþ224Þm3 n5 þð1043mþ392Þm4 n4 þ224m5 n3 ðA:16Þ f2 ðm;nÞ ¼ 224m2 n11 þð392mþ1043Þm2 n10 þð224m2 þ2086mþ1486Þm2 n9 þð1512m2 þ3021mþ672Þm2 n8 þð224m2 þ2344mþ1080Þm3 n7 þð424m2 þ646m272Þm3 n6 ð60mþ444Þm4 n5 ð88m2 þ222m16Þm4 n4 þ120ð1mÞm5 n3 þð329mÞm5 n2 þð6m2Þm6 nþm7 ðA:17Þ

References Amaouche, M., Mehidi, N., Amatousse, N., 2007. Linear stability of a two-layer film flow down an inclined channel: a second-order weighted residual approach. Phys. Fluids 19, 084106. Banner, M., Melville, W., 1976. On the separation of air flow over water waves. J. Fluid Mech. 77, 825–842. Barnea, D., Shioham, O., Taitel, Y., Dukler, A.E., 1982. Flow pattern transition for downward inclined two-phase flow: horizontal to vertical. Chem. Eng. Sci. 37, 735–740. Boomkamp, P.A.M., Miesen, R.H.M., 1996. Classification of instabilities in parallel two-phase flow. Int. J. Multiphase Flow 22, 67–88. Brauner, N., Moalem Maron, D., 1993. The role of interfacial shear modelling in predicting the stability of stratified two-phase flow. Chem. Eng. Sci. 48, 2867– 2879. Charru, F., Fabre, J., 1994. Long waves at the interface between two viscous fluids. Phys. Fluids 6, 1223–1235.

29

Hesla, T.I., Pranckh, F.R., Preziosi, L., 1986. Squire’s theorem for two stratified fluids. Phys. Fluids 29, 2808–2811. Hooper, A.P., Boyd, W.G.C., 1983. Shear-flow instability at the interface between two viscous fluids. J. Fluid Mech. 128, 507–528. Jeffreys, H., 1925. On the formation of water waves by wind. Proc. R. Soc. Lond. A 107, 189–206. Joseph, D.D., Renardy, Y.Y., 1993. Fundamentals of Two-Fluid Dynamics: Part I: Mathematical Theory and Applications. Springer, New-York. Kharif, C., Giovanangeli, J.-P., Touboul, C., Grade, L., Pelinovsky, E., 2008. Influence of wind on extreme wave events: experimental and numerical approaches. J. Fluid Mech. 594, 209–247. Kawai, S., 1982. Structure of air flow separation over wind wave crests. Bound.-Lay. Meteorol. 23, 503–521. Kuru, W.C., Sangalli, M., Uphold, D.D., Mcready, M.J., 1995. Linear stability of stratified channel flow. Int. J. Multiphase Flow 21, 733–753. Miles, J.W., 1962. On the generation of surface waves by shear flows. Part 4. J. Fluid Mech. 13, 433–448. Segal, V., 2008. Stability Analysis of long waves at the interface between two viscous fluids in an inclined channel. MSc thesis, Tel-Aviv University. Tilley, B.S., Davis, S.H., Bankoff, S.G., 1994. Linear stability of two-layer fluid flow in an inclined channel. Phys. Fluids 6, 3906–3922. Ullmann, A., Zamir, M., Ludmer, Z., Brauner, N., 2003a. Stratified laminar countercurrent flow of two liquid Phases in inclined tube. Int. J. Multiphase Flow 29, 1583–1604. Ullmann, A., Zamir, M., Gat, S., Brauner, N., 2003b. Multi-holdups in co-current stratified flow in inclined tubes. Int. J. Multiphase Flow 29, 1565–1581. Ullmann, A., Goldstien, A., Zamir, M., Brauner, N., 2004. Closure relations for the shear stresses in two-fluid models for stratified laminar flows. Int. J. Multiphase Flow 30, 877–900. Vempati, B., Oztekin, A., Neti, S., 2010. Stability of two-layered fluid flows in an inclined channel. Acta Mech. 209, 187–199. Yiantsios, S.G., Higgins, B.G., 1988. Linear stability of plane Poiseuille flow of two superposed fluids. Phys. Fluids 31, 3225–3238. Yih, C.S., 1967. Instability due to viscosity stratification. J. Fluid Mech. 27, 337–352. Zamir, M., 2003. An investigation into the flow phenomena in phase transition extraction based on solvent mixtures with critical point miscibility. PhD thesis, Tel-Aviv University.