The effect of transient viscoelastic properties on interfacial instabilities in superposed pressure driven channel flows

J. Non-Newtonian Fluid Mech. 80 (1999) 217 – 249

The effect of transient viscoelastic properties on interfacial instabilities in superposed pressure driven channel flows Herambh K. Ganpule, Bamin Khomami * Department of Chemical Engineering and Materials Research Laboratory, Washington Uni6ersity, St. Louis 63130, Missouri, USA Received 15 December 1997

Abstract In a recent study, Ganpule and Khomami (submitted to J. Non-Newtonian Fluid Mech.) have shown that in order to accurately describe the experimentally observed interfacial instability phenomenon in superposed channel flow of viscoelastic fluids, a constitutive equation that can accurately depict not only the steady viscometric properties of the experimental test fluids, but also their transient viscoelastic properties must be used in the analysis. In the present study, the effect of differences in transient viscoelastic properties which can arise either due to the differences in the predictive capabilities of various constitutive models or from the presence of multiple modes of relaxation on the interfacial instabilities of the superposed pressure driven channel flows has been investigated. Specifically, a linear stability analysis is performed using nonlinear constitutive equations which predict identical steady viscometric properties but different transient viscoelastic properties. It is shown that different nonlinear constitutive equations give rise to the same mechanism of interfacial instability, but the boundaries of the neutral stability contours and the magnitudes of the growth/decay rates, especially at intermediate and shortwaves, are shifted due to the overshoots in the transient viscoelastic responses predicted by the constitutive equations. In addition, the effect of the presence of multiple modes of relaxation on interfacial stability is studied using single and multiple mode upper convected Maxwell (UCM) fluids and it is shown that pronounced differences in the intermediate and shortwave linear stability predictions arise due to the fact that the increase in the number of modes gives rise to additional fast as well as slow relaxation modes and the presence of these additional relaxation modes gives rise to differences in the transient viscoelastic response of the fluids in the absence of any overshoots. The effect of fluid inertia on the interfacial stability of viscoelastic liquids is examined and it is shown that at longwaves, inertia has a pronounced effect on the stability of the interface, whereas at shortwaves, elastic and viscous effects dominate. Furthermore, the mechanism of viscoelastic interfacial instabilities is studied by a careful examination of disturbance eigenfunctions as well as performing a disturbance energy analysis. The results indicate that the mechanism of viscoelastic interfacial instabilities can be described in terms of interaction of mechanisms of purely viscous and purely elastic instabilities. However, since more than one mechanism for the instability is at work, the disturbance energy analysis can not clearly distinguish between them due to the fact that the eigenfunctions used in the energy analysis contain the information regarding both viscous and elastic effects. Hence, the mechanism of the instability must be determined by a careful examination of disturbance eigenfunctions. © 1999 Elsevier Science B.V. All rights reserved. * Corresponding author: Fax: +1 314 935 7211. 0377-0257/99/$ - see front matter © 1999 Elsevier Science B.V. All rights reserved. PII S0377-0257(98)00085-8

218

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

Keywords: Interfacial instability; Viscoelastic fluids; Linear stability analysis; Energy analysis

1. Introduction The need for the development of specifically engineered structures of polymers has resulted in the proliferation of products that are manufactured using various coating and multilayer extrusion technologies. Likewise, the study of interfacial instabilities in multilayer flow of viscoelastic liquids has received significant attention in recent years. Aside from their technological significance, multilayer flows can be effectively used to investigate the effect of fluid elasticity on stability of generic viscoelastic flows with interfaces. Specifically, the base flow of multilayer plane Poiseuille flows can be obtained fairly easily, either analytically or numerically. More importantly, the presence of an interface in the flow allows one to experimentally observe the onset of the instability when disturbances are introduced in an otherwise steady flow, and to measure the growth or decay rates of these disturbances. This in turn allows one to assess the constitutive complexity required to accurately describe the onset conditions and the dynamics of the interfacial instabilities. Hence, it can be established whether or not the correct physics required to predict viscoelastic instabilities in the presence of interfaces can be incorporated into existing constitutive models. The phenomenon of interfacial instability in the multilayer flow of immiscible Newtonian and polymeric fluids has been studied by a number of investigators. Most studies have been theoretical in nature and have employed the method of linear stability analysis, with the exception of a few attempts at weakly nonlinear analyses [1,2]. The earlier studies used longwave asymptotic techniques to examine the linear stability of the flow. Using this technique, Yih [3] considered the plane Couette flow (PCF) and plane Poiseuille flow (PPF) of two superposed Newtonian fluids and showed that even in the limit of vanishing Reynolds number, viscosity stratification is sufficient to cause instability. Waters [4] and Khomami [5] studied the longwave instability of two superposed (modified) power law fluids in PCF and PPF respectively and showed that shear thinning has a pronounced effect on the stability of the interface, particularly if it occurs in the less viscous layer. The earlier efforts to identify the effect of elasticity stratification on the interfacial stability of plane channel flows [6,7], yielded misleading results due to an incorrect shear stress interfacial boundary condition. In more recent studies, Renardy [8], Chen [9], Su and Khomami [10] have used the correct boundary condition and shown that even in the absence of viscosity and density mismatches, elasticity stratification is sufficient to cause interfacial instabilities at small or vanishing Reynolds numbers (i.e. a purely elastic instability). Moreover, they observed that elasticity stratification tends to destabilize the interface when the less elastic fluid is the majority component (i.e. occupies more than half the channel). Since the asymptotic methods are applicable only in the limit of longwave and shortwave disturbances, Yiantsios and Higgins [11] used the compound matrix numerical method to study the stability behavior of the two layer PPF of Newtonian fluids to disturbances of arbitrary wavelength and showed that longwave disturbances are stabilized by the presence of a thin layer of less viscous fluid (the ‘thin layer effect’) and intermediate and short-wave disturbances are

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

219

stabilized by interfacial tension. Su and Khomami [12–14] used spectral based (namely spectral-tau and pseudospectral) techniques to analyze the stability of superposed PPF of Newtonian and viscoelastic fluids by using a number of quasilinear constitutive equations such as the Upper convected Maxwell, Oldroyd-B, and modified Oldroyd-B models. They demonstrated that the relative contribution from elasticity stratification is at least on the same order of magnitude as the contribution from viscosity stratification [10]. They also showed that elasticity stratification can affect the interfacial stability at all disturbance wavelengths, however this effect is more pronounced at intermediate and short wavelengths. It was also shown that waves that have a dimensionless wavenumber (when length is nondimensionalized by the length of the more viscous/elastic layer) of around unity are the most dangerous waves. Most importantly, the authors have compared the predictions of their linear stability analyses with the experimentally determined neutral stability contours as well as growth/decay rate measurements of Wilson and Khomami [15 – 17]. These comparisons indicate that one-mode quasilinear viscoelastic models can predict the stability of the polymeric test systems at most qualitatively. In particular, they can accurately predict the critical layer depth ratio, and the most dangerous wavenumber quantitatively, but can only predict the neutral stability contour and the growth rate behavior qualitatively. Based on their studies, the authors provided ample evidence that to obtain a quantitative agreement between the experimental findings and theoretical analyses, one must use nonlinear and/or multimode constitutive models which can accurately predict the shear rate dependence of viscosity as well as elasticity exhibited by the polymeric fluids utilized in the experiments. In a recent theoretical study [18], we have conducted such an investigation. The constitutive equations used were Modified Oldroyd-B, single mode Giesekus, Modified Phan-Thien Tanner (MPTT), and the multimode Giesekus (MMG). All of the models can predict the shear rate dependence of viscosity as well as the elasticity of the test fluids used in the experiments of Wilson and Khomami [15,16] and Khomami and Ranjabaran [19] but with a different degree of flexibility and accuracy. The results of the theoretical studies were compared with the experimental studies of Wilson and Khomami [15,16] and Khomami and Ranjabaran [19] both of which used polypropylene (PP) and high density polyethylene (HDPE) melts as test fluids. Comparisons were made based on the neutral stability contours in the plane of dimensionless wavenumber (a) and layer depth ratio (e), and on growth rate versus dimensionless wavenumber plots. The important conclusions of this study are as follows. Below a certain depth ratio called the critical layer depth ratio, the interface is stable to disturbances of all wavelengths. All four models predict the critical layer depth ratio accurately. At subcritical layer depth ratios, viscosity stratification stabilizes the interface at longwaves (this is nothing but the ‘thin layer effect’) and elasticity stratification at intermediate waves. At supercritical layer depth ratios, viscosity stratification is destabilizing in the limit of longwaves. Elasticity stratification has a stabilizing effect at intermediate waves, especially at high values of the depth ratio. In this region, the 4-mode Giesekus model is found to produce quantitative agreement with the experimental findings whereas the MPTT model does not. Overall, the MMG model most accurately predicts the experimentally observed instability in terms of the neutral stability contour. The MPTT model, despite possessing accurate steady state viscometric properties, predicts the instability only qualitatively.

220

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

For most bulk flow instabilities, such a favorable comparison between experimental and theoretical (i.e. with the MMG model) findings in predicting the critical onset conditions would lead one to conclude that the constitutive equation can accurately describe the flow behavior of the polymeric fluids. However, as mentioned earlier, the presence of an interface allows one not only to actually observe the onset of the instability, but also to measure the growth and decay rates of the instability. The comparison of the theoretically and experimentally obtained growth rate plots provide an extremely stringent test of the ability of the constitutive equation to accurately depict the temporal and spatial feature of the instability. In the above study, such a comparison was performed. The first important observation is that all of the models correctly predict the most dangerous wavenumber ae1. Moreover, the constitutive equations which accurately depict the shear rate dependence of viscosity as well as elasticity (i.e. MPTT and MMG) do indeed produce a much more accurate description of the instability phenomenon as opposed to those which do not (i.e. Modified Oldroyd-B and single mode Giesekus). Hence, overall, the 4-mode Giesekus model accurately predicts the interfacial instability phenomenon in the terms of the critical layer depth ratio, the most dangerous wavenumber, the neutral stability contour, and the growth rate behavior. This comparison emphasizes two predictions of the linear stability analysis using the above models. Firstly, although the MPTT model can accurately depict the steady shear properties of the polymeric test fluids, it can not predict the instability phenomenon (especially the growth rate behavior) with quantitative accuracy. Secondly, in spite of having identical behavior in steady shear, the MMG and MPTT models give rise to different linear stability predictions, both in terms of the neutral stability contour and the growth rates. In light of these observations, in this study, we have examined the effect of transient viscoelastic properties of the fluids on the stability of the interface. In particular, we have examined the effect of inherent differences in the transient response of the Giesekus and MPTT constitutive equations on the linear stability predictions. Furthermore, we have also examined the effect of the presence of multiple modes of relaxation that are known to give rise to different transient rheological properties on the stability of the interface. Note that the presence of second normal stresses does not influence this 2-D analysis. However, the effect of these stresses which manifest themselves only in a 3-D analysis, will be a topic of a future study. Further, the difference in the elongational behavior of the two models is not expected to significantly influence the interfacial stability since we are considering the linear stability of a shearing flow. In the same previous study [18], we have also investigated the mechanisms of longwave and shortwave viscous and purely elastic interfacial instabilities using the method of disturbance-energy analysis as well as providing physical arguments based on the longwave and short-wave disturbance eigenfunction behaviors in the respective cases. The mechanism of purely viscous shortwave interfacial instability was found to be interfacial friction (i.e. interfacial velocity/vorticity mismatch) which is caused by viscosity stratification. The mechanism of purely viscous longwave instability was found to be due to the effect of inertia (represented by the Reynolds stress term of the disturbance-energy equation) coupled with the interfacial velocity/vorticity mismatch. The mechanism of purely elastic shortwave as well as longwave instability was found to be due to the coupling between the jump in the base flow normal stresses across the interface and the perturbation velocity field in the bulk. These mechanisms are valid for the pressure driven as well as drag driven flows. In this study, we have studied the effect of the presence of

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

221

viscosity as well as elasticity stratification (i.e. viscoelastic instabilities) on interfacial instabilities and their mechanisms. The objective is to check whether the physical arguments and energy analysis employed to explain the purely viscous and purely elastic instabilities can be extended to the viscoelastic flow instabilities. We have also examined the effect of dynamic viscoelastic properties on the mechanism of interfacial instabilities using the White-Metzner (WM) and single mode Giesekus constitutive models.

2. Problem formulation

2.1. Go6erning equations Fig. 1 depicts the flow geometry considered in this study. We consider the plane superposed pressure driven flow of two distinct polymer fluids through a parallel channel flow geometry. The governing equations are the equations of conservation of mass and momentum, and the constitutive equation. The equations of continuity and motion can be expressed as 9( · U( = 0,

(1)

DU( = − 9( P − 9 · t˜ +rg¯, Dt

r

(2)

in which U( , P( , g¯, and t˜ denote the velocity vector, isotropic pressure, gravitational acceleration and the deviatoric stress tensor, respectively. The constitutive equations used in this study are the single and multimode versions of the UCM, the WM, and the Giesekus model. In these equations, the total stress consists of a polymeric contribution and a solvent contribution t˜¯ = t˜¯ p + t˜¯ s

(3)

The solvent contribution is purely viscous, t˜¯ s = − hsg˜; ,

(4)

where g˜; = 9U( +9U( T

(5)

is the rate of strain tensor. In the WM equation, the polymeric part of the stress satisfies the following equation:

Fig. 1. Schematic of flow channel.

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

222

t˜¯ p + l1(g; )t˜¯ p(1) = − hp(g; )g˜; .

(6)

For the Giesekus equation, the polymeric stress satisfies t˜¯ p + l1,0t˜¯ p(1) −

al1,0 t˜¯ p · t˜¯ p = − hp,0g˜; . hp,0

(7)

Here, l1 is the dominant relaxation time of the polymer, a is the adjustable parameter of the models, hp is the shear rate dependent polymeric viscosity, hs is the solvent viscosity, subscript 1 denotes the upper-corrected derivative and the subscript 0 denotes the zero shear rate value of the parameters. For the Giesekus model, a is the mobility parameter associated with either an anisotropic Brownian motion or an anisotropic hydrodynamic drag on the polymer molecules [20]. For the WM equation, hp and l1 are shear rate dependent and are defined as: hpg; = hp,0G1(g; ), C1(g; ) = C1,0G2(g; ),

(8)

C1(g; ) = l1,0G3(g; ), 2hp(g; )

(9)

l1(g; ) = where

G1(g; ) =[1 + (k1g; )2](n1 − 1)/2, G2(g; ) = [1+(k2g; )2](n2 − 2)/2, G3(g; ) =

G2(g; ) . G1(g; )

C1 is the shear rate dependent first normal stress coefficient, and (g; ) is the second invariant of the rate of strain tensor defined as g; = (tr(g˜; · g˜; ))1/2. To study the effect of the presence of multiple modes of relaxation on the stability of the interface, we have used the UCM model which is a simplification of the Giesekus equation with hs = 0, and a = 0. In such multimode models, the total polymeric stress consists of contributions from each individual mode, M

tp = % t m p,

(10)

m=1

where M is the total number of modes. For the multimode UCM model, each of the t m p satisfies the Giesekus equation (Eq. (7)) with a = 0. The following set of dimensionless variables is introduced to nondimensionalize the above equations:

 

x1 x2 P( t¯ ij t( U0 U( k , y= , P= , tij = , t= , Uk = , h0U0 d1 d1 h0U0 d1 U0 d1 d1 hk rk (l1)k (l2)k d2 , bk = , o= , (Rm)k = , rk = , uk = d1 h1 r1 (l1)1 (l1)k h1g; int,1d1 (r2 − 1)gd1 Ca = , Fr= s U 20 rkU0d1 (l1,0)kU0 (l2,0)kU0 Rek = , Wek = , Lk = , h1 d1 d1 x=

(11)

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

223

where U0 is the interfacial velocity, g; int,1 is the base flow shear rate at the interface in layer 1, subscript k = 1, 2 where 1 and 2 denote the upper and lower layers respectively, s is the interfacial tension, g is the gravitational acceleration and l2 =l1(hs/hs +hp0) is the retardation time. For the base flow, the governing equations are solved numerically with the appropriate boundary conditions, namely no slip at the solid walls, and continuity of velocity and of shear stress at the interface. The numerical technique used is a multidomain spectral technique that is discussed briefly in the next section. To examine the interfacial stability of the flow, equations governing the evolution of small normal mode perturbations applied to the base flow are formed. The perturbation quantities are assumed to have an exponential time dependence and a periodic spatial dependence of the form: qˆ (x, y, u) = q(y) exp[ia(x − cu)],

(12)

where, q denotes an arbitrary perturbation variable, q(y) is the amplitude of the disturbance, a is the dimensionless wave number of the disturbance in the flow direction x, u is the dimensionless time, and c= cr + ici is the dimensionless complex wavespeed. The stability of the flow is governed by the sign of the imaginary part of c: if it is positive, the amplitude of the disturbance increases with time and the flow is unstable; if negative, the flow is stable; and if zero, the flow is neutrally stable. Then, the governing equations are linearized with respect to perturbation quantities and the equation of motion, continuity, and the constitutive equations are recast in terms of the stream function. Then, the pressure terms in the x and y components of the equation of motion are eliminated by cross differentiation and the stability governing equations are obtained. The boundary conditions for the perturbed flow problem are as follows: No slip at the solid walls: f1(1) = f%1(1)= 0,

(13)

f2( − e) =f%2(− e)= 0,

(14)

Continuity of the x and y components of velocity at the interface:

f = 0,

f% +

(15)

* *

dU f%1 = 0, c − U1 dy

(16)

balance of stresses normal to the interface:

Fxy =

iaf1

txx − tyy , c −U1

(17)

and balance of stresses tangential to the interface:

* 

ia Re (U − c)f% −

dU f dy

n*

+ ia Fxx −Fyy + F%xy =

ia(F+a 2)f1 , c−U1

(18)

224

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

where S = s/r1d1U 20 and F= (rk − 1)gd1/U 20 are the dimensionless groups expressing the effects of interfacial tension s and gravitational acceleration g, respectively, and A means the jump in quantity A across the interface. In the present study, we have neglected the effects of interfacial tension and density stratification, therefore S=F=0.

2.2. Solution procedure The base flow problem constitutes a nonhomogenous boundary value problem (BVP), the constant pressure gradient term being the nonhomogenous term. The use of nonlinear constitutive equations makes is impossible to obtain an analytical solution to the base flow problem for arbitrary values of the model parameters. To overcome this problem, we have used a multidomain spectral method to solve the base flow problem. Details of the method are given elsewhere [18], hence not being reproduced here. The perturbation flow problem constitutes an eigenvalue problem (EVP) with the complex wavespeed c as the eigenvalue and the disturbance amplitudes f, Fxy, Fxx, and Fyy as eigenfunctions. In this study, we first used longwave asymptotic analysis to examine the stability of the interface to longwave disturbances. This analysis not only provides guidelines about the critical parameters, but also provides an important check to the numerical predictions. However, as mentioned earlier, the dominant mode of interfacial instability can be due to disturbances of arbitrary wavelength. The EVP constituted by the stability governing equations cannot be solved asymptotically for dimensionless wavenumber 0(1), and hence needs to be solved numerically. The numerical method used is again the multidomain spectral tau method. The stability governing equations and the boundary conditions constitute a generalized EVP A−cB=0. The QR algorithm is used to solve for the eigenvalues. Again, the details of the method are given in our earlier publication [18], hence not being reproduced here. It is worth noting that the nonlinearity of the constitutive equations gives rise to steep shear and normal stress boundary layers, thus necessitating the use of multidomain spectral methods as opposed to singledomain spectral techniques used in most prior studies to obtain convergent solutions to the base flow as well as the perturbation flow problem. The numerical technique is discussed in much more detail in [18].

3. Results and discussions The effect of dynamic rheological properties on the stability of the interface is examined using two approaches. First, the effect of overshoots in rheological properties (such as the viscosity and first normal stress difference during the start-up of steady shear) is examined by performing a linear stability and disturbance-energy analysis using nonlinear constitutive equations such as the single mode Giesekus and WM. Secondly, the effect of the presence of multiple modes of relaxation on the stability of the interface is examined by studying the differences in their dynamic rheological properties using a quasilinear constitutive model, namely the UCM model.

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

225

Table 1 Parameters of the four systems GG1, GG2, GG3, and GG4 of Giesekus fluids System-layer

hp

l1

a

GG1-layer GG1-layer GG2-layer GG2-layer GG3-layer GG3-layer GG4-layer GG4-layer

1.0 0.5 1.0 0.5 1.0 0.5 1.0 0.5

0.4 0.1 0.4 0.1 0.6 0.15 0.2 0.09

0.1 0.05 0.4 0.3 0.4 0.3 0.4 0.3

1 2 1 2 1 2 1 2

3.1. Effect of o6ershoots in dynamic rheological properties on the stability of the interface To study the effect of differences in the constitutive behaviors of WM and Giesekus models on the stability of the interface, we have performed a linear stability analysis of the superposed flow of these fluids. We have studied four systems of a superposed flow of two Giesekus fluids. They are called GG1, GG2, GG3, and GG4, and the parameters of the Giesekus model used in these systems are shown in Table 1. Note that GG1, GG2, and GG3 are progressively more nonlinear in nature due to the increasing values of the l1 and/or a parameters. Also, note that layer 2 is the least viscous as well as being the least elastic layer. The Carreau type of shear rate dependence used for the viscosity and first normal stress coefficient (C1) of the WM model provides the flexibility to exactly fit the parameters of the WM model so that it gives a steady shear behavior identical to that of a Giesekus fluid. We have fitted the parameters of the WM fluid for each of the GG1, GG2, GG3, and GG4 systems. The corresponding systems of the superposed flow of WM fluids are denoted by WW1, WW2, WW3, and WW4 respectively. The parameters of the WM model used in these systems are shown in Table 2. It should be noted that GG1, GG2, GG3, and the corresponding WM systems are creeping flow systems (i.e. Re= 0), whereas GG4 and WW4 are not. Table 2 Parameters of the four systems WW1, WW2, WW3, and WW4 of WM fluids System-layer

hp, 0

k1

n1

C1, 0

k2

n2

WW1-layer WW1-layer WW2-layer WW2-layer WW3-layer WW3-layer WW4-layer WW4-layer

1.0 0.5 1.0 0.5 1.0 0.5 1.0 0.5

0.27 0.05 0.47 0.12 0.71 0.17 0.24 0.10

0.13 0.17 0.08 0.14 0.08 0.14 0.08 0.14

0.8 0.1 0.8 0.1 1.2 0.15 0.4 0.09

0.27 0.05 0.47 0.12 0.71 0.17 0.24 0.10

0.65 0.68 0.59 0.66 0.59 0.66 0.59 0.66

1 2 1 2 1 2 1 2

226

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

Fig. 2. Steady state viscometric properties of the Giesekus and matched WM fluids from the GG1, WW1, GG2, WW2, GG3, WW3, GG4, and WW4 systems. Solid lines correspond to the Giesekus fluids and the dashed ones to the WM.

The steady state viscosity and first normal stress coefficient behavior of the fluids in these systems is shown in Fig. 2. As seen from these figures, the steady viscometric properties of the WM fluids are virtually identical to those of the corresponding Giesekus fluids. Before comparing the results of the stability analysis, it was checked to ensure that the flow described by the Giesekus and WM constitutive models with identical steady rheological properties gives rise to identical base flow velocity and stress profiles. A representative base flow velocity, shear stress, and first normal stress difference (N1) profiles is shown in Fig. 3. As expected, the two

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

227

models do give rise to indistinguishable velocity and stress profiles. Note that the Giesekus model predicts a nonzero second normal stress difference in planar channel flows, whereas the WM model does not. However, the present study is only concerned with a 2-D stability analysis. Hence, the presence of second normal stresses does not affect the linear stability analysis [21]. To capture the effects of second normal stresses on interfacial stability, a full 3-D stability analysis is required. Such an analysis is underway and will be the subject of a future publication. Fig. 4 shows the results of the linear stability analysis in terms of a neutral stability contour in the plane of dimensionless wavenumber, a, and layer depth ratio, e for the GG1 and WW1 systems. As seen from this figure, the two neutral stability contours are qualitatively similar. The neutrally stable line of constant depth ratio corresponds to the critical depth ratio, ec. It is the depth ratio at which the shear rate across the interface is continuous (i.e. the maximum in the

Fig. 3. Base flow velocity, shear stress, and normal stress profiles for the GG2 (solid lines) and WW2 (dashed lines) systems. e =0.5.

228

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

Fig. 4. Neutral stability contour for the GG1 (solid lines) and WW1 (dashed lines) systems.

velocity profile occurs at the interface); when eBec the maximum in the velocity profile occurs in the more viscous layer, and when e\ ec it occurs in the less viscous layer. As is clear from this figure, both models predict the same value of the critical depth ratio. This is due to the fact that the value ec is governed by viscosity stratification, and since both models predict identical steady state viscosity behavior, they yield the same value of ec. In the intermediate and shortwave regime, the stability predictions of the two models are different, especially at subcritical layer depth ratios. This fact can be more clearly seen from the plots of growth/decay rate versus layer depth ratio in the long, intermediate, and shortwave regimes. This is shown in Fig. 5. Longwave linear stability predictions (Fig. 5(a)) of the Giesekus and WM models differ only slightly. This is because at longwaves, the longest relaxation time plays the most significant role in determining the stability of the interface. Hence, it is not surprising that both the models give rise to very similar predictions. In the intermediate and shortwave regimes, however, the growth/decay rates predicted by the two models are significantly different. In particular, in the limit of shortwaves, growth rates predicted by the Giesekus model are always smaller than those predicted by WM. This indicates that the nonlinearity of the Giesekus model results in a decreased driving force of the shortwave instability. The physical mechanism responsible for these differences will be discussed shortly. As mentioned earlier, the nonlinearity of the GG2 and GG3 systems is progressively more pronounced. Hence, the comparison of the linear stability predictions from these, and the corresponding WM systems (i.e. WW2 and WW3) shows the same trends as discussed above for GGl/WW1 systems; but the differences between the corresponding Giesekus and WM systems become more pronounced. Figs. 6 and 8 show the neutral stability contours for the GG2/WW2 and GG3/WW3 systems respectively, whereas Figs. 7 and 9 show the corresponding growth/decay rate versus layer depth ratio plots for long, intermediate, and shortwave disturbances. Clearly, increasing the severity of the nonlinearity leads to more and more pronounced differences in the GG2/WW2 and GG3/WW3 systems (as compared to GGl/WW1 system), both in terms of the boundaries of the neutral stability contours as well as the magnitudes of the growth/decay rates, especially at intermediate and shortwaves. Again, note that the longwave behavior is virtually identical due to the reason explained above.

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

229

Since the steady state rheological behavior of the two fluids in each case is identical, these differences in the linear stability predictions (i.e. at intermediate and shortwaves) must arise due to the differences in the transient rheological response of the fluids. These differences in stability predictions and their dependence on dynamic viscoelastic properties can be best understood by examining the physical mechanism of the shortwave instability. In the intermediate and shortwave regime, elasticity has a pronounced effect on the stability of the interface. This is due to the fact that in the normal stress interfacial boundary condition (Eq. (17)), the term representing the first normal stress difference is multiplied by the wavenumber. Further, as shown by Hinch [23], inertia is necessary to cause the instability due to viscosity stratification

Fig. 5. Growth rate versus dimensionless wavenumber for the GG1 (solid lines) and WW1 (dashed lines) systems for: (a) longwave, a=0.01; (b) intermediate wave, a= 8; and (c) shortwave disturbances, a = 100.

230

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

Fig. 6. Neutral stability contour for the GG2 (solid lines) and WW2 (dashed lines) systems.

alone. Since the flows under consideration are creeping flows, the instability is caused by the disturbances generated due to elasticity stratification, even though viscosity stratification does play an important role in the dynamics of the instability. The mechanism of purely elastic shortwave instability has been studied by Ganpule and Khomami [18] and is reproduced here to illustrate the effect of transient viscoelastic properties on interfacial instability. In the limit of shortwaves, the balance of normal stresses across the interface is shown in Fig. 10. The jump in the first normal stresses across the interface results in a net force being exerted on the interface along the flow direction. This causes the interface to move in the direction of the applied force. This movement of the interface sets up a perturbation velocity in the x-direction which is maximum at the interface and, in the limit of shortwaves, decays exponentially fast away from the interface. These disturbance velocities can be converted into disturbance vorticities. At a position midway between the crest and the trough (x= p/2a), vorticity in the less elastic layer (layer 1) is negative (clockwise) and that in the more elastic layer (layer 2) is positive (counter-clockwise). Thus the vorticity in layer 1 leads to depletion of fluid above the crest and is destabilizing whereas that in layer 2 leads to depletion of fluid below the crest and hence is stabilizing. The rate of decay of perturbation velocity, and hence the magnitude of the perturbation vorticity is larger in the less elastic layer due to the coupling of the equation of motion and the constitutive equation. Therefore, the destabilizing vorticity field set up in layer 1 dominates over the stabilizing vorticity field set up in layer 2, and the interface becomes unstable. Note that the growth rate of the instability is determined by the relative magnitudes of the stabilizing and destabilizing vorticities in the two layers. Ganpule and Khomami [18] have examined this mechanism for the case of purely elastic (i.e. in the absence of viscosity stratification) creeping flow shortwave instability. In the presence of viscosity and elasticity stratification and in the absence of inertia, the driving force behind the instability is still the jump in the base flow normal stresses across the interface, however, the relative magnitudes of the perturbation vorticities are governed by the elasticity as well as viscosity differences between the two fluids. This is true for the case in consideration, namely for the GG and WW systems. To examine the possibility that the differences in linear stability predictions arise due to the differences in the transient viscoelastic properties of the fluids, the

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

231

viscosity and first normal stress coefficients of the two fluids during the start up of steady shear flow were determined. Figs. 11–13 show the plots of the transient viscosity (h + ) and transient first normal stress coefficient C1+ as a function of time for the GG1/WW1, GG2/WW2, and GG3/WW3 systems respectively. As seen from these figures, the transient behaviors of the Giesekus and WM fluids are almost identical at low shear rates. However, at high enough shear rates, the Giesekus model shows an overshoot for h + as well as C1+ , whereas for the WM fluid, these parameters monotonically approach their respective steady state values. For all of the systems, the magnitude of the overshoots increases with increasing shear rate. The overall magnitudes of the overshoots (i.e. for all shear rates) increase from GG1 to GG2 to GG3 systems. At long time, the transient properties of the Giesekus and WM fluids attain their respective steady state values which are equal due to the matched steady state viscometric properties. In the limit of longwaves, the longest relaxation time plays the most significant role in determining the stability of the system. Since the viscosities and first normal stress coefficients of the WM and Giesekus fluids at long times are equal, it is not surprising that the two models predict almost identical longwave behavior. At shortwaves, however, the short time scale

Fig. 7. Growth rate versus dimensionless wavenumber for the GG2 (solid lines) and WW2 (dashed lines) systems for: (a) longwave, a=0.01; (b) intermediate wave, a= 8; and (c) shortwave disturbances, a = 100.

232

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

Fig. 8. Neutral stability contour for the GG3 (solid lines) and WW3 (dashed lines) systems.

behavior is dominant. As seen from Figs. 11 and 12, and 13, the short and intermediate time behaviors of the Giesekus and WM models show pronounced differences. In particular, due to the overshoots in the dynamic viscosity and normal stress, magnitudes of the perturbation vorticities in both the layers are reduced i.e. the magnitude of the stabilizing vorticity in the more elastic layer as well as that of the more significant destabilizing vorticity in the less elastic layer is reduced. Fig. 14 shows the perturbation vorticity profile in the channel for the GG3 and WW3 systems. Clearly, the magnitudes of the perturbation vorticities for the GG3 system are smaller than those for the WW3 system. But since the magnitude of shear rate in the less viscous layer is more than that in the more viscous layer, the overshoots in viscosity as well as normal stress are more pronounced in the less viscous and less elastic layer. Thus, the reduction in the more significant destabilizing vorticity in the less viscous/less elastic layer is more pronounced than that in the more viscous/more elastic layer. Hence the magnitude of the growth rates for the GG systems is smaller than those for the WW systems. As seen from Figs. 11 and 12and Fig. 13, as the nonlinearity of the Giesekus model becomes more pronounced, so do the overshoots at shorter time scales (i.e. at intermediate and shortwaves). Therefore, the differences in magnitudes of the growth rates predicted by the Giesekus and WM models increase from GG1/WW1 to GG2/WW2 to GG3/WW3 systems. These differences in the shortwave growth rate behavior are most significant at depth ratios significantly larger or smaller than the critical depth ratio. This is because at the critical depth ratio, the shear rate at the interface is zero in both the layers. At values of depth ratio close to ecrit, the shear rates at the interface are small in magnitude in both the layers. Since the transient viscoelastic properties of the Giesekus and WM fluids are similar at small values of shear rates (Fig. 13), the differences in their linear stability predictions are small when eecrit. Moreover, the differences in the growth rate behavior are more pronounced at subcritical depth ratios. This is because at these depth ratios, the less viscous fluid is present as a thin layer, and hence the shear rates in that fluid are larger in magnitude thus leading to larger overshoots in viscosity and first normal stress difference. Since the longwave stability of the flow is predominantly governed by the longest relaxation behavior time (i.e. the steady state behavior), the mechanism of the longwave viscoelastic instability is expected to be a combination of those proposed by Hinch et al. [24] for purely elastic instability

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

233

and by Ganpule and Khomami [18] for purely viscous instability. Again, since the flow under consideration is a creeping flow, the instability is caused by the disturbances generated due to elasticity stratification, even though viscosity stratification does play an important role in the dynamics of the instability. Transition from longwave to shortwave behavior, i.e. the behavior at intermediate wavelengths, is governed by the transient viscoelastic behavior at intermediate time scales and the manner in which it affects the shift from the longwave instability mechanism to the shortwave one. Fig. 15 shows the plot of growth/decay rates for the GG2/WW2 system at a depth ratio of 0.2. Since the differences in the transient properties of the Giesekus and WM fluids are largest at intermediate times (Fig. 13), the differences in the growth/decay rates are the largest at intermediate waves. This is of significant practical interest since the longwave disturbances can be stabilized using the ‘thin layer effect’, shortwave disturbances can be stabilized by interfacial tension, and the stability of the interface to intermediate wave disturbances is often the key in establishing a stable coextrusion processing window.

Fig. 9. Growth rate versus dimensionless wavenumber for the GG3 (solid lines) and WW3 (dashed lines) systems for: (a) longwave, a=0.01; (b) intermediate wave, a= 8; and (c) shortwave disturbances, a = 100.

234

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

Fig. 10. Schematic of the mechanism of shortwave purely elastic interfacial instability. Reproduced from [18].

3.2. Effect of inertia As mentioned earlier, it has been shown by Hinch [23] that inertia is necessary to cause interfacial instability due to viscosity stratification alone. Thus, if inertia is present in the above GG/WW systems, one expects to observe the instability which is driven by the jump in viscosity as well as normal stresses across the interface. To study the effect of inertia on the stability of the interface at shortwaves, we have performed a linear stability analysis of the inertial viscoelastic system of two Giesekus fluids, namely GG4, the corresponding system of two WM fluids, namely WW4 and a system of two inelastic fluids called VV4 which have the same viscosity behavior as the GG4 and WW4 systems. These systems were specifically selected to study the effect of inertia on interfacial instability by significantly reducing the jump in the base flow normal stresses across the interface. Therefore, the elastic forces in these systems are small in comparison to the previous GG/WVV systems, hence inertial effects are amplified. Figs. 16 and 17 show the plots of growth rate versus Reynolds number for the WW4, and VV4 systems for shortwave and longwave disturbances respectively. The linear stability predictions of the GG4 systems are not shown in these figures because they are exactly identical to those of the WW4 system. The reason for this is that the elasticities of the fluids were chosen to be small in order to observe the effect of inertia. Since the nonlinearity of the Giesekus model is governed by the value of a as well as l1, small elasticities imply that the nonlinearity of the Giesekus model – which gives rise to the overshoots of the viscosities and normal stresses during the start up of steady shear–is very weak. Therefore, the WW4 and GG4 systems produce the same linear stability predictions. As seen from Fig. 16, the growth rate of the shortwave instability is not very sensitive to changes in Reynolds number. This is not surprising because at shortwaves, the effective length scale of the flow problem is the disturbance wavelength which is infinitesimally small. Therefore, the effective Reynolds number is very small and viscosity and elasticity dominate inertia. As shown by Hinch [23] and Ganpule and Khomami [18], the mechanism of purely viscous as well as purely elastic shortwave instability requires the diffusion of disturbance vorticity. The plot of laplacian of the disturbance vorticity as a function of the Reynolds number for the shortwave

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

235

instability is shown in Fig. 18. As seen from this figure, the laplacian of the disturbance vorticity and hence the disturbance vorticity diffusion is not very sensitive to changes in Reynolds number. Therefore, Reynolds number does not exhibit a pronounced effect on the shortwave interfacial instability. As seen from Fig. 17, inertia has a significant stabilizing or destabilizing effect on the interface at longwaves depending on the depth ratio. In particular, inertia enhances the stabilizing or destabilizing effect of viscosity and/or elasticity stratification. This pronounced effect of inertia on the longwave stability is not surprising since at longwaves, the effective length scale of the flow problem is the channel size. Therefore, the effective Reynolds number is large in magnitude. The dependence of vorticity diffusion on inertia can be best understood by the examination of the dimensionless vorticity transport equation: 92v ¯ = Re





(v ¯ + u¯ · 9( v ¯ −v ¯ · 9( u¯ . (u

(19)

Fig. 11. Dynamic viscoelastic properties for the GG1 (solid lines) and WW1 (dashed lines) systems during the start up of steady shear flow. Dynamic viscosities normalized by the respective steady state values,h + (t, g; )/h(g; ) (a) Layer 1, (c) Layer 2; Dynamic first normal stress coefficients normalized by the respective steady state values, C1+ (t, g; )/ C1(g; ) (b) Layer 1, (d) Layer 2.

236

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

Fig. 12. Dynamic viscoelastic properties for the GG2 (solid lines) and WW2 (dashed lines) systems during the start up of steady shear flow. Dynamic viscosities normalized by the respective steady state values, h + (t, g; )/h(g; ) (a) Layer 1, (c) Layer 2; Dynamic first normal stress coefficients normalized by the respective steady state values, C1+ (t, g; )/ C1(g; ) (b) Layer 1, (d) Layer 2.

Fig. 19 shows the plot of laplacian of the disturbance vorticity as a function of the Reynolds number for longwave disturbances. Clearly, for longwave disturbances, disturbance vorticity diffusion shows a strong dependence on Reynolds number. Moreover, the vorticity diffusion in the less viscous layer has a stronger dependence on the Reynolds number. This is due to the fact that for a fixed depth ratio, the effective Reynolds number in the less viscous layer is larger due to the smaller viscosity. Ganpule and Khomami [18] have shown that this difference in the vorticity diffusion between the more viscous and less viscous layer is responsible for the longwave stability/instability and gives rise to the ‘thin layer effect’. Figs. 16 and 17 again bear out the fact that the effect of elasticity stratification is more pronounced at shortwaves. For the longwave stability/instability, the difference between the VV4 and WW4 predictions is very small; i.e. the addition of elasticity stratification (i.e. WW4 system) to the purely viscous system (i.e. VV4) leads to a very small change in the growth/decay rate behavior. On the other hand, for the shortwave instability the addition of elasticity stratification to the purely viscous system leads to a large change in the magnitudes of the growth rates of the instability.

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

237

3.3. Effect of the presence of multiple modes of relaxation on the stability of the interface As mentioned earlier, we have also analyzed the effect of the presence of multiple modes of relaxation on interfacial instability. Due to the polydispersity of the polymeric materials, the presence of multiple modes of relaxation and its effect on the stability of superposed multilayer flows is of much practical interest. To study the effect of the presence of multiple modes of relaxation on the stability of the interface, a linear stability analysis of two superposed UCM fluids was performed. In layer 2, a single mode UCM fluid is used, and in layer 1, various fluids with different numbers of modes are utilized. The steady shear properties of the fluids are summarized in Table 3. In determining the parameters of a given fluid for different number of modes, it is required that their steady shear properties remain unchanged even if the number of modes is changed [22]. This implies the presence of additional fast-relaxing and slow-relaxing modes when the number of modes is increased. Moreover, the steady shear viscosities of the fluids in layers 1 and 2 are identical. Hence, any instability in this system is only due to elasticity stratification.

Fig. 13. Dynamic viscoelastic properties for the GG3 (solid lines) and WW3 (dashed lines) systems during the start up of steady shear flow. Dynamic viscosities normalized by the respective steady state values, h + (t, g; )/h(g; ) (a) Layer 1, (c) Layer 2; Dynamic first normal stress coefficients normalized by the respective steady state values, C1+ (t, g; )/ C1(g; ) (b) Layer 1, (d) Layer 2.

238

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

Fig. 14. Disturbance vorticity profiles in the channel for the GG2 (solid line) and WW2 (dashed line) systems for shortwave disturbances, e=0.2, a= 80.

The neutral stability contours for all three cases are shown in Fig. 20. As expected, the presence of multiple modes of relaxation clearly changes the neutral stability contour. The effect is more clearly seen in the plots of growth/decay rates versus layer depth ratio at long, intermediate, and shortwaves, shown in Fig. 21. This effect of the presence of multiple modes of relaxation on the stability of the interface can be explained by the differences in the transient behaviors. The h + and C1+ data for one, two and four mode fluids used in layer 1 are shown

Fig. 15. Growth/decay rate versus dimensionless wavenumber for the GG1 (solid line) and WW1 (dashed line) systems, e =0.3.

Fig. 16. Growth/decay rate of the shortwave instability (a= 80) as a function of Reynolds number for the VV4 (solid lines) and WW4 (dashed lines) systems, (a) e=0.3, (b) e = 3.0.

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

239

Fig. 17. Growth/decay rate of the longwave instability (a= 0.01) as a function of Reynolds number for the VV4 (solid lines) and WW4 (dashed lines) systems, (a) e =0.3, (b) e= 3.0.

in Fig. 22. It can be seen in this figure that as the number of modes increases, the transient response of the fluid in the startup of steady shear flow becomes faster at small times and more gradual at long times. This is because the increase in the number of modes implies presence of additional fast-relaxing as well as slow-relaxing modes.

Fig. 18. Laplacian of the disturbance vorticity as a function of Reynolds number for the WW4 system for shortwave disturbances, a=80, e =0.3.

Fig. 19. Laplacian of the disturbance vorticity as a function of Reynolds number for the WW4 system for longwave disturbances, a=0.01, e=0.3.

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

240

Table 3 Parameters of the UCM fluids used to study the effect of the presence of multiple modes of relaxation on the stability of the interface

Layer 1

Layer 2

Number of modes

l1

hp

Shp, i

Shp, il1, i

1 2 4 1

0.1 1.0, 0.01 10.0, 0.46, 0.022, 0.001 0.001

1.0 0.09, 0.91 0.005, 0.076, 0.67, 0.25 1.0

1.0 1.0 1.0 1.0

0.1 0.1 0.1 0.001

Again, the longwave behavior of the three fluids is identical since their viscosities and first normal stress coefficients at long times are identical. The above proposed instability mechanism can be used to explain the shortwave linear stability behavior as follows. At small times, the fluids show a more rapid increase in, h + and C1+ as the number of modes is increased. This results in a progressive decrease in the magnitude of the stabilizing disturbance vorticity near the interface in the more elastic layer 1 for shortwave disturbances. This results in increased growth rates of the shortwave instability as the number of modes are increased as seen from Fig. 21. The behavior of the multimode fluid at intermediate time scales is quite complex as seen from Fig. 22 in that, the magnitudes of the dynamic viscosities and normal stresses increase monotonically but depend significantly on the time scale. As a result, the boundaries of the neutral stability contour and the magnitudes of growth/decay rates at intermediate waves do not shift monotonically as the number of modes are increased. The differences in the growth/decay rates are still however the largest at intermediate time scales since the differences in the dynamic behavior are most pronounced at intermediate time scales Fig. 23. Note that the differences in the transient viscoelastic behavior are sufficient to cause differences in the linear stability predictions. The differences in the transient viscoelastic behavior can arise either due to the constitutive nonlinearities leading to overshoots in the transient viscoelastic properties at high shear rates or due to changes in transient behavior due to the presence of multiple modes of relaxation or a combination of both.

Fig. 20. Neutral stability contours for the 1-mode (solid line), 2-mode (dashed line), and 4-mode (dotted line) UCM fluids.

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

241

Fig. 21. Growth/decay rates of the: (a) longwave; (b) intermediate wave; and (c) shortwave instability as a function of depth ratio for the 1-mode (solid lines), 2-mode (dashed lines), and 4-mode (dotted lines) systems of UCM fluids.

3.4. Energy analysis The mechanism of interfacial instability has been clearly ascertained for the purely viscous and purely elastic instabilities by Ganpule and Khomami [18] by examining the disturbance energy equation. This equation evaluates the balance of mechanical energy in the flow system to determine the stabilizing and destabilizing effects of the coupling of velocities and stresses from the base flow and perturbation flow. The method of energy analysis used by Ganpule and Khomami [18] and in this study is similar to that developed by Joo and Shaqfeh [25–27] for the analysis of viscoelastic Taylor-Couette and Taylor-Dean flows. For the purpose of studying the effect of combined viscosity and elasticity stratification and that of the transient viscoelastic properties of the fluids on the mechanism of interfacial instabilities, we have used the Giesekus and WM constitutive equations. The disturbance energy equation is obtained by multiplying the linearized perturbation momentum equation by the perturbation velocity and integrating it over a unit volume in space (i.e. V( : −eByB 1, 0 BxB 2p/a) and one periodic cycle in time (0BtB2p/)acr)). The temporal integration is necessary

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

242

Fig. 22. Transient viscoelastic properties of the one mode (solid lines), two mode (dashed lines), and four mode (dotted lines) UCM fluids during the startup of steady shear flow at g; = 1; normalized by the respective steady state values, (a) viscosity, (b) first normal stress coefficient.

because the instability is an overstable one (i.e. cr "0). The disturbance energy equation has the following form:



− Re



dKE dVD + frey = fpre − fvis − fel +fjump − +fp +fNLG, dt dt

(20)

where fp = fpv1 + fpv2 + fps + fpv4. (21) All of the terms in Eqs. (20) and (21) with the exception of fNLG have been described by Ganpule and Khomami [18]. The term fNLG represents the contribution of the quadratic stress tensor term in the Giesekus model to the mechanical energy of the system.

Fig. 23. Disturbance growth/decay rates as a function of dimensionless wavenumber for the 1-mode (solid lines), 2-mode (dashed lines), and 4-mode (dotted lines) systems of UCM fluids, (a) e =3.0, (b) e = 0.3.

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

243

Fig. 24. Terms of the disturbance-energy equation as a function of depth ratio for the viscoelastic longwave (a =0.01) instability of the (a) GG2 system, (c) WW2 systems. The interfacial and bulk contributions to fpv1 (b) GG2 system, (d) WW2 system.

&

fNLG = 9 · (t · t) · 6% dV( .

(22)

In Eq. (20), dKE/dt signifies the rate of change of kinetic energy with time, frey corresponds to the Reynolds stresses, fpre is the energy associated with perturbation pressure, and fvis is the viscous dissipation energy term. dVD/dt is the rate of change of energy associated with total viscosity dissipation (i.e. viscosity dissipation associated with solvent and polymer). In Eq. (21), fpv1, fpv2, fpv3 and fps are the energies associated with the coupling between base flow and perturbation flow velocities and stresses and their gradients as defined in [18]. In the presence of inertia, dKE/dt can be used as the term which indicates the stability or instability of the flow. If the flow is unstable, perturbation velocities grow with time. Therefore the perturbation kinetic energy = 6% · 6%dV( , which is positive definite, grows with time. Thus, when dKE/dt is positive, the flow is unstable, and when it is negative the flow is stable. As the flow goes from being stable to unstable (or vice versa) one can identify the mechanism of instability by identifying the term in the energy equation that tracks with the dKE/dt which goes from being negative to positive (or vice versa).

244

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

For creeping flow, the arguments that apply to dKE/dt can be applied to dVD/dt, the rate of change of total (i.e. solvent plus polymer) viscosity dissipation. When dVD/dt is positive, the flow is unstable, and when it is negative the flow is stable. As the flow goes from being stable to unstable (or vice versa) one can identify the mechanism of instability by identifying the term in the energy equation that tracks with the dVD/dt which goes from being negative to positive (or vice versa). As in [18], one would expect that the terms associated with the jump in physical properties across the interface would drive the instability. Therefore, the terms appearing in Eqs. (20) and (21) have been split so that the energy associated with the bulk and the interface could be separately studied. Clearly, the mechanism of viscoelastic instability for all three systems of Giesekus fluids is expected to be identical and the results of the energy analysis are expected to be qualitatively similar for all three GG systems. Similarly, one would expect that a consistent reproducible mechanism of instability could be identified from all three WW systems. Fig. 24 shows the results of the energy analysis for the GG2 and WW2 system at longwaves. From Fig. 24(a), fjump can be identified as the mechanism of instability for the GG2 system

Fig. 25. Terms of the disturbance-energy equation as a function of depth ratio for the viscoelastic shortwave (a = 80) instability of the (a) GG2 system, (c) WW2 systems. The interfacial and bulk contributions to fpv1 (b) GG2 system, (d) WW2 system.

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

245

Fig. 26. Terms of the disturbance-energy equation as a function of depth ratio for the viscoelastic longwave (a =0.01) instability of the (a) GG1 system, (b) WW1 systems, (c) GG3 system, (d) WW3 system.

since it tracks the rate of change of total viscous dissipation, dVD/dt. Fig. 24(b) shows the interfacial and bulk contributions to fpv1 for this system. The interfacial contribution to fpv1 tracks dVD/dt and is therefore also the mechanism of instability. From Fig. 24(c, d), fjump or the interfacial contribution to fpv1 can not be identified as the mechanism of instability for the WW2 system. Moreover, none of the terms track the rate of change of total viscous dissipation. Fig. 25 shows the results of the energy analysis for the GG2 and WW2 system at shortwaves. From Fig. 25(a), fpv2, fNL, fps, fvis can all be identified as the mechanism of instability for the GG2 system. Fig. 25(b) shows the interfacial and bulk contributions to fpv1 for this system neither of which track dVD/dt From Fig. 25(c, d), fps and the interfacial contribution to fpv1 can be identified as the mechanism of instability for the WW2 system. Fig. 26 shows the results of the longwave energy analysis for the GG1, GG3, WW1, and WW3 systems. It can be clearly seen from Fig. 26(a, c) that the mechanism identified above (namely fjump) for the longwave instability of the superposed Giesekus fluids can no longer be reproduced from the energy analysis of the GG1 and GG3 systems. As seen from Fig. 26(b, d), none of the terms of the disturbance energy equation track dVD/dt for the WW1 system, whereas for the WW3 system, fjump does track dVD/dt.

246

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

Fig. 27 shows the results of the shortwave energy analysis for the GG1, GG3, WW1, and WW3 systems. From Fig. 27(a), fpv2 and fps can both be identified as the mechanism of instability for the GG1 system, whereas for the GG3 system (Fig. 27(c)), fpv2, fNL, fps, and fvis can all be identified as the mechanism of shortwave instability. As seen from Fig. 27(b, d), none of the terms of the disturbance energy equation track dVD/dt for the WW1 system, where as for the WW3 system, fps does track dVD/dt. Thus, the energy analysis indicates that the three GG systems do not exhibit qualitatively similar disturbance energy behavior in that the mechanism identified from one GG system does not hold for another. Moreover, in some cases a single term of the disturbance energy equation can not be identified as the mechanism of instability. The same is true for the three WW systems. This seemingly surprising failure of the disturbance-energy analysis to identify the mechanism of instability is due to the fact that in the presence of viscosity and elasticity stratification, there is more than one mechanism of instability at work. In particular, Huang and Khomami [28] have shown by asymptotic analysis that at longwaves when the viscosity and elasticity ratios are of 0(1), the viscoelastic interactions as opposed to the purely viscous or purely elastic forces are the main driving force of the interfacial instabilities in pressure driven channel flows. They have also

Fig. 27. Terms of the disturbance-energy equation as a function of depth ratio for the viscoelastic shortwave (a = 80) instability of the (a) GG1 system, (b) WW1 systems, (c) GG3 system, (d) WW3 system.

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

247

shown that the driving force of the shortwave instability is always due to (i.e. at all viscosity and elasticity ratios) the viscoelastic interaction terms [28]. In case of the viscoelastic instabilities studied here, the structure of the eigenfunctions depends on both mechanisms of the instability (i.e. purely viscous and purely elastic). Since the disturbance energy equation is constructed utilizing the perturbation velocity and stress eigenfunctions which themselves manifest all of the mechanisms at work, hence any of the terms alone signifying a simple interaction between base flow and perturbation flow velocities and stresses is not necessarily responsible for the rate of change of total viscous dissipation. Therefore, this type of analysis can be effectively used to track the mechanism of instability when only one mechanism (i.e. for the case of purely viscous or purely elastic instability) is involved as shown by Joo and Shaqfeh [25–27] and Ganpule and Khomami [18]. However, when more than one mechanism is at work, it is not surprising that the energy analysis can not distinguish between the different mechanisms present.

4. Conclusions The effect of the differences in transient viscoelastic properties which can arise due to the differences in various constitutive equations or from the presence of multiple modes of relaxation on the linear stability of the superposed pressure driven channel flows has been studied. To study the effect of overshoots in transient viscoelastic properties on interfacial stability, three systems of two Giesekus fluids and three systems of two White-Metzner fluids with steady viscometric properties identical to those of the Giesekus fluids were used. It is shown that the same mechanism of instability is at work for both models, however, the boundaries of the neutral stability contour and the magnitudes of the growth/decay rates are shifted due to the differences in the transient rheological response exhibited by the two models. In particular, the Giesekus model exhibits overshoots in viscosity and first normal stress coefficient during the start up of steady shear whereas the WM model does not. The effect of this difference on the stability of the interface depends on the wavenumber of the disturbance. This is due to the fact that the viscoelastic behavior at long times (i.e. the steady state behavior) governs the longwave stability, whereas the behavior at short times governs the shortwave stability of the interface. Since the steady viscometric properties of the Giesekus and WM systems are identical, their longwave interfacial stability predictions are identical. However, at intermediate and shortwaves, the two constitutive models exhibit pronounced differences. This is due to the intermediate and short time behavior of the two models being significantly different. At shortwaves, the overshoots in the start up viscosities of the Giesekus fluid contribute to a decrease of the stabilizing disturbance vorticity in the more elastic/viscous layer and a more pronounced decrease in the more significant destabilizing disturbance vorticity in the less elastic/viscous fluid leading to a reduction in the magnitudes of shortwave growth rates. The stability of the interface to intermediate wave disturbances depends on the viscoelastic behavior of the fluids at intermediate times. Since the differences in the transient behaviors of the Giesekus and WM fluids are most pronounced at intermediate times, the differences in their linear stability predictions are most pronounced at intermediate waves. It is also shown that inertia has a pronounced effect on interfacial instability at longwaves. At shortwaves, viscosity and elasticity dominated inertia, hence inertial effects are not significant.

248

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

To study the effect of the presence of multiple modes of relaxation on the stability of the interface, a linear stability analysis was performed using various multimode UCM fluids. It is shown that the presence of multiple modes of relaxation has a significant effect on the stability of the interface. This is also due to the differences in the transient viscoelastic properties of the fluids having different number of modes of relaxation. Again, the longwave behavior of all the systems considered was identical since the steady viscometric behavior of all the systems is identical. The shortwave, and more significantly the intermediate wave stability behavior shows pronounced dependence on the presence of multiple modes of relaxation. Thus, the differences in transient viscoelastic behavior—which can arise either due to the overshoots in the transient viscoelastic properties or due to the presence of multiple modes of relaxation—are sufficient to produce differences in the linear stability predictions. Furthermore, the mechanism of viscoelastic interfacial instabilities is studied by performing a disturbance energy analysis. The results indicate that the mechanism of viscoelastic interfacial instabilities can be described in terms of interaction of mechanisms of purely viscous and purely elastic instabilities. However, since more than one mechanism for the instability is at work, the disturbance energy analysis can not clearly distinguish between them due to the eigenfunctions used in the energy analysis containing information regarding both viscous and elastic effects. Hence, the mechanism of the instability must be determined by a careful examination of disturbance eigenfunctions.

Acknowledgements This work was supported in part by the National Science Foundation under Grant No. CTS-9304035 and CTS-9840937.

References [1] R. Yuriko, Weekly nonlinear behavior of periodic disturbances in two-layer Pouette-Poiseuille flow, Phys. Fluids A 10 (1989) 1666. [2] A.P. Hooper, R. Grimshaw, Nonlinear instability at the interface between two viscous fluids, Phys. Fluids 28 (1) (1985) 37. [3] C.S. Yih, Instability due to viscosity stratification, J. Fluid Mech. 27 (1967) 337. [4] N.D. Waters, The stability of two stratified power law liquids in Couette flow, J. Non-Newtonian Fluid Mech. 12 (1983) 85. [5] B. Khomami, Interfacial stability and deformation of two stratified power law fluids in plane Poiseuille flow, Part I. stability analysis, J. Non-Newtonian Fluid Mech. 37 (1990) 19. [6] C.S. Li, Stability of two superposed elasticoviscous liquids in plane Couette flow, Phys. Fluids 12 (1969) 531. [7] N. D. Waters, A.M. Keeley, The stability of two stratified non-Newtonian liquids Couette flow, J. Non-Newtonian Fluid Mech. 24 (1987) 161. [8] Y. Renardy, Stability of the interface in two-layer Couette flow of upper convected Maxwell liquids, J. Non-Newtonian Fluid Mech. 28 (1988) 99. [9] K.P. Chen, Elastic instabilities of the interface in Couette flow of viscoelastic liquids, J. Non-Newtonian Fluid Mech. 40 (1991) 261.

H.K. Ganpule, B. Khomami / J. Non-Newtonian Fluid Mech. 80 (1999) 217–249

249

[10] Y.Y. Su, B. Khomami, Purely elastic interfacial instabilities in superposed flow of polymeric fluids, Rheol. Acta 31 (1992) 413. [11] S.G. Yiantsios, B.G. Higgins, Numerical solution of eigenvalue problems using the compound matrix method, Phys. Fluids 31 (1988) 3225. [12] Y.Y. Su, B. Khomami, Numerical method solution of eigenvalue problems using spectral techniques, J. Comput. Phys. 100 (2) (1992) 297. [13] Y.Y. Su, B. Khomami, Interfacial stability of multilayer viscoelastic fluids in slit and converging channel die geometries, J. Rheol. 36 (1992) 2. [14] Y.Y. Su, B. Khomami, Stability of multilayer power law and second order fluids in plane Poiseuille flow, Chem. Eng. Comm. 109 (1991) 209. [15] G.M. Wilson, B. Khomami, An experimental investigation of interfacial instabilities in the multilayer flow of viscoelastic fluids; part I. incompatible polymer systems, J. Non-Newtonian Fluid Mech. 41 (1992) 355. [16] G.M. Wilson, B. Khomami, An experimental investigation of interfacial instabilities in the multilayer flow of viscoelastic fluids; part II. Part elastic and nonlinear effects in incompatible polymer systems, J. Rheol. 37 (1993) 315. [17] G.M. Wilson, B. Khomami, An experimental investigation of interfacial instabilities in the multilayer flow of viscoelastic fluids; part III. compatible polymer systems, J. Rheol. 37 (1993) 341. [18] H.K. Ganpule, B. Khomami, An investigation of interfacial instabilities in the superposed channel flow of viscoelastic fluids, Submitted to J. Non-Newtonian Fluid Mech. [19] B. Khomami, M.M. Ranjbaran, Experimental studies of interfacial instabilities in multilayer pressure-driven flow of polymeric melts, Rheol. Acta 36 (1997) 1–22. [20] R.B. Bird, C.F. Curtiss, R.C. Armstrong, O. Hassager, Dynamics of Polymeric Liquids, vol. 2, Wiley, New York, 1987. [21] H.K. Ganpule, B. Khomami, A 3-D Stability Analysis of Multilayer Flow of Viscoelastic Fluids: Effect of Second Normal Stresses, Presented at the 68th Annual Meeting of the Society of Rheology, Galveston, Texas, February, 1997. [22] R.B. Bird, C.F. Curtiss, R.C. Armstrong, O. Hassager, Dynamics of Polymeric Liquids, vol. 1, Wiley, New York, 1987. [23] E.J. Hinch, A note on the mechanism of the instability at the interface between two shearing fluids, J. Fluid Mech. 144 (1984) 463. [24] E.J. Hinch, O.J. Harris, J.M. Rallison, The instability mechanism for two elastic liquids being co-extruded, J. Non-Newtonian Fluid Mech. 43 (1992) 311. [25] Y.L. Joo, E.S.G. Shaqfeh, Viscoelastic Poiseuille flow through a curved channel: a new elastic instability, Phys. Fluids 7 (3) (1991) 1691. [26] Y.L. Joo, E.S.G. Shaqfeh, A purely elastic instability in Dean and Taylor-Dean flow, Phys. Fluids 3 (4) (1992) 524. [27] Y.L. Joo, E.S.G. Shaqfeh, The effects of inertia on the viscoelastic Dean and Taylor-Couette flow instabilities with application to coating flows, Phys. Fluids 11 (4) (1992) 2415. [28] C.T. Huang, B. Khomami, Role of Fluid Elasticity on Stability of Multilayer Coating Flows, Presented at the 68th Annual Meeting of the Society of Rheology, Galveston, Texas, February, 1997.

.