Direct numerical simulations of three-layer viscosity-stratified flow.

Cao, Qing ; Sarkar, Kausik ; Prasad, Ajay K. 等

Department of Mechanical Engineering, University of Delaware, Newark, DE, U.S. 19716 Two-dimensional simulations of flow instability at the interface of a three-layer, density-matched, viscosity-stratified Poiseuille flow are performed using a front-tracking/finite difference method. This is an extension of the study for the stability of two-layer viscosity-stratified flow of Cao et al., Int. J. Multiphase Flow, 30, 1485-1508 (2004). We present results for large-amplitude non-linear evolution of the interface for varying viscosity ratio m, Weber number We, and phase difference between the perturbations of the two interfaces. Strong non-linear behaviour is observed for relatively large initial perturbation amplitude. The higher viscosity fluid is drawn out as a finger that penetrates into the lower viscosity layer. The finger originates at the crest of the perturbation at the interface. The simulated interface shape compares well with previously reported experiments. Increasing interfacial tension retards the growth rate of the interface as expected, whereas increasing the viscosity ratio enhances it. The sinuous instability appears to evolve faster than the varicose one. For certain flow parameters the high-viscosity finger displays a bulbous tip, which is also seen in our previously conducted experiments and two-layer results, although it is less pronounced. The low-viscosity intruding finger does not display this curious bulbous tip. Drop formation is precluded by the two-dimensional nature of the calculations.

On a effectue a l'aide d'une methode de differences finies avec capture de front des simulations bidimensionnelles de l'instabilite de l'ecoulement a l'interface d'un ecoulement de Poiseuille stratifie en viscosite mais de masse volumique a trois couches. Ceci fait suite a l'etude de la stabilite de l'ecoulement stratifie en viscosite a deux couches de Cao et al., Int. J. Multiphase Flow, 30, 1485-1508 (2004). Nous presentons des resultats pour l'evolution non lineaire a grande amplitude de l'interface pour un rapport de viscosite m, un nombre de Weber We et une difference de phases entre les perturbations des deux interfaces. Un comportement non lineaire marque est observe pour une amplitude de perturbation initiale relativement grande. Le fluide ayant la viscosite la plus elevee prend la forme d'un doigt penetrant la couche de viscosite inferieure. Le doigt prend naissance a la crete de la perturbation a l'interface. La forme simulee de l'interface se compare bien aux experiences decrites anterieurement. L'augmentation de la tension interfaciale retarde la vitesse de croissance de l'interface comme prevu, tandis que l'augmentation du rapport de viscosite la renforce. Il semble que l'instabilite sinueuse evolue plus vite que la variqueuse. Pour certains parametres d'ecoulement, le doigt de plus forte viscosite montre une pointe bulbeuse, ce qu'on a egalement observe dans nos experiences anterieures et resultats a deux couches, bien que de facon moins prononcee. Le doigt intrusif de faible viscosite ne montre pas cette curieuse pointe bulbeuse. La formation de gouttelettes est evitee par la nature bidimensionnelle des calculs.

Keywords: three-layer viscosity-stratified flow, direct numerical simulations, front tracking, finger formation

INTRODUCTION

This work is motivated by a problem encountered in certain mixing processes that employ a centre line injector upstream of a static mixer. Static or motionless mixers have been adopted for a large number of blending and dispersion operations because they offer attractive features such as closed-loop operation, and no moving parts, in contrast to continuously-stirred tank reactors. The two liquids to be mixed (for example, a polymer melt and an additive) are forced under high pressure through the mixer placed inside a tube, and the two liquids are cut and folded repeatedly as they negotiate the bends and openings within the mixer (see Ventresca et al. (2002) for a photograph of a static mixer). A typical inlet geometry to the static mixer is the centre line injector depicted schematically in Figure 1a. For certain flow parameters of the injected liquid and the co-flowing liquid (such as viscosity ratio, interfacial tension, and volume flux ratio) a wavy oscillation could develop immediately downstream of the injector as shown in Figure 1b. An extreme manifestation of this instability could lead to a breakup of the injected stream leading to axially segregated clumps of the injected liquid. Because static mixers are less efficient at removing axial concentration gradients than radial gradients, this instability could have a detrimental impact on the overall mixed product quality. The instability was studied experimentally using a mixer with axisymmetric geometry by our group and the results are available in Cao et al. (2003). The primary focus of the investigation was the viscosity stratification between the injected stream and the co-flowing stream, which is the primary cause of the instability.

[FIGURE 1 OMITTED]

Subsequently, we performed direct numerical simulations of viscosity-stratified two-layer plane Poiseuille flow employing a front-tracking/finite difference method (Cao et al., 2004). We demonstrated in that paper that the front-tracking method reproduced analytically predicted growth rates for the linear instability regime, and also provided new results for the nonlinear regime for large interface amplitudes, such as finger formation. These elongated fingers closely matched flow features observed in our experiments. This paper also contains an extensive review of two-layer viscosity-stratified studies, most of which are analytical.

The current paper is an extension of Cao et al. (2004) in that we now examine three-layer viscosity stratified flow, which better matches the centre line injector flow, albeit in a planar geometry. In contrast to Cao et al. (2004) our current focus is not the linear instability regime; instead we now concentrate exclusively on the large amplitude non-linear instability. The three-layer configuration encompasses a myriad of flow and geometrical parameters and it is impractical to explore all possible parameter combinations here; therefore, we have chosen to fix certain parameters and vary others as described subsequently. In particular, the three-layer configuration allows the exploration of sinuous and varicose interfacial modes.

Multi-layer viscosity-stratified flow has been a topic of many theoretical investigations, and a brief review is presented here. Yih (1967) used a long-wave perturbation analysis to show that two-layer, viscosity-stratified plane Poiseuille flow and plane Couette flow can be unstable for arbitrarily small Reynolds numbers. The growth rate was found to be proportional to [[alpha].sup.2]Re, where [alpha] is the dimensionless wave number and Re is the Reynolds number. Yih's long-wave perturbation method was applied to a number of multi-layer shearing flows. Among them, Li (1969) expanded Yih's (1967) work by considering plane Couette flow with three layers of fluids of different viscosities and densities. His principal finding was that while a single surface of viscosity discontinuity will cause instability, the presence of an additional surface of discontinuity may or may not be stabilizing depending on the values of the depth ratio and viscosity ratio at the additional interface. Hickox (1971) applied Yih's method to axisymmetric vertical pipe flow of two fluids with different densities wherein the core is less viscous than the annulus, and concluded that the primary flow was always unstable to either asymmetric or axisymmetric disturbances, and second, that the instability was primarily due to viscosity stratification.

Weinstein and Chen (1999) considered the stability of three-layer flow with vanishing Reynolds number and negligible surface tension down an inclined plane at finite wavelengths. They found that the largest growth rates were achieved when a low-viscosity and thin layer is centrally located in the film. The maximum growth rates were found at intermediate wavelengths.

Experimental studies of multi-layer viscosity stratified flow are not available in the literature. However, viscosity stratified flow in an axisymmetric geometry (core-annular flow) has been studied by a few researchers. For example, Rozen and Baldyga (2003) present photographs of the varicose instability forming down stream of a centre line injector for the case in which the core flow is less viscous than the annular flow (viscosity ratio is about 400). Our own experimental work (Cao et al. 2003) pertained to viscosity-stratified flow with vanishing interfacial tension using laser-induced fluorescence (LIF) and particle image velocimetry (PIV). A comparison of experimental stable cases and exact solutions revealed the existence of a thin interfacial layer that smoothes out the discontinuity of the velocity gradient at the interface. We also observed two kinds of unstable modes for the first time: (1) wavy core-flow with fissures; and (2) wavy core-flow with core breakup. Our results confirmed that viscosity stratification can cause instabilities even when the Reynolds number is O(1).

The typical procedure for most of the analytical studies listed above is to start with the Navier-Stokes equations, introduce a perturbation on the primary motion, derive the Orr-Sommerfeld equations, and solve them using a suitable numerical method. Direct numerical solutions (DNS) of the Navier-Stokes equations are difficult due to the unsteady evolution of the interface between dissimilar fluids, implying that the interface shape must be determined concurrently using equations that are coupled with the Navier-Stokes equations. Such a numerical solution is particularly attractive because it does not face constraints that are typical in analytical treatments, such as linearity and small-amplitude perturbations. For example, DNS calculations using the volume of fluid (VOF) method for two-layer Couette flow by Coward et al. (1997) showed that the interface evolves to form waves with a steep front for certain parameter regimes. Subsequently, finger formation at the interface in the non-linear regime was presented for the same problem by Li et al. (1998) and Renardy and Li (1999) also using the VOF method. More recently, Rozen and Baldyga (2003) demonstrated good agreement between their experimentally observed varicose instability in viscosity-stratified core-annular flow and their simulations using a commercial code; however, neither their experiments nor simulations extended to the regime where finger formation is evident.

The front-tracking/finite difference method was originally formulated by Tryggvason and co-workers (for example, see Unverdi and Tryggvason, 1992). Recently Sarkar and Schowalter (2001a, b) used it to study drop deformation in time dependent flows at finite Reynolds numbers. Zhang et al. (2002) applied the method to a gravity-driven two-layer fluid flow in an inclined channel. A growing finger was observed at the interface under certain conditions. The front-tracking method was also applied by Tauber et al. (2002) to conduct calculations on the Kelvin-Helmholtz stability of an interface between two immiscible fluids. Cao et al. (2004) studied the linear and non-linear stability of a pressure driven two-layer Poiseuille flow using a front-tracking method. The growth rates of linear stability were compared with analytical analysis, and the effect of viscosity ratio and surface tension on the propagation of the interface was presented. The appearance of fingers at the interface was investigated. The front-tracking method offers greater generality over boundary element methods, and is an attractive alternative to finite difference (with a body-fitted coordinate system) or finite element implementations.

In this paper, we extend our study of two-layer flow (Cao et al. 2004) to a three-layer viscosity-stratified Poiseuille channel flow. The next two sections describe the mathematical formulation of the problem and its numerical implementation. Results from a detailed two-dimensional computation are presented for different parameters in the fourth section. We apply the method to study non-linear stabilities and obtain new and interesting results for large interfacial amplitude evolution for varying viscosity ratio and interfacial tension. The relevant parameters are systematically varied and the effects collated and explained. We summarize our findings in the last section.

MATHEMATICAL FORMULATION

Governing Equations

The velocity field u and the pressure p satisfy the equation of momentum conservation:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (1)

in the entire domain consisting of the two continuous fluid domains [[OMEGA].sub.1] and [[OMEGA].sub.2] (Figure 2). Here [sigma] is the interfacial tension, [partial derivative]B is the interface (front) consisting of points [x.sub.B], [kappa] the local curvature, n the normal to the surface, [delta](x - [x.sub.B]) is the Dirac delta function (two-dimensional for the present problem), [mu] is the viscosity, and the superscript T represents the transpose of the velocity gradient [nabla]u. Interfacial tension produces a jump in the normal stress across the interface, and is represented as (singular) distributed body force, anticipating its numerical implementation to be described below. The flow field is incompressible with density [rho],

[FIGURE 2 OMITTED]

[nabla] x u = 0 (2)

Note that the velocity field satisfies a single equation in both phases with a spatially varying viscosity [mu](x,t) that satisfies

D[mu]/Dt [equivalent to] [partial derivative][mu]/[partial derivative]t + u x [nabla][mu] = 0 (3)

Moreover, by applying the momentum Equation (1) in a rectangular element of vanishing thickness straddling the front, one can recover velocity and shear stress continuity across the front, and the jump in the normal stress due to interfacial tension.

Geometry and Initial Conditions

Figure 2 depicts the two-dimensional flow geometry under consideration. We consider two liquids that are co-flowing in a channel. The inner liquid has a thickness [d.sub.1], with a volume flux [Q.sub.1] and viscosity [[mu].sub.1]. The outer liquid consists of two layers of fluid; the upper and lower liquid layers each have thickness [d.sub.2], volume flux [Q.sub.2], and viscosity [[mu].sub.2].

We use the velocity profile of the primary flow as the initial velocity field. The primary flow has only one nonzero velocity component, u(y). The corresponding steady-flow governing equation is:

[d.sup.2]u/d[y.sup.2] = 1/[mu] dp/dx (4)

The relevant boundary conditions are no-slip at the upper and lower walls, velocity continuity at the interfaces, and continuity of shear stress at the interfaces. Periodic boundary conditions are applied at the inlet and outlet region to simulate the spatially evolving flow. We define the thickness ratio as n = [d.sub.1]/([d.sub.1]+2[d.sub.2]) and the viscosity ratio as m = [[mu].sub.1]/[[mu].sub.2].

Equation (4) is easily solved to yield:

[u.sub.1]/[U.sub.0] = 1 + 1/[n.sup.2](m - 1) - m [(y/[d.sub.1]/2 + [d.sub.2]).sup.2] (5)

[u.sub.1]/[U.sub.0] = m/[n.sup.2](m - 1) - m [[(y/[d.sub.1]/2 + [d.sub.2]).sup.2] -1] (6)

Where [U.sub.0] is the centre line velocity given by:

[U.sub.0] = -[([d.sub.1]/2 + [d.sub.2]).sup.2]/2[[mu].sub.2] [[n.sup.2] m - 1/m - 1] dp/dx (7)

Interface

The velocity at a point on the interface u([x.sub.B]) is related to the field velocity using the property of the delta function:

u([x.sub.B])[[integral].sub.[OMEGA] dx[delta](x - [x.sub.B])u(x) (8)

As noted earlier, the interface conditions of stress and velocity continuities are automatically met by the governing equation with spatially varying viscosities and the distributed forces (due to interfacial tension) in the field equation.

For the simulations, we introduce a sinusoidal perturbation at the interface, y = [y.sub.0] + [a.sub.0] cos(2[pi]x/[lambda][phi].), where y is the vertical position of the interface, [y.sub.0] is the vertical position of the interface without perturbation, [a.sub.0] is the initial amplitude of the perturbation , x is the downstream distance, [lambda] is the wavelength, and [phi] is the phase. We consider two configurations for the phases of the two fronts, one is [[phi].sub.1] = [[phi].sub.2] = [pi], the other is [[phi].sub.1] = 0 and [[phi].sub.1] = [pi]; these represent the sinuous and varicose modes, respectively. Than et al. (1987) have shown in their three-layer stability analysis that it is sufficient to consider even (sinuous) and odd (varicose) modes in order to recover the full spectrum of the problem.

NUMERICAL IMPLEMENTATION

The physical domain is represented numerically as a box of size [L.sub.x] and [L.sub.y]. We have used [L.sub.x] = [L.sub.y] = [d.sub.1] + 2[d.sub.2] and a 256 x 256 grid. The interfaces between the two fluids are described by line elements created initially by placing points on the line. The movement of the element vertices describes the evolving shape of the interface. An adaptive regridding scheme is implemented to prevent the elements from distorting excessively. The scheme creates/destroys elements by insertion/removal of points on the existing front.

Front Tracking

The properties of the outer fluid (such as [[mu].sub.2] and [[mu].sub.2]) could be different from those in the inner fluid, although the current simulations are for density-matched liquids. The present method reduces the three layers to a single layer with spatially varying properties:

[mu(x) = [[mu].sub.2] + ([[mu].sub.1] - [[mu].sub.2])I (9)

The indicator function I(x) is 0 when x belongs to the outer fluid and 1 otherwise. The following equation for I(x) can be derived:

[[nabla].sup.2]I(x) = [[integral].sub.[partial derivative]B d[x.sub.B] [nabla] x n[delta](x - [x.sub.B] (10)

A smooth representation of the a-function is required for the numerical implementation of Equations (1), (8) and (10) (Sarkar and Schowalter, 2001a):

D(x - [x.sub.B] = [D.sup.1] (x - [x.sub.B])[D.sup.1](y - [y.sub.B] (11)

where

[D.sup.1] (x - [x.sub.B] = 1/4[DELTA]x [1 + cos [pi]/2[DELTA]x (x - [x.sub.B])] for |x - [x.sub.B] [less than or equal to]2[DELTA]x (12)

and [DELTA]x is the grid spacing. The approximation of the delta function is coupled with the discretization of the computational domain: as the discretization length [DELTA]x approaches zero, the approximant approaches infinity, as required of a family of regular functions approaching a delta function. Essentially, the sharp interface separating the phases is replaced by a region of sharp variation in properties that has a finite thickness of approximately 4[DELTA]x. This continuously stratified fluid interface is more or less analogous to real situations involving miscible fluids. On a cautionary note, Wilson and Rallison (1999) found for channel flow of elastic liquids that as the thickness of the layer over which the elastic properties vary is increased, the instability mechanism is countered by convective effects and the growth rate is reduced. A similar effect could apply in our simulations as well.

Finite Difference

This formulation leads to a system of partial differential equations with smooth spatially varying coefficients. The front is decoupled from the underlying flow equation, and is retained only as a means for computing properties at successive time steps. We use an MAC type operator splitting/projection finite difference method. The original method solves the system in two explicit steps. A detailed description of the method can be found in Sarkar and Schowalter (2001a). The numerical solutions were confirmed to be grid-independent by performing grid-refinement studies, and were also shown to be independent of the chosen time-step (Cao et al. 2004).

RESULTS

Relevant dimensionless parameters include the viscosity ratio m, the thickness ratio n, dimensionless velocities u * and v *, dimensionless wavelength e*, dimensionless amplitude of the perturbation a *, dimensionless time t *, and the phase difference between the two interfacial waves, A[phi]:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (13)

where u, v, [lambda], a, and t are the corresponding dimensional quantities. We consider two baseline cases in the paper (Table 1). In the first baseline case, the inner fluid is more viscous than the outer fluid with m = 10, n = 0.375, Reynolds number Re = 5, Weber number We = [infinity], and phase difference [DELTA][phi] = 0. Reynolds number is defined as Re = [U.sub.0][d.sub.1]/[v.sub.1]1, where [U.sub.0] is the centre line velocity, and d1 and [v.sub.1] are the thickness and kinematic viscosity, respectively, of the inner fluid. Weber number is defined as We = [rho][U.sub.0][d.sub.1], where sigma is the interfacial tension. For the second baseline case, we switch the viscosities of the two fluids (m = 0.1) with all the other parameters unchanged. We use these two parameter sets based on our original experiments (Cao et al. 2003) of instabilities during centre line injection of a lower viscosity liquid into a higher viscosity miscible co-flow in a pipe with relatively low Reynolds numbers (up to 13). In the current work Re is fixed at 5, and n is fixed at 0.375 to resemble the similar value of n in the experiments. For viscosity ratio m, we pick these two extreme baseline cases because we are interested primarily in the effect of viscosity ratio. However, m = 2 and m = 0.5 are also examined. The fourth parameter in Table 1 is the phase difference between two interfaces [DELTA][phi] = [[phi].sub.1] - [[phi].sub.2]. Based on our own experimental observations, [DELTA][phi] is set to 0 in baseline cases (sinuous instability). But we also study [DELTA][phi] [not equal to] 0 because Joseph and co-workers (Bai et al. 1992) have observed bamboo (or varicose) instabilities corresponding to [DELTA][phi] = [pi]. The last parameter in Table 1 is Weber number. Our original experiments were conducted with aqueous solutions of CMC in both streams, with negligible interfacial tension. So We = [infinity] is a natural choice for the baseline cases. Additional calculations were conducted for different Weber numbers. Two other parameters that are not listed in Table 1 are [lambda] * and [a.sub.0] *. In the current study, we fixed these two parameters as [lambda] * = 1 and [a.sub.0] * = 0.02, based on our experimental observations (Cao et al., 2003). It should be noted that although our experiments were conducted in an axisymmetric geometry, the interface exhibits 3-D features. Owing to the 2-D nature of our simulation, such features cannot be reproduced; however, we will demonstrate good qualitative agreement with most of the observed flow characteristics.

First, we investigate the flow with m = 10, We = [infinity], and [DELTA][phi] = 0. Figure 3 depicts the interface position in the computational domain for different times with a time increment of 5. The two fronts are slightly distorted at t * = 5. At t * = 10 steepening of the interface at the crest is evident, which is significantly different from t * = 5. By t * = 15, a finger of high-viscosity liquid is drawn into the surrounding low-viscosity co-flow. The finger elongates continuously with time and experiences a simultaneous thinning. Interestingly, the finger of high-viscosity liquid terminates in a bulbous tip. Similar shapes for two-dimensional fingers have been reported in Zhang et al. (2002), where a two-layer gravity driven fluid in an inclined channel was considered using a front-tracking method. As the high-viscosity finger grows in length, so does a corresponding finger of low-viscosity liquid penetrate into the high-viscosity layer; the low-viscosity intrusion does not display this curious bulbous tip. Our experiments (Figure 1b) provide validation for our computational results. In Figure 1b, "A" corresponds to the intrusion of more viscous fluid into less-viscous inner fluid and terminates with a bulbous tip. In contrast, "B" corresponds to the intrusion of less-viscous fluid into the more-viscous outer layer and displays a sharp tip. Several other instances of sharp and bulbous tips are evident in Figure 1b. The 2-D nature of the current calculations precludes the breakup of the elongated finger into drops because, unlike an axisymmetric filament, surface area for a finger represented by a two-dimensional sheet actually increases upon the superposition of a small sinusoidal undulation.

[FIGURE 3 OMITTED]

We plot in Figure 4 the acceleration field obtained by subtracting the flow field at t * = 25 from the flow field at an earlier time separated by [DELTA]t* = 0.25. Since the two fronts behave identically, we focus on the upper front. The acceleration field indicates motion everywhere in the channel with greater activity in the vicinity of the fingertips. This is not surprising given that finger penetration is the most prominent motion for the flow. This motion is characterized by two counter-rotating eddies riding with the tips. Interfacial tension has a very strong effect on the evolution of the interfaces as depicted in Figure 5a. We choose front length as a measure of interface evolution because the amplitude of the interface saturates after a point. For zero interfacial tension (We = [infinity]), the interface grows rapidly as already shown in Figure 3. When the interfacial tension is increased to We = 134, the interfacial growth is highly damped; when interfacial tension is increased further to We = 67, the wave begins to diminish in amplitude to the point that the interface eventually assumes a flat profile.

[FIGURES 4-5 OMITTED]

Figure 6 presents interface profile evolution for a smaller viscosity-stratification (m = 2) as compared to the baseline case 1 (Figure 3). Front steepening is evident at t * = 5 for the m = 2 case, with a high-viscosity finger already well formed by t * = 10. In contrast, these two features only develop by t * = 10 and 15 for the m = 10 case (Figure 3). Thereafter, the front length grows in a linear manner with time. The growth rate of the interface length with time for m = 2 is plotted in Figure 5b. Although the interface length for m = 2 increases more rapidly than m = 10 in the short-term, the longer-term behaviour is opposite: at about t * = 24, the m = 10 interface length overtakes the m = 2 front length. The interface shape exhibits distinct differences as m is reduced from 10 to 2. The finger of high-viscosity fluid grows in close proximity to the parent layer for m = 2, and is aligned in a highly parallel manner to it. In contrast, the central layer for m = 10 exhibits greater undulations with the low-viscosity fluid penetrating deeper into it. The bulbous fingertip is not evident for m = 2.

[FIGURE 6 OMITTED]

The acceleration field for m = 2 (Figure 7) indicates that the central layer is more mobile as compared to the more viscous central layer in Figure 3. Two counter-rotating vortices are clearly visible in each of the three layers. Vectors adjacent to the fingertips indicate that the high-viscosity finger is penetrating faster into the low-viscosity co-flow than vice-versa.

[FIGURE 7 OMITTED]

Figure 8 is a result for parameters identical to baseline case 1 (Figure 3) except that [DELTA][phi] = [pi]. The interface evolution is slightly delayed in comparison to the [DELTA][phi] = 0 case. Front steepening is not evident even at t * = 10. Thereafter, a high-viscosity finger originates at the crest of the interface and penetrates into the low-viscosity layer. Interestingly, the finger develops a pronounced bulbous tip at t * = 25 and beyond, similar to the [DELTA][phi] = 0 case. Such a bulbous tip is not visible at the tip of the low-viscosity finger. The [DELTA][phi] = [pi] case develops more slowly than the [DELTA][phi] = 0 case for all t * as depicted in Figure 5b. Thus, the sinuous instability appears to evolve faster than the varicose one. This result provides an explanation for the fact that we always observed the sinuous instability in our experiments (Cao et al, 2003).

[FIGURE 8 OMITTED]

The acceleration field is presented for the [DELTA][phi] = [pi] case in Figure 9. A pair of counter-rotating vortices is evident in the vicinity of the high-viscosity fingertip, again similar to baseline case 1 (Figure 4). The local flow field generated by these vortices is probably responsible for the accumulation and growth of fluid at the fingertip.

[FIGURE 9 OMITTED]

Now, we turn our attention to baseline case 2, wherein the central layer is less viscous than the outer layers (m = 0.1). The interface evolution with time is depicted in Figure 10. The calculation could not be carried as far as the previous cases due to the much smaller time-step for this case. The nominal velocity profile for this case resembles a "jet-flow", where the central layer travels rapidly downstream owing to its lower viscosity, whereas the outer layers are relatively sluggish. A low-viscosity finger originates at the crest and penetrates into the high-viscosity outer layer. The high-viscosity layer also penetrates into the low-viscosity layer; in this case, a bulbous tip again appears in the high-viscosity finger, although it is less pronounced.

[FIGURE 10 OMITTED]

Figure 11 displays the acceleration field for baseline case 2, at t * = 20. Due to its lower viscosity, the central layer appears more mobile than the outer layers. The motion is quite complex with a number of vortices distributed within the central layer. The high-viscosity finger experiences a shear stress from the lower viscosity layer that tends to elongate it over much of its length. However, the vectors are directed opposite to the high-viscosity finger at its tip, which may explain the slight thickening of the fingertip.

[FIGURE 11 OMITTED]

The effect of interfacial tension on the evolution of the interfaces for m = 0.1 case is depicted in Figure 12a. For zero interfacial tension (We = [infinity]), the interface grows rapidly as already shown in Figure 10. When the interfacial tension is increased to We = 1.34, the interfacial growth is highly damped and the wave decreases in amplitude to the point that the interface eventually assumes a flat profile; when interfacial tension is increased further to We = 0.67, the wave begins to diminish in amplitude more rapidly than We = 1.34 case.

[FIGURE 12 OMITTED]

The interface for the m = 0.5 case (Figure 13) evolves in a qualitatively similar manner to the m = 2 case (Figure 6), except that the former develops more rapidly. This is confirmed by plotting front length with time (see Figures 5b and 12b). Additionally, the m = 0.5 case develops more rapidly than the m = 0.1 case in the short-term, but is overtaken by the latter at about t * = 11 (Figure 12b). Thus, the long-term growth rate of the front length is more rapid for the baseline cases of m = 10, and m = 0.1. Here again, the implication is, greater the viscosity stratification, greater the rate of evolution of the interface. Figure 14 shows the acceleration field for the m = 0.5 case, at t * = 30. A complicated distribution of several vortices at the central layer characterizes the flow.

[FIGURES 13-14 OMITTED]

Interface evolution for m = 0.1, [DELTA][phi] = [pi] is plotted in Figure 15. The interface develops more rapidly when compared with the m = 10, [DELTA][phi] = [pi] case (Figure 8; compare also Figures 5b and 12b). The [DELTA][phi] = [pi] case develops more slowly than the [DELTA][phi] = 0 case for all t * as depicted in Figure 12b. Thus, the sinuous instability appears to evolve faster than the varicose one for both m = 0.1 and m = 10 cases. Interestingly, the finger forms at the trough of the wave. A close-up of the acceleration field in the neighbourhood of the finger clearly indicates that the flow in the central layer is responsible for drawing the finger out. The total acceleration field for this case is shown in Figure 16 and confirms that the flow responsible for finger elongation is part of a counter-rotating vortex pair situated near the finger.

[FIGURES 15-16 OMITTED]

CONCLUSIONS

The stability of three-layer density-matched viscosity-stratified Poiseuille flow using a front-tracking/finite difference method was examined. The simulations have reproduced features observed experimentally, such as the sharp curvature of the interface and the fingers penetrating from the more viscous fluid into the less-viscous layer. We examine the non-linear evolution of the interface for n = 0.375, and Re = 5 with varying m, We, and [DELTA][phi]. Interfacial tension impedes finger development and hinders the growth of the front length greatly. The front length grows more rapidly for larger m when m >1, and for smaller m when m < 1. The finger appears to initiate at a location near the crest of the interface. The varicose interface mode develops more slowly than the sinuous one, which probably explains why we only observed the latter in our experiments. The high-viscosity finger displays a bulbous tip for certain flow parameter conditions; this is apparently in response to the local acceleration field in the fluid adjacent to the fingertip.

REFERENCES

Bai, R., K. Chen and D. D. Joseph, "Lubricated Pipelining: Stability of Core-Annular Flow. Part V: Experiments and Comparison with Theory," J. Fluid Mech. 240, 97-132 (1992).

Cao, Q., K. Sarkar and A. K. Prasad, "Direct Numerical Simulations of Two-Layer Viscosity-Stratified Flow," Int. J. Multiphase Flow 30, 1485-1508 (2004).

Cao, Q., A. L. Ventresca, K. R. Sreenivas and A. K. Prasad, "Instability Due to Viscosity Stratification Downstream of a Centre line Injector," Can. J. Chem. Eng. 81, 913-922 (2003).

Coward, A. V., Y. Y. Renardy, M. Renardy and J. R. Richards, "Temporal Evolution of Periodic Disturbances in Two-Layer Couette Flow," J. Comp. Phys. 132, 346-361 (1997).

Hickox, C. E., "Instability Due to Viscosity and Density Stratification in Axisymmetric Pipe Flow," Phys. Fluids 14, 251-262 (1971).

Li, C.-H., "Instability of Three-Layer Viscous Stratified Flow," Phys. Fluids 12, 2473-2481 (1969).

Li, J., Y. Y. Renardy and M. Renardy, "A Numerical Study of Periodic Disturbances on Two-Layer Couette Flow," Phys. Fluids 10, 3056-3071 (1998).

Renardy, Y. Y. and J. Li, "Comment on 'A Numerical Study of Periodic Disturbances on Two-Layer Couette Flow' [Phys. Fluids 10, 3056 (1998)]," Phys. Fluids 11, 3189-3190 (1999).

Rozen, A. and J. Baldyga, "Influence of Viscosity Difference on the Instability of the Core-Annular Flow," Task Quarterly (Poland) 7, 425-438 (2003).

Sarkar, K. and W. R. Schowalter, "Deformation of A Two-Dimensional Drop at Non-Zero Reynolds Number in Time-Periodic Extensional Flows: Numerical Simulation," J. Fluid Mech. 436, 177-206 (2001a).

Sarkar, K. and W. R. Schowalter, "Deformation of a Two-Dimensional Viscous Drop in Time-Periodic Extensional Flows: Analytical Treatment," J. Fluid Mech. 436, 207-230 (2001b).

Tauber, W., S. O. Unverdi and G. Tryggvason, "The Nonlinear Behavior of a Sheared Immiscible Fluid Interface," Phys. Fluids 14, 2871-2885 (2002).

Than, P. T., F. Rosso and D. D. Joseph, "Instability of Poiseuille Flow Two Immiscible Liquids with Different Viscosities in a Channel," Int. J. Eng. Sci. 25, 189-204 (1987).

Unverdi, S. O. and G. Tryggvason, "A Front-Tracking Method for Viscous, Incompressible Multi-Fluid Flows," J. Comput. Phys. 100, 25-37 (1992).

Ventresca, A., Q. Cao and A. K. Prasad, "The Influence of Viscosity Ratio on Mixing Effectiveness in a Two-Fluid Laminar Motionless Mixer," Can. J. Chem. Eng. 80, 614-621 (2002).

Weinstein, S. J. and K. P. Chen, "Large Growth Rate Instabilities in Three-Layer Flow Down an Incline in the Limit of Zero Reynolds Number," Phys. Fluids 11, 3270-3282 (1999).

Wilson, H. J. and J. M. Rallison, "Instability of Channel Flows of Elastic Liquids Having Continuously Stratified Properties," J. Non-Newtonian Fluid Mech. 85, 273-298 (1999).

Yih, C.-S., "Instability Due to Viscosity Stratification," J. Fluid Mech. 26, 337-352 (1967).

Zhang, J., M. J. Miksis, S. G. Bankoff and G. Tryggvason, "Nonlinear Dynamics of an Interface in an Inclined Channel," Phys. Fluids 14, 1877-1885 (2002).

* Author to whom correspondence may be addressed. E-mail address: prasad@me.udel.edu

Manuscript received October 21, 2005; revised manuscript received June 1, 2006; accepted for publication June 1, 2006.

Qing Cao, Kausik Sarkar and Ajay K. Prasad *

Table 1. Parameters for two baseline cases

Name Re n m [DELTA][phi] We

Baseline case 1 5 0.375 10 0 [infinity]

Baseline case 2 5 0.375 0.1 0 [infinity]

NOMENCLATURE

a amplitude of perturbation (m)
[a.sub.0] initial amplitude of perturbation (m)
d thickness of layer (m)
I indicator function (-)
[L.sub.x] horizontal size of calculation box (m)
[L.sub.y] vertical size of calculation box (m)
m viscosity ratio (-)
n thickness ratio (-)
p static pressure (Pa)
Q volume flux ([m.sup.3]/s)
Re Reynolds number (-)
t time (s)
u axial velocity (m/s)
[U.sub.0] centre line velocity (m/s)
We Weber number (-)
x downstream distance (m)
[x.sub.B] points on the boundary (m)
y vertical position of interface (m)
[y.sub.B] vertical position of a point on the interface (m)
[y.sub.0] vertical position of the unperturbed interface (m)

Greek symbols

[alpha] dimensionless wave number (-)
[delta] Dirac delta function ([m.sup.-2])
[DELTA]t time difference (s)
[kappa] local curvature (1/m)
[lambda] wavelength of perturbation (m)
[micro] dynamic viscosity (Pa x s)
[upsilon] kinematic viscosity (m2/s)
[rho] density (kg/[m.sup.3])
[delta] interfacial tension (N/m)
[phi] phase of perturbation (-)
[omega] fluid domain (-)

Subscripts

1 lower flow
2 upper flow

Superscripts

T symbol of transpose (-)
* dimensionless symbol (-)