Linear stability analysis of flow of an Oldroyd-B fluid through a linear array of cylinders

J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

Linear stability analysis of flow of an Oldroyd-B fluid through a linear array of cylinders M.D. Smith1 , Y.L. Joo2 , R.C. Armstrong∗ , R.A. Brown Department of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA Received 24 September 2001; received in revised form 28 August 2002

Abstract The linear stability of the flow of an Oldroyd-B fluid through a linear array of cylinders confined in a channel is analyzed by computing both the steady-state flows and the linear stability of these states by finite element analysis. The linear stability of two-dimensional base flows to three-dimensional perturbations is computed both by time integration of the linearized evolution equations for infinitesimal perturbations and by an iterative eigenvalue analysis based on the Arnoldi method. These flows are unstable to three-dimensional perturbations at a critical value of the Weissenberg number that is in very good agreement with the experimental observations by Liu [Viscoelastic Flow of Polymer Solutions around Arrays of Cylinders: Comparison of Experiment and Theory, Ph.D. dissertation, Massachusetts Institute of Technology, Cambridge, MA, 1997] for flow of a polyisobutylene Boger fluid through a linear periodic array of cylinders. The wave number of the disturbance in the neutral or spanwise direction also is in good agreement with experimental measurements. For closely spaced cylinders, the mechanism for the calculated instability is similar to the mechanism proposed by Joo and Shaqfeh [J. Fluid Mech. 262 (1994) 27] for the non-axisymmetric instability observed in viscoelastic Couette flow, where perturbations to the shearing component of the velocity gradient interact with the polymer stresses in the base flow. However, when the cylinder spacing is increased, we discover a new mechanism for instability that involves the coupling between the elongational component of the velocity gradient and the streamwise normal stress in the base state. In addition, the most unstable disturbance in the Couette geometry is over-stable (leading to time-periodic states) whereas the instability calculated for closely spaced cylinders grows monotonically with time. This apparent inconsistency is resolved by the observation that in a complex flow, the evolution of the perturbations to the base flow can be transient in a Lagrangian frame of reference, while steady in an Eulerian frame. © 2002 Elsevier Science B.V. All rights reserved. Keywords: Stability; Cylinder; Three-dimensional perturbation; Protean coordinates

∗ Corresponding author. Tel.: +1-617-253-4581. E-mail address: [email protected] (R.C. Armstrong). 1 Present address: KLA-Tencor, Austin, TX 78759, USA. 2 Present address: School of Chemical and Biomolecular Engineering, Cornell University, NY 14850, USA.

0377-0257/02/$ – see front matter © 2002 Elsevier Science B.V. All rights reserved. PII: S 0 3 7 7 - 0 2 5 7 ( 0 2 ) 0 0 1 6 2 - 3

14

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

1. Introduction Instabilities in viscoelastic flows are emerging as a rich class of fluid mechanical phenomena that are of much importance in polymer processing and coating operations. These instabilities often cannot be explained by consideration of only the more common mechanisms, such as inertia or surface tension; and a complete understanding of these flow transitions requires consideration of the non-Newtonian behavior of complex fluids such as polymer melts and solutions. Because these instabilities are much less well understood than flow transitions observed in Newtonian fluids, it is useful to examine instabilities where the transition is driven solely by the interaction between large normal stresses in the fluid and curved streamlines. This mechanism for viscoelastic flow transitions has been demonstrated in simple model geometries such as circular Couette flow and Dean flow [2–5]. For both of these geometries, the instability takes the form of toroidial vortices superimposed on the one-dimensional steady base flow, and a linear stability analysis with an Oldroyd-B model can qualitatively predict the onset and structure of the transitions observed in experiments performed with Boger fluids [2–4]. However, although the curved streamline instability mechanism can also qualitatively explain experimentally observed transitions in viscoelastic flows in more complex geometries, such as lid-driven cavity flow [6] and flow around a cylinder confined in a channel [7], most of the progress in the analysis of the stability of viscoelastic flows has been restricted to simple rectilinear flows. The analysis of viscoelastic flows in which the base flow is one-dimensional is straightforward, because the equations governing the linear stability of the base flow are separable. As a result, the linear stability problem can be written as a generalized eigenvalue problem [5] or as a set of ordinary differential equations [2] with a relatively small number of degrees of freedom. These linear problems can be solved on a desktop computer by using standard numerical analysis. By contrast, the methods which have been used to analyze one-dimensional viscoelastic flows are difficult to apply to viscoelastic flows in non-trivial geometries where the base state is not known in closed form and spatial discretization of the flow geometry leads to a large number of degrees of freedom. This paper is an extension of our earlier work [8,9] that demonstrated numerical methods suitable for determining the linear stability of steady, two-dimensional, viscoelastic flows to two- and three-dimensional perturbations. In this method, the two-dimensional, steady base flow of an Oldroyd-B model is calculated by using the DEVSS-G/SUPG finite element discretization. The linear stability of the base state is then determined by time integration of the evolution equations for the perturbations generated by linearization of the equations of motion and of the constitutive equation about the base flow. The long-time growth rate of the perturbations to the kinematics and to the polymer contribution to the stress gives an estimate for the real part of the leading eigenvalue, and the form of the perturbations at long times gives an approximation to corresponding eigenfunction. We also developed [8,9] an iterative eigenvalue calculation based on the Arnoldi method to compute a subset of the eigenspectrum corresponding to the eigenvalues associated with the most dangerous flow disturbance. Note that when we refer to the “most dangerous” eigenvalue, we mean the eigenvalue with the largest growth rate (i.e. the greatest real part), not the eigenvalue with the largest magnitude; the former will have a real part equal to zero at neutral stability, and the latter is related to the spectral radius of the operator and will likely have a very large, negative real part, which will be strongly stabilizing. In this study, we apply these numerical methods to examine the linear stability of the flow of an Oldroyd-B fluid through a linear array of cylinders confined in a channel, shown in Fig. 1. This flow of

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

15

Fig. 1. Flow past a linear array of cylinders confined in a channel. The dimensionless cylinder spacing is given by L, and the dimensionless channel height is 2H. Both dimensions are scaled by the cylinder radius, R.

a polyisobutylene Boger fluid has been studied experimentally by Liu [1]. Liu found a transition from a two-dimensional flow to a three-dimensional, time-dependent state at a critical value of the Weissenberg number, Wecrit , that depends on the spacing of cylinders. However, the exact nature of this interaction is complicated, because the kinematics is a combination of elongation and shear, and the streamwise curvature changes sign as the fluid moves over the top of the cylinder and into the cylinder wake. The goal of this study is to calculate the onset of instability of the two-dimensional base flow to three-dimensional disturbances and to analyze the form of the predicted instability to determine if the mechanisms for the circular Couette and Dean flows can be applied. We review the experimental results of Liu [1] in order to make comparisons with our numerical results. The flow is described in a Cartesian coordinate system where x is the streamwise direction, y is the cross-channel direction, and z is the neutral direction. The origin is at the center of one of the cylinders, with the cylinder radius, R, chosen as the characteristic length scale. The geometry is then described by the dimensionless channel half-height, H, and the dimensionless distance between the cylinders, L. In both the experiments by Liu [1] and the calculations presented here, the channel half-height is H = 2. Liu performed experiments for L = 2.5, L = 6, and L = ∞ (an isolated cylinder). The characteristic velocity of the flow is the cross-channel average velocity 
H  1 L V ≡ (1) vx − , y, 0 dy, 2H −H 2 where vx (x, y, z) is the streamwise velocity. The definition of the Weissenberg number is We ≡

λV , R

(2)

where the time constant of the fluid, λ, is defined as λ≡

Ψ1 . 2(η0 − ηs )

(3)

In Eq. (3), Ψ 1 is the first normal stress coefficient, η0 is the viscosity of the fluid, and ηs is the viscosity of the solvent. The definitions for the time constant and for the Weissenberg number are the same as those used by Liu, except that for the polyisobutylene Boger fluid, the first normal stress coefficient and the viscosity are shear rate dependent, so that the time constant must be evaluated at the characteristic shear rate of the flow, V/R. In the zero-shear-rate limit, the ratio of the solvent viscosity to the total viscosity of the fluid, β ≡ ηs /η0 , is equal to 0.67, the value used in the experiments [1]. In the experiments for closely spaced cylinders (L = 2.5) and small values of We, the flow is steady, and a recirculating region is observed between the cylinders. Laser Doppler velocimetry measurements

16

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

for We of 0.68, 1.32, and 1.47 indicate that the flow is essentially two-dimensional in the central region away from the sidewalls of the channel. A transition to a three-dimensional flow is detected and the flow becomes time-periodic at a critical value of the Weissenberg number, Wecrit = 1.53. The disturbance in the streamwise velocity is localized near the recirculation region with a dimensionless wave number of kz,x = 3.3 ± 0.7 while the disturbance to the velocity in the neutral direction vz has a wave number of kz,x = 2.7 ± 0.6. The dimensionless time period (scaled with R/V) of the new time-periodic state is measured to be T ∼ = 1950 ± 730. The transition between the three-dimensional, time-periodic state and the two-dimensional, steady-state occurs at Wecrit = 1.53 both as the flow rate is increased or decreased, indicative of a supercritical bifurcation from the steady base flow to the time-periodic state. Liu has also partially determined the dependence of this instability on cylinder spacing, L. Experiments were performed for L = 2.5, L = 6, and for an isolated cylinder (L → ∞), all with the dimensionless channel half-height of H = 2. Most interesting is the observation that the critical value, Wecrit , decreases with increasing L. For a single cylinder in a channel (L → ∞, H = 2) Wecrit was measured to be 0.53. A detailed comparison between the mechanism for the viscoelastic instability observed in the flow around the linear array of cylinders and the mechanisms for the viscoelastic transitions observed in the circular Couette flow and the Dean flow is reported here. Hence, it is useful to review both of these mechanisms. For both geometries, the viscoelastic transition is driven by perturbations to the radial velocity, vˆ r , that ultimately lead to perturbations in the “hoop” stress, τˆpθθ , where we adopt the notation that perturbation and steady-state variables are designated with a carat and the superscript “ss”, respectively. The hoop stress interacts with the curvature in the streamlines, given by 1/r, through the term τˆpθθ /r in the r-momentum equation, reinforcing the perturbation to the radial component of the velocity. However, the mechanism by which the perturbation vˆ r leads to the perturbation τˆpθθ are different for the Couette flow and for the Dean flow. In the circular Couette flow, the perturbations vˆ r are followed by a coupling between perturbations to the polymer stress τˆ p and the base state velocity gradient ∇v ss [2]. For the axisymmetric mode, perturbations to the radial extensional flow ∂ˆvr /∂r lead to perturbations in τˆprr . Next, the perturbation τˆprr couples ˆprθ that also interacts with with the base state shearing flow ∂vss θ /∂r, which gives rise to the perturbation τ the base shearing flow ∂vss /∂r to finally give rise to the perturbations to the hoop stress τˆpθθ [4]. For θ the non-axisymmetric mode, perturbations to ∂vr /∂θ lead to perturbations in the shear stress τˆprθ . Next, the perturbations τˆprθ couple with the base state shearing flow ∂vss θ /∂r to drive perturbations to the hoop stress τˆpθθ [2]. For both the axisymmetric and the non-axisymmetric Couette instabilities, these multi-step mechanisms lead to time-periodic states, i.e. the leading set of eigenvalues in the linear stability analysis is a complex conjugate pair. In the viscoelastic Dean flow, both the outer and inner cylinders are stationary and the base flow is driven by an azimuthal pressure gradient. This leads to a radial variation in the stress across the gap between the cylinders. The mechanism for the transition in Dean flow is described by the coupling of the ss perturbation velocity field vˆ r to the polymeric stress gradient in the base flow dτpθθ /dr that directly drives perturbations to the “hoop” stress, τˆpθθ [3]. Joo and Shaqfeh [3] note that the condition ss dτpθθ

dr

< 0,

(4)

is a necessary condition for this type of instability mode. The mechanism for the elastic Dean instability is quite distinct from the Couette mechanism, because this coupling is negligible in the Couette flow

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

17

where the gradients in the base state stress are small. Furthermore, this single-step mechanism leads to a transition to a steady-state flow, i.e. the leading eigenvalue in the linear stability analysis has zero imaginary part. In Section 2, the governing equations are given, and the numerical methods developed in [9] are briefly outlined, along with several recent improvements that lead to greater computational efficiency. A Protean coordinate system is defined in Section 3 that makes explicit the interactions between the normal stress and the curvature of the streamlines. This coordinate system is used to analyze the linear stability results. In Section 4, linear stability results are presented for flow through the cylinder array, and we compare these results with experimental observations in Section 5. In Section 6, the mechanism for the transition from the two-dimensional base flow to a three-dimensional flow for an Oldroyd-B model is analyzed. Finally, in Section 7, we summarize our results and offer our conclusions.

2. Governing equations and numerical methods We model the rheology of the MIT Boger fluid used in the experiments of Liu [1] with the Oldroyd-B constitutive equation. The Oldroyd-B model may be derived from the kinetic theory of polymeric liquids by modeling the fluid as a dilute solution of bead-spring dumbbells with a linear spring force. The governing equations for the DEVSS-G finite element formulation in dimensionless form are written as   ∂τ p τ p + We (5) + v · ∇τ p − {GT · τ p + τ p · G} = −(1 − β)[G + GT ], ∂t −∇ 2 v + (1 − β)∇ · [G + GT ] + ∇ · τ p + ∇p = 0,

(6)

∇ · v = 0,

(7)

G − ∇v = 0,

(8)

where v is the dimensionless velocity vector, p and τ p are the dimensionless pressure and polymeric contribution to the extra stress tensor, and G is the velocity gradient. The time constant in the Oldroyd-B model is consistent with the definition given by Eq. (3). The viscosity and the first normal stress coefficient predicted by the Oldroyd-B model are independent of shear rate. The dimensionless solvent viscosity is selected as β = 0.67 in order to match the zero-shear-rate properties of the MIT Boger fluid [1]. The numerical methods used to analyze the viscoelastic flow are only described briefly; more detail is available in [9]. Steady-state, two-dimensional base flows are calculated by using the DEVSS-G finite element discretization of Eqs. (6)–(8), combined with the streamwise upwind Petrov–Galerkin (SUPG) discretization of the hyperbolic constitutive equation (Eq. (5)). The finite element representation is constructed from piecewise-continuous, low-order polynomials based on a standard mixed-method for viscous flows. Lagrangian biquadratic approximations are used for the velocity components and bilinear approximations are applied for the pressure, polymer stress, and velocity gradient interpolants. ss Transitions from the steady-state, two-dimensional base flow, denoted as [τ ss p (x1 , x2 ), v (x1 , x2 ), ss ss p (x1 , x2 ), G (x1 , x2 )], to a three-dimensional, possibly time-dependent state are computed by considering the linear stability of the base flow to a set of small, random perturbations. Each variable is

18

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

written as a sum of the base state and a small amplitude disturbance as     ss    ˆ τ (x, y, t) τ p (x, y) τ p (x, y, z, t)   p       v(x, y, z, t)   v ss (x, y)    ˆ v (x, y, t)    ikz z    , = e   + ε Re  ˆ y, t)     p(x, y, z, t)   pss (x, y)   p(x,       ˆ G(x, y, z, t) Gss (x, y) G(x, y, t)

(9)

ˆ is the disturbance and ε  1. The form of the disturbance in the neutral direction ˆ G) where (τˆ p , vˆ , p, is written in Fourier form with wave number kz , and the two-dimensional (x, y)-dependence of the ˆ is approximated with the same finite element representation as used to represent perturbation (τˆ p , vˆ , p, ˆ G) the steady base state. Time dependence of the perturbation is included in the unknown coefficients in the finite element representation of each variable. The equation set (5)–(8) is linearized about the base state solution for ε  1, and the resulting equation set is discretized by the DEVSS-G/SUPG method as described in [9]. This procedure results in a set of linear differential-algebraic equations (DAEs) for the finite element ˆ which are represented as coefficients in the approximations for the variables (τˆ p , vˆ , p, ˆ G) M·

du = A · u, dt

(10)

ˆ and M ≡ diag[We M, ˜ 0, 0, 0] is the mass matrix where u is the vector of unknowns, u ≡ [τˆ p , vˆ , p, ˆ G], ˜ is the mass matrix for the constitutive equation, formed by taking the inner for Eq. (10). The matrix M product of the test functions for the constitutive equation with the approximating functions for the stress tensor. The DAE set Eq. (10) is integrated in time by using the θ-method operator-splitting technique [8,10,11]. Each θ-method step consists of three substeps that are concisely written in terms of two linear operators A1 and A2 (A ≡ A1 + A2 ) as un+θ − un un+1−θ − un+θ = A1 (un+θ ) + A2 (un ), M· = A1 (un+θ ) + A2 (un+1−θ ), θ t (1 − 2θ)t un+1 − un+1−θ = A1 (un+1 ) + A2 (un+1−θ ), (11) M· θ t where A1 and A2 are given by M·



−ωτˆ p



  ˆ +G ˆ T)   −∇ 2 vˆ + ∇ · τˆ p + ∇ pˆ + (1 − β)∇ · (G , A1 ≡    ∇ · vˆ   ˆ − ∇ vˆ G   ss ˆ ˆ T v · ∇τ ss ˆ +G ˆ T )] ˆ p Gss }−{τˆ p Gss }T ) − (1−β)(G [−(1−ω)τˆ p − We(v ss · ∇ τˆ p −{τ ss p · G} − {τ p · G} +ˆ p −{τ   0   A2 ≡  .   0 0

√ The parameters θ and ω are chosen as θ ≡ (1 − 2)/2 and ω ≡ (1 − 2θ)/(1 − θ) following [9].

(12)

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

19

The θ-method is an efficient time integration scheme because it decouples at each time step the stress ˆ The first substep described by Eqs. (11) update from the update for the kinematic variables (ˆv , p, ˆ G). and (12) consists of an explicit update for the polymer stress (τˆ p ) followed by the solution of a modified ˆ in which the polymer stress is treated as data. The second Stokes problem for the variables (ˆv , p, ˆ G), substep requires an implicit update of τˆ p , and the third substep is identical to the first. Each of these steps involves the solution of a large, sparse set of linear algebraic equations that is simply represented ˆ In our previous paper [8] and for the results presented in [9], we have solved each of the ˆ · uˆ = b. as A linear subproblems with the frontal elimination subroutine MA42 from the Harwell Subroutine Library ˆ =L ˆ · Uˆ , is written to disk for reuse with other [12]. The LU-decomposition for each subproblem, i.e. A ˆ right-side vectors b. Although the θ-method is much more efficient than a fully implicit time integration method, the modified Stokes problem that typically accounts for more than 90% of the computational cost (CPU time) for each θ-method step becomes very large, especially for the full channel simulations of a widely spaced array of cylinders. Furthermore, although the operation count for the solution of the triangular linear systems may be relatively small compared with the cost of the initial LU-decomposition, the computational time can be unacceptably large due to the access time required to read the LU-decomposition of A from storage. For these reasons, an efficient, parallel time integration scheme for calculating viscoelastic flows developed by Caola et al. [13] is applied: the linear system that results from the generalized Stokes problem is solved by using an algorithm that combines a parallel preconditioner with a Krylov subspace method. This block complement and additive levels method (BCALM) preconditioner is based on treating pressure unknowns separately from the velocity and gradient variables, and the resulting iterative method is demonstrated to have high parallel efficiency. Robustness and good performance result from using geometric partitioning of the unknowns by the nested bisection routines from CHACO [14]. The implicit update for τˆ p was calculated by using preconditioned BiCGStab. Because the stress problem is not block singular, two-level additive Schwarz was utilized without the block complement operation. This combination of the θ-method operator-splitting technique with the iterative, parallel solution of the modified Stokes problem produces an efficient time integration algorithm. The linear stability of the calculated base flow can be determined by long-time integration of the linearized equations for the perturbations starting from random initial conditions for the disturbance of the polymer stress. This leads to a growth rate of the perturbations equal to the real part of the most dangerous eigenvalue and a structure of the long-time solutions that gives the corresponding eigenfunction. We used this approach to analyze the stability of planar Couette flow of an Oldroyd-B fluid to two-dimensional perturbations [8] and the stability of circular Couette flow of an Oldroyd-B fluid to three-dimensional perturbations [9]. For both test problems, the long-time growth rates were in good agreement with the real parts of the leading eigenvalues obtained by using other spectral (Chebyshev) methods of analysis [5,8]. A Krylov subspace method also is used for computing a subset of the eigenvalues of the generalized eigenvalue problem that results if the long-time response of the perturbations is assumed to have the form eσ t . Time integration of the linearized equations for the perturbations can be written as a linear transformation, denoted as TN , that advances the variables N time steps equal to a period of time t ≡ N t. A second eigenvalue problem is defined for the operator TN as TN s = µs,

(13)

20

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

where s is a vector containing the coefficients for the finite element representation of τˆ p and µ is an eigenvalue. The eigenvalues of the linear operator TN are approximately related to the eigenvalues of the original linear stability problem, σ, as µi ≈ e(N t)σi .

(14)

We apply the Arnoldi method [15,16] to the eigenvalue problem described by Eq. (13) by using the public domain code ARPACK [17]. The combination of time integration, as represented by the linear operator TN , with the Arnoldi algorithm is an efficient means to calculate accurate approximations for the leading eigenvalues and eigenvectors. This approach has been described in [9], and also has been used by Ramanan et al. [18] to determine the stability of a time-periodic viscoelastic flow. We have demonstrated [9] for calculations of the linear stability of circular Couette flow of a viscoelastic fluid that the hybrid Arnoldi/time integration approach leads to eigenvalues and eigenvectors that are in good agreement with other results in the literature [19]. 3. Protean coordinate system For all of the intercylinder spacings investigated in this study, the instability is driven by an interaction between the normal stresses and the curvature of the streamlines. The streamlines in the base flow demonstrate large curvature in two distinct regions of the flow: the curvature is negative as the fluid flows over the top of the cylinder, and the curvature is positive in the gap between the cylinders. Whereas the calculations are performed in a Cartesian coordinate system, the analyses of the stability results and the base states are simplified if the stress is decomposed by using basis vectors that are the unit tangents and the unit normals to the streamlines for the two-dimensional base state—the so-called Protean coordinates [20]. The unit vectors describing this coordinate system are vss vss y x e1 ≡ ss ex + ss ey , ||v || ||v ||

e2 ≡

−vss y ||v ss ||

ex +

vss x ey , ||v ss ||

(15)

where e1 points in the streamwise and e2 in the cross-stream directions. These coordinates are used to ss ss , and the shear component of the stress τ12 , with define the streamwise normal component of the stress τ11 ss ss respect to the material surfaces of the fluid. The polymer contributions for the base flow (τ11 , τ12 ) are shown in Fig. 2 for L = 2.5 and We = 1.46. A single streamline is shown in Fig. 2 that is very close to the boundary of the recirculation found in the gap between the cylinders. Both the shear and normal components of the polymer stress relax as the fluid moves across the top of the recirculation region. For this value of We, the relaxation time of the polymer is sufficiently long relative to the time required for a fluid particle to traverse the gap between the cylinders that stress is convected from one cylinder to the adjacent downstream cylinder. The constitutive equation and momentum equation are written in the Protean coordinate system with basis vectors given by Eq. (15). The coordinate surfaces are defined by the partial differential equations Dq1 = ||v ss ||, Dt ∂q2 = −vss y, ∂x

(16) ∂q2 = vss x, ∂y

(17)

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

21

Fig. 2. Contours of the 11-component (left) and 12-component (right) of the polymer stress for L = 2.5 and We = 1.46 calculated with mesh M2. The heavy black streamline shows the position of the recirculating region between the cylinders. On the left, 15 equally spaced contours between −150 and −10 are shown, and 15 equally spaced contours between −4.6 and 3 are shown on the right.

where q1 is the arc length along a streamline and q2 is the streamfunction. The scale factors for the new coordinate system are h1 = 1,

h2 =

1 , ||v ss ||

(18)

and the gradient operator is ∇ = e1

∂ ∂ + e2 ||v ss || . ∂q1 ∂q2

(19)

The derivatives of the basis vectors with respect to q1 are given by the Frenet formula as ∂e1 = κe2 , ∂q1

∂e2 = −κe1 , ∂q1

(20)

where κ is the curvature and is given by κ≡

2 ss ss ss ss ss 2 ss (vss x ) Gxy − 2vx vy Gxx − (vy ) Gyx

||v ss ||3

=

e1 · Gss · e2 , ||v ss ||

(21)

where the velocity gradient tensor Gss is given by Eq. (8). The derivatives of the basis vectors with respect to q2 are calculated from ∇e1 and ∇e2 in Cartesian coordinates by rewriting these in terms of e1 and e2 with the aid of Eq. (15) and then subtracting the first term in Eq. (19), which is given by Eq. (20). This calculation yields ∂e1 = −Γ e2 , ∂q2

∂e2 = Γ e1 , ∂q2

(22)

22

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

where Γ ≡

e2 · Gss · e2 e1 · Gss · e1 = − . ||v ss ||2 ||v ss ||2

(23)

The momentum equation is written in Protean components as   ∂G11 ∂G12 ∂G1z ∂p −β − vΓ(G11 − G22 ) + v − 2κG12 + ∂q1 ∂q1 ∂q2 ∂z +

∂τp11 ∂τp12 ∂τp1z = 0, − vΓ(τp11 − τp22 ) + v − 2κτp12 + ∂q1 ∂q2 ∂z

(24)

  ∂G12 ∂G22 ∂G2z ∂p −β − 2vΓG12 + κ(G11 − G22 ) + v + ∂q2 ∂q1 ∂q2 ∂z +

∂τp12 ∂τp22 ∂τp2z − 2vΓτp12 + κ(τp11 − τp22 ) + v + = 0, ∂q1 ∂q2 ∂z

(25)

  ∂p ∂G1z ∂G2z ∂Gzz −β − vΓG1z + v − κG2z + ∂z ∂q1 ∂q2 ∂z ∂τp2z ∂τp1z ∂τpzz − vΓτp1z + v + − κτp2z + = 0, ∂q1 ∂q2 ∂z

(26)

for the 1, 2, and transverse or z-components, respectively. The equations describing the Oldroyd-B model also are written in the Protean coordinate system. When this is done the constitutive equation becomes a set of ordinary differential equations, with the only derivative in the 1-direction. To write these equations, we need the convected derivative of the perturbations to the stress tensor. Because the stress tensor is symmetric, the six components of the tensor (v ss · ∇ τˆ p ) are {v ss · ∇ τˆ p }11 = ||v ss ||

∂ˆτp11 − 2||v ss ||κˆτp12 , ∂q1

{v ss · ∇ τˆ p }12 = ||v ss ||

∂ˆτp12 + ||v ss ||κ(ˆτp11 − τˆp22 ), ∂q1

{v ss · ∇ τˆ p }1z = ||v ss ||

∂ˆτp1z − ||v ss ||κˆτp2z , ∂q1

{v ss · ∇ τˆ p }22 = ||v ss ||

{v ss · ∇ τˆ p }2z = ||v ss ||

∂ˆτp2z + ||v ss ||κˆτp1z , ∂q1

{v ss · ∇ τˆ p }zz = ||v ss ||

∂ˆτp22 + 2||v ss ||κˆτp12 , ∂q1

∂ˆτpzz . ∂q1

(27)

The relationships outlined above make explicit the interactions between the different components of the stress and the curvature κ, as well as the interactions between the different components of the stress and the variable Γ . These equations prove as useful for describing the interaction between curved streamlines and normal stress for general flows with a two-dimensional base state as cylindrical coordinates are in

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

23

describing the circular Couette and Dean flows. We will use this coordinate system in the analysis of the results in the following sections. 4. Stability of flow around a linear array of cylinders 4.1. Cylinder spacing of L = 2.5, symmetric case The stability of the flow of an Oldroyd-B fluid through a linear array is analyzed for an intercylinder spacing of L = 2.5, a geometry studied in the experiments of Liu [1]. The disturbances are assumed to be reflectively symmetric about the centerplane of the flow (y = 0). Results for additional cylinder spacings and non-symmetric perturbations are presented in Section 4.2. Three finite element meshes have been used in this study (labeled M1, M2, M3). Mesh M2 is shown in Fig. 3. Each mesh in this sequence contains 50% more refinement than the previous mesh in each coordinate direction, with the additional elements being uniformly distributed across the entire domain. The important dimensions of each mesh are listed in Table 1. The size of the element directly on top of the cylinder is important for resolving the dynamics of the shear flow that develops between the cylinder and the channel wall, whereas the element at (x, y) = (1.25, 0.40) is located along the top of the recirculating region of the flow between adjacent cylinders where the instability is observed experimentally.

Fig. 3. Finite element mesh M2 for calculation of the flow of an Oldroyd-B fluid through a linear, periodic array of cylinders confined in a channel with H = 2 and L = 2.5. Table 1 Mesh dimensions for meshes (M1–M3) for the study of the flow through a linear array of cylinders with L = 2.5. Critical values of the Weissenberg number for three-dimensional perturbations with kz = 3.125 and β = 0.67 Mesh

x at (0, 1)

y at (0, 1)

x at (1.25, 0.40)

y at (1.25, 0.40)

d.f. in modified Stokes problem

Wecrit

M1 M2 M3

0.1065 0.0784 0.0523

0.0165 0.00917 0.00444

0.0391 0.0311 0.0198

0.0606 0.0424 0.0295

38,838 77,568 171,372

1.395 1.434 1.475

24

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

Fig. 4. Evolution of the norms of the perturbations to the stress (filled symbols) and to the velocity (open symbols) for flow through the periodic array with L = 2.5 and We = 1.46. The wave numbers kz of the perturbations in the neutral direction are 0 (䉬), 1.25 (䊏), 2.5 (䊉), 3.125 (
), 3.75 (), 5 (䉱), and 7.5 (䉲).

It is interesting at this point to summarize the computational cost of these calculations. Most of the calculations were performed on a cluster of 500 MHz Pentium III Xeon PCs. When the frontal solver was used to perform time integration with mesh M2, a single time step required 53 s, whereas a single time step with the parallel time integration scheme required 76 s on a single dual-processor machine (using both processors), and 20 s on seven dual-processor machines. Thus, a typical calculation, such as required for one of the results shown in Fig. 4, required about 15 h with a single processor using the frontal method, or about 6 h on seven dual-processor machines. It should also be noted that for larger calculations, the parallel method became much more efficient, as shown in [13]. In addition, multiple processors were required in order to have enough memory for the full channel results presented in Sections 4.2 and 4.3. For more information regarding efficiency and implementation of the parallel time integration scheme, the interested reader is referred to [13]. Steady, two-dimensional base states were calculated with meshes (M1–M3) by using the DEVSS-G/ SUPG discretization combined with Newton’s method. The linear stability of the two-dimensional base flow was analyzed by calculation of the leading eigenvalue by using the time integration and Krylov methods described in Section 2. For mesh M2, time integration of the linear evolution equations for the perturbations was performed for several transverse wave numbers, kz . Sample dynamics for We = 1.46 for the norms of the perturbations to the stress and the velocity as functions of time are shown in Fig. 4. The norm used here is the Euclidean norm of the vector of coefficients in the finite element representation for the perturbations to the polymer stress and to the velocity. For all wave numbers, numerically stable calculations were performed with a dimensionless time step of t = 0.0375. The growth rates of the perturbations for each wave number are shown in Fig. 5 as calculated from the long-time response. The points shown in Fig. 5 are the average of the growth rates of the stress and velocity norms. The error bars show the actual values of the growth rate for the stress norm and for the velocity norm. If only one eigenfunction were present, then the growth rate of the stress and velocity norms would be the same, so that differences between the two growth rates give an indication of the possible error in using the growth rate of a measure of the flow (such as shown in Fig. 4) as an estimate for the real part of the dominant eigenvalue.

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

25

Fig. 5. Growth rate for the perturbations to the stress and the velocity for L = 2.5 and We = 1.46 as a function of kz , obtained by fitting exponential curves to the norms in Fig. 4.

The results in Fig. 5 clearly indicate that disturbances with a range of kz are unstable for We = 1.46 and L = 2.5. The perturbations with kz = 3.125 are the fastest growing, and the sharpness of the peak in the growth rate as a function of kz demonstrates relatively strong wave number selection. Because simulations were carried out at intervals of 0.625 in kz , the wave number for the maximum growth rate is estimated as kz = 3.1 ± 0.3. Time integration for perturbations with kz = 3.125 was repeated with all three finite element meshes and for a range of We. Interpolation of the long-time growth rates yielded estimates of the critical Weissenberg number, Wecrit , for which Re(σ) = 0. These results are shown in Table 1 for the three meshes. The values of Wecrit are reasonably consistent for each mesh. The consistency of the spatial structure of the eigenfunctions with mesh refinement also is demonstrated in Fig. 6, where the trace of the perturbations to the polymer stress is shown along the line (x, y, z) = (L/2, y, 0) that is midway between the cylinders. Gibbs’ phenomenon was observed in the calculation near y = 0.4 where the magnitude of the stress in the base state changes rapidly with y; however, the region over which the oscillations occur shrinks as the mesh is refined. The best estimate of Wecrit was obtained by extrapolation of the critical values from each mesh to the limit of zero element size. This leads to an estimate of the critical Weissenberg number of Wecrit = 1.55. A subset of the leading eigenvalues in the linear stability problem also were calculated by using the Arnoldi algorithm to demonstrate that the eigenvalue associated with the instability computed by time integration is distinct from the continuous part of the eigenspectrum expected for a shear flow. The continuous part of the eigenspectrum arises in all shearing flows of viscoelastic fluids [21] and appears in the complex plane as a line from −1/We to −(1/We) − iα where α is the streamwise wave number of the disturbance. It is difficult to estimate accurately the eigenvalues associated with the continuous

26

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

Fig. 6. Trace of the perturbation to the polymer stress along the line x = L/2 and z = 0 for meshes M1 (䊐), M2 (䊊), and M3 () for L = 2.5 and We = 1.46. Also shown is Tr(τ ss p ) for mesh M2 (䊉).

part of the spectrum, because the associated eigenfunctions have singular spatial structure and converge slowly with mesh refinement [22]. We have demonstrated [8,9] that if the mesh is not properly refined, a fictitious, numerical instability can arise. The relationship between the most dangerous eigenvalue in the stability of the flow through the cylinder array and the continuous part of the eigenspectrum was examined through two sets of calculations. First, estimates for the nine eigenvalues with the largest real part (one real eigenvalue and eight complex conjugate pairs) are calculated by the Arnoldi method. A large tolerance was used in these calculations because only rough estimates for the eigenvalues from the continuous spectrum are desired; an eigenvalue was considered converged if the Ritz estimate for the error [17] in the eigenvalue was smaller than 1×10−2 . This calculation was repeated with a stricter tolerance as a stopping criterion for the Arnoldi method (Ritz estimates smaller than 1 × 10−4 ). For this stricter tolerance, only the leading eigenvalue converged. Both calculations were performed on mesh M2. The most dangerous eigenvalue had zero imaginary part Im(σ) = 0 with a real part Re(σ) in good agreement with the long-time growth rate obtained by time integration. The eigenfunctions corresponding to the subdominant modes were shear-dominated, as demonstrated by small, mesh-sized cellular structures along the channel wall and the cylinder surface. The structures of the perturbations to the velocity and stress in the subdominant modes were identical to those observed in [9] where we examined two-dimensional perturbations to the flow through linear arrays of cylinders in which only members of the continuous spectrum were found. Similar structures in the perturbations to the stress and velocity were observed by Brown et al. [23] along the surface of a sphere falling in a tube. The separation of the most dangerous eigenvalue from the continuous part of the spectrum is shown in Fig. 7 where the eigenspectra are presented for several values of We, as calculated with mesh M2. The separation of a discrete mode from the continuous spectra is very different from the eigenspectra calculated for the stability of the cylinder array to two-dimensional perturbations, where only a continuous

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

27

Fig. 7. Eigenvalue estimates for flow through the cylinder array with L = 2.5 are shown on the left for We = 1.46. On the right are eigenvalue estimates for We = 1.37 (䊐), We = 1.46 (), and We = 1.55 (䊊).

spectrum was observed [9]. In addition to the results presented here, calculations also were performed with a Runge–Kutta time integrator combined with the discontinuous Galerkin discretization of the polymer contribution to the stress. The predictions of stability and spatial structures of the perturbations at long times were very similar to the results presented here. We believe that the instability calculated here is an accurate reflection of the dynamics of the Oldroyd-B model and not a numerical artifact. This conclusion is supported by the facts that the calculations are reproduced on several different meshes, with two different spatial discretizations and two different time integrators, and the unstable mode is discrete from the continuous portion of the spectrum. Throughout the rest of this paper, results are not presented for mesh refinement studies that have been performed for other cylinder spacings. Instead, the presentation focuses on the calculated spatial structure of the most dangerous perturbation and the mechanism that drives the flow instability. As shown by the plots in Figs. 8 and 9, the intensities of the perturbations to the polymer contribution to the stress and to the velocity for the most dangerous disturbance are concentrated directly above the recirculation region and along the fore and aft surfaces of the cylinder. The components of the perturbations to the polymer stress with the largest relative magnitude (referred to there as the dominant components) are the streamwise normal stress τˆp11 , the shear stress in the xy-plane τˆp12 , and the shear component τˆp1z . The phase of the perturbations in the neutral direction is portrayed in Fig. 10. The measure Tr(τˆ p ) is in phase with the perturbations to the streamwise component vˆ 1 and cross-stream component vˆ 2 of the velocity, but out of phase with vˆ z . The phase of the perturbations also is shown in Fig. 9, where contour plots are shown in the xy-plane for a value of z where the perturbations have the largest magnitude. This figure demonstrates that the shear components (ˆτp1z , τˆp2z ) are in phase with vˆ z and out of phase with the other components of τˆ p .

28

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

Fig. 8. Magnitude of the perturbations to the velocity, ||ˆv ||, in the plane z = 0 (left) and in the plane kz z = π/2 (right) for mesh M3 with kz = 3.125 and We = 1.56. In the plane z = 0, the perturbations to vˆ z are zero, whereas in the plane kz z = π/2, the perturbations to vˆ 1 and vˆ 2 are zero. Regions with small perturbations are dark, and regions with large perturbations are light.

Along the top of the recirculation boundary and near kz z = π, the positive perturbations to the normal stress decrease the magnitude of the large, negative normal stresses in the base flow. Because the streamlines are curved, τˆp11 appears in the linearized momentum balance for the cross-stream velocity through the term κˆτp11 in the linearized form given by Eq. (25), and this drives the flow away from the center of curvature, which is observed in Fig. 10 as the negative value of vˆ 2 at kz z = π. This downward motion decreases the size of the recirculation region, and to conserve mass, fluid flows away from kz z = π in the positive and negative z-directions both above and below the recirculation boundary (see Fig. 8). The flow in the z-direction is concentrated above the rear separation point and in the recirculation region. Because the velocity is zero on the surface of the cylinder, shear stresses (ˆτp1z , τˆp2z ) develop on the surface of the cylinder that are in phase with vˆ z . As shown in Fig. 2, the normal stresses in the flow are generated along the top of the cylinder and convected into the cylinder wake. The decrease in the size of the recirculation region for kz z = π increases the path length along a perturbed streamline from the rear of one cylinder to the front of the next cylinder, and this reduces the amount of stress convected across the gap. Thus, the fluid motion reinforces the positive perturbation to the large negative normal stresses along the top of the recirculation boundary found in the base state. An estimate of the three-dimensional stress and velocity fields was obtained by multiplying the most dangerous eigenfunction calculated with the Arnoldi method by a small amplitude (ε), and adding the result to the steady two-dimensional base flow. An isosurface of the trace of the stress tensor obtained by this method is shown in Fig. 11. The periodic variation with the wave number kz = 3.125 is obvious. Estimates for the streamlines were calculated by a similar procedure and are shown in Fig. 12. The recirculating region has grown larger for values of z where the streamlines are converging, and the recirculating region has become smaller for values of z with diverging streamlines. It is interesting to note that the converging/diverging streamlines of the three-dimensional flow have a structure that is very similar to that depicted in Fig. 26 of the study by McKinley et al. [7], which describes the three-dimensional state observed for the flow of a Boger fluid around an isolated cylinder confined in a channel.

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

29

Fig. 9. Perturbations to the polymer contribution to the stress for mesh M3 with kz = 3.125 and We = 1.56. Each contour plot is for a value of z where the stress component is largest in the magnitude. For the top two rows, z = 0; and for the bottom row, kz z = π/2.

30

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

Fig. 10. The perturbations to the trace of the polymer stress (䊐) and to the velocity, vˆ 1 (䉱), vˆ 2 (䉬), and vˆ z (䊉) as a function of kz z along the line (1.25, 0.488, z).

4.2. Widely spaced array of cylinders: with and without symmetry Additional calculations were performed to determine the scaling of the critical Weissenberg number, Wecrit , with the cylinder spacing L. All calculations in this section were performed with a mesh similar to mesh M2 shown in Fig. 3, except that additional elements were placed along the boundaries x = ±L/2

Fig. 11. Isosurface for the trace of the polymer contribution to the stress tensor for We = 1.46 obtained by summation of the dominant eigenfunction with the two-dimensional base state. Flow is from upper left to lower right. The portion of the isosurface along the channel wall that obscured the view of the cylinders has been removed.

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

31

Fig. 12. First-order approximation to the streamlines in the three-dimensional flow formed by summation of the two-dimensional base state with the eigenfunction of the leading eigenvalue. Flow is from left to right, and three cylinders are shown viewed from above. Streamlines are placed at equal intervals in the z-direction a distance 1.1 radii above the axis of the leftmost cylinder. The Tr(τ p ) isosurface in Fig. 11 is also shown.

to maintain a level of mesh refinement with element dimensions similar to those presented in Table 1 for mesh M2. In order to compare the flow structure and the onset point of the instability as the intercylinder spacing L is increased, it is convenient to define the Deborah number as λV De ≡ , (28) RL which scales the time constant of the fluid λ with the time required to travel from one cylinder to the adjacent cylinder downstream. Keep in mind that L is dimensionless, and R is the (dimensional) cylinder radius. For De = 0.5, the recirculating region shrinks as L is increased until L = 3, where the single recirculation splits into two smaller circulating regions attached the fore and aft of the cylinder along the symmetry plane. For L = 3.25, only a single recirculation is found in front of the cylinder. The size of the recirculation region shrinks to less that 5% of the radius of the cylinder for L = 3.5. As the spacing between the cylinders is increased, the onset of the flow instability described in the previous section is calculated to occur at larger values of We, as shown in Fig. 13. A linear regression on the data in Fig. 13 gives Wecrit = 0.27L + 0.72,

(29)

or Decrit = 0.27 +

0.72 . L

(30)

Also, for L > 2.5, the dominant wave number shifts to kz = 2.5. The forms of the instability for all of the cylinder spacings in the range 2.5 ≤ L ≤ 3.5 are shown in Fig. 14, along with the streamline that approximately delineates the recirculating regions in the base state. The form of the perturbations to the stress for L = 4.0 is very similar to that shown in Fig. 14 for L = 3.5.

32

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

Fig. 13. The critical value of the Weissenberg number as a function of the cylinder spacing, L, for symmetric and off-centered systems. The shear-rate-dependent Wecrit obtained from the experiments by Liu is also shown for comparison.

Fig. 14. Tr(τˆ p ) for supercritical values of We for all of the cylinder spacings. Plots are shown for L = 2.5 with We = 1.46 (top left), L = 2.75 with We = 1.58 (top right), L = 3 with We = 1.64 (bottom left), and L = 3.25 with We = 1.68 (bottom right). Two base state streamlines are shown in each plot, one delineating recirculating regions and the other tracing a path over the cylinder and into the region in the wake where Tr(τ ss p ) has the largest magnitude.

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

33

Fig. 15. Perturbations to the polymer contribution to the stress, τˆp11 and τˆp12 for L = 3.5 with We = 1.80. Each contour plot is for a value of z where the stress component is largest in magnitude (z = 0).

For all values of L, the portion of the perturbation with the largest magnitude appears in the region between the cylinders. For L > 2.5 the perturbations are not concentrated on the recirculation boundary. The fact that the recirculation itself does not play a critical role in the mechanism for the flow transition is emphasized by the presence of the flow instability for L = 3.5 where the very small recirculating region has little or no influence on either the base state or the perturbations to the base state as shown in Fig. 15. For all cylinder spacings, the largest perturbation to τ p begins from midway between the cylinders and extends into the region where the fluid moves up the face of the downstream cylinder. Linear stability calculations also were performed without the assumption of symmetry about the centerplane (y = 0). Meshes used were similar to M2 and had approximately twice as many degrees of freedom than those with symmetry. Time integration of general, random perturbations to the polymer contribution to the stress for L = 2.5 and L = 3.5 led to growth rates that were very similar to those computed for the symmetric case. The form of the perturbations for t  1 appeared to be a mixture eigenfunctions for several closely spaced eigenvalues. In order to distinguish between these eigenmodes, the stability calculation for L = 2.5 was performed with the Arnoldi method. A single, three-fold degenerate discrete eigenvalue was computed that corresponded to two symmetric eigenfunctions and one non-symmetric eigenfunction, as shown in Fig. 16. We note that the flow structure without the assumption of symmetry after a long period of time integration is very similar to that with the symmetry condition at the centerplane (y = 0). 4.3. Effect of geometrical imperfection on the stability characteristics Whereas the flow structure after the transition is experimentally observed to be anti-symmetric about the centerplane, our linear stability analysis with the Oldroyd-B model does not predict preferential selection of an asymmetric mode. In this section, we present linear stability results for the geometry in which the cylinder is positioned slightly away from the centerplane. Again, all calculations in this section are performed with a mesh similar to mesh M2 (Fig. 3), except that the cylinders are displaced 5% off-center from the centerplane, so that the full channel must be considered.

34

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

Fig. 16. Shown in the top row are eigenfunctions for Tr(τˆ p ) calculated for the full channel with a centrally located cylinder for L = 2.5, We = 1.46, and kz = 3.125. The left plot shows the non-symmetric eigenfunction and the other plots in the top row show the two degenerate, symmetric eigenfunctions. Shown in the bottom row is the non-symmetric eigenfunction for the full channel with an off-centered cylinder for L = 2.5, We = 1.41, and kz = 2.5.

Our calculations reveal two important results. First, the growth rate increases when the cylinders are positioned off-center; and thus the flow becomes more unstable, as shown in Fig. 13. The critical We decreases roughly 10% regardless of the cylinder spacing. The critical wave number also shifts to lower values. In addition, the flow structure at a slightly supercritical condition shows that the most unstable mode at long times is an asymmetric one as seen in Fig. 16. It is found that for closely spaced cylinders (L = 2.5), the mode in the wider half-channel is strongly excited, whereas that in the narrower half-channel is substantially suppressed. For widely spaced cylinders, two asymmetric modes are located in the wider half-channel as shown in Fig. 17. These results demonstrate that a geometrical imperfection can lead to asymmetric flow structures that are close to those observed experimentally.

5. Comparison with experiments The predictions of the linear stability analysis reported here for an Oldroyd-B fluid can be compared, at least qualitatively, with the experimental observations of Liu [1]. The most extensive experimental

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

35

Fig. 17. Perturbations to the polymer contribution to the stress, τˆp11 , for centered (top row) and off-centered (bottom row) placement of the cylinders with L = 3.5 and kz = 2.5. For the symmetric geometry, We = 1.75; and We = 1.7 for the off-centered cylinder placement. Each contour plot is at a value of z or x for which the stress component is largest in magnitude (z = 0 for the xy-plane, and x = 1.25 for the yz-plane).

data were collected for the parameters L = 2.5, H = 2. The predictions of the calculations for the most dangerous wave number, shown in Fig. 5 as kz = 3.1 ± 0.3, are in reasonable agreement with the experiments. Liu measured kz,x = 3.3 ± 0.7 for the disturbance to the streamwise velocity component, vx , and kz,z = 2.7 ± 0.6 for the disturbance to the neutral velocity component vz . There is also remarkable agreement between the critical value Wecrit = 1.55 computed by extrapolation of the values in Table 1 and the experimental measurement of the shear rate-dependent value of Wecrit = 1.53. One potential discrepancy between the calculated linear flow transition and the experimentally observed three-dimensional state is that the calculated structure is a steady flow, whereas Liu [1] observed time-periodic fluctuations with a dimensionless period of T ∼ = 1950 ± 730. As shown in Fig. 4, an accurate estimate for the growth rate requires time integration only to t ∼ = 50. Simulation times two orders-of-magnitude longer would be required to estimate T accurately for the disturbance if the calculated

36

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

transition were oscillatory with the experimentally observed period. If calculations were performed until t ∼ = 2T with t = 0.0375, about 105 time steps would be required. As noted at the beginning of Section 4.1, if we used seven dual-processor workstations, each time step would take about 20 s to perform for mesh M2, and the entire calculation would require about 24 days. In addition, further studies would be required with different values of t in order to ensure that the result had converged with time step. Because of the enormous computation time, these calculations were not performed. The most obvious discrepancy between Liu’s experiments and the calculations with the Oldroyd-B fluid model is in the prediction of the dependence of Wecrit on the cylinder spacing L. As described by the correlation Eq. (29), the analysis with the Oldroyd-B fluid model predicts that Wecrit increases with increasing L, whereas the experiments suggest that Wecrit decreases with L. The discrepancy is consistent with the calculations reported by Smith et al. [9] in which the Oldroyd-B model failed to predict an instability expected for the isolated cylinder. 6. Mechanism of instability The mechanism for the flow instability is studied by detailed analysis of the results for the intercylinder spacing L = 2.5, and then these results are extended to account for the effect of intercylinder spacing. Here we use an energy analysis similar to the one described by Joo and Shaqfeh [3] that has been adapted for the analysis of the finite element results. First, the linearized constitutive equation for τˆ p is rewritten as   ∂τˆ p ss ss T ss ss ss ˆ ss ˆ T + v · ∇ τˆ p + vˆ · ∇τ p − {τ p · G} − {τ p · G} − {τˆ p · G } − {τˆ p · G } τˆ p = −We ∂t T

ˆ +G ˆ ), − (1 − β)(G

(31)

and is substituted into the momentum equation to give ˆ +G ˆ T ) − We∇ ∇ pˆ − ∇ · (G   ∂τˆ p ss ss T ss ss ss ˆ ss ˆ T · + v · ∇ τˆ p + vˆ · ∇τ p − {τ p · G} − {τ p · G} − {τˆ p · G } − {τˆ p · G } = 0. ∂t

(32)

A mechanical energy balance is formed by taking the inner product of the perturbation to the velocity vˆ with Eq. (32) and integrating the resulting equation over one period of the domain. Integration by parts yields
ˆ : ∂τˆ p dV = φτ p v1 + φτ p v2 + φτ p G1 + φτ p G2 + φvis , (33) G ∂t in which

 ˆ : (v ss · ∇ τˆ p ) dV, φτ p v1 ≡ − G  ˆ : (ˆv · ∇τ ss φτ p v2 ≡ − G p ) dV,  ss ˆ ˆ ˆ T φτ p G1 ≡ G : ({τ p · G} + {τ ss p · G} ) dV,  ˆ : ({τˆ p · Gss } + {τp · Gss }T ) dV, φτ p G2 ≡ G
1 ˆ : (G ˆ +G ˆ T ) dV. φvis ≡ − G We

(34)

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

37

The terms (φτ p v1 , φτ p v2 , φτ p G1 , φτ p G2 ) on the right side of Eq. (33) represent the transfer of energy from the base flow to the perturbations. The term φvis is the viscous dissipation found in the Reynolds–Orr energy equation [24]. It is always negative and stabilizes the flow. If the sum of all of these terms is positive, then the energy of the perturbations increases with time, and the flow is unstable. It is important to note that the Reynolds–Orr equation describes the rate of change of the kinetic energy of the perturbations to the flow, whereas Eq. (33) describes the exchange of energy between the base state and the perturbations in an inertialess flow. As in [2,3], the disturbance power created or dissipated by the perturbations to the polymer stresses is defined as
1 ˆ : τˆ p dV. E≡ G (35) 2 Because the energy analysis is performed on an eigenfunction obtained from either time integration or the Arnoldi method, the corresponding eigenvalue σ is known, and the time dependence of both τˆ p and ˆ will be the same. This implies that G   ˆ : ∂τˆ p = ∂ G ˆ : τˆ p = 1 ∂ (G ˆ : τˆ p ), G (36) ∂t ∂t 2 ∂t which leads to dE = φτ p v1 + φτ p v2 + φτ p G1 + φτ p G2 + φvis . dt

(37)

This formulation of the energy analysis is very similar to the result presented in [3], but here we have ˆ tensor, which is very formulated the terms on the right side of Eq. (37) as integrals involving the G convenient for analyzing results calculated with the DEVSS-G finite element formulation. The only remaining choice is the normalization of the eigenfunctions, which we choose as
ˆ : τˆ p dV = 1. G (38) This choice makes the integral term on the left side of Eq. (37) equal to the eigenvalue, i.e.

1∂ ˆ dE ˆ : τˆ p ) dV = σ. = (G : τˆ p ) dV = σ (G dt 2 ∂t

(39)

The energy balance given by Eq. (37) can be used to determine the critical couplings that cause the onset of the viscoelastic instability by examination of each of the terms on the right side of the equations at subcritical and supercritical values of We. It is important to note that the most important couplings are not the terms with the largest magnitudes, but rather the terms that change most rapidly with We, because these are the terms that will cause the sum in the energy balance to change from negative (stable) to positive (unstable). Joo and Shaqfeh [3] examined the magnitudes of the coupling terms given by Eq. (34) at subcritical and supercritical values of We to determine the most important interactions between the perturbations and the base flow for the instabilities observed in viscoelastic Couette and Dean flows. For axisymmetric Couette flow, the energy analysis showed that the instability is driven by interactions between perturbations to the stress and the base state velocity gradient. That is, the term φτ p G2 is positive and increases as We increases [3]. The term φτ p G1 is relatively small and stabilizing. Their energy analysis also demonstrates that the

38

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

ss coupling between τpθθ and ∂ˆvr /∂θ (described by the term φτ p G1 ) is dissipative and stabilizes the flow and that the dominant term that governs the onset of the transition is still φτ p G2 [2]. In viscoelastic Dean flow, one of the dominant terms in the energy analysis is the term that governs the onset of the Couette instability, φτ p G2 . However, although φτ p G2 is positive, it decreases with increasing We; and a second term, φτ p v2 , governs the onset of the transition, because it is both positive and increases with increasing We. The term φτ p v2 describes the coupling between gradients in the base state stress and perturbations to the velocity, so that the importance of this term in the energy analysis provides strong evidence supporting the mechanism in the viscoelastic Dean instability. The energy method is applied to the flow of an Oldroyd-B fluid for the two different cylinder spacings shown in Fig. 18. The dominant balance is between the terms φτ p G1 and φτ p G2 for both cylinder spacings. For L = 3.5, the term φτ p v1 is also important. The term φτ p G2 is destabilizing, because it has a relatively large magnitude and increases with increasing We. This term is also the dominant term in the energy analysis of the Couette flow. The term φτ p G1 , which couples the base state polymer stress with perturbations to the velocity gradient, stabilizes the flow. Similar results were obtained for other values of L (2.25, 3.0, and 4.0). Although the dominant term in the energy analysis indicates that the mechanism for the transition found in flow past the linear array of cylinders is similar to the mechanism for viscoelastic Couette flow, the transition in the linear array geometry is to a three-dimensional steady-state flow. The multi-step mechanism that drives the transition to time-periodic states in the Couette flow may drive a stationary

Fig. 18. The values of the terms in the energy balance calculated by using mesh M3 for two different cylinder spacings. The different symbols represent the value of dE/dt (䊐), φτ p v1 (), φτ p v2 (), φτ p G1 (), φτ p G2 () and φvis (䊊) for a subcritical and a supercritical value of We.

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

39

transition in the linear array because the contributions to the energy equation are highly non-localized as shown in Fig. 19 for L = 2.5. The energy production term φτ p G2 and the total energy rate-of-increase term dE/dt are highly localized near the recirculation boundary. Although the term φτ p v1 does not control the onset of the instability for L = 2.5, it becomes more important for larger cylinder spacings; and even for the L = 2.5 spacing, the magnitude of this coupling in the region between the cylinders indicates that convection of the perturbations to τ p by the base state velocity field is important. As L increases, the region of the energy production moves to the centerplane; and it is convected onto the upstream face of the next cylinder as seen in Fig. 20. The multi-step, time-dependent mechanism for the instability in circular Couette flow is very similar to the mechanism for the transition in the linear array, except that each step in the mechanism for the transition in the linear array occurs sequentially in space as a fluid particle moves along a streamline instead of sequentially in time. This is most easily demonstrated by extracting streamlines from the finite element results and examining the stresses (τ p11 , τ p12 , τ p22 ) in Protean coordinates in the plane kz z = π. The stress and curvature extracted along a streamline passing through the point (0, 1.1, π/kz ) are shown in Fig. 21. The perturbations to τ p12 and τ p22 are mostly confined to the region between the cylinders where the curvature of the streamline in the two-dimensional base state is positive. The component τˆp11 is non-zero along the entire streamline. The convected interactions between the components of the stress are demonstrated by the fact that the largest perturbations to τ p12 and τ p22 occur upstream of the maximum value of τˆp11 . The mechanism for the flow instability is obtained by writing the linearized constitutive equation for the polymer stress in Protean coordinates and integrating the components of this equation along a streamline. The base state stress and kinematics and the perturbations to the kinematics are introduced as data, resulting in a set of linear ordinary differential equations for the stress components. The importance of each term in this set of ordinary differential equations for (ˆτp11 , τˆp12 , τˆp22 ) is determined by plotting the magnitude of each of the forcing terms as a function of position along the streamline. The two dominant forcing terms in the equation for τˆp22 are shown in Fig. 22. Keeping only these terms gives the following linearized evolution equation for τˆp22 :   τp22 1 ss dˆ ss ˆ G12 + 2ˆτp22 Gss =− (40) v + σ τˆp22 + 2τp21 22 . dq1 We However, the evolution of τˆp22 is not critical to the mechanism for the instability, because the dominant ss ˆ 12 . The dominant forcing terms term in the evolution equation for τˆp12 is the coupling between τp11 and G for τˆp12 are shown in Fig. 23. Keeping only the terms shown in Fig. 23, the evolution equation for τˆp12 is vss

dˆτp12 ss ˆ ˆp22 Gss = −{ˆv · ∇τpss }12 − vss κ(ˆτp11 − τˆp22 ) + τˆp11 Gss 12 + τp11 G12 + τ 21 dq1 ss ˆ ss = −{ˆv · ∇τpss }12 + τp11 G12 + τˆp22 (Gss 12 + G21 ),

(41)

ss where κ ≡ Gss 12 /v . The mechanism for growth of τˆp12 is a combination of the mechanisms for the axisymmetric (ˆτp22 Gss 21 ) ss ˆ and the non-axisymmetric (τp11 G12 ) transitions in viscoelastic Couette flow. One interesting difference is that the time derivative terms in the Couette flow have been replaced by the convected derivative term on the left side of Eq. (41). An additional convection term is given by the first term on the right side of Eq. (41), which accounts for convection of the base state polymer stress by the perturbations to the

40

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

Fig. 19. Contour plots of the integrands in the energy analysis (Eqs. (33) and (34)) for the dominant disturbance at We = 1.56. Calculations were performed for L = 2.5 with mesh M3.

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

41

Fig. 20. Contour plots of the integrands in the energy analysis (Eqs. (33) and (34)) for the dominant disturbance at We = 1.80. Calculations were performed for L = 3.5 with mesh M3.

42

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

Fig. 21. Perturbations to the polymer stress, τˆp11 (䊉), τˆp12 (䉱), and τˆp22 (䊏), along a streamline passing through (0, 1.1, π/kz ) for L = 2.5 (top) and L = 3.5 (bottom). The dotted curve shows the curvature of the streamline, and the vertical dashed lines indicate where the curvature changes signs.

velocity. This term is negative as the fluid moves over the top of the cylinder, so that convection of the base state in this region stabilizes the flow. Because the Protean coordinate system is defined in terms of the base state kinematics, linearization of the convected derivative of the stress with respect to perturbations to the velocity leads to terms in addition to those in Eq. (27). For the 12-component, the Protean coordinate

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

43

Fig. 22. Forcing terms in the differential equation for τˆp22 along the streamline passing through the point (0,1.1, π/kz ). Only the ss ˆ ss ˆ G12 (䊏) and τp22 G22 (䉱) are shown. The vertical dashed lines indicate where the curvature changes sign. leading terms τp21

Fig. 23. Forcing terms in the differential equation for τˆp12 along the streamline passing through the point (0, 1.1, π/kz ). Only ss ˆ ˆp22 Gss the leading terms τˆp22 Gss v · ∇ τˆ ss 12 (䊏), −{ˆ 21 (䊉) are shown. The vertical dashed lines indicate p }12 (䉱), τp11 G12 (䉬), and τ where the curvature changes sign.

44

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

ss ss ss ss Fig. 24. Contributions to the term {ˆv · ∇ τˆ ss ˆ 1 κ(τp11 − τp22 ) and −ˆv2 vss Γ(τp11 − τp22 ). The term {ˆv · ∇ τˆ ss p }12 from v p }12 is shown as (䊏), the term scaled with κ is shown as (- - -), and the term scaled with Γ is shown as (· · · ). The solid curve shows the sum ss ss (ˆv1 κ − vˆ 2 vss Γ)(τp11 − τp22 ).

form is {ˆv · ∇ τˆ ss ˆ1 p }12 = v

ss ∂τp12

∂q1

+ vˆ 2 vss

ss ∂τp12

∂q2

ss ss + (ˆv1 κ − vˆ 2 vss Γ)(τp11 − τp22 ).

(42)

As shown in Fig. 24, the third term on the right side of Eq. (42) accounts for the majority of the contribution to the convection term {ˆv · ∇ τˆ ss p }12 forcing term in Eq. (41). The data were extracted from a base state streamline in the plane kz z = π, so that vˆ 1 is positive. The product of vˆ 1 with the large, negative value of ss τp11 on the surface of the cylinder and the negative value of the curvature κ leads to a positive value of {ˆv · ∇ τˆ ss p }12 , which stabilizes Eq. (41). As in both the viscoelastic Couette flow and Dean flow instabilities, the key step in the mechanism is the evolution of τˆp11 . The important forcing terms in the differential equation for τˆp11 are shown in Fig. 25. For this equation we can neglect fewer terms, so that the “simplified” evolution equation is   dˆτp11 1 ss ˆ =− + σ τˆp11 + 2vss κˆτp12 + 2τp11 vss G11 + 2ˆτp11 Gss τp12 Gss (43) 11 + 2ˆ 21 . dq1 We The term 2ˆτp12 Gss 21 in Eq. (43) is found in both the axisymmetric and non-axisymmetric viscoelastic Couette mechanisms, but the terms involving G11 are new. The “stretching” or “streamwise” terms ss ˆ 2τp11 G11 and 2ˆτp11 Gss 11 are destabilizing in the region of the flow where the curvature is negative, i.e. where the fluid is moving over the top of the cylinder. The Couette-like term 2ˆτp12 Gss 21 is destabilizing only in the region between the cylinders and acts in concert with the perturbations to τ p12 and τ p22 . It should be noted that the Couette-like term 2ˆτp12 Gss 21 becomes unimportant when the intercylinder spacing is increased, whereas the new coupling between the elongational component of the velocity gradient in

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

45

Fig. 25. Forcing terms in the differential equation for τˆp11 along a streamline passing through the point (0, 1.1, π/kz ). Only the ss ˆ leading terms 2vss κˆτp12 (䊏), 2τp11 G11 (䉱), 2ˆτp11 Gss τp12 Gss 11 (䉬), and 2ˆ 21 (䊉) are shown. The vertical dashed lines indicate where the curvature of the streamline changes sign.

the base state and the perturbations to the normal stress is substantially enhanced. Therefore, this coupling represents a new mechanism that is not found in the shear-dominated circular Couette and Dean flows. Eqs. (40), (41) and (43) describe a mechanism that is similar to the mechanism for the transitions observed in viscoelastic Couette flow. The magnitudes of the terms in Eq. (41) indicate that the mechanism is most similar to the non-axisymmetric Couette mechanism, due to the key role of the perturbation to

46

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

Fig. 26. Comparison of the perturbations to the stress obtained by integration of the simplified equations in the text (solid curves) with the finite element data, which were calculated with the full set of equations. Results are shown for τˆp22 (䊏), τˆp12 (䉱), and τˆp11 (䊉).

G12 in generating large perturbations to τ p12 , which in turn excite large perturbations to τ p11 . Integration of Eqs. (40), (41) and (43) using the kinematics of the base flow as data shows good agreement with the values of the perturbations to the stress calculated from the linear stability analysis. Sample results are shown in Fig. 26. The importance of the term 2ˆτp12 Gss 21 is further emphasized by the fact that small errors in τˆp12 during the integration of Eq. (41) can lead to large errors in the integration of Eq. (43) for τˆp11 . For the results shown in Fig. 26, the value of τˆp12 in the term 2ˆτp12 Gss 21 is taken from the finite element results as data, whereas all other values of τˆp are obtained by integration of the simplified equations for the stress perturbations. The last step in the instability mechanism is the interaction between the streamwise normal stress τˆp12 and the momentum equation. The 11-component of the polymer stress enters the 2-component of the momentum equation through the curvature term κτ p11 . After linearization, this interaction has the form ss κˆτp11 + κˆ τp11 , where the perturbation to the curvature of the streamlines is given by κˆ =

ˆ 12 (ˆv1 + vˆ 2 )Gss G 12 − . ss ss 2 v (v )

(44)

Eq. (44) holds only for values of kz z where vˆ z = 0, such as kz z = π. The linearized contributions from ss κτp11 by the perturbations are shown in Fig. 27. The sum of these two terms, κˆτp11 + κˆ τp11 , is shown in Fig. 28 along with the perturbation to the normal component of the velocity vˆ 2 . The gradients in vˆ 2 ss along the streamline are largest where κˆτp11 + κˆ τp11 is large. This fact demonstrates that the perturbations ˆ 12 through the 2-component of the to the curvature and the normal stress reinforce perturbations to G ss ˆ ˆ momentum equation. The perturbation to G12 provides feedback into the dominant term, τp11 G12 , in the evolution equation for τˆp12 , Eq. (41).

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

47

Fig. 27. Contributions to the 2-component of the momentum equation through the linearized form of the term κτp11 . The term ss κˆτp11 is shown by (䊏), and the term κˆ τp11 is shown by (䉱). The vertical dashed lines indicate where the curvature in the base state changes signs.

ss Fig. 28. The interaction between the linearized curvature term κˆτp11 + κˆ τp11 , shown as (䊏), and the normal component of the perturbation to the velocity, vˆ 2 , shown as (䉱). The vertical dashed lines indicate where the curvature in the base state changes signs.

48

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

7. Conclusions A transition in the flow of an Oldroyd-B fluid through a linear, periodic array of cylinders from a two-dimensional steady-state to a three-dimensional steady-state has been computed by using a linear stability analysis based on the DEVSS-G/SUPG finite element discretization. The calculation of a transition in a complex, two-dimensional geometry demonstrates the applicability of the suite of numerical methods presented in [9] for analysis of important flow problems. The stability calculations for the Oldroyd-B model reasonably predict the experimentally observed flow transition of a polyisobutylene Boger fluid when the cylinders are closely spaced. The critical value of the Weissenberg number and the spatial form of the calculated three-dimensional state are in good agreement with the experimental observations of Liu [1], who detected the onset of a time-periodic instability at Wecrit (γ) ˙ = 1.53 with a spatial wave number in the neutral direction of about kz = 3. The calculations predict Wecrit = 1.55 and kz = 3.1 for these values. The spatial structure of the instability converges as the mesh is refined; and the dominant eigenvalue is shown to be a discrete mode, clearly separated from the continuous eigenspectrum. The experiments of Liu show a flow transition to a time-periodic state with an extremely long time period (T ∼ = 1950 ± 730). The calculations with the Oldroyd-B fluid model seem to predict the transition to a steady-state, three-dimensional flow structure, although the time integrations were not performed for long enough times to see such a large period oscillation. The analysis with the Oldroyd-B fluid fails, however, to predict the correct dependence of the flow instability upon cylinder spacing L, much as it fails to predict the critical value of We for the instability about a single cylinder [9]. We believe this shortcoming is caused by a change in physical mechanism for the instability from the shear-like mode for closely spaced cylinders to a mode dominated by the transient extensional rheology in the wake of the cylinder as seen by McKinley et al. [7] and Byars [25]. The Oldroyd-B model does not capture this rheological behavior adequately. Another study [26] shows that the closed form of the adaptive length scale constitutive equation (ALS-C), which was developed to capture more accurately transient rheological behavior [27], qualitatively predicts the instability for an isolated cylinder and the dependence of Wecrit upon L for the periodic array. Energy analysis of the calculated transition demonstrates that the dominant contributions to the energy balance involve couplings between the velocity gradient and the polymer stress. These couplings also are found in the analysis of the transition observed in viscoelastic Couette flow [2]. However, the transition found for the linear array is steady, whereas the transition found for the Couette flow is time-periodic. This apparent inconsistency may be resolved by the observations that for the linear array the contributions to the energy balance are highly localized, and that the evolution of the perturbations to the polymer stress, although steady in an Eulerian frame, is unsteady in a Lagrangian frame of reference. A Protean coordinate system is used to analyze the finite element stability results in the context of the curved streamline instability mechanisms. The analysis shows that the coupling between the fluctuation in ss the G12 component of the velocity gradient of the perturbation and the base state normal stress τp11 leads to growth of τˆp12 . These perturbations to τp12 excite growth of τˆp11 through a coupling with the base state velocity gradient Gss 21 . These two steps are precisely the mechanism for generating large perturbations to τ p11 in the transition observed in viscoelastic, non-axisymmetric Couette flow [2]. However, this sequence is localized to the regions between the cylinders. Additional interactions not found in viscoelastic Couette flow also are observed in the flow where ˆ 11 the curvature changes signs and the fluid moves over the top of the cylinder. Here new terms τ ss G p11

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

49

and τˆp11 Gss 11 involving elongational components of the kinematics are significant, and these interactions dominate the evolution of the perturbations to τ p11 for larger cylinder spacings where the velocity in the wake of the cylinder can approach the maximum velocity in the open channel. These terms should be especially interesting in the stability analysis of the viscoelastic flow around an isolated cylinder in a channel where the elongational flow in the wake is very strong. However, as shown in [9], the large polymer stress that results from convection of stress into the wake and the elongational kinematics along the centerplane have made calculation of base states for the isolated cylinder very difficult for We ≥ 0.75. Acknowledgements This work was supported in part by the ERC Program of the National Science Foundation under Award Number EEC-9731680. References [1] A.W. Liu, Viscoelastic Flow of Polymer Solutions around Arrays of Cylinders: Comparison of Experiment and Theory, Ph.D. dissertation, Massachusetts Institute of Technology, Cambridge, MA, 1997. [2] Y.L. Joo, E.S.G. Shaqfeh, Observations of purely elastic instabilities in the Taylor–Dean flow of a Boger fluid, J. Fluid Mech. 262 (1994) 27–73. [3] Y.L. Joo, E.S.G. Shaqfeh, A purely elastic instability in Dean and Taylor–Dean flow, Phys. Fluids 4 (1992) 524–543. [4] R.G. Larson, E.S.G. Shaqfeh, S.J. Muller, A purely elastic instability in Taylor–Couette flow, J. Fluid Mech. 218 (1990) 573–600. [5] M. Avgousti, A.N. Beris, Non-axisymmetric modes in viscoelastic Taylor–Couette flow, J. Non-Newtonian Fluid Mech. 50 (1993) 225–251. [6] P. Pakdel, G.H. McKinley, Elastic instability and curved streamlines, Phys. Rev. Lett. 77 (1996) 2459–2462. [7] G.H. McKinley, R.C. Armstrong, R.A. Brown, The wake instability in viscoelastic flow past confined circular cylinders, Phil. Trans. R. Soc. Lond. A 344 (1993) 265–304. [8] R. Sureshkumar, M.D. Smith, R.C. Armstrong, R.A. Brown, Linear stability and dynamics of viscoelastic flows using time-dependent simulations, J. Non-Newtonian Fluid Mech. 82 (1999) 57–104. [9] M.D. Smith, R.C. Armstrong, R.A. Brown, R. Sureshkumar, Finite element analysis of two-dimensional viscoelastic flows to three-dimensional perturbations, J. Non-Newtonian Fluid Mech. 93 (2000) 203–244. [10] R. Glowinski, O. Pironneau, Finite element methods for Navier–Stokes equations, Ann. Rev. Fluid Mech. 24 (1989) 167– 204. [11] P. Saramito, A new ϑ-scheme algorithm and incompressible FEM for viscoelastic fluid flows, Math. Model. Numer. Anal. 28 (1994) 1–35. [12] Harwell Subroutine Library Release 11, available from AEA Technology, Building 424.4, Harwell, Oxfordshire OX11 0RA, UK. [13] A.E. Caola, Y.L. Joo, R.C. Armstrong, R.A. Brown, Highly parallel time integration of viscoelastic flows, J. Non-Newtonian Fluid Mech. 100 (2001) 191–216. [14] B. Hendrickson, R. Leland, The CHACO User’s Guide Version 2.0, 1995 [15] Y. Saad, Numerical Methods for Large Eigenvalue Problems, Manchester University Press, Manchester, UK, 1992. [16] L.N. Trefethen, D. Bau, Numerical Linear Algebra, SIAM, Philadelphia, 1997. [17] R.B. Lehoucq, D.C. Sorensen, C. Yang, ARPACK Users Guide: Solution of Large Scale Eigenvalue Problems by Implicitly Restarted Arnoldi Methods, ftp://ftp.caam.rice.edu/pub/software/ARPACK. [18] V.V. Ramanan, K.A. Kumar, M.D. Graham, Stability of viscoelastic shear flows subjected to steady or oscillatory transverse flow, J. Fluid Mech. 379 (1999) 255–277. [19] R. Sureshkumar, A.N. Beris, M. Avgousti, Non-axisymmetric subcritical bifurcations in viscoelastic Taylor–Couette flow, Proc. R. Soc. Lond. A 447 (1994) 135–153.

50

M.D. Smith et al. / J. Non-Newtonian Fluid Mech. 109 (2003) 13–50

[20] K. Adachi, Calculation of strain histories in Protean coordinate systems, Rheol. Acta 22 (1983) 326–335. [21] H.J. Wilson, M. Renardy, Y. Renardy, Structure of the spectrum in zero Reynolds shear flow of the UCM and Oldroyd-B liquids, J. Non-Newtonian Fluid Mech. 80 (1999) 251–268. [22] M.D. Graham, Effect of axial flow on viscoelastic Taylor–Couette instability, J. Fluid Mech. 360 (1998) 341–374. [23] R.A. Brown, M.J. Szady, P.J. Northey, R.C. Armstrong, On the numerical stability of mixed finite-element methods for viscoelastic flows governed by differential constitutive equations, Theor. Comput. Fluid Dyn. 5 (1993) 77–106. [24] P.G. Drazin, W.H. Reid, Hydrodynamic Stability, Cambridge University Press, Cambridge, UK, 1981. [25] J.A. Byars, Experimental Characterization of Viscoelastic Flow Instabilities, Ph.D. Dissertation, Massachusetts Institute of Technology, Cambridge, MA, 1996. [26] Y.L. Joo, I. Ghosh, R.C. Armstrong, R.A. Brown, Use of a New Adaptive Length Scale Model in the Simulation of Flow of a Viscoelastic Fluid around a Linear Array of Cylinders in a Channel, in preparation. [27] I. Ghosh, Y.L. Joo, G.H. McKinley, R.A. Brown, R.C. Armstrong, A new model for dilute polymer solutions in flows with strong extensional components, J. Rheol. 46 (2002) 1057–1090.