The stability of two immiscible fluids with different densities and viscosities is examined for channel flow. A multi-domain Chebyshev collocation spectral method is used for solving the coupled Orr-Sommerfeld stability equations for the entire spectrum of eigenvalues and associated eigenfunctions. Numerical solution of the eigenvalue problem is obtained with the QZ eigenvalue solver and is validated with analytical results derived in the long and short wave limits. A parametric study is carried out to investigate the spectral characteristics and eigenfunction structures related to the shear and interfacial modes of instability. The interactions and mode switching between the two instability modes are investigated. Under certain conditions, the two modes are found to become unstable simultaneously. Mode coalescence can occur in either the stable or the unstable part of the spectrum. In general, the eigenfunction of the most dangerous mode is observed to vary sharply in the neighborhood of the interface and the critical points.

The stability of parallel two-phase flows has received much attention in the literature because of their importance in engineering applications, such as coating, polymer extrusion and oil transportation as well as for their significance in natural phenomenon.1,2 The study of the general stability problem for parallel two-phase flow is thus important and started with the pioneering works of Miles,3 Yih,4,5 and others. Miles3 investigated the stability of a liquid film sheared by a turbulent boundary layer and discussed various physical mechanisms with regards to the transfer of energy from the wind to the interfacial waves. Charles and Lilleleht6 experimentally investigated interfacial waves generated by a mean-flow instability. The appearance of waves in their experiments apparently coincided with the transition to turbulence of the less viscous fluid. Cohen and Hanratty7 determined the critical value of the gas velocity at which the waves first start to appear. Craik8 considered the instability of relatively thinner fluid films and pointed out the existence of two different instability modes as a function of liquid layer thickness.

A long-wave linear stability analysis by Yih5 for channel flow, illustrated the origin of instability associated with the difference in the viscosity of the two fluids. The author reported that apart from the Tollmien-Schliching (shear) mode of instability, there exists a second instability mode (interfacial), which is invoked by the viscosity difference between the two fluids. Thus, the viscosity contrast can give rise to an interfacial mode that is unstable even for arbitrarily small values of the Reynolds number. In general, the interfacial mode can coexist with the instability of the shear mode. Wave generation in two-fluid systems was subsequently studied by considering hydrodynamic stability1,9 in both of the fluids. Yiantsios and Higgins10 carried out a numerical study of the linear stability of plane Poiseuille flow of two superposed fluids and confirmed previous findings. The studies by Yih4,5 and Hooper11 focused on the interfacial instability mode for relatively small values of the Reynolds number. For large values of the Reynolds number, Özgen et al.12 found that shear modes dominate the instability behavior in the case of gas boundary layer over a newtonian liquid.

Hooper and Boyd13 examined the flow of two superposed fluids in linear shearing motion bounded by a wall. For large values of the Reynolds number, waves longer than the thickness of the lower fluid arise, that can be longer or shorter than the thickness of the viscous sublayer in the upper fluid, depending on whether the bottom layer is more or less viscous. Timoshin14 considered the same problem using the triple-deck theory and recovered the interfacial modes predicted by Hooper and Boyd.13 Timoshin and Hooper15 extended the triple-deck analysis to examine a broader range of parameters, including variations in film thickness in particular. In these studies, the two instability modes related to Tollmien-Schlichting and interfacial waves were shown to have points of mode coalescence located in either the stable or unstable part of the spectrum. Ozgen16 also studied mode coalescence for the two-fluid boundary-layer problem. He showed that coalescence occurs in the unstable part of the domain when the sheared fluid thickness is sufficiently large. He also reported that for a specific set of parameters, mode coalescence occurs in both the stable and unstable parts of the spectrum.

To date, a detailed characterization of shear and interfacial modes of instability has been carried out mostly for two-phase boundary layer and shear flow problems. The corresponding behavior of two-phase channel flow remains relatively unexplored. For example, the structure of the spectrum and the eigenfunction profiles of the two instability modes are not well studied. Moreover, the effect of the interaction and coalescence of the modes on instability mechanisms is also not well known. The domains of the existence of shear and interfacial modes of instabilities are thus not well studied. It is therefore also not clear for example how the two-phase boundary layer problem compares with the two-phase channel flow problem.

The present work is concerned with the linear stability analysis of two-phase channel flow. The shear, and interfacial modes are identified and the effects of specific parameters on these modes are investigated. The objective is to understand the nature of competing stability mechanisms. The primary contributions of this study are the characterization of the spectrum and the identification of coalescence for the shear and interfacial modes of instability. Simultaneous instability of these modes, as well as their coalescences in both the stable and the unstable parts of the spectrum is shown to exist. It is found that the eigenfunctions associated with the most dangerous modes exhibit characteristic extrema in the neighborhood of the interface and near the critical points. Such information is particularly useful for the initialization of direct numerical simulation of two-phase flow for small amplitude perturbations.

This article is organized as follows: In Sec. II, we present the geometry and governing equations. The linearized stability equations are presented in Sec. III and the asymptotic analysis in the long and short wave limits is derived in Sec. IV. Then, in Sec. V, we present the numerical method and in Sec. VI, we discuss the stability results. A discussion regarding the effect of the base state and comparison with experimental results is presented in Sec. VII. We summarize our conclusions in Sec. VIII.

We consider the flow of two superposed, parallel fluid layers of different viscosities and densities in a channel. The upper layer (vapor phase for example) is a fluid of viscosity μ1, density ρ1, and thickness d1 while the lower layer (liquid phase) has viscosity μ2, density ρ2, and thickness d2.

Squire’s transformation17 for the parallel flow of a uniform fluid and later extensions by Yih4,5 for stratified fluids indicate that the stability of a three-dimensional disturbance can be determined from that of an equivalent two-dimensional disturbance at a lower Reynolds number.1,9 It is therefore sufficient to consider only two-dimensional disturbances. In the following, we adopt a cartesian coordinate system (X, Y) such that X represents the streamwise direction and Y is the direction normal to the direction of motion. The undisturbed flow is parallel to the X-axis with origin at the interface. The schematic representation of the two-layer fluid system in a channel is shown in Figure 1.

FIG. 1.

Schematic of two-layer fluid flow in a channel.

FIG. 1.

Schematic of two-layer fluid flow in a channel.

Close modal

The dimensionless parameters used for the linear stability analysis include the ratios of viscosity, layer thickness, and density, denoted, respectively, by

m = μ 2 μ 1 , n = d 2 d 1 , r = ρ 2 ρ 1 .

The Reynolds number, Weber number, and the Froude number are, respectively, of the forms

Re = ρ 1 d 1 U 0 μ 1 , We = ρ 1 d 1 U 0 2 γ , Fr = U 0 g d 1 ,

where g is the gravitational acceleration and γ is the interfacial tension. The velocity at Y = 0 is denoted by U0.

The motion of each Newtonian fluid is governed by the conservation of mass and momentum. The dimensionless governing equations for incompressible flow in the upper layer are

. u 1 = 0 ,
(1)
u 1 t + ( u 1 . ) u 1 = 1 Re Δ u 1 p 1 1 Fr 2 e y .
(2)

In the lower liquid layer, the equations read

. u 2 = 0 ,
(3)
u 2 t + ( u 2 . ) u 2 = m r Re Δ u 2 1 r p 2 1 Fr 2 e y ,
(4)

where u1 and u2 are the dimensionless velocities in the upper and lower layers, respectively, and p denotes the pressure. The above equations are made dimensionless by scaling velocity with U0, time with d1/U0, and pressure with ρ 1 U 0 2 .

The basic flow is steady and parallel to the x-axis. The basic velocity profiles in the upper and lower layers, scaled with U0, are given by

U 1 = 1 + α 1 y + β 1 y 2 , U 2 = 1 + α 2 y + β 2 y 2 ,
(5)

where, x = X/d1 and y = Y/d1 and

α 1 = m n 2 n ( 1 + n ) , β 1 = ( m + n ) n ( 1 + n ) ,
α 2 = m n 2 m n ( 1 + n ) , β 2 = ( m + n ) m n ( 1 + n ) .

To understand the mechanisms governing the stability of two-layer channel flow, we carry out a linear stability analysis in the domain −n < y < 1. Velocity, pressure, and interfacial position are decomposed, respectively, as

u i = U i + u i , v i = v i , p i = P i + p i , h = h 0 + η ,
(6)

where Ui and Pi are the dimensionless velocities and pressure of the basic flow, respectively. Subscript i denotes the phase and primes denote perturbations. The linearized problem in the x-y plane can be expressed in terms of the streamfunctions ψi, where

u i = ψ i / y , v i = ψ i / x .

The linearized equations and boundary conditions suggest a normal mode analysis where the disturbance quantities can be expressed as

( ψ i , p , η ) = ( ϕ i ( y ) , f i ( y ) , h ) exp [ i α ( x c t ) ] , i = 1 , 2 ,
(7)

where ϕi(y) and fi(y) are the disturbance amplitudes, α is the dimensionless wavenumber defined by 2πd1/λ, with wavelength λ, c = cr + icI is the dimensionless wave velocity, and the most rapidly growing wave is the one for which the rate of amplification αcI is a maximum.

Upon substitution of the expressions (6) and (7) into the governing equations (1)(4) and by neglecting quadratic terms, we obtain the linearized vorticity equations for disturbance quantities, which lead to the following Orr-Sommerfeld equations in each fluid:

ϕ 1 I V 2 α 2 ϕ 1 + α 4 ϕ 1 = i α Re [ ( U 1 c ) ( ϕ 1 α 2 ϕ 1 ) U 1 ϕ 1 ] , 0 < y < 1 ,
(8)
ϕ 2 I V 2 α 2 ϕ 2 + α 4 ϕ 2 = i α m r n Re [ ( U 2 c ) ( ϕ 2 α 2 ϕ 2 ) U 2 ϕ 2 ] , n < y < 0 .
(9)

The above equations will be solved subject to eight boundary conditions, two at the wall and two kinematic and four dynamics conditions at the interface. At the walls, y = − n and y = 1, the no slip condition gives

ϕ 1 ( 1 ) = 0 , ϕ 1 ( 1 ) = 0 ,
(10)
ϕ 2 ( n ) = 0 , ϕ 2 ( n ) = 0 .
(11)

At the interface, y = η(x, t), the velocity and tangential shear stress are continuous and the jump in the normal stress is balanced by surface tension. The kinematic free surface condition also holds at the interface. The linearized boundary conditions satisfy

ϕ 1 ( 0 ) = ϕ 2 ( 0 ) ,
(12)
ϕ 1 ( 0 ) ϕ 2 ( 0 ) = ϕ 1 ( 0 ) c ̃ ( α 2 α 1 ) ,
(13)
ϕ 1 ( 0 ) + α 2 ϕ 1 ( 0 ) = m ( ϕ 2 ( 0 ) + α 2 ϕ 2 ( 0 ) ) + 2 ( m β 2 β 1 ) ϕ 1 ( 0 ) c ̃
(14)
i α Re ( c ̃ ϕ 1 ( 0 ) + a 1 ϕ 1 ( 0 ) ) ϕ 1 + 3 α 2 ϕ 1 ( 0 ) + i r α Re ( c ̃ ϕ 2 ( 0 ) + α 2 ϕ 2 ( 0 ) ) + m ( ϕ 2 ( 0 ) 3 α 2 ϕ 2 ( 0 ) ) = i α Re ( F + α 2 / W e ) ϕ 1 ( 0 ) / c ̃ ,
(15)

where c ̃ = c U 1 ( 0 ) and F = ( r 1 ) g d 1 / U 0 2 is the ratio of the gravitational force to the inertial force. In Eqs. (13)(15), the kinematic boundary condition h = ϕ 1 ( 0 ) / c ̃ was incorporated. Equations (8) and (9) together with boundary conditions (10)(15) constitute the generalized eigenvalue problem for eigenvalue, c = c(α, Re, We, Fr, m, n, r) from which instability is determined by ci > 0. Neutral stability is associated with ci = 0. The case r = 1,  n = 1, corresponds to a single fluid system, which was treated by Yih5 and Higgins and Yiantsios.10 

The asymptotic analysis of the eigenvalue for long waves was carried out by Yih5 for two-layer Couette and Poiseuille flows in the case of equal densities. Here, we consider the stability problem in the long wave limit including the effect of density and viscosity difference between the two fluids. We shall use the method of non-singular perturbation around α = 0 to derive the analytic solutions. The method proceeds along the lines of Higgins and Yiantsios10 analysis. In long-wavelength asymptotics, one expands the eigenfunctions and the eigenvalue as a regular perturbation series in α,

ϕ 1 ϕ 10 + α ϕ 11 + α 2 ϕ 12 + ,
ϕ 2 ϕ 20 + α ϕ 21 + α 2 ϕ 22 + ,
(16)
c c 0 + α c 1 + α 2 c 2 + .

When the expansions defined by Eq. (16) are substituted into Eqs. (8)(13) and terms of like order in α are grouped together, a simplified set of fourth-order equations and boundary conditions for ϕ1 and ϕ2 result. These equations may be solved sequentially to determine ϕ1 and ϕ2 and an eigenvalue relation for ci at each order of the approximation.

In the first approximation, all terms containing α in the differential system are ignored. For α = 0, Eqs. (8) and (9) become

ϕ 10 I V = 0 , ϕ 20 I V = 0 .
(17)

At the leading order, the boundary conditions satisfy

ϕ 10 ( 1 ) = 0 , ϕ 10 ( 1 ) = 0 ,
ϕ 20 ( n ) = 0 , ϕ 20 ( n ) = 0 ,
ϕ 10 ( 0 ) = ϕ 20 ( 0 ) ,
ϕ 10 ( 0 ) ϕ 20 ( 0 ) = ϕ 10 ( 0 ) c 1 ( α 2 α 1 ) ,
(18)
ϕ 10 ( 0 ) m ϕ 20 ( 0 ) = 0 ,
ϕ 10 ( 0 ) m ϕ 20 ( 0 ) = 0 .

The eigenfunctions ϕ10 and ϕ20, at zero order, are third-order polynomials of the forms

ϕ 10 = A 10 + B 10 y + C 10 y 2 + D 10 y 3 ,
ϕ 20 = A 20 + B 20 y + C 20 y 2 + D 20 y 3 .
(19)

From boundary conditions (18), we find the following expressions for the constants in (19):

B 10 = m + 3 n 2 + 4 n 3 2 n 2 ( 1 + n ) , B 20 = n 3 + 3 m n + 4 m 2 m n ( 1 + n ) ,
C 10 = m C 20 , C 20 = m + n 3 m n 2 ( 1 + n ) ,
(20)
D 10 = m D 20 , D 20 = n 2 m 2 m n 2 ( 1 + n ) .

Without loss of generality, we have set the arbitrary constants to unity, i.e., A10 = A20 = 1. At this order of approximation, the eigenvalue c0 is real and may be expressed in terms of the thickness ratio n and viscosity ratio m,

c 0 = 1 + 2 n ( m n 2 ) ( m 1 ) m 2 + 2 m n ( 2 + 3 n + 2 n 2 ) + n 4 .
(21)

Note that this equation is independent of r and the Eqs. (19)-(21) are identical to Yih’s5 Equations (33) and (34).

In the second approximation, all terms containing α2 and higher orders of α are ignored, and the equations to be solved are

ϕ 11 I V = i α Re [ ( U 1 c 0 ) ϕ 10 U 1 ϕ 10 ] ,
ϕ 21 I V = i α Re m 1 n r [ ( U 1 c 0 ) ϕ 20 U 2 ϕ 20 ] .
(22)

The solutions of these equations are again polynomials, though of higher order

ϕ 11 = B 11 y + C 11 y 2 + D 11 y 3 + i α Re h 1 ( m , n , y ) ,
ϕ 21 = B 21 y + C 21 y 2 + D 21 y 3 + i α Re m 1 r h 2 ( m , n , r , y ) .
(23)

The functions h1(m, n, y) and h2(m, n, y) result from the inhomogeneous terms in perturbation equations (22). Further details are given in the Appendix. In Eq. (23), the first three terms constitute the particular solution. The term of zero degree in y is taken to be zero, as demanded by condition (12).

Boundary conditions (10), (11), and (14) assume the forms

B 11 + C 11 + D 11 + i α Re h 1 ( m , n , 1 ) = 0 ,
(24)
B 11 + 2 C 11 + 3 D 11 + i α Re h 1 ( m , n , 1 ) = 0 ,
(25)
B 21 n + C 21 n 2 D 21 n 3 + i α Re m 1 r h 2 ( m , n , y = n ) = 0 ,
(26)
B 21 n 2 C 21 n + 3 D 21 n 2 + i α Re m 1 r h 2 ( m , n , y = n ) = 0 ,
(27)
m C 21 = C 11 .
(28)

We are then left with conditions (13) and (15) to contend with. To the present order of approximation O(α), Eq. (15) assumes the following form at y = 0:

m ϕ 21 ϕ 11 = i α Re [ ( F ϕ 10 / c ̃ ) r ( c ̃ ϕ 20 + α 2 ϕ 20 ) + ( c ̃ ϕ 10 + α 1 ϕ 10 ) ] .
(29)

Conditions (12) and (13), when applied to the first approximation, must satisfy

ϕ 10 ( 0 ) ϕ 20 ( 0 ) = ϕ 10 ( 0 ) ( α 2 α 1 ) / c ̃ and ϕ 10 ( 0 ) = ϕ 20 ( 0 ) .
(30)

Hence, (29) can be written further as

m ϕ 21 ϕ 11 = i α Re [ ( F ϕ 10 / c ̃ ) ( r 1 ) ( c ̃ ϕ 10 + α 1 ϕ 10 ) ] , at y = 0 .
(31)

Thus,

6 m D 21 6 D 11 = i α Re [ F / c ̃ ( r 1 ) ( c ̃ B 10 + α 1 ) ] .
(32)

As for (13), its form for the second approximation takes some care because c ̃ also suffers a perturbation. With this in mind, and remembering that both ϕ11(0) and ϕ21(0) are zero, (13) becomes

ϕ 21 ( 0 ) ϕ 11 ( 0 ) = c 1 ( α 2 α 1 ) ϕ 10 ( 0 ) / c ̃ 2 ,
(33)

in which c1 is the change in c and is of course identical to the change in c ̃ = c 1 . Hence,

( B 21 B 11 ) c ̃ 2 = c 1 ( α 2 α 1 ) .
(34)

By solving Eqs. (24)-(28) and (32), the six coefficients B11, C11, D11, B21, C21, and D21 can be found. Then c1 is determined from (34). The result is

c 1 = i c I , c I = 8 α Re H 1 ( m , n , r , F ) ,
(35)

in which H1 depends on the viscosity, thickness ratios m and n, the density ratio r, and the constant F as shown in the Appendix. Note that this asymptotic analysis is applicable to fluids of different densities as well as different viscosities. After a brief calculation, we can show that (35) gives for n = r = 1,

c I = O ( ( m 1 ) 2 ) ,

which confirms Yih’s results in the long wave limit and shows that the instability is due to viscosity variation, m≠1.

For n ≠ 1, we can see that cI vanishes as m − 1 as m approaches 1. Thus, it is interesting to find that cI vanishes as either m − 1 or (m − 1)2 as m approaches unity depending upon the layer thickness. The mode considered here is different from the long wave modes in the usual theory and thus, does not contradict the known result that there are no neutral modes for long waves at low Reynolds numbers.

The asymptotic analysis in the short wave limit proceeds along the same lines as the analysis of Yiantsios and Higgins10 for the plane Poiseuille flow and Hooper and Boyd13 for the plane Couette flow. The flow coordinates (x, y) are scaled with α, i.e.,

x 1 = α x , y 1 = α y ,

and the following asymptotic expansions are assumed for ϕ1, ϕ2, and c:

ϕ 1 ϕ 10 + α 1 ϕ 11 + α 2 ϕ 12 + ,
ϕ 2 ϕ 20 + α 1 ϕ 21 + α 2 ϕ 22 + ,
(36)
c c 0 + α 1 c 1 + α 2 c 2 + .

The perturbation analysis leads to the following expression for the wave speed:

c = 1 + 5 8 α 1 Re 2 ( α S ) ( 1 m ) ( 1 + m ) 2 m α 3 + i [ R 1 ( α S ) 2 ( 1 + m ) α 1 + ( R 1 α 1 2 ( α S ) ( 1 m ) 2 2 ( 1 + m ) m 2 3 8 Re 3 ( α S ) 2 ( 1 + m ) 3 ) α 3 ] .
(37)

The above formula was derived under the assumption that αS = O(1) and r = 1. This result shows that interfacial tension and viscosity stratification play an important role in the short wave limit. The stabilizing effect of the interfacial tension enters at O(α−1) in the wave speed expansion, whereas the destabilizing effect of viscosity stratification enters at O(α−3).

Upon neglecting the effects of gravity and interfacial tension, one obtains a relatively simple expression

c = 1 + i Re α 1 2 ( 1 m ) 2 2 ( 1 + m ) m 2 α 3 .
(38)

This result shows that short waves are unstable unless r > m2 > 1 or 1 > m2 > r.

It is important to calculate the eigenvalues and eigenfunctions accurately for the entire spectrum of unstable modes by numerically solving the generalized eigenvalue problem. In this study, we use a Chebyshev collocation spectral method for this purpose. This choice is appropriate for non-periodic, bounded domains.18 The basic unknowns are defined in terms of Chebyshev polynomials. These polynomials are orthogonal on [ − 1, 1]. Therefore, we transform Orr-Sommerfeld equations (8) and (9) to the interval [ − 1, 1] by a change of the independent variable y. This is easily achieved by means of the linear transformations,

z = 2 n y + 1 , n < y < 0 ,
(39)
z = 2 y 1 , 0 < y < 1 .
(40)

After transforming the Orr-Sommerfeld equations in this manner, we approximate the eigenfunctions ϕ1(z) and ϕ2(z) by the truncated Chebyshev expansions,

ϕ 1 ( i ) ( z ) = n = 0 N 1 a n T n ( i ) ( z ) and ϕ 2 ( i ) ( z ) = n = 0 N 2 b n T n ( i ) ( z ) ,
(41)

where Tn(z) are the Chebyshev polynomials, i denotes the ith derivative with respect to z, and an and bn are the discrete Chebyshev expansion coefficients.

To achieve the highest possible accuracy, it is sometimes important to control the number of points in a specified subsection of the physical domain. Derivatives can be obtained by differentiating the Chebyshev polynomials in (41). This is different from Trefethen18 where orthogonality conditions and recurrence relations for Chebyshev polynomials were employed to arrive at expressions without derivatives.

Upon substitution of (41) into the discretized Orr-Sommerfeld equations and boundary conditions, we obtain a generalized eigenvalue problem of the form

A Ψ = c B Ψ ,
(42)

where the (N1 + N2 + 2) components’ vector Ψ is given by

Ψ = ( a 0 , a 1 , a 2 , , a N 1 , b 0 , b 1 , b 2 , , b N 2 ) T .
(43)

The matrix B is singular, because some of the boundary conditions lead to empty rows. Due to this, the QZ algorithm may produce singular eigenvalues, which can potentially interfere with finite eigenvalues. We have chosen to use the first, second, N1th, N1 + 1th, N1 + 2th, N1 + 3th, last, and next-to-last row of B to implement the eight boundary conditions. The same rows in the matrix A can be chosen as a complex multiple of the corresponding rows in B. By carefully selecting this complex multiple, the spurious modes associated with the boundary conditions can be mapped to an arbitrary location in the complex plane. We thus have a N1 + N2 + 2 matrix eigenvalue problem which is solved using the QZ algorithm.

In this section, we present the stability behavior of two-phase flow in a channel. A wide range of new stability features is discussed. These include the structure of the spectra, the behaviors of the eigenfunctions, and mode coalescence among unstable modes. Eigenvalues are calculated numerically and compared with the asymptotic solution for validation. The number of Chebyshev modes used in the computation varies. We carefully checked that we used a sufficient number of modes to obtain converged results.

Before carrying out the parametric study, we first validated our numerical method by comparing our results with the numerical results of Higgins and Yiantsios in Ref. 10. As shown in Table I, the numerical values of the complex wave speed of the least stable modes predicted by our numerical method agree very well with those predicted by Higgins and Yiantsios in Ref. 10.

TABLE I.

Asymptotic and numerical values for the complex wave speed c (m = 5, n = 1, r = 1).

α Re αS Asymptotic solution Numerical solution of Higgins and Yiantsios10  Numerical solution (present work)
10  0.999 97−0.008 122i  0.999 98−0.008 199i  0.999 982−0.008 199i 
10  0.999 30−0.040 182i  0.999 56−0.041 184i  0.999 561−0.041 184i 
10  0.999 94−0.016 460i  0.999 96−0.016 537i  0.999 964−0.016 537i 
10  0.998 61−0.083 134i  0.999 07−0.083 491i  0.999 073−0.083 491i 
20  0.999 997−0.004 140i  0.999 997−0.004 145i  0.999 997−0.004 145i 
20  0.999 913−0.020 727i  0.999 929−0.020 75i  0.999 929−0.020 75i 
20  0.999 993−0.008 307i  0.999 994−0.008 312i  0.999 994−0.008 312i 
20  0.999 826−0.041 641i  0.999 855−0.041 664i  0.999 855−0.041 664i 
α Re αS Asymptotic solution Numerical solution of Higgins and Yiantsios10  Numerical solution (present work)
10  0.999 97−0.008 122i  0.999 98−0.008 199i  0.999 982−0.008 199i 
10  0.999 30−0.040 182i  0.999 56−0.041 184i  0.999 561−0.041 184i 
10  0.999 94−0.016 460i  0.999 96−0.016 537i  0.999 964−0.016 537i 
10  0.998 61−0.083 134i  0.999 07−0.083 491i  0.999 073−0.083 491i 
20  0.999 997−0.004 140i  0.999 997−0.004 145i  0.999 997−0.004 145i 
20  0.999 913−0.020 727i  0.999 929−0.020 75i  0.999 929−0.020 75i 
20  0.999 993−0.008 307i  0.999 994−0.008 312i  0.999 994−0.008 312i 
20  0.999 826−0.041 641i  0.999 855−0.041 664i  0.999 855−0.041 664i 

We compare our results with the asymptotic solution and find good agreements as shown in Figures 2(a)-2(b). In the long wave limit, the imaginary part of an unstable mode is plotted as a function of αRe for m = 30, n = 0.1, r = 100, Re = 5000, We = 80, and Fr = 36. As shown in Figure 2(a), the numerical solution agrees well with asymptotic expression (35). In Figure 2(b), we also compare asymptotic formula (37) with the numerical solution of the eigenvalue problem, for m = 30, n = 1, r = 1, Re = 5000, We = ∞, and Fr = 36. For α > 100, the numerical values for cI are in good agreement with asymptotic results. This helps to develop confidence in the numerical solution for intermediate wavenumbers that is needed for identifying the most dangerous mode. We have further checked our numerical solution with the asymptotic results of Yih in Ref. 5 and Hooper and Boyd in Refs. 13 and 19 as well as the numerical results of Boomkamp.20,21

FIG. 2.

Imaginary part of the wave speed ci, for an unstable mode, in the long and short wave limits, respectively. Flow parameters: (a) m = 30, n = 0.1, r = 100, Re = 5000, We = 80, and Fr = 36. (b) m = 30, n = 1, r = 1, Re = 5000, We = ∞, and Fr = 36; the solid lines represent the analytic solutions given by Eq. (37); the dashed lines are the numerical solutions.

FIG. 2.

Imaginary part of the wave speed ci, for an unstable mode, in the long and short wave limits, respectively. Flow parameters: (a) m = 30, n = 0.1, r = 100, Re = 5000, We = 80, and Fr = 36. (b) m = 30, n = 1, r = 1, Re = 5000, We = ∞, and Fr = 36; the solid lines represent the analytic solutions given by Eq. (37); the dashed lines are the numerical solutions.

Close modal

Some of the stability characteristics of plane Poiseuille flow have been studied in the literature. To extend these previous studies, we first reproduce the results reported by Drazin,1 Schmid,9 Orszag,22 and Dongarra23 for single fluid Poiseuille flow for Re = 10 000, α = 1, m = 1, n = 1, and r = 1. We also show that in addition to the shear mode found in the classical hydrodynamic stability theory, additional interfacial modes may occur due to viscosity and density stratifications in two-fluid flows as reported in previous works of Yih5 and Higgins and Yiantsios.10 

When we introduce a viscosity contrast between the two-layers, several features emerge. As we increase the viscosity ratio, the lower fluid spanning the region −1 ≤ y ≤ 0 becomes more viscous with respect to the layer above, and a substantial change in the structure of the spectra and in the shapes of the eigenfunctions is observed. To understand why, it is instructive to consider the structures of the spectra and to examine how the eigenfunctions behave near the interface and the critical points as we increase the viscosity ratio. This is also of interest for understanding the nature of the competitive instabilities arising in two-phase flows and for identifying the shear and interfacial modes.

In the absence of viscosity contrast (m = 1), we obtain the classical Y-shaped eigenspectrum. Figure 3 shows a portion of the set of discrete eigenvalues cn for nonzero α. The eigenvalues are located on three main branches which have been labeled A(cr → 0), P(cr → 1), and S(cr → 2/3) by Mack.24 The least stable eigenvalues on the A and P branches are related to the shear and interfacial modes, respectively. As example, the least stable mode which corresponds to the eigenvalue c = 0.237 526 + 0.003 739 6i is located on the A branch and is associated with the shear mode.

FIG. 3.

Orr-Sommerfeld spectrum for plane Poiseuille flow with Re = 10 000, wave numbers α = 1, N1 = N2 = 100, m = 1, n = 1, and r = 1.

FIG. 3.

Orr-Sommerfeld spectrum for plane Poiseuille flow with Re = 10 000, wave numbers α = 1, N1 = N2 = 100, m = 1, n = 1, and r = 1.

Close modal

As shown in Figure 4, the spectra for m = 1/2 and m = 10 have structures similar to the one described for m = 1, although there are some noteworthy differences. Most importantly, the S-branch splits into two new branches (denoted S1 and S2) and two new branches P1 and P2 appear to intersect the real axis at values close to the interface velocity (cr ≈ 1). Branches S1 and S2 converge to the branch S as m → 1; and the branches A1, A2, P1, and P2 likewise converge to the A and P branches.

FIG. 4.

Orr-Sommerfeld spectrum for Re = 10 000, wave numbers α = 1, N1 = N2 = 100, β = 0, m = 1, n = 1, and r = 1. (a) m = 1/2, (b) m = 10.

FIG. 4.

Orr-Sommerfeld spectrum for Re = 10 000, wave numbers α = 1, N1 = N2 = 100, β = 0, m = 1, n = 1, and r = 1. (a) m = 1/2, (b) m = 10.

Close modal

For m = 1/2, the shear mode c = 0.1713 − 0.000 49i is observed in Figure 4(a) to be stable and located on the A branch, and the least stable interfacial mode c = 0.9999 + 10−9i is located on the P branch. The highly damped S1 and S2 modes occur with phase speeds cr → 0.64 and cr → 0.74 for m = 1/2 and with phase speeds cr → 0.61 and cr → 1.42 for m = 10.

For the viscosity ratio m = 10, another unstable mode c = 1.453 + 0.2205i appears in the upper plane ci > 0, attributed to an interfacial instability mode. This eigenvalue is insensitive to an increase in the Reynolds number and hence, remains unchanged at very high Reynolds numbers as illustrated in Table II. Consequently, the instability of this mode is not primarily due to inertial effects. It is also noted that the real part of this eigenvalue remains constant (cr = 1.454) at any Reynolds number Re > 13 000 and the imaginary part of the eigenvalue ci → 0.2236, which has not been noted in the past. The exact limiting values at other wavenumbers are difficult to determine analytically except in the long and short wave limits. Note that in the inviscid case, the asymptotic results in free surface shear flow, with the numerical solution of the eigenvalue problem for plane Poiseuille flow, are derived recently by Kaffel and Renardy.25,26 In their study, it is shown that, in the short wave limit, there are two unstable eigenmodes localized near the free surfaces which can be combined into even and odd eigenfunctions. The derivation of these asymptotics showed that the difference between the two eigenvalues tends to zero exponentially.

TABLE II.

Numerical values for the complex wave speed c for α = 1, m = 10, n = 1, r = 1, and We = ∞.

Re c
10 000  1.453+0.2205i 
12 000  1.453+0.2213i 
15 000  1.454+0.2216i 
30 000  1.454+0.2222i 
50 000  1.454+0.2229i 
100 000  1.454+0.2230i 
500 000  1.454+0.2233i 
1 000 000  1.454+0.2236i 
100 000 000  1.454+0.2236i 
Re c
10 000  1.453+0.2205i 
12 000  1.453+0.2213i 
15 000  1.454+0.2216i 
30 000  1.454+0.2222i 
50 000  1.454+0.2229i 
100 000  1.454+0.2230i 
500 000  1.454+0.2233i 
1 000 000  1.454+0.2236i 
100 000 000  1.454+0.2236i 

Typical shapes of eigenfunctions, given in term of normal velocities v, are plotted for Re = 10 000, α = 1, m = 1, n = 1, and r = 1 as shown in Figure 5. The solid and dashed curves are associated with the real and imaginary parts of v, denoted as v(r) and v(i), respectively, and all eigenfunctions are normalized such that v(r) tends to unity at y = 0. These eigenfunctions are associated with three eigenvalues, one from each of the three branches. It is apparent from the shape of the eigenfunctions why one usually designates the A modes as wall modes (Figure 5(a)) and the P modes as center modes (Figure 5(b)). The A modes, which have their largest variation close to the wall, have rather small phase velocities, whereas the P modes, which have their maxima close to the center of the channel (y = 0), have much higher phase speeds. In the absence of viscosity and density stratifications, the S modes, which are highly damped, have a phase speed that is nearly equal to 2/3 as shown by Schmid and Henningson9 and have their maxima close to y = ± 1/2 (Figure 5(c)).

FIG. 5.

Orr-Sommerfeld eigenfunctions for plane Poiseuille flow with Re = 10 000, wave numbers α = 1, m = 1, n = 1, and r = 1. (a) A-branch; (b) P-branch; and (c) S-branch. The solid and dashed lines represent the real and imaginary parts of the normal velocity.

FIG. 5.

Orr-Sommerfeld eigenfunctions for plane Poiseuille flow with Re = 10 000, wave numbers α = 1, m = 1, n = 1, and r = 1. (a) A-branch; (b) P-branch; and (c) S-branch. The solid and dashed lines represent the real and imaginary parts of the normal velocity.

Close modal

The eigenfunctions in Figure 5(a) are associated with the first even shear mode of the A family that corresponds to the slightly unstable eigenvalue c = 0.237 526 + 0.003 739 6i, whereas the eigenfunctions in Figure 5(b) are associated with the interfacial mode of the P family. The dashed curve in Figure 5(a) shows that the critical points (ci = 0) for the flow occur at y = ± 0.8732, and the imaginary part of the eigenfunction exhibits characteristic extrema in the neighborhood of the critical points. Another important feature is that, in the absence of viscosity contrast (m = 1), all the eigenfunctions in Figure 5 are parabolic and symmetric with respect to y = 0.

As we change the viscosity ratio, distinct features emerge. In Figures 6-7, the eigenfunctions, associated with the A, P, S1, and S2 modes, are plotted for m = 1/2 and m = 10, respectively. In Figures 6(a) and 7(a), we show that the eigenfunctions associated with the shear mode (located on the A branch) are parabolic but not symmetric. It is noted that the imaginary part of the eigenfunctions v(i) undergoes a rapid change in its derivative near the walls. At the interface, y = 0, there is also a dramatic change in the shapes and magnitudes of these functions for m = 10, as shown in Figure 7(a).

FIG. 6.

Orr-Sommerfeld eigenfunctions for stratified flow with Re = 10 000, wave numbers α = 1, m = 1/2, n = 1, and r = 1. (a) A-branch; (b) P-branch; (c) S1-branch; and (d) S2-branch. The solid and dashed lines represent the real and imaginary parts of the normal velocity.

FIG. 6.

Orr-Sommerfeld eigenfunctions for stratified flow with Re = 10 000, wave numbers α = 1, m = 1/2, n = 1, and r = 1. (a) A-branch; (b) P-branch; (c) S1-branch; and (d) S2-branch. The solid and dashed lines represent the real and imaginary parts of the normal velocity.

Close modal
FIG. 7.

Orr-Sommerfeld eigenfunctions for stratified flow with Re = 10 000, wave numbers α = 1, m = 10, n = 1, and r = 1. (a) A-branch; (b) P-branch; and (c) S1-branch; (d) S2-branch. The solid and dashed lines represent the real and imaginary parts of the normal velocity.

FIG. 7.

Orr-Sommerfeld eigenfunctions for stratified flow with Re = 10 000, wave numbers α = 1, m = 10, n = 1, and r = 1. (a) A-branch; (b) P-branch; and (c) S1-branch; (d) S2-branch. The solid and dashed lines represent the real and imaginary parts of the normal velocity.

Close modal

The eigenfunctions associated with the interfacial P mode have their maxima close to the center of the channel, but shifted to the left when m < 1 (Figure 6(b) for m = 1/2), and shifted to the right when m > 1 (Figure 7(b) for m = 10). Those associated with the S1 modes are highly damped in the less viscous layer (i.e., 0 ≤ y ≤ 1 for m = 1/2 and −1 ≤ y ≤ 0 for m = 10) as shown in Figures 7(c) and 6(c). Opposite trends are found for the eigenfunctions associated with the S2 modes (Figures 7(d) and 6(d)), which are highly damped in the layer which is more viscous, unlike the eigenfunctions of the S modes for m = 1, which exhibit translational symmetry as shown in Figure 5(c).

As we increase the viscosity ratio, the eigenfunctions associated with both the interfacial or the shear modes undergo a rapid change at the interface. This behavior is shown in Figure 8 for m = 2, m = 5 and m = 10, and m = 20. It is also noted that the eigenfunctions associated with the interfacial modes have their maxima close to the center of the channel and have much higher amplitudes for moderate values of the viscosity ratios. Whereas those associated with the shear modes have their largest variation close to the wall. Moreover, both types of eigenfunctions are damped in the lower layer and become centered around y = 1/2 for large viscosity ratios.

FIG. 8.

Orr-Sommerfeld eigenfunctions v(r) for stratified flow with Re = 10 000, wave numbers α = 1, n = 1, and r = 1. (a) Interfacial P-branch mode and (b) shear A-branch mode. Solid lines (m = 2), dashed-dashed lines (m = 5), dashed-dotted lines (m = 10), and dotted-dotted lines (m = 20).

FIG. 8.

Orr-Sommerfeld eigenfunctions v(r) for stratified flow with Re = 10 000, wave numbers α = 1, n = 1, and r = 1. (a) Interfacial P-branch mode and (b) shear A-branch mode. Solid lines (m = 2), dashed-dashed lines (m = 5), dashed-dotted lines (m = 10), and dotted-dotted lines (m = 20).

Close modal

In the limit m → ∞ and for large values of the viscosity ratio, the lower layer tends to become passive with regard to any shear mode that is excited in the upper fluid. Consequently, in this limit, the behavior of the eigenfunction corresponding to the unstable shear mode will be similar to that of the eigenfunction for a single fluid except it now will be centered around y = 1/2 as shown in Figure 9. It is also observed that the real and imaginary parts of eigenfunctions associated with the P mode have their extrema close to the center of the channel and shifted to the right for m > 1 as shown in Figure 7(b), and will be centered around y = 1/2 as m → ∞ (not shown here), instead of y = 0 as is the case of the P modes eigenfunctions for m = 1, n = 1, r = 1, and Re = 10 000 (Figure 5(b)).

FIG. 9.

Orr-Sommerfeld eigenfunctions for stratified flow of the least stable A-branch mode with Re = 10 000, wave numbers α = 1, n = 1, and r = 1. (a) Real part of the normal velocity and (b) imaginary part of the normal velocity. Solid lines (m = 103), dashed-dashed lines (m = 104), and dashed-dotted lines (m = 105).

FIG. 9.

Orr-Sommerfeld eigenfunctions for stratified flow of the least stable A-branch mode with Re = 10 000, wave numbers α = 1, n = 1, and r = 1. (a) Real part of the normal velocity and (b) imaginary part of the normal velocity. Solid lines (m = 103), dashed-dashed lines (m = 104), and dashed-dotted lines (m = 105).

Close modal

In Figure 10(a), we also plot the neutral stability curves for the shear mode in the plane of the Reynolds number and the wavenumber for the representative values of the viscosity ratio m = 1, m = 10, and m = 30. In the calculations, we have set n = 1, r = 1, and We = ∞. The curve for m = 1 is the well known result for plane Poiseuille flow in a channel. The critical values of the Reynolds number and the wavenumber are Rec = 5772 and αc = 1.02, respectively. The shapes of these neutral stability curves remain qualitatively the same as the viscosity ratio varies and the critical Reynolds numbers depend on the viscosity ratio. The neutral curves for m = 10 and m = 30 show different values of the Reynolds numbers and the wavenumbers. The critical Reynolds numbers for the shear mode are Rec = 10 228 and 4425.2, and the wavenumbers at criticality are αc = 1.274 and 1.453, respectively, for m = 10 and m = 30. Similar features are found for m = 1/10, m = 1/3, and m = 1/2 as shown in Figure 10(b).

FIG. 10.

Neutral stability curves for the least stable modes (a) m = 1, m = 10, and m = 30; (b) m = 1/10, m = 1/3, and m = 1/2. Flow parameters: We = ∞, r = 1, and n = 1.

FIG. 10.

Neutral stability curves for the least stable modes (a) m = 1, m = 10, and m = 30; (b) m = 1/10, m = 1/3, and m = 1/2. Flow parameters: We = ∞, r = 1, and n = 1.

Close modal

In Figure 11, we plot the real parts of the y-velocity disturbance function v(r) for both the shear mode and the interfacial mode at criticality. It is shown that the eigenfunctions associated with the interfacial modes have their maxima close to the center of the channel and have much higher amplitudes for moderate values of the viscosity ratios. Whereas those associated with the shear modes have their largest variation close to the wall. Moreover, both types of eigenfunctions damped in the lower layer become centered around y = 1/2 for large viscosity ratios. This means that, even at criticality, all these eigenfunctions have the same shapes and similar behaviors as shown previously in Figure 8.

FIG. 11.

Orr-Sommerfeld eigenfunctions v(r) for stratified flow at criticality. (a) Interfacial P-branch mode and (b) shear A-branch mode. Solid lines (m = 10, n = 1, r = 1, Rec = 10 228, We = ∞, and αc = 1.274), dashed-dashed lines (m = 30, n = 1, r = 1, Rec = 4425.2, We = ∞, and αc = 1.453).

FIG. 11.

Orr-Sommerfeld eigenfunctions v(r) for stratified flow at criticality. (a) Interfacial P-branch mode and (b) shear A-branch mode. Solid lines (m = 10, n = 1, r = 1, Rec = 10 228, We = ∞, and αc = 1.274), dashed-dashed lines (m = 30, n = 1, r = 1, Rec = 4425.2, We = ∞, and αc = 1.453).

Close modal

For certain parameter ranges, it is not possible to classify a mode as a shear mode or an interfacial mode. The problem arises when the modes occur close together and intersect. The aim of this section is to explore the general features of mode coalescence at low, and large Reynolds numbers for few representative values of the viscosity ratios, density ratios, and thickness ratios. We will discuss in detail several instances of mode coalescence between the shear mode and the interfacial mode which will take place either in the unstable part and/or in the stable part of the spectrum.

In Figure 12, we plot the temporal growth rates versus α for the unstable interfacial modes at low Reynolds number (Re = 10) with equal densities and thickness of layers (r = 1 and n = 1) for representatives values of the viscosity ratio. An increase of the viscosity ratio may have a non monotonic destabilizing effect on the growth rate for the choice of small values of viscosity ratios (m = 1/3, m = 1/2, m = 4/5). This effect become destabilizing at large values of viscosity ratio (m = 10, m = 30, m = 100). Moreover, no unstable shear mode is found and the interfacial modes become the dominant unstable modes at high viscosity ratio. Thus, to quantify the unstable wave patterns and to examine the interactions between the interfacial and shear modes at low Reynolds number, we consider specific parameters combinations with the choice of a small viscosity ratio (m = 1/20).

FIG. 12.

Temporal growth rates for disturbances of the unstable interfacial mode at low Reynolds R = 10 with We = ∞, n = 1, and r = 1. (a) m = 1/3, m = 1/2, and m = 3/4 and (b) m = 10, m = 30, and m = 100.

FIG. 12.

Temporal growth rates for disturbances of the unstable interfacial mode at low Reynolds R = 10 with We = ∞, n = 1, and r = 1. (a) m = 1/3, m = 1/2, and m = 3/4 and (b) m = 10, m = 30, and m = 100.

Close modal

One of the important features of the flow concerns the growth rate distributions and the coalescence in the range of most unstable waves. One example calculation is carried out to compute the growth rates for the unstable interfacial and the unstable shear modes based on the flow parameters for the experimental setup of Kao and Park27Re = 230, We = 63.81, Fr = 1.1974, m = 1/20, n = 1, and r = 1.16. Figure 13 shows the growth rates for both interfacial and shear modes. It is noted that no crossover between modes occurs and the maximum growth rate for the interfacial mode is about two orders of magnitude larger than that for the shear mode.

FIG. 13.

Temporal growth rates for disturbances. (a) Unstable interfacial mode and (b) unstable shear mode. Flow parameters of Kao and Park:27Re = 230, We = 63.811, Fr = 1.197, m = 1/20, n = 1, and r = 1.16.

FIG. 13.

Temporal growth rates for disturbances. (a) Unstable interfacial mode and (b) unstable shear mode. Flow parameters of Kao and Park:27Re = 230, We = 63.811, Fr = 1.197, m = 1/20, n = 1, and r = 1.16.

Close modal

For a specific choice of the viscosity ratio, the thickness ratio and the density ratio, the flow may become unstable to more than one unstable mode and coalescences between different modes may occur in the stable or the unstable region of the spectrum. These possibilities are displayed in Figure 14 where the growth rates are plotted against wave number α for few representative cases. The flow parameters are the same as those used in the calculations of Figure 14. It is worth emphasizing that the coalescence is very sensitive to the choice of the governing parameters. For the parameters’ combination m = 1/20, n = 1, r = 1.5, Re = 10, and We = ∞, a crossover between the interfacial mode and the shear mode takes place in the stable range of the wave number as shown in Figure 14(a). When either the density ratio or the thickness ratio is decreased slightly, only interfacial modes are unstable and no coalescence between unstable modes occurs as shown in Figure 14(b) for r = 0.9 and in Figure 14(c) for n = 0.8. On the other hand, there are instances when the growth rates exhibit complicated nonmonotonic behavior showing multiple crossings between the least stable modes. As illustrated in Figure 15, when the flow parameters are kept fixed (Re = 10, We = ∞, G = 0.6, m = 1/20, and n = 1) and the density ratio varies (up to 5 and 30 in this example), one can see that the coalescence occurs either in the stable region of the spectrum (Figure 15(a)) for a moderate density ratio of r = 5 or in the unstable region of the spectrum (Figure 15(b)) at higher density ratio (r = 30). Let us consider now the coalescence between modes for larger Reynolds numbers (Re = 4000) for specific parameters combinations. The results reveal there are instances when the flow becomes unstable to more than one unstable mode and coalescences between different modes may occur in the stable and/or in the unstable region of the spectrum. At small viscosity ratio m = 1/20 with n = 1 and r = 1, coalescence between modes may occur in the stable region of the spectrum (Figure 16(a)). For a greater density ratio (r = 2), the crossovers take place in the unstable region of the spectrum (Figure 16(b)).

FIG. 14.

Temporal growth rates for disturbances at Re = 10 with We = ∞, G = 0.6 and Fr = ((r − 1)/G)1/2. (a) m = 1/20, n = 1, r = 1.5; (b) m = 1/20, n = 1, r = 0.9; (b) m = 1/20, n = 0.8, r = 1.5. The solid and dashed lines are associated with the interfacial modes and the shear modes, respectively.

FIG. 14.

Temporal growth rates for disturbances at Re = 10 with We = ∞, G = 0.6 and Fr = ((r − 1)/G)1/2. (a) m = 1/20, n = 1, r = 1.5; (b) m = 1/20, n = 1, r = 0.9; (b) m = 1/20, n = 0.8, r = 1.5. The solid and dashed lines are associated with the interfacial modes and the shear modes, respectively.

Close modal
FIG. 15.

Temporal growth rates for disturbances. (a) r = 5; (b) r = 30. Flow parameters : Re = 10, We = ∞, Fr = ((r − 1)/G)1/2, m = 1/20 and n = 1.

FIG. 15.

Temporal growth rates for disturbances. (a) r = 5; (b) r = 30. Flow parameters : Re = 10, We = ∞, Fr = ((r − 1)/G)1/2, m = 1/20 and n = 1.

Close modal
FIG. 16.

Temporal growth rates for disturbances. (a) m = 1/20, n = 1, r = 1; (b) m = 1/20, n = 1, r = 2; (c) m = 10, n = 1, r = 1.5; (d) m = 100, n = 1, r = 1.5; (e) m = 10, n = 1, r = 10; and (f) m = 10, n = 0.7, r = 1. Flow parameters: Re = 4000, We = ∞, Fr = ((r − 1)/G)1/2, m = 1/20, and n = 1.

FIG. 16.

Temporal growth rates for disturbances. (a) m = 1/20, n = 1, r = 1; (b) m = 1/20, n = 1, r = 2; (c) m = 10, n = 1, r = 1.5; (d) m = 100, n = 1, r = 1.5; (e) m = 10, n = 1, r = 10; and (f) m = 10, n = 0.7, r = 1. Flow parameters: Re = 4000, We = ∞, Fr = ((r − 1)/G)1/2, m = 1/20, and n = 1.

Close modal

At higher viscosity ratio and using flow parameters of Figures 16–18 one sees that for m = 10, the curves of the shear and interfacial modes coalesce at a point located in the stable region of the spectrum. An increase of the viscosity ratio results in a coalescence between the most unstable waves occurs in the stable and the unstable regions of the spectrum as shown in Figure 16(c) and Figure 16(d), respectively, for m = 10 and m = 100. However, no coalescence occurs with an increase of the density ratio (r = 10) as shown in Figure 16(e). It is also noted that with a decrease of the thickness (n = 0.7), the crossing between the shear and interfacial modes is shifted in the stable region associated with higher wave numbers (Figure 16(f)).

FIG. 17.

Temporal growth rates for disturbances at Re = 10 000 and We = ∞. (a) m = 1/20, n = 1, r = 1.5; (b) m = 1/20, n = 0.8, r = 1.5; and (c) m = 1/20, n = 1, r = 0.9.

FIG. 17.

Temporal growth rates for disturbances at Re = 10 000 and We = ∞. (a) m = 1/20, n = 1, r = 1.5; (b) m = 1/20, n = 0.8, r = 1.5; and (c) m = 1/20, n = 1, r = 0.9.

Close modal
FIG. 18.

Temporal growth rates for disturbances at Re = 10 000 with We = ∞, G = 0.6 and Fr = ((r − 1)/G)1/2. (a) m = 10, n = 1, r = 1; (b) m = 50, n = 1, r = 1; (c) m = 10, n = 1, r = 20; (d) m = 10, n = 1/2, r=1.

FIG. 18.

Temporal growth rates for disturbances at Re = 10 000 with We = ∞, G = 0.6 and Fr = ((r − 1)/G)1/2. (a) m = 10, n = 1, r = 1; (b) m = 50, n = 1, r = 1; (c) m = 10, n = 1, r = 20; (d) m = 10, n = 1/2, r=1.

Close modal

In order to discuss the relevance of the stability analysis presented above, we repeat the stability analysis for different base states as well as compare our prediction of the most dangerous mode with the disturbance wavelengths observed in an experimental study reported in the literature.

The linear stability analysis presented above is carried out on the basis of the standard approach with a fully developed laminar velocity profile as the base state. In practice, this condition would exist downstream of the entrance where the two viscous boundary layers merge to form the fully developed velocity profile. The spatial location at which this occurs is a function of the flow rate and the viscosity ratio. For such cases, the stability analysis presented in this study would apply downstream of the point of fully developed flow. For problems where the domain length is insufficient, the velocity profile will be in the developing phase and will vary spatially. The stability analysis of such problems has not been investigated in the past. Here, we consider the linear stability behavior with respect to partially developed base states to obtain a basic understanding of the role of the base state profiles.

Figure 19(a) shows the velocity profiles for the partially developed two-fluid flow at various stages of development. These profiles are obtained with direct numerical simulation of 2-D two-phase flow.28 The location of these profiles downstream of the entrance depends on the Reynolds number. However, the profile itself is independent of the Reynolds number for a specific value of the viscosity ratio. In other words, the same profile would be obtained at different spatial locations for different values of Re. The profiles in Fig. 19(a) show various stages of development towards the fully developed state. The least developed profile is quite flat in the middle with steep wall gradients. The corresponding dispersion curves are shown in Fig. 19(b). The maximum growth rate can be observed to occur for the least developed profile U1. It decreases for U2 but then increases for U3 and Us profiles. The wavenumber with the highest growth rate also varies from 0.5 to about 3.5. Results in Fig. 19(b) indicate that both the maximum growth rate and the corresponding most dangerous modes depend on the base state. Hence, the analysis based on the fully developed base state should be most valid in practice when the domain is sufficiently long so as to allow the flow to achieve a fully developed state. Characterization of the location where the fully developed two-fluid flow downstream of the entrance would exist, which needs to be quantified in future studies as a function of the viscosity ratio and the layer thickness.

FIG. 19.

(la) Developing and developed velocity profiles. (b) Growth rate, σ, versus the wavenumber, α, for the most unstable mode for m = 10, n = 1, r = 1 and Re = 10.

FIG. 19.

(la) Developing and developed velocity profiles. (b) Growth rate, σ, versus the wavenumber, α, for the most unstable mode for m = 10, n = 1, r = 1 and Re = 10.

Close modal

One of the important contemporary applications of two-fluid flow stability analysis is related to the cooling of electronic devices via two-phase heat transfer in thin gap channels.29,30 This is a promising technology for removing large amounts of heat by the evaporation of liquid layers. The deformations and emerging patterns observed at the liquid vapor interface are of particular interest as they influence temperature fluctuations within the layer, which in turn govern the heat flux at the heated wall and the rate of evaporation at the interface. The present work takes the first step in examining the structure of the waves and instabilities that occur in adiabatic two-phase flows in such systems. For adiabatic flows, it has been observed that the fully developed nonlinear state is stable and does not transition to turbulence. This is because of relatively small liquid Reynolds numbers and finite channel lengths in such applications.

We compare the predictions of the most dangerous mode from our analysis with the experimental observation of wave patterns in the experiments of Kabov et al.31 We have computed growth rates for the interfacial modes and the shear modes based on the flow parameters used by the authors. The wavelength obtained from the linear stability analysis, expressed in dimensional form, is given by

λ ̄ = 2 π d 1 α c ,

where αc is the most unstable wavenumber related to the maximum growth rate, σm.

Figure 20 shows the growth rate σ versus the wavenumber α for the most unstable mode. The predicted wavelength, αc = 1.9 cm, appears to agree well with the experimental observation of the measured wavelength of interfacial waves, suggesting that the linear stability approach accurately captures the wavelength of interfacial waves. Note that in the absence of the actual flow profile, we have used a fully developed base state for this comparison. The actual velocity profile may be different due to nonuniform entrance flow and spanwise interfacial features. The relatively good agreement thus needs to be further substantiated when more experimental results become available. In a separate study,32 we have carried out a more detailed comparison with experimental observations of two-phase annular flow.

FIG. 20.

Growth rate σ versus the wavenumber α for the most unstable mode for the experiment of Kabov et al.31 Flow parameters: m = 3.57; n = 1; r = 124.65; We = 1.26, and Fr = 10.

FIG. 20.

Growth rate σ versus the wavenumber α for the most unstable mode for the experiment of Kabov et al.31 Flow parameters: m = 3.57; n = 1; r = 124.65; We = 1.26, and Fr = 10.

Close modal

In this study, a linear stability analysis of two-phase channel flow is performed. The linearized system is solved numerically using a multi-domain Chebychev collocation spectral method. The algorithm has proven efficient in accurately computing the full eigenspectrum. Numerical solutions are validated against asymptotic results in the long and short wave limits. The spectrum and eigenfunction profiles associated with A, P, and S modes are described. The least stable modes located on the A and P branches are identified. These are shown to be associated with the shear and interfacial modes, respectively. The effects of viscosity ratio and Reynolds number on the behavior of the two instability modes are investigated. For m > 10, an unstable interfacial mode appears, which is shown to be asymptotically independent of the Reynolds number. The asymptotic value depends on the viscosity ratio.

Observations of the eigenfunction profiles reveal sharp variations around the interface for large values of the viscosity ratio. The profiles however transform progressively to parabolic shapes when the viscosity ratio is increased. The maximum amplitude of the eigenfunction occurs close to the center of the channel, shifted to the right or to the left depending upon whether the viscosity ratio is greater or less than one, respectively. The maximum amplitude asymptotes to the channel center for large viscosity ratios. A parametric study is performed to understand the mode coalescence phenomenon. It is shown that for certain parameter combinations, coalescence may occur in both the stable and unstable flow regimes.

This study was supported by the Office of Naval Research, Grant No. N00014-11-1-0468 monitored by Dr. Mark Spector.

The functions h1(m, n, y) and h2(m, n, y) in (22) are given by

h 1 = ( y 4 p 1 ( m , n , y ) ) / ( 840 n 3 ( n + 1 ) 2 ( m 2 + 4 m n 3 + 6 m n 2 + 4 m n + n 4 ) ) ,

with

p 1 = d 0 + d 1 y + d 2 y 2 + d 3 y 3 ,
d 0 = 70 n 2 ( n + 1 ) ( m 3 + 2 m 2 n 3 + 8 m 2 n 2 + 5 m 2 n + 2 m 2 + 2 m n 5 + 5 m n 4 + 8 m n 3 + 2 m n 2 n 5 ) ,
d 1 = 7 ( m 4 + 8 m 3 n 3 + 7 3 n 2 + 3 m 3 n 8 m 2 n 6 46 m 2 n 5 57 m 2 n 4 35 m 2 n 3 10 m 2 n 2 8 m n 8 24 m n 7 43 m n 6 23 m n 5 2 n 9 4 n 8 9 n 7 6 n 6 ) ,
d 2 = 7 ( m n 2 ) 2 ( m 2 + 4 m n 3 + 6 m n 2 + 4 m n + n 4 ) ,
d 3 = 2 ( m n 2 ) ( m + n ) ( m 2 + 4 m n 3 + 6 m n 2 + 4 m n + n 4 ) ,
h 2 = ( y 4 p 2 ( m , n , y ) ) ( 840 m 2 n 3 ( n + 1 ) 2 ( m 2 + 4 m n 3 + 6 m n 2 + 4 m n + n 4 ) ) ,

with

p 2 = e 0 + e 1 y + e 2 y 2 + e 3 y 3 ,
e 0 = 70 m 2 ( n + 1 ) ( m 3 2 m 2 n 3 8 m 2 n 2 5 m 2 n 2 m 2 2 m n 5 5 m n 4 8 m n 3 2 m n 2 + n 5 ) ,
e 1 = 7 ( 6 m 4 n 3 + 9 m 4 n 2 + 4 m 4 n + 2 m 4 + 23 m 3 n 4 + 43 m 3 n 3 + 24 m 3 n 2 + 8 m 3 n + 10 m 2 n 7 + 35 m 2 n 6 + 57 m 2 n 5 + 46 m 2 n 4 + 8 m 2 n 3 3 m n 8 7 m n 7 8 m n 6 n 9 ) ,
e 2 = 7 ( n m 2 ) 2 ( m 2 + 4 m n 3 + 6 m n 2 + 4 m n + n 4 ) ,
e 3 = 2 ( m n 2 ) ( m + n ) ( m 2 + 4 m n 3 + 6 m n 2 + 4 m n + n 4 ) .

The function H1 in (22) is given by

H 1 = k = 1 k = 6 w k / q k ,

with

w 1 = n 5 ( n + 1 ) 4 ( 16 m + 132 m n + 525 m n 2 + 1350 m n 3 + 2493 m n 4 + 3384 m n 5 + 3312 m n 6 + 2112 m n 7 + 768 m n 8 + 13 n 4 r + 54 n 5 r + 189 n 6 r + 204 n 7 r + 192 n 8 r + 4 n 3 + 27 n 4 + 90 n 5 + 171 n 6 + 288 n 7 + 240 n 8 + 192 n 9 + 315 m n 2 r + 1098 m n 3 r + 2187 m n 4 r + 2700 m n 5 r + 1968 m n 6 r + 768 m n 7 r + 40 m n r ) ,
w 2 = n 4 F / 24 m r / 21 4 n r / 21 m n / 240 m n 2 / 15 463 m n 3 / 1680 43 m n 4 / 84 58 m n 5 / 105 83 n 2 r / 420 53 n 3 r / 560 + 13 n 4 r / 21 + 629 n 5 r / 336 + 317 n 6 r / 105 + 8 n 7 r / 5 + n 2 / 35 + 179 n 3 / 560 + 191 n 4 / 140 + 5653 n 5 / 1680 + 527 n 6 / 105 + 183 n 7 / 35 + 236 n 8 / 105 11 m n 2 r / 420 71 m n 3 r / 1680 m n 4 r / 5 + 37 m n r / 1680 + m n 3 F / 24 ,
w 3 = r ( n 4 + 20 n 3 + 100 n 2 + 246 n + 160 ) ,
w 4 = n 3 ( n + 1 ) 2 ( 28 n 19 m 179 m n 691 m n 2 1583 m n 3 2419 m n 4 2312 m n 5 1280 m n 6 + 53 n 2 r + 603 n 3 r + 2103 n 4 r + 4471 n 5 r + 5355 n 6 r + 4116 n 7 r + 1344 n 8 r + 245 n 2 + 1012 n 3 + 2656 n 4 + 4898 n 5 + 6832 n 6 + 6484 n 7 + 4368 n 8 + 1344 n 9 346 m n 2 r 1082 m n 3 r 1457 m n 4 r 1008 m n 5 r 48 m n r ) ,
w 5 = 20 n m 6 m n + 100 n 2 + 246 n 3 + 160 n 4 + 1 , w 6 = n 3 r ( n + 6 ) ,
q 1 = 15 ( m ( 4 n 3 + 6 n 2 + 4 n ) + m 2 + n 4 ) 3 , q 2 = m ( 4 n 3 + 6 n 2 + 4 n ) m 2 n 4 ,
q 3 = 3360 m ( n + 1 ) 2 , q 4 = 105 ( m ( 4 n 3 + 6 n 2 + 4 n ) + m 2 + n 4 ) 2 ,
q 5 = 3360 ( n + 1 ) 2 , q 6 = 3360 m 2 ( n + 1 ) 2 .
1.
P.
Drazin
,
Introduction to Hydrodynamic Stability
(
Cambridge University Press
,
2002
).
2.
D. D.
Joseph
and
Y. Y.
Renardy
, “
Fundamentals of two-fluid dynamics. Pt. II: Lubricated transport, drops and miscible liquids
,” in
Fundamentals of Two-Fluid Dynamics. Pt. II: Lubricated Transport, Drops and Miscible Liquids
,
Interdisciplinary Applied Mathematics
Vol.
4
(
Springer-Verlag
,
1993
), p.
459
.
3.
J. W.
Miles
, “
The hydrodynamic stability of a thin film of liquid in uniform shearing motion
,”
J. Fluid Mech.
8
,
593
610
(
1960
).
4.
C.-S.
Yih
, “
Stability of liquid flow down an inclined plane
,”
Phys. Fluids
6
,
321
334
(
1963
).
5.
C.-S.
Yih
, “
Instability due to viscosity stratification
,”
J. Fluid Mech.
27
,
337
352
(
1967
).
6.
M.
Charles
and
L.
Lilleleht
, “
An experimental investigation of stability and interfacial waves in co-current flow of two liquids
,”
J. Fluid Mech.
22
,
217
224
(
1965
).
7.
L. S.
Cohen
and
T. J.
Hanratty
, “
Generation of waves in the concurrent flow of air and a liquid
,”
AIChE J.
11
,
138
144
(
1965
).
8.
A. D.
Craik
, “
Wind-generated waves in thin liquid films
,”
J. Fluid Mech.
26
,
369
392
(
1966
).
9.
P. J.
Schmid
and
D. S.
Henningson
,
Stability and Transition in Shear Flows
(
Springer Science & Business Media
,
2001
), Vol.
142
.
10.
S. G.
Yiantsios
and
B. G.
Higgins
, “
Linear stability of plane poiseuille flow of two superposed fluids
,”
Physics of Fluids
31
,
3225
(
1988
).
11.
A.
Hooper
, “
Long-wave instability at the interface between two viscous fluids: Thin layer effects
,”
Physics of Fluids (1958-1988)
28
,
1613
1618
(
1985
).
12.
S.
Özgen
,
M.
Carbonaro
, and
G. S. R.
Sarma
, “
Experimental study of wave characteristics on a thin layer of de/anti-icing fluid
,”
Phys. Fluids
14
(
10
),
3391
3402
(
2002
).
13.
A.
Hooper
and
W.
Boyd
, “
Shear-flow instability due to a wall and a viscosity discontinuity at the interface
,”
J. Fluid Mech.
179
,
201
225
(
1987
).
14.
S.
Timoshin
, “
Instabilities in a high-reynolds-number boundary layer on a film-coated surface
,”
J. Fluid Mech.
353
,
163
195
(
1997
).
15.
S.
Timoshin
and
A.
Hooper
, “
Mode coalescence in a two-fluid boundary-layer stability problem
,”
Phys. Fluids
12
,
1969
1978
(
2000
).
16.
S.
Özgen
, “
Coalescence of Tollmien–Schlichting and interfacial modes of instability in two-fluid flows
,”
Phys. Fluids
20
,
044108
(
2008
).
17.
H.
Squire
, “
On the stability for three-dimensional disturbances of viscous fluid flow between parallel walls
,”
Proc. R. Soc. London, Ser. A
142
,
621
628
(
1933
).
18.
L. N.
Trefethen
,
Spectral Methods in MATLAB
(
SIAM
,
2000
), Vol.
10
.
19.
A.
Hooper
and
W.
Boyd
, “
Shear-flow instability at the interface between two viscous fluids
,”
J. Fluid Mech.
128
,
507
528
(
1983
).
20.
P.
Boomkamp
and
R.
Miesen
, “
Classification of instabilities in parallel two-phase flow
,”
Int. J. Multiphase Flow
22
,
67
88
(
1996
).
21.
P.
Boomkamp
and
R.
Miesen
, “
Nonaxisymmetric waves in core-annular flow with a small viscosity ratio
,”
Phys. Fluids A
4
,
1627
1636
(
1992
).
22.
S. A.
Orszag
, “
Accurate solution of the Orr–Sommerfeld stability equation
,”
J. Fluid Mech.
50
,
689
703
(
1971
).
23.
J.
Dongarra
,
B.
Straughan
, and
D.
Walker
, “
Chebyshev tau-qz algorithm methods for calculating spectra of hydrodynamic stability problems
,”
Appl. Numer. Math.
22
,
399
434
(
1996
).
24.
L. M.
Mack
, “
A numerical study of the temporal eigenvalue spectrum of the blasius boundary layer
,”
J. Fluid Mech.
73
,
497
520
(
1976
).
25.
A.
Kaffel
and
M.
Renardy
, “
On the stability of plane parallel viscoelastic shear flows in the limit of infinite weissenberg and Reynolds numbers
,”
J. Non-Newtonian Fluid Mech.
165
,
1670
1676
(
2010
).
26.
A.
Kaffel
and
M.
Renardy
, “
Surface modes in inviscid free surface shear flows
,”
ZAMM-J. Appl. Math. Mech./Z. Angew. Math. Mech.
91
,
649
652
(
2011
).
27.
T. W.
Kao
and
C.
Park
, “
Experimental investigations of the stability of channel flows. Part 2. Two-layered co-current flow in a rectangular channel
,”
J. Fluid Mech.
52
,
401
423
(
1972
).
28.
Z.
Qin
,
K.
Delaney
,
A.
Riaz
, and
E.
Balaras
, “
Topology preserving advection of implicit interfaces on Cartesian grids
,”
J. Comput. Phys.
290
,
219
(
2015
).
29.
A.
Bar-Cohen
and
E.
Rahim
, “
Modeling and prediction of two-phase microgap channel heat transfer characteristics
,”
Heat Transfer Eng.
30
,
601
625
(
2009
).
30.
A.
Bar-Cohen
,
C.
Hollaway
, and
D.
Sharar
, “
Liquid film wave patterns and dryout in microgap channel annular flow
,” in
Electronic Equipment Cooling, Proceedings of the 15th International Heat Transfer Conference, Kyoto, Japan, August 10–15
(
Begell House
,
2014
).
31.
O.
Kabov
,
D.
Zaitsev
,
V.
Cheverda
, and
A.
Bar-Cohen
, “
Evaporation and flow dynamics of thin, shear-driven liquid films in microgap channels
,”
Exp. Therm. Fluid Sci.
35
,
825
831
(
2011
).
32.
A.
Bar-Cohen
,
C. A.
Holloway
,
A.
Kaffel
, and
A.
Riaz
, “
Waves and instabilities in high quality adiabatic flow in microgap channels
,”
Int. J. Multi Phase Flow
(submitted).