## Abstract

The current work deals with the numerical analysis of linear stability problems in a stratified plain Poiseuille flow of air over water with equal layer heights. The interaction and branch exchange between Tollmien–Schlichting instability in air and interfacial instability is discovered and investigated. This effect is shown to stabilize disturbances with wavelengths of the order of channel height for interfacial waves and to produce a closed stable region inside the neutral curve of the interfacial mode. The behaviour of three unstable modes in the problem, corresponding to Tollmien–Schlichting type instability in air and water layers and interfacial instability respectively, has been studied in detail. Neutral conditions for all three modes and the stable region have been calculated.

## 1. Introduction

This paper revisits the classical linear stability problem for a two-fluid Poiseuille flow in a planar channel. The problem has been approached by many authors from different angles. A long-wave asymptotic analysis performed by Yih (1967) illustrates the origin of interfacial instability in the flow. In general, interfacial modes can coexist with instability of a Tollmien–Schlichting variety, which is more common for a channel flow. This was studied computationally by Blennerhassett (1980) where the focus was on linear and nonlinear waves in a two-layer flow of air over water with equal layer depths. Yiantsios & Higgins (1988) performed accurate computations for linear instabilities in a Poiseuille flow with two superimposed fluids. Their work was centred on the behaviour of the interfacial wave mode and its dependence on the relative depths and kinematic viscosities of fluid layers. Hooper (1989) studied instabilities in a Poiseuille–Couette channel flow of two viscous fluids. An extensive review of the stability problem for a two-layer plane channel flow can be found in Joseph & Renardy (1993).

The formulation in this paper follows the early work by Blennerhassett (1980). We consider a two layer air–water configuration, with the majority of the results obtained for the case of equal layer depths. The focus will be on two issues. First, in a search for novel instability modes, we extend the results of Blennerhassett to higher Reynolds numbers. The second, and perhaps more interesting issue, stems from the following remarkable property. Blennerhassett's neutral stability curves in fig. 2 of his paper (Blennerhassett 1980), reproduced in our figure 1*a* (where *k* and *Re* denote the wavenumber and Reynolds number) from our computations, clearly show a broadband interfacial instability in the flow starting from the Reynolds number 3848.24.

In addition, a narrow band of waves marked by a neutral curve starting at the Reynolds number of 8641.96, indicates the appearance of a second instability mode. Blennerhassett attributed this second instability to a Tollmien–Schlichting mechanism. However, we found that the narrow neutral curve in figure 1*a* was in fact a loop closed down at Reynolds numbers of order 20 000, see figure 1*b*. As we shall show in this paper, the region inside the loop corresponds to decaying waves, contrary to what might be inspected in the case of Tollmien–Schlichting instability. From a physical standpoint, of major significance is the fact (entirely omitted in Blennerhassett's computations) that neutral-wave solutions for the points on the loop can be connected through a wavenumber continuation with the solutions on the broad-band interfacial neutral curve. This observation (illustrated in detail below) raises the question of physical identity for various instability modes in the flow. We will show that the neutral curves in figure 1 (and further neutral curves discovered at higher Reynolds numbers) are a product of a complicated interaction and intersection of several key modes associated with three more or less distinct instability mechanisms, namely the interfacial instability and two Tollmien–Schlichting instabilities, for the flow in the air and water layers. In particular, Tollmien–Schlichting instabilities are found to arise as two initially distinct modes at Reynolds numbers of order of 15 000, which is beyond the Reynolds number range covered in Blennerhassett (1980).

It is not uncommon for a neutral curve of a layered flow to have several branches. Yiantsios & Higgins (1988) were probably the first to note that the interfacial waves can be stabilized in a certain range of wavenumbers. Later, the same effect was noted by Kuru *et al*. (1995), Miesen & Boersma (1995) and Sangalli *et al*. (1997). It appears that under certain conditions the curve of the growth rate as a function of the wavenumber, *c*_{i} (*k*), for the interfacial mode can develop a local minimum that produces a stable region with *c*_{i} (*k*)<0. The formation of this stable region can also be seen in the study of interfacial waves behaviour in oil–water channel flows presented by McCready *et al*. (1997). Another reason for the appearance of complicated neutral curve diagrams in two-fluid flow is a mode crossing between interfacial and Tollmien–Schlichting waves considered for a boundary layer with high Reynolds numbers in Ozgen *et al*. (1998) and Timoshin & Hooper (2000). In our study, we find both the interfacial mode splitting and interfacial/Tollmien–Schlichting wave interaction are present in certain ranges of the wavenumber. In what follows we refer to the modes of interfacial instability, Tollmien–Schlichting instability due to air flow and Tollmien–Schlichting instability in water (liquid) as the I-mode, TSA-mode and TSL-mode, respectively.

## 2. Problem formulation

Assuming a flow in a horizontal channel, we use the upper (air) layer density and viscosity as reference quantities. The velocity is nondimensionalized on the integral mean value, *U*, of the undisturbed velocity across both layers in order to keep the formulation identical with that of Blennerhassett (1980). The coordinates are nondimensionalized on the air layer thickness *L*, with the origin at the undisturbed interface. For the majority of this study the base water layer is assumed to have the same thickness as the layer of undisturbed air. The media properties are given in table 1. Surface tension coefficient, *S*, was taken to be 0.0735 N m^{−1} to allow comparison with Blennerhassett (1980). The flow is laminar and entirely two-dimensional.

Subscripts ‘a’ and ‘l’ are used to denote air and liquid, respectively. Ratios of viscosity coefficients and densities are defined as *μ*=*μ*_{l}/*μ*_{a} and *ρ*=*ρ*_{l}/*ρ*_{a}, correspondingly.

Double subscripts denote dimensionless ratios of the corresponding values in the two layers, for example, the ratio of the viscosity coefficients is denoted as *μ*_{la}=*μ*_{l}/*μ*_{a}. Superscript ‘o’ is used for the base flow solution.

The horizontal velocity component in the base flow can be taken in the form(2.1)with the constant parameters defined as(2.2)

The linear stability problem then reduces to a pair of Orr–Sommerfeld equations written in terms of the perturbation in the vertical velocity *v*_{a,l}(*y*)exp(i*k*(*x*−*ct*)) (in each layer for a wavenumber *k* and phase speed *c*) in the following form:(2.3)

The required boundary conditions follow from the velocity and stress continuity at the interface *y*=0. Tangential stress continuity yields(2.4)

Normal stress condition is given by(2.5)

Also, the continuities of normal and tangential velocity leads to(2.6)and(2.7)respectively. At the channel walls (*y*=±1) no-slip conditions are used:(2.8)

Here, the governing parameters are defined as *We*=ρ_{l}*U*^{2}*L*/*S*, *Fr*=*gL*/*U*^{2}, *Re*=*ρ*_{a}*UL*/*μ*_{a}, for the Weber number, Froude number and Reynolds number, respectively. The dimensional surface tension coefficient is *S* and *g* stands for the acceleration due to gravity. The amplitude of the interface perturbation, *f*, derived from the kinematic condition at the interface, is given by

Then, following Blennerhassett (1980), we note that *Fr*, *Re* and *We* are not independent for fixed fluids as *F*=*FrRe*^{2}=*gL*^{3}/*ν*_{a}^{2} and *t*=*Re*^{2}/*We*=*SL*/(*ρ*_{l}*ν*_{a}^{2}) depend on *L* only. Hence, we introduce a dimensionless parameter,

representing the non-dimensional wavenumber of the inviscid surface wave with the minimum phase speed (see Blennerhassett 1980). All the computations presented in this paper have been performed for *k*_{s}=29.18 to allow comparisons with the results in Blennerhassett (1980).

## 3. Numerical scheme and solution method

We used a compact finite-difference scheme of the fourth order of accuracy derived on the principles described by Hirsh (1975) and Li (1998). The Orr–Sommerfeld equations are written in the form(3.1)

Let *v*^{II}=d^{2}*v*/d*y*^{2} and *v*^{IV}=d^{4}*v*/d*y*^{4}. The following relations are true to *O*(Δ*y*^{4}) on a uniform grid:(3.2)

Substituting expression for *v*^{IV} from Orr–Sommerfeld equations into these difference equalities we obtain a system of linear equations for discretized variables. Let the computational grid contain *N*_{l} points in water layer and *N*_{a} points in air layer. Then the combined linear system (with boundary conditions included) will contain 2(*N*_{l}+*N*_{a}) equations. Orr–Sommerfeld equations in water and air are written at points *n*=1 … *N*_{l}−2, *N*_{l}+1, … *N*_{a}+*N*_{l}−1 (all indices start from 0), boundary conditions at the upper and lower walls are written at points *n*=0 and *n*=*N*_{a}+*N*_{l}−1 and boundary conditions at the air–water interface are written at points *n*=*N*_{l}−1 and *n*=*N*_{l}, respectively.

In order to maintain fourth order accuracy, derivatives in boundary conditions are discretized with fourth order accuracy using a five point one-side approximation. The zero perturbation condition for normal and tangential velocity at the upper wall are then given by(3.3)and(3.4)respectively.

On the air–water boundary, normal and tangential velocity conditions are given by(3.5)and(3.6)respectively. The tangential stress condition becomes(3.7)

Finally, the normal stress continuity condition is given by(3.8)

Also, velocity perturbations must vanish at the lower channel wall; therefore(3.9)and(3.10)

Thus, we obtain a system of linear equations with a system matrix that differs from a tridiagonal one only in lines containing boundary conditions. To build the dispersion relation we use the algorithm described applied by Thomas (1953) and Blennerhassett (1980). Gaussian elimination with partial pivoting is performed (taking into account the matrix structure both in storage and algorithm) in order to bring the system matrix to the block-triangular form so that in the matrix of the linear system, *a*_{i,i} is the only non-zero element in the *i*th row and *i*th column. Hence, for the system to have a non-zero solution, *a*_{i,i} must be equal to zero. Therefore, we can regard *a*_{i,i} as a function of (*k*, *Re*, *c*) and apply an iterative method (we used a modified Newton–Raffson method with regularization multiplier described, for example, in Kogan *et al*. (2000) to find roots of *a*_{i,i}(*c*)=0 for given (*k*, *Re*). The code built based on the proposed scheme was tested against results published by Orzag (1971), Blennerhassett (1980) and Yiantsios & Higgins (1988) and demonstrated good accuracy and convergence characteristics. In our opinion, the benefits of compact differencing combined with the generally available computational power practically eliminate the need for alternative (e.g. spectral) methods in multi-layer stability calculations. A detailed description of the scheme derivation and numerical tests can be found in the Ph.D. thesis by Shapiro (2004).

## 4. Results

Interfacial instability gives rise to the broadband neutral curve shown in figure 1*a*. Our results for this neutral curve are in good agreement with the numerical and asymptotic results presented in Blennerhassett (1980). However, as was mentioned in the introduction, the second neutral curve shown in fig. 2 of Blennerhassett (1980) and again in our figure 1*a*, is in fact a closed loop which does not belong to a second, distinct, mode of instability in the flow. In order to understand the events leading to the formation of the neutral loop, we consider the real and imaginary parts of the phase speed, *c*_{r}(*k*) and *c*_{i}(*k*) respectively, as a function of the wavenumber for a fixed Reynolds number in the flow. An example is shown in figure 2 for *Re*=10 000, which corresponds to the middle part of the loop.

As follows from figure 2, the loop is formed by two different modes, with the region inside the loop corresponding to stable disturbances. Later, we will show that the modes involved can be connected with the interfacial instability (I-mode) and Tollmien–Schlichting instability in air (TSA-mode), as labelled in the diagram. However, one should bear in mind that the mode identification is a matter of some subtlety, in this case because, in fact, the I-mode in figure 2 is connected to the lower branch of the broadband neutral curve whereas the TSA-mode is obtained by continuation from the upper branch of the broadband neutral curve. In this case, to be able to trace mode behaviour in a consistent manner, we have chosen to identify I and TSA modes by their longwave behaviour.

At lower Reynolds numbers the interfacial waves neutral curve is formed by one mode (figure 3); therefore the mode structure shown in figure 2 has no direct continuation into lower Reynolds numbers.

Let us consider the variation of the two modes with the Reynolds number in more detail. To this end, the values of the phase speed, *c*, for the two modes are first calculated at a fixed point in the wavenumber-Reynolds number plane (*k*=1.4, *Re*=10 000 in our computations) and then extended parametrically to other Reynolds numbers in as wide a range of wavenumbers as necessary. The results can be summarized as follows.

For a sufficiently low Reynolds number, the I-mode forms the broadband neutral curve whereas the TSA-mode remains stable (see figure 3). As the Reynolds number increases, the distance between the mode curves in (*k*, *Re*) parameter plane decreases. The imaginary part of the phase speed for the I-mode develops a local minimum, which at gives rise to a narrow range of stable waves with *c*_{i}<0, hence indicating the beginning of the stable loop (figure 4*a*). At the same time, the imaginary part of the phase speed for the TSA-mode develops a spike nearly touching the neutral level with *c*_{i}=0 (figure 4*b*). For the Reynolds numbers in figure 4, the I-mode alone is responsible for the broadband interfacial instability (for both short and long waves) and for the neutral point of the loop and a range of stable waves inside the loop.

At the Reynolds number *Re*=9800 (*ca*), the I-mode and TSA mode intersect and exchange branches as shown in figure 5. Now the loop is formed by two different modes (as in figure 2 earlier) with the mode identity essentially lost. A branch of the previously distinct TSA-mode is now connected with the interfacial I-mode for longer waves, the other branch of the TSA-mode is connected to the short-wave I-mode.

At still higher Reynolds numbers, the distance between the modes first increases, but then it begins to decrease. The TSA-mode develops a local maximum in *c*_{i} and the modes exchange branches once again between *Re*=12 400 and *Re*=12 450 (see figure 6). After the second mode crossing, the mode configuration returns to the same arrangement as before mode intersection, with both the loop and the broadband neutral curve formed by the I-mode alone. At still higher Reynolds numbers, the local range of stable waves of the I-mode narrows down and eventually vanishes at approximately *Re*=20 553. This point marks the closure of the loop.

The two intersections between the I-mode and TSA-modes described above take place at points representing stable waves. Before the first and immediately after the second mode crossing, i.e. at *Re*<9600 and , the TSA-mode is stable. However, at somewhat higher Reynolds numbers, starting from approximately *Re*=13 767, the TSA mode becomes unstable and forms its own neutral curve as can be seen in figure 7.

In our search for the starting points for the modes at *k*=1.6 *Re*=10 000 we came across a root with *c*=−0.136 440 2−0.000 788 2i. At higher Reynolds numbers, this root leads to instability associated with a Tollmien–Schlichting mechanism acting in the liquid layer (water) or TSL-mode, for short. According to our computations this mode is well separated from the other two modes.

The overall patterns of the neutral curves for the three modes described above are shown in figure 8. The critical conditions are summarized in the table 2.

At still higher Reynolds numbers, intersections between the I-mode and TSA-mode reappear again, however, accurate computations for these events become very demanding, hence only a brief summary of our results will be presented here. Just before the next mode crossing, the *c*_{i}(*k*) distribution for the I-mode again develops a local minimum that produces a stable region starting from approximately *Re*=50 450, while the TSA-mode neutral curve is shifted towards lower wave numbers (figure 9). The mode crossing itself occurs between *Re*=83 000 and *Re*=83 400 (figure 10). Afterwards, the local minimum in the graph of *c*_{i}(*k*) disappears and the modes become well separated.

In our description above, we assigned the TSA and TSL abbreviation to two distinct modes with reference to Tollmien–Schlichting instability mechanisms operational in air and in water separately. We see that in our case the Tollmien–Schlichting instability mechanism appears in distinct forms in the two media. This can probably be attributed to the large viscosity and density ratios for the two fluids involved. Nevertheless, a number of indirect observations and further parametric studies allow the modes to be linked with a particular physical mechanism. Consider, for example, the phase speed of neutral TSA and TSL waves. In our computations above, the base flow speed at the interface is *u*_{a}^{o}(0)=0.3076. As follows from our computations, the phase speed for the unstable TSL-waves is lower than that at the interface and it tends to zero as *Re*→∞, whereas the phase speed for unstable I-waves and TSA-waves exceeds the flow speed at the interface which is in agreement with results reported by Sangalli *et al*. (1997).

An alternative argument helping to identify the two Tollmien–Schlichting modes can be based on the variation of the relative thickness of the liquid layer. Consider the case of decreasing water layer thickness. In the limit of a thin water film, the velocity distribution in the base channel flow approaches a parabolic Poiseuille flow profile in air and a Couette flow linear profile in water. In this limit, the TSA mode should turn into the standard Tollmien–Schlichting instability in a channel filled with air. At the same time, since the Couette flow is stable to linear perturbations, one could expect some form of degeneration for the TSL-mode. In order to illustrate these considerations, we performed a set of computations for a fixed thickness of the air layer and a variable relative thickness of the liquid layer denoted as *n* (in the computations above, *n*=1). The base flow profile specified in §2 was altered accordingly. As expected, with the water layer thickness decreasing, the neutral curve for the TSA-mode approaches the neutral curve of the Tollmien–Schlichting mode in a single fluid case. This is illustrated in figure 11*a* for the neutral curves and in figure 12 for the critical parameters. For the results shown in figure 11*a*, the Reynolds number and wavenumber were renormalized to comply with the limit case of the Poiseuille air flow by taking the half-width of the air layer as a characteristic length and mean velocity by profile in the air flow as a characteristic velocity. The behaviour of the TSL-mode for a decreasing liquid layer thickness is shown in figure 11*b*. For these computations the Reynolds number and wavenumber were renormalized by taking the half-width of the water layer as a characteristic length and mean velocity by profile for *n*=1 in both layers. The critical Reynolds number for this mode increases and tends to infinity for a film of vanishing thickness.

## 5. Discussion

The flow configuration considered in this paper is a classical test bed for the stability theory of layered flows. The complicated pattern of the stability diagram described above came as a surprise to us. One might question the practical relevance of the linear stability results at Reynolds numbers above the first critical value where the most interesting events seem to be taking place in our work. We believe, however, that many features discovered in our computations are transferable to other, more realistic, regimes for flows in different geometries or with different combinations of fluids, with even more complex forms of mode interactions to be expected in flows with three or more distinct layers. In all such cases a clear understanding of the physical mechanisms behind observed or predicted instability is important. Not less important is a realization of the fact that the physics of instability in layered flows can easily be obscured by the mode crossing phenomena.

Based on the results obtained, we arrive at the following picture of the linear stability in this case. The TSL-mode was clearly identified as Tollmien–Schlichting mode by its asymptotic behaviour as *Re*→∞, wave speed on the neutral curve tends to zero (checked to *Re*∼10^{6}). The wavespeed on the neutral curve also tends to zero and is always lower than the interface speed *O*(10^{−2}). I-Mode can be definitely attributed to the surface waves following the analysis of Blennerhassett (1980). Calculations made for the point of intersection between the lower branch of the neutral curve and the *Re* axis were in good comparison with the results from the asymptotics of Blennerhassett (1980). The nature of the TSA-mode is clear from its asymptotic behaviour as water film thickness tends to zero. Therefore, the TSA-mode can be clearly attributed to the shear mode in air and TSL-mode—to the liquid shear mode.

The wavespeed of the neutrally stable solutions for Tollmien–Schlichting mode in air is faster than the velocity of the free boundary and this mode interacts with surface waves. The wavespeed of neutrally stable solutions for Tollmien–Schlichting mode in water is lower than the speed of the interface. That mode does not interact with two other modes and produces a neutral curve typical for Poiseuille–Couette flows (see, for example, Cowley & Smith 1985).

The interaction between the gas shear mode and interfacial mode occurs for the wavenumbers *k*∈[1, 2], which corresponds to wavelengths of the order of the channel height. The minimum initially developed by the *c*_{i}(*k*) for the interfacial mode is located approximately at the same *k* as the maxima of the stable gas shear mode. Also, the wave speeds of the interfacial mode and shear mode are close enough in the domain of interaction in *k*. It is possible that what we observe is in fact energy transfer from the interfacial waves to the gas shear mode, which takes place for the disturbances of *O*(1) wavelengths. This leads to the excitation of the stable gas shear mode and stabilization of the unstable interfacial mode. Apparently, the conditions that allow for this energy transfer to take place can be periodical in *Re* as the interaction and branch exchange is repeated at higher Reynolds numbers. Therefore, it is possible that the stabilization of the interfacial waves can be attributed to the energy transfer between interfacial and shear mechanisms of instability and energy transfer across the spectrum of the interfacial waves, rather than the interaction between different parts of the interfacial instability spectrum alone.

Another example of the observed interaction between gas shear and interfacial modes can be found in the results obtained by Timoshin & Hooper (2000), who have shown that Tollmien–Schlichting mode and interfacial waves mode can have points of coalescence. The results presented in our paper extend these observations to the case of channel flows. We can also note that the effect of interfacial wave stabilization and interaction between interfacial and gas shear modes has been observed by different authors for a wide variety of problems, which suggests that it possesses certain structural stability.

## Acknowledgements

The authors are grateful to Professor B. J. Boersma for pointing out the significance of the velocity profile curvature at the interface. Many useful suggestions made by an anonymous referee are also greatly appreciated.

## Footnotes

- Received May 30, 2004.
- Accepted November 17, 2004.

- © 2005 The Royal Society