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.

## I. INTRODUCTION

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. Miles^{3} 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 Lilleleht^{6} 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 Hanratty^{7} determined the critical value of the gas velocity at which the waves first start to appear. Craik^{8} 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 Yih^{5} 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 stability^{1,9} in both of the fluids. Yiantsios and Higgins^{10} carried out a numerical study of the linear stability of plane Poiseuille flow of two superposed fluids and confirmed previous findings. The studies by Yih^{4,5} and Hooper^{11} 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 Boyd^{13} 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. Timoshin^{14} considered the same problem using the triple-deck theory and recovered the interfacial modes predicted by Hooper and Boyd.^{13} Timoshin and Hooper^{15} 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. Ozgen^{16} 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.

## II. GOVERNING EQUATIONS AND BASIC FLOW

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 *d*_{1} while the lower layer (liquid phase) has viscosity *μ*_{2}, density *ρ*_{2}, and thickness *d*_{2}.

Squire’s transformation^{17} for the parallel flow of a uniform fluid and later extensions by Yih^{4,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.

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

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

where *g* is the gravitational acceleration and *γ* is the interfacial tension. The velocity at *Y* = 0 is denoted by *U*_{0}.

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

In the lower liquid layer, the equations read

where **u**_{1} and **u**_{2} 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 *U*_{0}, time with *d*_{1}/*U*_{0}, and pressure with $ \rho 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 *U*_{0}, are given by

where, *x* = *X*/*d*_{1} and *y* = *Y*/*d*_{1} and

## III. THE LINEARIZED STABILITY PROBLEM

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

where *U _{i}* and

*P*are the dimensionless velocities and pressure of the basic flow, respectively. Subscript

_{i}*i*denotes the phase and primes denote perturbations. The linearized problem in the

*x*-

*y*plane can be expressed in terms of the streamfunctions

*ψ*, where

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

where *ϕ _{i}*(

*y*) and

*f*(

_{i}*y*) are the disturbance amplitudes,

*α*is the dimensionless wavenumber defined by 2

*πd*

_{1}/λ, with wavelength λ,

*c*=

*c*+ i

_{r}*c*is the dimensionless wave velocity, and the most rapidly growing wave is the one for which the rate of amplification

_{I}*αc*is a maximum.

_{I}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:

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

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

where $ c \u0303 =c\u2212 U 1 ( 0 ) $ and $F= ( r \u2212 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= \varphi 1 ( 0 ) / c \u0303 $ 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 *c _{i}* > 0. Neutral stability is associated with

*c*= 0. The case

_{i}*r*= 1,

*n*= 1, corresponds to a single fluid system, which was treated by Yih

^{5}and Higgins and Yiantsios.

^{10}

## IV. ASYMPTOTIC ANALYSIS FOR LONG AND SHORT WAVES

### A. Instability in the long wave limit

The asymptotic analysis of the eigenvalue for long waves was carried out by Yih^{5} 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 Yiantsios^{10} analysis. In long-wavelength asymptotics, one expands the eigenfunctions and the eigenvalue as a regular perturbation series in *α*,

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 *c _{i}* 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

At the leading order, the boundary conditions satisfy

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

Without loss of generality, we have set the arbitrary constants to unity, i.e., *A*_{10} = *A*_{20} = 1. At this order of approximation, the eigenvalue *c*_{0} is real and may be expressed in terms of the thickness ratio *n* and viscosity ratio *m*,

Note that this equation is independent of *r* and the Eqs. (19)-(21) are identical to Yih’s^{5} 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

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

The functions *h*_{1}(*m*, *n*, *y*) and *h*_{2}(*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).

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:

Hence, (29) can be written further as

Thus,

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

in which *c*_{1} is the change in *c* and is of course identical to the change in $ c \u0303 =c\u22121$. Hence,

By solving Eqs. (24)-(28) and (32), the six coefficients *B*_{11}, *C*_{11}, *D*_{11}, *B*_{21}, *C*_{21}, and *D*_{21} can be found. Then *c*_{1} is determined from (34). The result is

in which *H*_{1} 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,

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 *c _{I}* vanishes as

*m*− 1 as

*m*approaches 1. Thus, it is interesting to find that

*c*vanishes as either

_{I}*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.

### B. Instability in the short wave limit

The asymptotic analysis in the short wave limit proceeds along the same lines as the analysis of Yiantsios and Higgins^{10} for the plane Poiseuille flow and Hooper and Boyd^{13} for the plane Couette flow. The flow coordinates (*x*, *y*) are scaled with *α*, i.e.,

and the following asymptotic expansions are assumed for *ϕ*_{1}, *ϕ*_{2}, and *c*:

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

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

This result shows that short waves are unstable unless *r* > *m*^{2} > 1 or 1 > *m*^{2} > *r*.

## V. NUMERICAL SOLUTION

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,

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

where *T _{n}*(

*z*) are the Chebyshev polynomials,

*i*denotes the

*i*th derivative with respect to

*z*, and

*a*and

_{n}*b*are the discrete Chebyshev expansion coefficients.

_{n}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 Trefethen^{18} 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

where the (*N*_{1} + *N*_{2} + 2) components’ vector Ψ is given by

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, *N*_{1}th, *N*_{1} + 1th, *N*_{1} + 2th, *N*_{1} + 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 *N*_{1} + *N*_{2} + 2 matrix eigenvalue problem which is solved using the QZ algorithm.

## VI. STABILITY RESULTS

### A. Spectra and eigenfunctions

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.

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 *c _{I}* 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}

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 Dongarra^{23} 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 Yih^{5} 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 *c _{n}* for nonzero

*α*. The eigenvalues are located on three main branches which have been labeled

*A*(

*c*→ 0),

_{r}*P*(

*c*→ 1), and

_{r}*S*(

*c*→ 2/3) by Mack.

_{r}^{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 6

*i*is located on the A branch and is associated with the shear mode.

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 *S*_{1} and *S*_{2}) and two new branches *P*_{1} and *P*_{2} appear to intersect the real axis at values close to the interface velocity (*c _{r}* ≈ 1). Branches

*S*

_{1}and

*S*

_{2}converge to the branch

*S*as

*m*→ 1; and the branches

*A*

_{1},

*A*

_{2},

*P*

_{1}, and

*P*

_{2}likewise converge to the

*A*and

*P*branches.

For *m* = 1/2, the shear mode *c* = 0.1713 − 0.000 49*i* 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^{−9}*i* is located on the *P* branch. The highly damped *S*_{1} and *S*_{2} modes occur with phase speeds *c _{r}* → 0.64 and

*c*→ 0.74 for

_{r}*m*= 1/2 and with phase speeds

*c*→ 0.61 and

_{r}*c*→ 1.42 for

_{r}*m*= 10.

For the viscosity ratio *m* = 10, another unstable mode *c* = 1.453 + 0.2205*i* appears in the upper plane *c _{i}* > 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 (

*c*= 1.454) at any Reynolds number

_{r}*Re*> 13 000 and the imaginary part of the eigenvalue

*c*→ 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.

_{i}^{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.

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 Henningson^{9} and have their maxima close to *y* = ± 1/2 (Figure 5(c)).

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 6*i*, 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 (*c _{i}* = 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*, *S*_{1}, and *S*_{2} 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).

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 *S*_{1} 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 *S*_{2} 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.

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)).

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 *Re _{c}* = 5772 and

*α*= 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

_{c}*m*= 10 and

*m*= 30 show different values of the Reynolds numbers and the wavenumbers. The critical Reynolds numbers for the shear mode are

*Re*= 10 228 and 4425.2, and the wavenumbers at criticality are

_{c}*α*= 1.274 and 1.453, respectively, for

_{c}*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).

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.

### B. Coalescence

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).

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 Park^{27} *Re* = 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.

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)).

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)).

## VII. DISCUSSION

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.

### A. Influence of the base state

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 *U*_{1}. It decreases for *U*_{2} but then increases for *U*_{3} and *U _{s}* 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.

### B. Comparison with experiments

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

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.

## VIII. SUMMARY AND CONCLUSIONS

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.

## Acknowledgments

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

### APPENDIX: FUNCTIONS *h*_{1} AND *h*_{2}

The functions *h*_{1}(*m*, *n*, *y*) and *h*_{2}(*m*, *n*, *y*) in (22) are given by

with

with

The function *H*_{1} in (22) is given by

with