Stability of visoelastic flow between eccentric rotating cylinders

~nmlor Nm-Nemoahm Fiaid Meelmuics ELSEVIER

J. Non-Newtonian Fluid Mech., 63 (1996)97-120

Stability of visoelastic flow between eccentric rotating cylinders Amit Chawda, Marios Avgousti* Department of Chemistry and Chemical Engineering,,, Stevens Institute o! Technology, Hoboken, NJ 07030, U~S'A Received 15 July 1995: in revised form 20 December 1995

Abstract

The viscoelastic flow of an Upper Convected Maxwell fluid, confined between two infinitely long eccentric rotating cylinders, is investigated. The two-dimensional steady-state flow as well as the stability of the flow against true three-dimensional (azimuthally and axially periodic and radially non-periodic) disturbances is analyzed numerically using pseudospectral methods. In this numerical scheme, variables are expressed as Fourier series in the periodic direction and as Chebyshev polynomials in the radial direction. The linear stability analysis shows that the critical Reynolds number, corresponding to the onset of flow instability, increases with eccentricity for a Newtonian flow. The critical wavenumber in the axial direction is found to remain nearly constant in the eccentricity range between 0 and 0.5. Addition of small flow elasticity in high Re flows is found to destabilize the system, causing the critical wavenumbers to slightly increase. For a purely elastic flow, it is found that the critical Deborah number decreases with an increase in eccentricity and the critical wavenumber also decreases for the parameters examined in our study.

Keywords: Eccentric cylinder; Linear stability analysis; Pseudospectral methods; Viscoelastic flow

I. Introduction

The flow between eccentric rotating cylinders continues to attract considerable interest because of its immense theoretical and technical importance. The practical significance of this flow arises from the fact that it occurs in journal bearings which are widely used for hydrodynamic lubrication. While several analyses have been performed for Newtonian flows, the investigation of the viscoelastic flow, other than two-dimensional steady-state solutions, is still almost non-existent. In this study, we have undertaken a systematic analysis of the stability of viscoelastic flow between eccentric rotating cylinders.

* Corresponding author. Present address: Union Carbide Corporation, Polymers R&D, Weston Canal Center, P.O. Box 450, Somerset, NJ 08875, USA. 0377-0257/96/$15.00 © 1996 - Elsevier Science B.V. All rights reserved PII S0377-0257(96)01425-5

98

A. Chawda, M. Avgousti / J. Non-Newtonian Fluid Mech. 63 (1996) 97 120

The problem of two-dimensional flow of a Newtonian fluid in a journal bearing was first addressed by Reynolds [1]. He investigated the problem from the standpoint of lubrication theory, making the assumption that the radii of the cylinders are nearly equal (limit of small gap assumption). Later on, large numbers of other researchers also studied this problem. Complex variable techniques were used [2] to solve the biharmonic equation corresponding to the exact geometry, satisfied by the streamfunction in the Stokes flow. Attempts were also made to analyze flows involving inertial effects. It appears that the analysis of Ref. [3], which involved perturbation of the non-inertial (Stokes) solution, was the first to estimate the effects of inertia. It was later on pointed out [4], with the use of first-order initial perturbation theory, that the inertial correction of Ref. [3] is incorrect. Kulinski and Ostrach [5] performed a small perturbation analysis with the eccentricity ratio as a parameter to consider inertial effects. By neglecting curvature effects and using a perturbation series in the clearance ratio DiPrima and Stuart [6] obtained the linearized inertial correction to the Stokes approximation. Most of these earlier publications treated the problem by placing some restrictions on the geometry. Without making any assumptions on the geometry, San Andres and Szeri [7] analyzed the problem numerically using Galerkin's procedure with B-spline test functions. An important result of this study was that the points of separation and reattachment are Reynolds number dependent and also that the separation point moves in the direction of rotation upon increasing the Reynolds number. Another attempt at treating this problem numerically was made by Roberts et al. [8], who used a three-dimensional spectral algorithm for the solution of Stokes flow between eccentrically rotating cylinders. They calculated the load and couple on the inner cylinder and found their results to be in good agreement with those available from lubrication theory. All the works mentioned so far were limited to the study of Newtonian flow between eccentric rotating cylinders. Perhaps, the most complete analysis to date of two-dimensional steady-state purely elastic flow is by Beris et al. [9]. A spectral/finite element calculation was used to simulate an Upper Convected Maxwell (UCM) modeled flow between eccentric rotating cylinders. This was the first work to present numerically convergent results for large elasticity (Deborah) values in geometries with arbitrary eccentricities, thus overcoming the so-called high Weissenberg flow problem. This success was primarily attributed to the use of high order spectral approximations. Compared to the two-dimensional flow, the study of the dynamics of viscoelastic flows is much more difficult. Though a lot of effort has been devoted to the study of flow instabilities in the Taylor-Couette flow, analysis of flow instabilities corresponding to the case of eccentric rotating cylinders has been very limited. To the author's knowledge, the first theoretical treatment of the stability of flow between eccentric cylinders was by DiPrima [10]. A local stability analysis of Newtonian flows was performed for small gaps and the least stable position was found to be in the maximum gap region. The local analysis also showed that for a considerably wide range of values of eccentricity, the flow between eccentric rotating cylinders is less stable than the corresponding Taylor-Couette flow. These results were shown to be incorrect by Vohr [11], who performed an experimental study. Results showed that the eccentricity has a stabilizing effect over its whole range and the maximum vortex intensity occurs not at the maximum gap region, but at a position approximately 50 ° downstream from it.

A. Chawda, M. Avgousti / J. Non-Newtonian Fluid M~'ch. 63 (1996) 97-- 120

99

A non-local approach, taking the whole flow field into account, was followed by DiPrima and Stuart [12]. The most important result of this global analysis was that the eccentricity was found to have a stabilizing effect, which is in agreement with the experiments of Ref. [1 I]. However, results for large values of eccentricity do not agree with the experimental results of Ref. [i 1], since calculations are based on the assumption that both gap and eccentricity are small. A numerically assisted linear stability analysis involving pseudospectral methods was performed by Oikawa et al. [13]. Results of their analysis showed that the critical Reynolds number increases with eccentricity and the critical axial wavenumber is nearly constant for small eccentricities, while for larger eccentricities ( > 0.7) it increases. In addition, secondary flows for slightly supercritical Reynolds numbers were investigated and torque values tbr these states were reported. Recently, Dai et al. [14] also analyzed the problem of stability and bifurcation numerically. They achieved this by projecting Fourier components of the governing system of equations onto a polynomial sub-space with a B-spline basis, and locating the critical points by parametric continuation. With the help of this method, they were able to extend the results of Refs. [6] and [12] to higher eccentricity values and also found agreement with experimental results [11]. All the above mentioned works have been limited to the study of stability of only Newtonian flows. An investigation of stability of elastic flows in the eccentric rotating cylinder geometry is still lacking. Hence, in this work, we have attempted to extend the previous stability analysis of Newtonian flows to viscoelastic flows. The overall objective of the present study is to numerically simulate the two-dimensional steady-state viscoelastic flow between eccentric rotating cylinders and to examine the stability, of this flow against true three-dimensional (azimuthally and axially periodic and radially non-periodic) disturbances. The viscoelastic flow in our study is modeled by an Upper Convected Maxwell fluid. This is the simplest differential viscoelastic constitutive relation which is capable of predicting a fading memory effect and has been used in the past for studying the stability of the viscoelastic Taylor-Couette flow [15,16]. The numerically assisted stability analysis is performed by following the classical linear stability approach, which involves infinitesimal perturbation of the base flow equations and evaluation of the eignevalues associated with the linearized equations of motion and stress. Stability of the flow is determined by the sign of the largest real part of the eigenvalue. In the present work, we have employed a pseudospectral collocation method to perform both the stability analysis as well as to compute the two-dimensional steady-state flow. We organise the paper in the following fashion. In Section 2 we describe the two-dimensional steady-state problem corresponding to a UCM fluid confined between infinitely long eccentric rotating cylinders. We introduce the governing system of equations and explain the implementation of our numerical scheme. Finally, we validate our method and present some results of the steady-state analysis. In Section 3, a numerically assisted linear stability analysis (LSA) of the base flow is performed. We describe the implementation of the linear stability analysis toward three-dimensional disturbances and present its results for a selected range of geometric and flow parameters. Finally, in Section 4 we state our conclusions of both the steady-state analysis as well as the stability analysis and also present thoughts about future work.

100

A. Chawda, M. Avgousti / J. Non-Newtonian Fluid Mech. 63 (1996) 97-120

2. Two-dimensional steady-state viscoelastic flow analysis We consider two infinitely long eccentric cylinders with axes separated by a distance e, as shown in Fig. 1. The average gap between the cylinders is given by cS= R - Ro, where R and R0 are the radii of the outer and inner cylinders, respectively. A UCM modeled fluid is confined in the annulus between the two cylinders, the inner one of which is rotated at an angular velocity while the outer one is kept stationary. The two-dimensional steady-state flow problem is solved in terms of bipolar cylindrical coordinates• The azimuthal angle d~ and the radial coordinate ~ are related to the rectangular Cartesian coordinates (x, y) by equations: a sinh ~ x = cosh ~ + cos d~'

a sin qb Y = cosh ~ + cos ~ '

(2.1)

where a is the separation between the pole and the origin and is given by the equation: a

(1 - - E2) 1/2

- .o

6

F

/1 - -

E 2\

L'+, +

71/2 .

In the above equation, e and /L are the dimensionless eccentricity and the dimensionless gap, respectively. Eccentricity and gap are made dimensionless in the following manner: E = elc~,

It = 6 l R o

= (R -

(2•3)

Ro)lRo.

R

>

Fig. 1. Eccentric rotating cylinder geometry.

A. Chawda, M. Avgousti / J. Non-Newtmmm HuM Mech. 63 (1996) 97-120

101

The inner and outer boundaries of the annulus are given by: ~t ~ arcsin h(a)/R(),

,~2 - arcsin h(a)/R).

(2.4)

The continuity equation, the steady-state equation of motion and the constitutive equation for a U C M model are written in dimensionless form as V'v=0,

(2.5)

Re(v "Vv) = (V" T) - Vp,

(2.6)

T + DeTj = 9,

(2.7)

where v is the velocity vector, ~ is the strain rate tensor, T is the elastic part of the extra stress tensor and p is the hydrodynamic pressure. These equations have been made dimensionless by scaling position with R0, time with Ro/V, velocity with f~R0 and stress with r/V/Ro, where r/is the (constant) polymer viscosity, Reynolds and Deborah numbers in the above equations are defined as

Re -

RoVp 11

,

De = 2V/Ro,

(2.8)

where 2 is the relaxation time for the U C M fluid. In the U C M constitutive equation, T(t) denotes the upper convected derivative o f stress which is expressed as D T(,) = D---tT - [(Vv) t. T + T.(Vv)].

(2.9)

The components of the equation of motion are expressed in the bipolar coordinate system as

4

)

Z ~v~ _ + v~vep sin qb + Z_v, + - - sinh Re v~ a _~£ a a -~ a _

z ~T:~ a Og

Re t,¢ Za ~v, ~4 -

a

~(

a

~----~- +

,2

a

~v¢

v ~ sin qb + 2' v, a a 5qb + - -+ a Oqb

a

sinh

4 + 2

roy 4

a

a

sin qb -

(2.10)

- --

a a g'

sinh (

sin qb - 2

J a

sinh ~ - - - a Oqb'

(2.11 )

where Z = cosh ~ + cos qb. In order to eliminate pressure, Eqs. (2.10) and (2.11) are cross -differentiated and subtracted. Thus the following equation is obtained:

[,~2

~2 )(?) + b7 + 1

~2 { . T ~ _ T * * ) = [ T e r m l _ T e r m 2 ] + z aS z

'

(2.12)

where T e r m l and Term2 represent the derivatives of the convective terms with respect to 0 in Eq. (2.10) and with respect to ~ in Eq. (2.11), respectively.

A. Chawda, M, Avgousti / J. Non-Newtonian Fluid Mech. 63 (1996) 97-120

102

Defining a modified, dimensionless stress tensor T* [9] as

and gubstituting it in Eq. (2.11) leads to

(~

2

0~ ~22 + 1)

-~2 ,

T;,)=

Z a [Terml _ Term2]"

(2.14)

Also, the velocities ve and v, are expressed in terms of the stream function 9 as

e9

a9

and the definition of 9 is inserted in Eq. (2.14) and in the UCM constitutive equation. The components of the UCM constitutive equation are

T*e+De* +

5(~

+

sinqb - - ~

~"\ eqb sinh ~ +

+

-X-~

sin qb - 2;( ~ - ~ j

+

sinh~

- 2T** Z - - ~

~qb + 4" cos qb

= 2 e----~ e~ ' (2.16)

T'~4,+De*

Z~

+

sinqb ~

+

-Z~

+

sinh~

0qb

+ T~¢,(~* sinh ~ + ~ * sin qb) - T~,(Z 029.-~ + 4" cos ¢b)

029* 29"

= 9* + ~qb2

T~o,+De*

29"

~2 ,

(2.17)

Z--~ +4*sinqb

+ T;,

sinh ~ +

~

+

-Z

sin 4) + 27. ~ - - ~ ]

+9*sinh~

6T;,

- 2T~, - Z ~ - ~

+ 4" cosh

a2* * -

2 ~{ ~qb'

(2.18)

where De* is the modified Deborah number defined as De* = (Ro/a)De. The boundary conditions include the no-slip, no-penetration conditions on the cylinder walls and are given as: v,[~=¢1 = - 1; v,l¢= ¢2= 0; v~[~=~l = 0; v~[~=¢1 = 0. In terms of the stream function, these conditions are written as 9"(~1, qb) = Q'z, 9"(~2, qb) = 0,

~ ~

=~, = 1 + Q, sinh ~,, = ~2 = 0

where Q* is the dimensionless flow rate per unit length of the cylinders scaled with a V.

(2.19) (2.20)

A. Chawda, M. Avgousti / J. NoniNewtonian Fluid Mech. 63 (1996) 97 120

103

For a complete set of equations, a pressure periodicity condition has to be imposed [9]. The periodicity of p(.& ¢) is obtained by integrating Eq. (2.11) with respect to ¢ at radial position ~j and setting it to zero: 2~-

0

--Re v~ ~

Z~ s i n ¢ + v ¢

8¢

Z sinh~

~-.,=0

(2.21)

The coupled system of Eqs. (2.14) and (2.16)-(2.18) along with the boundary conditions and the periodicity condition imposed on the pressure are solved using pseudospectral collocation methods by expanding the modified stream function q, (c, ¢) and the components of the modified extra stress tensor T*(4, ¢) in terms of the Fourier series along the periodic qb (N+ nodes) direction and as Chebyshev polynomials in the non-periodic _~(N~ nodes) direction. In the implementation of the boundary conditions, we consider only Eq. (2.19). The other boundary condition shown by Eq. (2.20) is automatically satisfied by using the following definition. $* = (~ - ~2)2$°, (2.22) where $o is the new unknown function to be determined by the numerical solution (this was first advocated by W. Heinrichs (Spectral multigrid techniques for the Navier-Stokes equations, Comp. Meth. Appl. Mech. Eng., 106 (1993) 297-314) as indicated to us by one of the referees). The reason for using the above definition is that, in order to achieve numerical stability, an at least two-order higher polynomial expansion should be used for the approximation of the stream function than for the vorticity. The above procedure finally results in a set of 4 * N + 1 discretized algebraic equations, where N = N. × N¢. This set of equations is solved by Newton's method. Once the solution (i.e., the values of 4" and the components of T*) is calculated, the velocities and the pressure are calculated in a post-processing step. The velocities are directly recovered from the modified stream function and its derivatives as

q,* t,~= ~ -

+ --Z s i n e ,

(2.23)

~$* + -4" - sinh ~.

(2.24)

The pressure is recovered by integrating the ~-momentum Equation (2.10) and C-momentum, equation (2.11), starting from a reference position (_~o,¢o) which is arbitrarily set equal to zero. This position is taken to be at the surface of the inner cylinder, (~o -- 4 ~) and for ¢ = ¢o = 0. The integration is performed numerically. Once the pressure is obtained, the x- and y-components of the force per unit length on the inner cylinder made dimensionless with #o V are calculated for the case of Re = 0 as follows [9]: 2re

Z

- T~¢ +

Z

0

de, ~=41

(2.25a)

A. Chawda, M. Avgousti / J. Non-Newtonian FluM Mech. 63 (1996) 97 120

104

400

350 e~

.,=~

300 -

250 0

O 200

150 0

I

I

I

I

1

2

3

4

Deborah Number Fig. 2. M a g n i t u d e o f t h e f o r c e p e r u n i t l e n g t h o n t h e i n n e r c y l i n d e r f o r • = 0.1, p = 0.096 a n d Re = O.

2~

FY= f I l + c ° s q b c ° S h ~ T * c c ~ Z

s i n q b s i n h ¢ ( T ~ c ¢ - P ) ] z ~=~, dqb.

(2.25b)

0

To obtain force values around the outer cylinder, the above integration is performed at = ¢2. The integrals are evaluated numerically in the same manner as is done in the case of Eq. (2.21). The direction of the load is given by the angle 0: 0 = t a n - l(Fy/Fx).

(2.26)

In order to evaluate the accuracy of our computational scheme we performed calculations for a U C M fluid with E = 0.1, ~t = 0.1 and Deborah numbers ranging from 0 to 5 and compared our results with those presented in Ref. [12]. We report the force on the inner cylinder and its direction as given by Eqs. (2.25) and (2.26), respectively. These are shown in Figs. 2 and 3, respectively and they compare very well with the results of Ref. [9]. In addition, we determined the convergence of our calculations by comparing the results of the dimensionless flow rate per unit length, Q, for different mesh sizes at two fixed values of De. Table 1 shows the convergence characteristics of the numerical solution for E = 0.1 and/~ = 0.1. It can be seen that at low Deborah number, De = 1, five significant digits of accuracy can be obtained using just 8 nodal points in the periodic direction and 9 in the Chebyshev direction. However at higher Deborah number, De = 5, a mesh of 8 x 17 nodal points would be required to obtain results accurate to four significant digits. Furthermore, we have extended the results of Ref. [12] by including the inertial terms in the m o m e n t u m equations. Our work is the first to show results for the viscoelastic case with

A. Chawda, M. Avgousti / J. Non-Newtonian H u M Mech. 63 (1996) 97 120

105

-10 -20 -30 ID

-40

o

-50

o

-60

.o

-70

.,..~

-80 -90 -100 0

I

I

I

I

1

2

3

4

5

Deborah Number Fig. 3. Direction o f the force on the inner cylinder for ~ = 0.1,/~ = 0.096 and Re = 0.

non-zero Reynolds numbers. We performed calculations for small flow elasticities, De = 0.1 and R e = 0, 50 and 100. Figs. 4(a)-4(d) show the radial profiles of vr and Too (in cylindrical coordinates, transformed from the bipolar coordinates) at ~b = 3rr/2 corresponding to /~ = 1.0 and E--0.5 at D e = 0.1 and 0.0. It can be seen that the presence of a small amount of elasticity in the flow has a small effect on velocity profiles (maximum values of vr corresponding to D e = 0 and D e = 0 . 1 for R e = 100 only differ by 2%). However, the hoop stress profile changes significantly with D e since the normal force and extensional response of the material are altered.

3. Linear stability analysis of three-dimensional disturbances For the linear stability analysis of arbitrary non-axisymmetric flow disturbances, a truncated series of non-axisymmetric infinitesimal disturbances is superimposed onto the primary flow solution (Vss, Tss, Pss) as shown in Eqs. (3.1.)-(3.3). Table 1 C o m p a r i s o n o f flow rate values, Q, for different meshes at E = 0.1, p = 0.1 and Re = 0 Mesh N , × N~

Flow rate Q at De = 1

Flow rate Q at De = 5

8x 8x 16 x 16 x

0.4853592950225018 0.4853589211146384 0.4853589199867829 0.4853589201324097

0.4897233210194052 0.4895396471438791 0.4895443689076385 0.4895464637738026

9 17 17 33

106

A. Chawda, M . Avgousti / J. Non-Newtonian Fluid Mech. 63 (1996) 9 7 - 1 2 0

0.025

0.025

0.020

¢,q

Re~-O

............

Re=50

.......

Re=100 . . . . .

I

0.015 o O

Re=0

, . ~ . , o .~

"'-,,,

.,2

,--, 0 . 0 2 0 ¢'4 k~

",~

I

0.015-

•~0

0.010

0.010

0.005

0.005

0.000 0.0

I

I

I

I

0.2

0.4

0.6

0.8

Normalizedradialdistance

(a)

0.000 0.0

1.0

I

I

I

I

0.2

0.4

0.6

0.8

1.0

Normalizedradialdistance

(b)

2.5

0.30 Re---0

2 ~.~

............

Re=50 . . . . . . . .

0.20 -

Re=100 . . . . . . ¢~

............

............

Re=50

........

; , /

Re=100 . . . . . .

1.5

0.10

. ' ~ ' ~ ~

I

Re=0

-

/.~;# s S o'

t

1

0.00

7 ~ ~t

-

j s SS sos~ S

0.5-

-0.10

0 (d)

0.2

I

I

I

0.4

0.6

0.8

N o r m a l i z e d radial distance

-0.20 0.0

I

I

I

I

0.2

0.4

0.6

0.8

1.0

Normalizedradialdistance

(c)

Fig. 4. Radial profiles at ~ = - 1r/2 for • = 0.5, # = 1.0 and R e = 0, 50, 100. (a) Radial velocity (absolute value) ur for D e = 0 . 1 ; (b) radial velocity (absolute value) Vr for D e = 0 . 0 ; (c) azimuthal normal stress Too for De = 0.1; (d) azimuthal normal stress Too for De = 0.0.

v = Vss(~, 4 0 + ~

Ifl = (

Y2) -- 1 ~

'~l,(4) e ('~'+~-'+~*~

1

,

(3.1)

1,

(3.2)

L # = - N4,/2

T= Tss(~,~b)+ 9~F# = (,-. Y2) - 1 ~#(#) L #=

-

N,b/2

e (2t+ic~z+i#~b)

A. Chawda, M. Avgousti / J Non-Newtonian Fluid Mech. 63 (1996) 97 120

p = pss(4, 4') + ~

[fl=(Ni2)-L

P

=

-

I

::(~)

e ~;' + ~=

107

+ iflqb)],

1'3.3)

N,,'2

where 9t[.] represents the real part, ~/~(~), /~/~(~),fi/:(~) represent complex infinitesimal functions of the radial coordinate, N, is the number of nodal points in the qb direction, ~ is the wavenumber of the assumed periodicity in the axial direction z, and fl is the azimuthal wavenumber which can only assume integral values. The eigenfunctions ~/:(g), i'll(d), 15/~(~) are subsequently approximated using truncated series of Chebyshev polynomials. This means that the perturbed variables have the same double Fourier-Chebyshev expansion as the steady-state variables. For fl = n c N,, the perturbation involves a sinusoidal variation of the variables in the 4' direction of frequency n (nth harmonic). The dimensionless wavenumber ~ is defined as c~= 2rcRo/I~L, where L represents the wavenumber of the assumed periodicity of the perturbation in the axial direction a n d / ; is the dimensionless gap given by Eq. (2.3). This is a classical non-local normal mode form perturbation analysis in which the various directions are decoupled from each other (standard separation of variables). The summation is required to account for the inherent asymmetry of the solution in the case of eccentric cylinder geometry, thus taking into consideration the contribution of all the harmonics which can possibly be excited at a particular radial distance (5). Substituting the expressions given by Eqs. (3.1), (3.2) and (3.3) into the original equations (2.5), (2.6) and (2.7), and neglecting all terms of order higher than linear with respect to the infinitesimal functions, leads to a system of ten first-order homogenous ordinary differential equations: fl-(No/2)--

E

/] /1 -

1

a

-

H

(3.4)

1

[N4/2) - 1

E /1 ~

v¢::+ (zifl + sin 4') t,,, + i:w:/, e{i/:'b)= O,

_

N,5/2

Re(2ve/:) e(ip+)

N m/2

©v}

.

i s)v*/: a

Z--~- + v S s i n 4 , + Z fly, /; =

+

N+,2

X ~P/: + Z ~T~/, a ~5. a ~

v} sin 4' + Z - - ~ + 2v~, sinh

+ (zifl + 2 sin 4') ~ a

+ sinh ~ T,,:a

+

T ~ cs/:

sinh ~ ~ - -

a

(3.5)

i~r:_.} ed/;+)

INdo,,2)- 1

Re(;tv+/~) e(i/~*) fl

=

- (N,b/2)

Z [] =

-Re

v)--~- + Z - - ~ - 2 v )

sin 4, - v~ sinh ~) Vena

N 4~ '2

+ zv;ifl+Z-~

--v~sinh~

- 2 s i n h : T~,~a + ( z i f l + sin 4 , ) - ~

-- a

+i:~T,:/:}@ i/'*),

a

Z ~T4*/: a ~5 (3.6)

A. Chawda, M. Avgousti / J. Non-Newtonian Fluid Mech. 63 (1996)97-120

108 fl=(N¢/2)--

1

Re(Av~/j) e ~i/~ fl = - (N,/2)

E

-- R e

a= - N,/2

v ¢ -~

+ - iflv~

a

- kzp/~ +

a 5~

sinh 4 - -

a

+(xifl + sin ~b) T¢:/,a + i~T__:/,] e ~i/~¢~, fl = (N¢,/2)

[1

=

(3.7)

-- 1

- De(2T¢¢p) e 0fl¢~

2

_ (N~b/2)

fl = (N,/2)

-- 1 ~

r'~ ,'r,s

~] (-2Z-(l+ve.¢¢)-~ = _ N,/2 a

+De

[z ST~¢

,, 81)~[::

Z

- +De-v~ a

8T~[

:l

5--T

TSc, ]v~/~

8----~--+ 2(sin qb - zifl) T

2 sina qb )v¢/J + ( De ~8T~¢ 8---dp--2De sin qbT~___A¢ a

[

+ l+O(ifl-Zsinfl)"¢De-ZDeZ--~JT¢¢,+ a a

2v) sind~-z78vsc)T "(eOfl¢) 8dp] ¢¢/iJ , (3.8)

[]

= (N~b/2)

-- 1

-- De(2T¢¢a) e (ifl¢) fl = --Nob~2

fl=(Neo/2)--12

--

De Z. TS¢, 8v¢~8~

Xa(1 + DeT~¢) - 8v¢~ 8~ + DeZ-vSea --8~ 5T¢¢13

fl = _ N ~b/2

5T~¢ + Dez ~ + DeT~¢¢sinh ~ + (DeT~¢ + 1)(sin qb 8TS¢¢ + Dez 8dp -De

v~

gift)] v¢___~a

DeT~¢(sin 0 + zifl) - (DeT~¢¢+ 1)sinh ~ 1 v~_

sinh ~ + Z --@--; a

+ [De(zi fl - sin ~) vs* - De 7. -5vs~ k a a 84

+De(v) sin qb-Z 8qbJ

D e Z_ a -8v~ ~

s sinh 4 + 1] T,,/~ + De V~a

(3.9)

A. Chawda, M. Avgousti /J. Non-Newtonian

Fluid Mcch. 63 (1996) 97-120

109

$aT+b,i CI at V<,i 1 1

y +DeLv,+ 2(DeT&

+ 1) sinh 5

-

a

aT”,, -De-2Desinh~T&,-2X$‘(DeT&,+l) X a+

+

:

-2De x

De”igv$-2De-a /i = (N+;2)

-

-

c

/i=(N+:2)--

-

De(3,T,,,,) eci”@ I

(

y

I

c

- De(;l TbzB)eciB@)

/I= --(N+,;Z)

,

-De%T&,z

p = (N@,2) -

(3.10)

a w

I

/j = ~~(N+!2)

p = (N,+,‘Z)

av;

au_

X

+DeavF

aT+=, at

- iccv+,j- i/l $DeT&, a

+ 1)2;,,,

I

(3.13) where for simplicity, the hats on p, u and T have been omitted and the steady-state values of the variables are represented by a superscript s. The mathematical problem of stability is formally completed (boundary value problem) only when the above described system of equations is accompanied by the boundary conditions. These are six no-slip and no-penetration conditions on the cylinder walls and these are given as (3.14) 5 = t1 and 5 = t2 6.. = 6, = I?,= 0, at both

A. Chawda, M. Avgousti / J. Non-Newtonian Fluid Mech. 63 (1996) 97-120

110

The above system of linearized equations constitutes a homogenous boundary value problem. These equations are discretized using a mesh identical to that used for the calculation of the base flow. This means that we use Chebyshev spectral collocation to discretize the variables in the radial ~ direction and Fourier series expansion in the periodic qb direction. Once all the equations are discretized, the corresponding linearized boundary value problem results in a generalized complex eigenvalue problem which can be represented as A*U=2B*U,

(3.15)

where the vector U is generated by the values of the unknown variables at the grid points. Both A and B are complex-valued square matrices of dimension 10(N) x 10(N) where N is the total number of grid points, as mentioned before. This system is solved using standard generalized eigenvalue routines. The main problem associated with computing the full eigenvalue spectrum is that it is computationally very expensive, requiring a lot of CPU time and large memory allocation. Furthermore, for the linear stability analysis, we are interested only in eigenvalues with the largest real part and so there is no need to compute all the eigenvalues. One method that exploits this particular characteristic to perform eigenvalue computations efficiently is the Arnoldi method. In this method, the generalized problem is transformed to a system: (A -

crB) - 1BZ = ~Z,

(3.16)

whose eigenvalues ~ are related to those of the original system Eq. (3.15) by £ = 1/(2 - o-),

(3.17)

where o-, a complex number, is called the shift such that in this scheme eigenvalues located close to the shift convege the fastest. A restarted Arnoldi algorithm [17-20] is then used to solve Eq. (3.16). The first step involves forming the matrix (A - o-B) and computing the product of matrix B with the vector Z. The left hand side of Eq. (3.16) is then solved by performing the LU factorization of the matrix (A -aB) followed by forward and backward solves. Of course, the value of vector Z is not known a priori and therefore an iteration strategy, known as shift-and-invert process, is adopted. Then, an upper Hessenberg matrix of much lower order than the original matrices is formed and subsequently its eigenvalues which are approximations to the leading eigenvalues of the modified problem given by Eq. (3.16) are obtained. Finally, Eq. (3.17) is used to find the leading eigenvalues coresponding to the original system, Eq. (3.15). The stability of the flow is determined from the eigenvalues of Eq. (3.15). The eigenvalue characterizing the stability of the flow is the one with the largest real part. If this eigenvalue has Re(2) < 0, then the corresponding eigenmodes will decay and the flow is stable. At some critical value of Reynolds or Deborah number, the leading eigenvalue (i.e., the eigenvalue with the maximum real part) crosses the imaginary axis into the right half plane, so that the corresponding eigenmode becomes linearly unstable. If the imaginary part of this leading eigenvalue is zero, then we have a simple bifurcation to a steady secondary solution. On the other hand, if the imaginary part is non-zero, then the phenomenon is referred to as overstability and we have a Hopf bifurcation to a time periodic secondary solution, whose frequency at the onset of instability is determined by the magnitude of the imaginary part of the eigenvalue at this crossing. For this case, the stability characteristic eigenvalues appear in complex conjugate pairs.

A. Chawda, M. Avgousti / J. Non-Newtonian Huid Mech. 63 (1996) 97 120

111

Table 2 Comparison of Re'~ for ~ = 0.001,/t = 0.096 and e = 3.17 Source

Re'~,

Ref. [11], experimental Ref. [12] Ref. [14] This work

140+ 1.5% 139.94 140.0 140.1

3. I. Lmear stability analysis results In this section, we report the linear stability analysis results for the eccentric cylinder problem. Because of the large parameter space corresponding to the viscoelastic eccentric cylinder problem, we chose to confine ourselves to geometries previously used by others in either the same geometry [11,12,14] or the concentric cylinder case [15,16]. The results have been organized into three main subsections: (a) Newtonian flow, (b) Purely elastic flow, and (c) Viscoelastic flow. Because o f the large number of terms involved in the stability equations, it is necessary to establish the validity of the code. Hence, in the following subsections we also describe the validation procedure along with the results. 3. 1.1. Newtonian flow (De = O) As a first step to check the accuracy of our numerical scheme, we performed calculations at the asymptotic limit of low eccentricity, e =0.001, and compared these results with those available in the literature for the concentric cylinder case [11,12,14]. Comparison of the values of critical Reynolds n u m b e r calculated by us with the literature values for the T a y l o r - C o u e t t e case is shown in Table 2. As seen, there is an excellent agreement with the published results. It is important to note that the modified Reynolds number, Re', given here is defined as R e = d V p / q = Re/t~ where d = R - R0 and Re is Reynolds n u m b e r defined by Eq. (2.8). We have used the above definition throughout this section in order to readily make the comparisons with the published results. The superscript prime on Re has been omitted for the remainder of this section. In the next step of the validation procedure, we compared our results with those obtained by Dai et al. [14] for a large range of eccentricities. The results, presented in Fig. 5. show the variation of critical Re with eccentricity for axial wavenumber c~= 3.17 and gap /1 = 0.096 (which corresponds to Ro/R = 0.912 used in Ref. [14]). Our results compare very well with those of Ref. [14] for 0 < • < 0.4. For ~ = 0.5, however, we obtained Recr = 198.5 which lies in between the values obtained by Dai et al. [14] and those determined experimentally by Vohr [11]. This is not surprising since the numerical accuracy of the solutions deteriorates appreciably with increasing eccentricity. The value of e is determined from the minimum of the Reynolds number at the onset of instability plotted as a function of the wavenumber ~. This curve, usually referred to as the neutral stability curve, is shown in Figs. 6(a) and 6(b) for • = 0.2 and E = 0.4, respectively. As seen, the critical value for ,~ is fairly insensitive to eccentricity (~rit = 3.14 for • = 0.2 and 0~crit 3.2 for • = 0.4). This is in agreement with the findings of Ref. [14]. =

112

A. Chawda, M. Avgousti / J. Non-Newtonian Fluid Mech. 63 (1996) 97-120

As an additional step to verify the numerical scheme, we compared our results with those of Oikawa et al. [13]. These results were obtained using a 32 x 9 grid which represents the most dense grid we could compute due to computer memory limitations and are shown in Table 3. They are presented in the form of real and imaginary parts of the three most unstable eigenvalues corresponding to E = 0.7,/~ = 0.1, 7 = 4.126 and Recr= 307.59. It is important to note that in order to make comparisons, we have modified the eigenvalue results of Ref. [13] by dividing these with (/~Re). This is required since our definition of characteristic time (t = Ro/V) is not consistent with the definition (t = d2/v, where v is the kinematic viscosity) used in Ref. [13]. As seen from the table above, our results compare quite well with those in the literature.

3.1.2. Purely elastic flow (Re = 0, De ¢ 0) Due to lack of data available in the literature for the purely elastic viscoelastic flow in eccentric cylinder geometry, we validated our numerical scheme by comparing our results with these corresponding to the asymptotic limit of the concentric cylinder case. To check the accuracy of our results, we used the numerical scheme of Ref. [16] which investigated the stability of T a y l o r - C o u e t t e flow against non-axisymmetric disturbances. To make comparisons, we first varied the azimuthal wavenumber for the T a y l o r - C o u e t t e problem and obtained the corresponding critical Deborah numbers which are listed in Table 4 below. The case fl = 0 corresponds to the axisymmetric case. The implementation of the case is straightforward for the eccentric cylinder since this only requires reducing the number of grid points in the qb direction to one. Comparison of our result with the concentric cylinder solution for the axisymmetric case is shown in Table 5. In the next step of the validation procedure, we compared our results with the stability results for the non-axisymmetric case of the T a y l o r - C o u 200

190

180

170

160

150

140

I

I

I

I

O.l

0.2

0.3

0.4

0.5

Fig. 5. Variation o f critical Reynolds mmlber with eccentricity e for/~ = 0.096, De = 0 and ~t = 3.17.

A. Chawda, M. Avgousti / J. Non-Newtonian Fluid Mech. 63 (1996) 97 120

113

149

148.5

148

147.5

147 2.8

I 2.9

i 3

, 3.1

, 3.2

I 3.3

I 3.4

3.5

12 (a)

73.9

73.8

73.7

73.6

73.5

73.4

73.3

73.2 3

I

I

I

3.1

3.2

3.3

3.4

12

(b) Fig. 6. N e u t r a l s t a b i l i t y c u r v e f o r inertially d o m i n a t e d i n s t a b i l i t i e s f o r it ~ 0.096 a n d D e = 0; (a) • = 0.2: (b) • = 0.4.

A. Chawda, M. Avgousti / J. Non-Newtonian Fluid Mech. 63 (1996) 97-120

114

Table 3 C o m p a r i s o n of real a n d imaginary parts of the eigenvalues a for • = 0.7, ;~ = 0 . 1 , ~ = 4.126 a n d Recr = 307.59 This work

Ref. [13]

6 r

0 -4.480x10 -0.0310

3

O"1

0-~

O"1

±0.3728 0 0.7168

0 -3.566x10 -0.0329

±0.3732 0 0.7194

3

Table 4 Critical D e b o r a h n u m b e r c o r r e s p o n d i n g to different a z i m u t h a l w a v e n u m b e r s fl for Taylor C o u r e t t e flow;/~ = 0.0526, Re = 0, a n d c~= 7.8 A z i m u t h a l w a v e n u m b e r , fl

Critical D e b o r a h , Deer

0 1 2 3

1.563 1.184 1.116 1.163

Table 5 C o m p a r i s o n of critical eigenvalue for ;t = 0.0526, c~= 7.8 a n d Re = 0 A z i m u t h a l w a v e n u m b e r fl

0 2

Critical D e b o r a h n u m b e r Decr

1.563 1.116

O'imag

Taylor Couette

This work

0.70058 0.09538

0.70062 0.09541

ette flow. As seen from Table 4, the most unstable disturbance is the one which corresponds to fl = 2. Implementation of this case for the eccentric cylinder would require changing the number of grid points in the qb direction to N , = 4. Comparison with the T a y l o r - C o u r e t t e result is shown in Table 5. An excellent agreement can be seen between the results of our work and those for the T a y l o r - C o u e t t e flow which validates the correctness o f our numerical scheme. We now present some results for the U C M fluid at other values of eccentricity. All the presented results were performed with a mesh corresponding to N , = 16 and Ne = 17. This was

A. Chawda, M. Avgousti / J. Non-Newtonian Fluid Mech. 63 (1996) 97 120

115

determined to be the smallest mesh able to provide converged results, indicating the non-local nature of the instability. In Figs. 7(a) and 7(b), we show neutral stability curves for eccentricity E = 0.001 and E = 0.1, respectively. It can be seen that the critical wavenumber shows a significant variation with eccentricity. At ~ = 0.001, the critical wavenumber is ~ = 7.8 but for E = 0.1, ~ = decreases to about 5.6. Next, we show the effect of eccentricity on the stability of the flow in the absence of inertial terms. This is illustrated in Fig. 8 which is a plot of critical Deborah number vs. eccentricity. Although the axial wavenumber :~ shows a variation with eccentricity ~ as mentioned before, a constant value of :~ = 7.8 has been used to obtain results shown in Fig. 8. F r o m this figure, it can be concluded that, for the case of inertialess flow, eccentricity has a destabilizing effect. This is in contrast to the case of Newtonian flow where eccentricity stabilizes the flow. This indicates that the elongational character of the fluid plays an important role in the morphogenesis of this particular mode of instability. This is in contrast to the Taylor Couette shear-driven instabilities which are attributed to the normal t\)rce behavior (hoop stress distribution) of the fluid. Regarding the nature of instability, it is found that except for a very low eccentricity of E _< 0.05 (which is similar to the T a y l o r - C o u e t t e flow), for all other values of eccentricity in the range 0.05 < ~ < 0.3, the critical eigenvalue has a zero imaginary part. Hence an exchange of stability to a steady secondary flow structure occurs in this case.

3.1.3. Viscoelastic flow (high Re, low De) The accuracy of the numerical scheme with respect to either inertial or elastic terms (terms multiplied by Re and De) has already been established independently. In this subsection we directly present some new results corresponding to the presence of both inertia and elasticity terms so that both Re and De assume non-zero values. First, the effect of introducing a small a m o u n t of elasticity (small De) in the Newtonian flow is considered. The results of our calculation are presented in Fig. 9 which is a plot of critical Reynolds number vs. eccentricity corresponding to /L = 0.096, ~ = 3.17 and De = 0.05 and 0.1. As expected, an increase of Deborah number results in a decrease of the critical Reynolds number at any given value of eccentricity. It can be seen that the decrease in critical Re with eccentricity corresponding to two different values of De is approximately constant. This means that, under these conditions, elongational effects are not so important and are not the driwag force for this mode of instability. In Fig. 10, we show the neutral stability curve t\w E :--0.1 and De = 0.1. As shown, the critical wavelength value of 3.37 is obtained which is slightly higher than the value of 3.17 obtained for the Newtonian case. This means that elasticity has a destabilizing effect on high Reynolds n u m b e r flows and also decreases the critical wavelength.

4. Conclusions

First, we make a c o m m e n t on the numerical methods used. The use of matrix methods coupled with high order pseudospectral approximations in the investigation of stability problems is very powerful. The ordinary differential equation system resulting from the linear stability of the base flow is transformed into a generalized eigenvalue problem after spectral discretizations.

116

A. Chawda, M. Avgousti / J. Non-Newtonian Fluid Mech. 63 (1996) 97-120 1.16

1.15

1.14

1.13

1.12

1.11 7.4

I

I

I

I

I

I

7.5

7.6

7.7

7.8

7.9

8

8.

(a)

0.85

0.84 -

0.83 0.82

0.81 0.8 0.79 0.78 4.5

I

I

I

I

5

5.5

6

6.5

7

O~ (b) Fig. 7. N e u t r a l stability curve for elastically d o m i n a t e d instabilities for; (a) E = 0.001; (b) e = 0.1. O t h e r p a r a m e t e r s a r e / ~ = 0.0526, c~ = 3.17 a n d Re = O.

A. Chawda, M. Avgousti / J. Non-Newtonian Fluid Mech. 63 (1996) 97-120

117

1.2 ],l

10.90.80.70.60.5 0.4 0

I

I

I

I

I

0.05

0.1

0.15

0.2

0.25

0.3

E Fig. 8. Variation of critical Deborah number with eccentricity e for/z = 0.0526, ~ = 7.8 and Re = 0.

Thus, the solution of large eigenvalue problems is also critical to the success of these methods. We have found that the Arnoldi method calculates the less stable eignevalues efficiently and accurately. Its use has been instrumental in our investigations which involve prohibitively large matrices for the calculation of their entire eigenspectrum. 155

150

t" 145

140 135

~

o,'-

~

De=0.1

oOS

ooo~ooo°°°" 130

"'"

125

I

0

0.1

I

0.2

0.3

Fig. 9. Variation of critical Reynolds number with eccentricity ~ for /~ = 0.096, c~= 3.17 and De = 0.05, 0. [.

118

A. Chawda, M. Avgousti / J. Non-Newtonian Fluid Mech. 63 (1996) 97-120 132

131.8 -

131.6 -

131.4-

131.2 -

131

130.8

I

I

I

I

I

3.1

3.2

3.3

3.4

3.5

3.6

12 Fig. 10. N e u t r a l stability c u r v e f o r U p p e r C o n v e c t e d M a x w e l l fluid f o r e = 0.1, p = 0.096 a n d De = 0.1.

Several important conclusions can also be drawn from the linear stability analysis of the flow against true three-dimensional disturbances. It was found that in the case of Newtonian flows, if the eccentricity is increased, the corresponding critical Reynolds number also increases and therefore the flow is stabilized. This may be attributed to the fact that eccentricity adds a streamwise extension rate in the flow which has a stabilizing effect. At higher eccentricities, non-local streamwise extensional effects become an important factor in stabilization of flows and hence any study based upon "local" linear stability analysis would fail to predict correctly the critical conditions for the onset of instability. In this sense, the linear stability analysis which we followed is more correct than perturbation or local expansion-based methods since it is "global" in nature and does not involve any restriction on the geometry. The accuracy of our analysis is evident from the good agreement seen between our results of Recr and other available experimental as well as theoretical results. In addition, in this work we have successfully extended the linear stability analysis to viscoelastic flows. Hence, this becomes the first work to investigate the stability of elastic flows in the eccentric cylinder geometry. For purely elastic flows, we found that, at least for the parameter space searched here, as the eccentricity is increased the critical Deborah number decreases. We, therefore, conclude that for purely elastic flows eccentricity has a substantial destabilizing effect which is contrary to the effect observed for the case of Newtonian flows. We also analyzed the flows involving both the viscous and elastic terms and found that introducing elasticity in the flow lowers the critical Reynolds number. This is to be expected since elasticity is known to act as a destabilizing agent due to the addition of normal stres behavior.

A. Chawda, M. Avgousti / J. Non-Newtonian Fhdd Mech. 63 (1996) 97 120

119

Another important result of our analysis is that for purely elastic flows, the critical wavenumber is affected significantly by the eccentricity. The critical wavenumber was determined to be :~ = 7.8 for the very low eccentricity value of E = 0.001 and it dropped to ~ = 5.6 for e = 0.1, thus showing that the wavelength associated with the cellular structure increases with increasing eccentricity. On the other hand, for the Newtonian case the critical wavenumber was affected only slightly (small increase) by the eccentricity in the range 0-0.5 for a dimensionless gap of approximately 0.1. Some key observations regarding the nature of instability can also be made. In the case of Newtonian flows, the least stable eigenvalues are real for all values of eccentricities in the range 0 < e _< (t.5 with gap/L = 0.096, thus showing that a simple steady bifurcation occurs at the point o f instability. Only for the case of E = 0.7 and # = 0.1, was the critical eigenvalue found to have a non-zero imaginary part, showing that in this case a H o p f bifurcation to time periodic secondary solution occurs. For the purely elastic flows corresponding to a small gap (/L = 0.05), it was found that, except for a very low eccentricity of e_< 0.05 (which is more like the Taylor- Couette flow), for all other values of eccentricity in the range 0.05 < e < 0.3, the critical eigenvalue has a zero imaginary part. Hence an exchange of stability to a steady secondary flow structure occurs in this case. It is worth mentioning that this is a different type of bifurcating family from the T a y l o r - C o u e t t e one (the least stable real eigenvalue was not lbund in the Taylor Couette eigenspectrum). As with the T a y l o r - C o u e t t e problem, we can distinguish between two different mechanisms, the inertial one responsible for Newtonian instabilities and the elastic one responsible for viscoelastic instabilities. The governing physical phenomenon behind the inertially driven instabilities is the distribution of angular m o m e n t u m . In the presence of eccentricity, this is disturbed in such a way as to delay the onset of instabilities. On the other hand, when the instabilities are driven by the flow elasticity, what primarily matters is the magnitude of the stresses and most importantly for shear-dominated flows, the hoop stresses. Their magnitude considerably increases with the introducion of eccentricity through a combination of extensional and higher shear flow rate characteristics. This may help explain the fact that. in the inertialess viscoelastic case, eccentricity has a destabilizing effect.

References [1] O. Reynolds, On the theory of lubrication and its application to Mr. Beachamp Tower's experiments. Philos. Trans. R. Soc. London, 177 (1886) 157 234. [2] G. Wannier, A contribution to the hydrodynamics of lubrication, Q. Appl. Math., 8 (1950) 1 32. [3] M.M. Kamal, Separation in the flow between eccentric rotating cylinders, ASME J. Basic Eng., 88 (1966) 717. [4] B.Y. Ballal and R.S. Rivlin, Flow of a Newtonain fluid between eccentric rotating cylinders: inertial effects. Arch. Rational Mech. Anal., 62(1976) 237-294. [5] E. Kulinski and S. Ostrach, Journal bearing velocity profiles for small eccentricity and moderate modified Reynolds number, ASME J. Appl. Mech,, 34 (1967) 16-22. [6] R.C. DiPrima and J.T. Stuart, Flow between eccentric rotating cylinders, ASME J. Lubrication Technol., 94 (1974) 266 274. [7] A. San Andres and A.Z. Szeri, Flow between eccentric rotating cylinders, ASME J. Appl. Mech.. 51 (19841 869 878.

120

A. Chawda, M. Avgousti / J. Non-Newtonian Fluid Mech. 63 (1996)97-120

[8] G.W. Roberts, A.R. Davies and T.N. Phillips, Three-dimensional spectral approximations to Stokes flow between eccentrically rotating cylinders, Int. J. Numer. Method Fluids, 13 (1991) 217-233. [9] A.N. Beris, R.C. Armstrong and R.A. Brown, Spectral/finite-element calculations of the flow of a Maxwell fluid between eccentric rotating cylinders, J. Non-Newtonian Fluid Mech., 22 (1987) 129-167. [10] R.C. DiPrima, A note on the stability of the flow in loaded journal bearings, Am. Soc. Lubrication Engineers Trans., 6 (1963) 249-253. [11] J.H. Vohr, An experimental study of Taylor vortices and turbulence in flow between eccentric rotating cylinders, J. Lubrication Technol., Trans. ASME, 90 (1968) 285 296. [12] R.C. DiPrima and J.T. Stuart, Non-local effects in the stability of the flow between eccentric rotating cylinders, J. Fluid Mech., 54 (1972) 393-415. [13] M. Oikawa, T. Karasudani and M. Funakoshi, Stability of flow between eccentric rotating cylinders, J. Phys. Soc. Jpn., 58 (1989) 2355-2364. [14] Rui-Xiu Dai, Q. Dong and A.Z. Szeri, Flow between eccentric rotating cylinders: bifurcation and stability, Int. J. Eng. Sci., 30 (1992) 1323-1340. [15] M. Avgousti and A.N. Beris, Viscoelastic Taylor-Couette flow: Bifurcation analysis in the presence of symmetries, Proc. R. Soc. London, Ser. A, 443 (1993) 17 37. [16] M. Avgousti and A.N. Beris, Non-axisymmetric modes in viscoelastic Taylor-Couette flow, J. Non-Newtonian Fluid Mech., 50 (1993) 225-251. [17] Y. Saad, Numerical solution of large nonsymmetric eigenvalue problems, Comput. Phys. Commun., 53 (1989) 71. [18] K.N. Christodoulou and L.E. Scriven, Finding leading modes of a viscous free surface flow: an asymmetric generalized eigenproblem, J. Sci. Comput., 3 (1988) 355-406. [19] R. Natarajan, An Arnoldi°based iterative scheme for non-symmetric matrix pencils arising in rising element stability problems, J. Comput. Phys., 100 (1992) 128-142. [20] R. Sureshkumar and A.N. Beris, Linear stability analysis of viscoelastic Poiseuille flow using an Arnoldi-based orthogonalization algorithm, J. Non-Newtonian Fluid Mech., 56 (1995) 151-182.