A numerical simulation of the appearance of chaos in finite-length Taylor-Couette flow

Applied Numerical North-Holland

Mathematics

41

7 (1991) 41-71

A numerical simulation of the appearance of chaos in finite-length Taylor-Couette flow * C.L.

Streett

NASA Lungley Reseurch Center, Hompton, VA 23665, USA

M.Y.

Hussaini

Institute for Computer Applicutrons 23665. USA Revised

December

in Science und Engineering,

NASA

Lunglq

Reseurch Center, Hampton,

VA

1988

Abstruct Streett, C.L. and M.Y. Hussaini, Couette flow, Applied Numerical

A numerical Mathematics

simulation of the appearance 7 (1991) 41-71.

of chaos

in finite-length

Taylor

Taylor-Couette flow, the shear-driven flow between concentric cylinders, exhibits a wide variety of instabilities and modal changes, especially for the case of finite length to gap ratio. The numerical simulations presented here capture many of the experimentally observed features. including the moderately high Reynolds number regime in which temporally aperiodic behavior is seen. The exponential decay of the temporal frequency spectrum of these modes in the simulations indicate such flows possess a low-order chaotic character. In this paper, the spectral collocation methods used in this study are described, select axisymmetric simulations are reviewed, and initial results from three-dimensional unsteady simulations are presented.

1. Introduction The Taylor experiment on Couette flow between coaxial circular cylinders has been the subject of numerous theoretical and experimental studies [15]. This flow is rich in complex phenomena; so rich, in fact, that they are still being discovered [30], and our understanding of them is far from complete. In a typical execution of the experiment, the inner cylinder rotates with a constant angular velocity, while the outer cylinder, along with the top and bottom walls, is kept at rest. The relevant geometric parameters are the radius ratio of the inner and outer cylinders, and the aspect ratio, which is the ratio of the length of the annulus to its width. The dynamic parameter is the Reynolds number based on the angular velocity of the inner cylinder and the annulus width or gap. The Taylor-Couette flow is strongly dependent on all of these parameters. The theory for the infinite aspect ratio case (which neglects end-wall effects) and its correspondence to the experiments in cylinders of necessarily finite but large aspect ratio are * Presented

at SAE Aerotech

0168-9274/91/$03.50

‘88, Anaheim,

CA, October

0 1991 - Elsevier Science

Publishers

3-7, 1988. B.V. (North-Holland)

42

C. L. Street& M. Y. Hussaini

/ Chaos in Taylor-Couette

flow

reviewed by Di Prima and Swinney [15] and Di Prima [14]. Benjamin [5] has developed a rigorous qualitative theory for the existence of nonunique steady states for confined flows and their stability with particular reference to finite-length Taylor-Couette flow. His predictions have been qualitatively confirmed in a series of experiments [6] in cylinders of small aspect ratio with fixed end walls. They have been further confirmed by the numerical results of Cliffe and Mullin [13] who discretized the steady Navier-Stokes equations by a Galerkin finite-element method and solved the resulting nonlinear algebraic equations by the pseudo-arclength continuation method of Keller [21]. Among other numerical studies of the finite-length Taylor-Couette flow problem, those of Alziary de Roquefort and Grillaud [2] and Neitzel [28] are worth mentioning. Both investigations were based on the time-dependent vorticity/ stream-function formulation along with an equation for the azimuthal velocity. They used a finite-difference method to solve for the steady state by a time asymptotic approach. Their numerical results show the axial structure for small Reynolds number, the smooth development of the flow with rapid increase in vortex activity for Reynolds numbers in the quasicritical range, and multiple states for sufficiently large Reynolds number. The aspect ratio being large, their results are only pertinent to the exchange phenomena beyond the two-cell/ four-cell interactions examined by Benjamin. Experimental studies of the so-called “wavy” vortex flows which appear at higher Reynolds numbers show how richly varied the phenomena observable in this experiment are. Two such studies are those of Bust, Dornblaster and Koschmieder [ll], and Gorman and Swinney [18], in which the structures of the early periodic and quasiperiodic wavy vortex flows are examined in detail for high aspect ratios. For smaller aspect ratios pertinent to this work, the investigations of Ross and Hussain [31], Mullin and Benjamin [27], and Lorenzen, Pfister and Mullin [23] show the effect of finite aspect ratio on the critical Reynolds numbers for onset of periodic and quasiperiodic wavy modes. It is now well-known from experiments that Taylor-Couette flow at higher Reynolds numbers exhibits aperiodic behavior which is nonetheless deterministic and low-dimensional in the dynamical systems theory sense; that is, the flow is chaotic. This is to be distinguished from a fully turbulent flow, which is stochastic and of high enough dimension that the concept of an attractor is meaningless. Studies of the onset and nature of chaotic Taylor-Couette flow include those of Brandstater and Swinney [8], Fenstermacher, Swinney and Golub [16], Andereck, Liu and Swinney [3], and Pfister [30]. The basic tool for analyzing these flows is temporal frequency spectra from histories of a measured velocity; the characteristics of the broad-band noise in the spectrum can be used to distinguish between chaotic and turbulent behavior [8], as well as determining the nature of discrete oscillatory components of the flow. Phase-space portraits, Poincare sections, and calculations of Lyapunov exponents from time histories are used to analyze the dynamical behavior of the flow, in an attempt to relate the aperiodic transition of the system to theoretically developed scenarios of routes to chaos. Numerical investigations of wavy-vortex flow states were carried out by Jones [20], who computed the stability of the axisymmetric Taylor-Couette flow to non-axisymmetric disturbances for finite gap width; and by Walgraef, Borckmans and Dewel [37], who included finite-length effects in their stability calculations. Full numerical simulations of periodic wavyvortex flows were carried out by Moser, Moin and Leonard [26], and by Marcus [22]. Although the discretizations in the former were too coarse to produce meaningful results, the computations of the latter were carried out carefully and were in outstanding agreement with experiments for

C.L. Streett, M. Y. Hussaini / Chaos in Taylor-Couette

flow

43

the critical Reynolds number and wave propagation speed of the first wavy mode in the infinite-length case. The purpose of our continuing research effort is to solve the unsteady Navier-Stokes equations by highly accurate spectral collocation methods with a view to elucidate the underlying processes leading to laminar-turbulent transition in the Taylor-Couette flow. In a method development and verification phase of this work, we have performed a number of studies of axisymmetric state bifurcations which take place for aspect ratios in the range 1 < r < 6. These simulations are discussed in detail in [34,36]. The numerical schemes presented in [34] were extended for the computation of time-accurate solutions of the three-dimensional Navier-Stokes equations, to study wavy and chaotic Taylor-Couette flow in the finite aspect ratio case. Some preliminary results of this study were presented in [35]; analysis of the chaotic flow simulations is continuing, and will be presented in the future. 2. Numerical methods In this section, the two algorithms used in these studies are outlined. The first is a time-split scheme, implemented in both axisymmetric and three-dimensional time-accurate forms. The second algorithm employs a version of the staggered-mesh discretization of Bernardi and Maday and time-accurate form for axisymmetric [7]. This scheme was used in both a steady-state calculations. 2.1. Time-split scheme incompressible

The time-dependent conservation form:

Navier-Stokes

equations

in cylindrical

coordinates

are, in

2

--

g+g+

0)

(4

aw

a(uw>

at+-F-+Yae 1 d(i-24)+Lk+aM:=o -r

i3r

r

ae

aZ

(3) (4)

7

where

and U, u and w are the r, 8 and z components of velocity, respectively. The axisymmetric case is recovered by eliminating the partial derivatives with respect to 0. For the case where the outer cylinder and end walls are stationary and the inner cylinder rotates, the boundary conditions are: u=u=w=o on the outer cylinder, r=R 0’

u=u=w=o

at the top and bottom

u=w=o

on the inner cylinder,

wall,

z = *Il. r=R,,

44

C.L. Streett, M. Y. Hussaini

In the present satisfies:

calculations,

the azimuthal

/ Chaos in Taylor-Couette

velocity

flow

is split into two parts,

u = ut, + 6, where u,,

V2Ub = 0, Ub =

0

at r=R,

~Ri

at r=Ri.

andat

z=

kz,,

The quantity ub is computed and stored at the start of the calculation, which proceeds with the computation of u, u” and w at each time step. Note that these three velocity components satisfy homogeneous boundary conditions. A splitting method is employed to advance the solution from time t” to t”+l. Writing the Navier-Stokes equations in vector notation, with IL representing the velocity (u, C, w), we have: u,iu~vu=

-VPfvV211

(5a)

and V*u=O in the annulus,

(5b)

and:

u = 0

on the boundary

The split scheme first advances

F. un to an intermediate

solution

u* by solving:

u: + u* ‘vu*=u~2u*, U

*_

(6)

on F.

-g*

The intermediate solution is advanced

boundary condition u * = g * will be discussed from u* to un+’ via:

subsequently.

Finally,

the

n+l _ - -vpn+l, ut v

. un+l

A*ll

= 0,

‘+I=0

(7)

oriF,

where i? is the unit normal to the boundary F. Note that the final “pressure-correction” step by itself is a set of inviscid equations. It is well-posed when boundary conditions on the normal component of the velocity only are enforced. At the end of the full step there exists a nonzero tangential component of velocity on the boundary. The magnitude of this slip velocity can be reduced by a proper choice of the intermediate boundary condition on u *. Marcus [25] has described the difficulties which arise from the use of u* = 0 as the intermediate boundary condition. The conditions used here are based on the work of Fortin, Peyret and Temam [17]. Using backward Euler time discretization for (7) yields: U

n+l

=

IL* -

and the slip velocity 4.U

At

vP”+l,

on the boundary

“+‘(r=7^.(~*Ir-Aat

(8) is given by: vPn+‘jI.).

If 7^.u*= 0 on the boundary, then 7^. un+’ 1,‘ = 0( At) w h ere 7^ is the unit tangent boundary F. However, if VP n+l is expanded in a Taylor series about t = t”, i.e., VP”+’ = VP” + At VP*” + 0( At’),

(9) to the

C. L. Streett, M. I: Hussaini / Chaos in Tqlor-Couette

and the second

term is approximated

At VP,” = (VP” - VP’-‘)

flow

45

by + O(At’),

then (8) becomes $.“n+’

(,.=F(u*(,.-At(2vP”-VP”-‘)]

Thus, the slip velocity ?.u*

may be reduced

+O(At3).

to 0( At’) through

the intermediate

boundary

condition:

(,=FAt(2vP”-VP”-‘).

00)

Of course, the boundary condition ii . u = 0 is retained. The pressure step is actually carried out in two parts.

v2P n+l

First, the divergence

of (8) yields:

-&p.ld*, _

where v . u”+’ = 0 is enforced. Then, the velocities are updated using (8). Note that this formulation requires a boundary condition for the pressure. This poses a problem, since there is no natural boundary condition for pressure. The use of a condition derived from enforcing the normal momentum equation at the boundary was attempted, but was apparently inconsistent as it resulted in explosive instability. Fortunately, analysis of the split scheme itself yields a self-consistent pressure condition: ~~vP”+‘~~=o,

(12)

are both zero on the boundary. The error involved in this since both A.u* and ii.u*+’ specification is, we believe, related to the overall splitting error of the scheme. It is also known that the error due to imposition of an inconsistent Neumann pressure boundary condition is isolated to a thin “boundary layer” [29]. Zang and Hussaini [38] have extensively investigated a related split scheme in which two coordinate directions were periodic and the third employed general boundary conditions. Comparison between split and unsplit codes using the same discretization yielded agreement to five decimal places. However, they utilized a staggered grid for pressure in the nonperiodic direction, obviating the need for a pressure boundary condition. The unstaggered scheme used here requires a pressure boundary condition. Actually. the consistent pressure equation derived by Zang and Hussaini yields exactly the same condition on pressure as used here. No unusual instabilities were ever encountered with the present scheme. 2.1. I. Discretization and solution scheme Since no-slip boundary conditions are enforced in both the r and z directions. Chebyshev spectral representation of the Gauss-Lobatto points is appropriate in both directions. Azimuthal periodicity naturally allows the use of Fourier series in that direction. Collocation is used for a number of reasons: straightforward treatment of nonlinear terms and boundary conditions, capability to include coordinate stretchings, and ability to solve the resulting discrete equations rapidly. For further discussion of this form of discretization. see [12.32]. A coordinate stretching is employed in the radial direction to resolve the large gradients near the inner cylinder. The form of this stretching is: r=

[l+bex~(-a)l(K.-Ri)~ [l +b

exp(-ar,)]

+R

’

I’

(13)

46

CL. Streett, M. Y. Hussaini

/ Chaos in Taylor-Couette

flow

where yc is the radial coordinate in the computational space. Values of a = 2 and h = 2,. . . ,20 were typical in this work. The z direction was not stretched. Time discretization of the first step in the split scheme involves a low-storage mixed Runge-Kutta/Crank-Nicholson scheme [12]. Writing the semidiscrete equation for the first step symbolically for all velocity components as: u,=A(u)

+I+),

04)

where A(u) and D(u) represent advection and diffusion terms, respectively, the mixed scheme advances from time step t” to t * using a third-order Runge-Kutta scheme for the advection terms and Crank-Nicholson for the diffusion terms:

u” =

u(P),

H’=AtA(u’),

u1= u” + ;H’ + ; H*

=

At

=

At

gH2+

& At(D(u')

+D(u2)),

(15)

A( u’) - sH2,

u3= u2+ u( tn+‘)

D( u”) + D( u’)),

A( u’) - $H',

u2= u1+ H3

At(

=

AH3

+ + At( D( u’) + D( u’)),

u3.

The second step of the split scheme uses backward Euler time discretization, as mentioned in the previous section. The above scheme is stable to O(1) Courant numbers, which involves time steps many times larger than those desired for accuracy. Typically, a time step which corresponded to Courant numbers of 0.01 to 0.1 based on the smallest physical mesh spacing was used. The slip velocity resulting from the choice of time step was normally 8 to 10 orders of magnitude smaller than the maximum velocity in a given direction. Note that the above scheme reduces the problem to a sequence of uncoupled Helmholtz-Poisson equations to advance the discrete solution. Since one time step requires the solution of nine positive-definite Helmholtz equations with Dirichlet boundary conditions and one Poisson equation with pure Neumann boundary conditions (per Fourier azimuthal mode), a computationally efficient technique to solve such equations had to be employed if this study was to be feasible. A straightforward tensor-product decomposition technique is used here, in which the eigenvectors of the one-dimensional spectral operators are used to diagonalize the multidimensional operator [12]. Solution is carried out in the Fourier-transform space, since the Poisson operators for each Fourier wave number decouple. The eigenvectors are computed in a preprocessing step, and stored for use each time step; the resulting scheme requires less than 0.18 seconds on the NAS Cray 2 to solve a Helmholtz-Poisson equation with Dirichlet boundary conditions on a mesh with a Chebyshev-Chebyshev-Fourier discretization of 65 X 45 X 32 points, including arbitrary physical-space stretchings in both Chebyshev directions. The pure Neumann problem of the pressure discretization is handled through the influence matrix method; see [32] for details. For moderate to large aspect ratios (F - 3-20), a large number of points are required in the axial direction relative to the radial direction. Since the Courant number for a Chebyshev mesh is

C. L. Streett, M. Y. Hussaini / Chaos in Taj*lor-Couette

flow

47

proportional to 1/N2, this can produce a severe time step limitation for such cases. To reduce’ the impact of high resolution in the axial direction, a spectral multi-domain technique is employed, in which the axial extent is broken into two or more subdomains, each with independent collocations, and patched together along an interface line. The global flux-balance technique of Macaraeg and Streett [24] is used to ensure spectral accuracy of the overall discretization. 2.2. Staggered-mesh

scheme

2.2. I. Discretiratian In the incompressible Navier-Stokes equations, the pressure loses its identity as a thermodynamic quantity and becomes a kinematic variable whose purpose is to allow a solenoidal velocity field to satisfy the continuity equations. Since there are no physical boundary conditions on the pressure and the continuity equation does not require any additional boundary conditions on the velocity over those employed in the momentum equations, it is natural to relate the pressure in the interior of the domain with satisfaction of the continuity equation in the interior. A finite-element technique for implementing this idea which has appeared widely in the literature involves bilinear velocity elements (locating velocities at the corners of an element) and constant yields a correct number of unknowns or number of pressure elements. This arrangement equations and boundary conditions balance, but in theory and in practice is subject to a “checkerboard” oscillation of the pressure in the solution. This is known as a “spurious mode”; it does not affect the gradient of the pressure in the approximation, only the pressure itself. It has been argued that his oscillation does not affect the accuracy of the solution and may be filtered out. While this may be true in the finite-element of finite-difference context, filtering of a spectral solution to eliminate a spurious mode will almost surely destroy exponential-order error convergence. A finite-difference staggered-mesh discretization without spurious modes is used in the marker-and-cell method [19]. In this technique, the pressure is located at the center of the cell;

j+l

“I-+,j+1 pi,j+l V.

i-1 .j+k

vi A’4 +

j

pi-l,j

"i-4.j

pfj

ul+Js,j

pl+l,j

‘i ,j-)j j-l

p1+1 .j-1

i-l

Fig. 1. Marker-and-cell

i

1'1

staggered

mesh.

48

C.L. Street& M. Y. Nussuini / Chaos in Tuylor-Couette

flow

the velocity components are located in the center of the cell face which is normal to the direction of that velocity component (see Fig. 1). Although free of spurious modes, the no-slip boundary condition cannot be enforced directly in this method; it is set approximately using “ghost point” variables. The staggered-mesh spectral discretization of Bernardi and Maday [7] was developed theoretically to eliminate the appearance of spurious modes for the Stokes equations. The collocation implementation of this discretization yields a grid structure reminiscent of the marker-and-cell method. Before explaining the Bernardi-Maday discretization, a few basic concepts of spectral methods are required. The computational-space grid in a spectral collocation method may not be chosen arbitrarily, but is a consequence of the choice of basis functions and the quadrature formula used to derive a discrete approximation series from the analytic series. See [12] for further details. For spectral solution of differential equations, the most useful set of collocation points is that based on Gauss-Lobatto quadrature. This form yields collocation points at the extrema of the last retained basis function in the truncated series. Another useful collocation point set comes from Gauss quadrature, giving collocation points at the zeros of the first truncated basis function. The Gauss-Lobatto points (GL) include the endpoints of the domain, whereas the Gauss points (G) do not; the Gauss points from a series with N - 1 basis functions retained (GNPI) interleave with the Gauss-Lobatto point from a series with N basis functions (GL,). It is also possible to define basis functions for a spectral approximation series which yields collocation points which fall on the Gauss points, but which also include the endpoints of the domain. We refer to this set as G+. The G+ spectral series has theoretically similar accuracy as the GL series, except for nonsmooth functions; such a situation is, however, not an issue here. In the Bernardi-Maday staggered discretization in two dimensions, the pressure is located on the collocation points from a (G,_, x G,_,) representation; the u and u velocities (in the x and y directions, respectively) are on (GL, x GL,,) and (Gi+i x GL,) meshes, respectively. This yields a staggered mesh as shown in Fig. 2, with the pressure collocated only in the interior and both velocities have collocation points in their sets which fall on all at “cell centers”; boundaries. The continuity equation is enforced at the (G,_, X G,,,_,) points (i.e., the pressure collocation points); and the x and y momentum equations are enforced on the u and u

__

u, x-momentum

(GL x G’)

-----

v, y-momentum

(G’x

---

P, continuity

Fig. 2. Spectral staggered mesh.

(G x G)

GL)

C. L. Streett, M. Y. Hussainr / Chaos in Taylor-Couette

49

Jon

collocation point sets, respectively. The interpolations required by these locations of variables and equations enforcement are carried out through the use of the appropriate basis function series. For example, the two-dimensional steady Stokes equations may be written: x-momentum ZJ( D&_, D;+)

. (1, 1)~ - ( DcL, 0). ( IzL,

I)P

=f;;

y-momentum v

( D;+, D&)

. (1, 1)~ - (0,

D,,) . (1,

IzL)p

=f,;

continuity

differentiation where DG,_ and DG+ are the collocation boundary-extended Gauss points, respectively; and I:: from grid “Gl” to grid “G2”. 2.2.2. Solution scheme for steady state An iterative technique is used to solve the discretization [7]. The equations to be solved are: V2U-

operators on the Gauss-Lobatto and is the spectral interpolation operator

steady-state

vP+f,

Stokes

operator

based

on

this

(I6a)

V*U=O,

(16b)

where f is the momentum equation source term. In the present algorithm, the nonlinear convection terms are satisfied by a Newton iteration, with each step involving the solution of the above Stokes operator; thus, f contains the convection terms evaluated from the velocities at the previous Newton iteration. Spectral interpolation is used to evaluate the factors in these No difficulties were ever noted in this nonlinear terms at the appropriate mesh locations. interpolation procedure. Rewrite (16) symbolically as: LU-GP=f,

(174

DU=

(17b)

0.

Solving (17a) for U and substituting (DL?G)P=

-(DL-‘)f.

into (17b) yields: 08)

The solution of (18) for P is carried out by observing that D and G are first-order operators (divergence and gradient, respectively), and L is second-order (Laplacian). The operator DL- ‘G, thus, is essentially zeroth-order, aside from the effect of boundary conditions which are enforced in the operator L. Conjugate gradient iteration is therefore used for solution of (18), since it yields rapid convergence for problems with limited eigenvalue spread and requires only the application of the operator itself. This algorithm is generally known as the Uzawa algorithm [lo]. However, the eigenvalues of the pressure Stokes operator with spectral collocation discretization are apparently grouped far better than for the algorithm’s original finite-element application, resulting in much faster and more reliable convergence [4]. The typical spectral radius observed

C. L. Streett, M. Y. Hussaini / Chaos in Taylor-Couette froti

50

for the conjugate gradient iterations in this and other applications of the algorithm is between 0.05 and 0.15, with only an extremely weak degradation with mesh refinement up to very fine levels. Unfortunately, this advantageous eigenvalue grouping does not hold for the equivalent unsteady algorithm. Symbolically, we have: (L-X1)u-GP=f, Du=O,

which yields D(L-hl)-‘GP=

-DL?f,

where X - l/At from the time discretization. For the case of moderate to small time steps (At -=Kl), the pressure operator is essentially second-order. (It is interesting to relate this operator with the pressure step of the split scheme described in the previous section.) Its eigenvalues are thus poorly grouped and have a rapid (0( N4)) dependency on mesh refinement; the conjugate gradient solution algorithm typically fails to converge for grids finer than about 5 x 5. A multi-grid technique was developed to alleviate this difficulty [33]. Although the pressure operator is linear, it is expensive to generate explicitly. Therefore, the full-approximation scheme (FAS) multi-grid form of Brandt is used [9], and is reviewed here. The discrete problem we wish to solve is written:

A’Q’ = Ff. After a few iterations of the relaxation solver (conjugate gradient in this case), an approximation qf has been computed for which the high-frequency content of the error q f - Q f is small. For the coarse grid, on which the remaining error appears as high frequency, the problem to be solved is: A”Q” = F”,

where F” = R [ Ff - AfqF] + AcRqf.

The restriction operator R interpolates a function from the fine grid to the coarse grid. After an solution is obtained by application of the adequate approximation qc to the coarse-grid relaxation solver (or by direct solution, if the coarse grid is coarse enough), the fine-grid approximation is corrected via: qf=qf+P[qc-Rqf].

The prolongation operator P interpolates a function from the coarse grid to the fine grid. With the conjugate gradient scheme used as the relaxation solver, this multi-grid algorithm requires only applications of the pressure Stokes operator. Two relaxation iterations are typically performed at each level of the fixed V-cycle schedule; relaxations are performed only during the grid-coarsening phase of the cycle. Natural spectral interpolation operators are used for restriction and prolongation, with factor-of-two coarsening used. On the coarsest mesh, typically 17 x 17 points in a two-dimensional problem, a direct solution of the pressure Stokes operator is used, the operator having been generated, inverted and stored in a preprocessing step.

C. L. Streett, M. Y. Hussaini

/ Chaos in Taylor-Couette

flow

51

The multi-grid/conjugate-gradient algorithm is quite robust and relatively insensitive to choice of grid cycle scheduling. For solution of the steady problem, the multi-grid algorithm requires about the same amount of machine time as the single-grid conjugate-gradient algorithm; the improvement in convergence rate just offsets the overhead of the algorithm. For the unsteady problem, a small degradation occurs in the convergence rate with decreasing At relative to the steady case, but the machine time to reach a fixed residual level is roughly constant over a wide range of time steps. This is because the initial residual levels are lower for the case of smaller At; note that the residual level is a meaningful measure of the error in the continuity equation. The staggered-mesh multi-grid/conjugate-gradient coupled time-stepping scheme is considerably more expensive per time step than the nonstaggered time-split algorithm. On equivalent grids, the coupled scheme requires about 50-100 times the CPU time as the split scheme per time step; this is because a time step in the split scheme requires only the solution of 11 scalar uncoupled Helmholtz-Poisson equations per azimuthal plane, which can be carried out very rapidly using the tensor product diagonalization method. However, the coupled scheme is somewhat more accurate than the split scheme, due to the splitting errors and imposed pressure boundary conditions of the latter. To reduce the magnitude of these errors in the split scheme, the time step chosen must be between 10 and 100 times smaller than that required purely by accuracy of the temporal discretization. Thus neither scheme is clearly superior; comparison of these schemes for other transition simulations is ongoing.

3. Results 3.1. Axisymmetric

calculations

The first set of results presented here pertain to the axisymmetric two-cell/one-cell bifurcation, which occurs when the Taylor apparatus has an aspect ratio up to about 1.5. The form of the bifurcation is sensitive to this parameter; experiments show that this transition can change from supercritical to subcritical with variations in the aspect ratio of as little as 8% [l]. (See [34] for a more complete set of results than is presented here.) We investigate the geometry of Aitta, Ahlers and Cannel [l], where the radius ratio is 0.5. Their experimental results relate to three aspect ratios: 1.129, 1.266 and 1.281; for which the two-cell/ one-cell bifurcation is supercritical, transcritical, and subcritical, respectively. The order parameter used by Aitta, Ahlers and Cannel to quantify the asymmetry of the mode is:

+‘=

Jt”w(r, -=I I+“]

w(r,

z) dz z) 1 dz ’

J-z

where r = R, + 0.14( R, - Ri). The comparison of steady-state solutions with experiment in terms of order parameter versus Reynolds number for the three aspect ratios (1.129, 1.266 and 1.281) is shown in Figs. 3-5. Note that the agreement is generally good, with the critical Reynolds number predicted well in all three cases. The steepening of the #’ versus R curve near the critical point as r is increased is

52

C.L. Street& M. Y. Hussaini 1.0 -

V

flow

140

150

ss

0

.a-

/ Chaos in Taylor-Couette

exp

-

.6 .4 -

120

130 R

Fig. 3. Order parameter

of one-cell

mode versus

0

1.0 , ,a

w

Reynolds number, ment; T = 1.129.

comparison

of steady-state

results

with experi-

ss

exp

-

.6 .4 .2 0. I 130

140

150

160

R

Fig.

4. Order

parameter

for one-cell

0.

3

140

mode

versus Reynolds number, experiment; r = 1.266.

150

comparison

160

of steady-state

results

with

170

R

Fig. 5. Order

parameter

for one-cell

mode versus Reynolds number, comparison results with experiment: r = 1.281.

of steady-state

and time-accurate

also qualitatively captured. Figure 5 also shows results from time-accurate computations, which method also indicated the compare well with the steady-state solutions. The time-accurate presence of nonunique states and hysteresis, characteristic of this subcritical bifurcation. The time-accurate method produced both symmetric and asymmetric solutions at R = 145, depending on whether the Reynolds number was being raised or lowered, respectively; a very slow growth of the one-cell mode was observed at R = 147, indicated in Fig. 5 by the symbol with the upward arrow. This corresponds well with the point observed in the experiment at which the

C.L. Streett, M. Y. Hussaini

/ Chaos in Taylor-Couetteflow

53

time-accurate Fig. 6. Comparison

of cross-section

streamlines,

time-accurate

and steady-state

methods;

r = 1.281, R = 150.

symmetric mode lost stability altogether and the flow jumped to the one-cell mode. A comparison of cross-section streamlines at R = 150from the time-accurate and steady-state methods is shown in Fig. 6, indicating good pointwise agreement between the solutions. The results in Figs. 3-5 indicate a systematic difference between experimental and numerical results in terms of the level of IJ’, especially at higher Reynolds numbers. As described above, the order parameter used here is based on the axial velocity measured along a line at a particular radial position. Examination of the numerical results indicates that the order parameter is highly sensitive to this location. On the streamline plot from the steady-state method in Fig. 6, the measurement location is denoted by the leftmost arrow. When the measurement location is changed by just 0.02 gap (denoted by the rightmost arrow in Fig. 6) agreement with the experimental value of 4’ is obtained. Thus, the relatively small disagreement between computation and experiment may be attributed to a very small error in measurement location or to very small differences in the flow field in this region. Computations of axisymmetric bifurcation phenomena were also carried out for the geometry of Benjamin [6], for the range of aspect ratios in which a complicated set of two-cell/four-cell bifurcations occur. The numerical results uncovered an extreme sensitivity of bifurcation points to dynamic effects; specifically, to the finite rate of change of Reynolds number used to locate the bifurcation points. These effects, we feel, manifested in the experiments as a hysteresis in a modal change, where the steady-state calculations showed none; time-accurate calculations with a finite (but very small) rate of change of Reynolds number reproduced the experimental results. (See [35,36] for details.) Most of the axisymmetric bifurcations studied in this work concerned so-called primary and secondary modes; that is, modes which could be realized by quasistatic processes along

54

C.L. Street& M. Y. Hussaini Primary

Anomalous

I

Fig. 7. Cross-section

streamlines,

primary

/ Chaos in Taylor-Couette Mode

Mode

flow

(2.cell)

(3.cell)

I

two-cell mode and anomalous three-cell rate method; r = 3.0, R = 315.

mode, by multi-domain

continuous paths in the R-T space [5]. We also investigated “anomalous” typically produced experimentally by violent changes in the Reynolds number. three-cell asymmetric configuration, is shown in Fig. 7 for r = 3.0; also shown the two-cell primary mode for the same geometry and Reynolds number. discretization was used for these calculations; note that the cross-sectional smooth and continuous through the multi-domain interface line in the center. 3.2. Three-dimensional

time-accu-

modes, which are One such mode, a for comparison is The multi-domain streamlines are

calculations

In this section we describe results from our simulations of unsteady three-dimensional Taylor-Couette flow. A relatively low aspect ratio of 2.4 was chosen for this study to ease axial resolution requirements; since a large percentage of the early experimental investigations of higher Reynolds number Taylor-Couette flow with low aspect ratios was carried out by Benjamin, Mullin and coworkers, the radius ratio of their apparatus (0.615) was selected. Calculations of the very weak first wavy modes used a 33 x 21 X 16 mesh (z, Y, ~9) with

C.L. Streett, M. Y. Hussaini

/ Chaos in Taylor-Couette

Fig. 8. First wavy mode axial velocity in r - 0 plane at mid-height;

flow

55

T = 2.4, R = 1790.

subsequent refinement to 45 X 33 X 32 and finally to 65 X 45 X 32 as the Reynolds number was increased. Again, a very small time step, equivalent to a maximum Courant number of about 0.05, was required for accuracy. Data interpretation for these flows has been quite difficult, and is still ongoing. The complicated structure and inherently unsteady nature of these flows renders simple planar contour plots inadequate to convey all but the simplest flow characteristics. Visualization via animated color graphics has been quite useful; a 15-minute movie depicting the unsteady modes described below has been produced which partially captures some of the details of these flows. The primary axisymmetric flow state for r = 2.4 at low Reynolds numbers is two-cell; there are no secondary modes, either one-cell or four-cell, and no anomalous modes have been discovered experimentally at this aspect ratio. Thus, the onset of azimuthal waviness is free of uncertainty due to basic state selection. The basic state was seeded with a small-magnitude (- lo-‘) perturbation in the axial velocity; the perturbation was random in the u and 0 directions and constant in the z direction, so the perturbed velocity field satisfied the continuity equation. The axisymmetric state became unstable to the perturbation and energy in higher azimuthal modes began to grow, at about R = 1500;no attempt was made to determine this critical point precisely. This point is however in rough agreement with experiment [27]. The growth rate of the azimuthal modes is very low at Reynolds numbers just above critical. At R = 1790, the flow settles into a wavy pattern dominated by the m = 5 azimuthal wave number. The wave structure is highly swept, as can be seen in the contours of axial velocity shown in Fig. 8. The contours are taken on the mid-height plane as indicated in the sketch; the axial velocity is zero on this plane for the basic state, so variations are purely the result of

56

C. L. Streett, M. Y. Hussaini / Chaos in Taylor-Couette

.5 gap

.3 gap

~3gap

.7 gap Fig. 9. First wavy mode axial velocity in unwrapped

flow

z-8

plane;

r = 2.4, R = 1790.

cylinder at the bottom, outer at top, waviness. The r-8 plane is also shown “unwrapped’‘-inner periodic to left and right-to show the swept structure more clearly. As can be seen, most of the waviness is this mode is near mid-gap. This is further exemplified in Fig. 9, in which axial velocity contours on unwrapped z-6’ “cylinders” are shown for various radial positions. The waviness affects most of the axial extent of the flow near the inner cylinder, whereas it is quite confined to near the mid-height plane near the outer cylinder. This can also be seen in the plots of cross-sectional stream function shown in Fig. 10; waviness distorts the cells in a triangularshaped region from the inner cylinder top and bottom to the outer cylinder mid-height point. The two cross-sections plotted, taken one-half cycle out of phase of each other, are almost identical after top-to-bottom reflection. To help classify the oscillatory flow, time histories of the axial velocity were taken at several radial positions on the mid-height plane. The history taken at 0.10 gap for this first wavy mode is shown in Fig. 11, and is seen to be quite regular. The temporal frequency spectrum of this history is shown in Fig. 12, in log-log form. The dominant frequency of oscillation and two higher harmonics are clearly visible, the dominant frequency being about 1; the frequencies are nondimensionalized by the period of the inner cylinder rotation. Two very much weaker components appear, one at the lower frequency of about 0.2, the other at about 3.2. From this weakly oscillatory mode, the Reynolds number was raised rapidly to 2065. There the flow commenced a long and complicated transient period, showing a number of distinct

C. L. Street& M. I’. Hussaini

Streamlines

/ Chaos in Taylor-Couette

in azimuthal

57

flow

plane

1

t

Fig. 10. First wavy mode cross-section

stream function;

r = 2.4, R = 1790.

.06

.02

W -.02

-.06

-.,ot 0

)

I

.l

I

I .2

1

I .3

1 / .4

1 .5

I

I

.6

.7

I

I

.a

/

I .9

1

J 1.0

time

Fig. 11. First wavy mode time history of axial velocity at 0.10 gap; T= 2.4, R = 1790.

58

C.L. Streett, M. Y. Hussaini

-7t11, -2.0

-1.5

/ Chaos in Taylor-Couette

I1

1,

-1.0

-.5

It 0

I .5

I

I 1.0

/

flow

11

I 2.0

1.5

log Weq)

Fig. 12. First wavy mode frequency

spectrum

of axial velocity history at 0.10 gap; r = 2.4, R = 1790.

r!YY C-

-\

__--’

,

--

t=o

t = .26

t = .39

Fig. 13. Second wavy mode axial velocity in unwrapped r-13 plane at mid-height, rotations; r = 2.4, R = 2065.

sequence

spacing 0.13 inner cylinder

CL. Streett, M. Y. Hussaini

/ Chaos in Taylor-Couette

Fig. 14. Second wavy mode axial velocity in r- 0 plane at mid-height,

flow

first frame of Figs. 13; r = 2.4, R = 2065.

Fig. 15. Third wavy mode axial velocity in r ~ 0 plane at mid-height;

r = 2.4. R = 2065

59

60

C. L. Streett. M. Y. Hussaini

/ Chaos in Taylor-Couette

flow

,3 !YP

5

7 gap

3 gap

Fig. 16. Third wavy mode axial velocity in unwrapped

z - 13plane;

gap

r = 2.4, R = 2065.

modal structures. The first mode of this sequence (which we refer to overall as the second mode) is somewhat difficult to present via static displays. The essential feature is the appearance near the outer cylinder of another m = 5 mode, which moves slower than the original m = 5 wavy structure. The original swept wave slides past this new structure, which rotates at a speed about one-fourth of the speed of the original wave. Another feature which appears is a wave on the original swept structure; an extrema appears on the swept wave near the inner cylinder, and moves radially outward along the wave. These features can be observed in Figs. 13, which display a sequence in time of contours of axial velocity on the mid-height plane; the disk has been unwrapped (r-0) for clarity. The interval between displays is about 0.13 rotations of the inner cylinder. For illustration, the first frame of Figs. 13 is shown in disk form in Fig. 14. The more slowly moving structure near the outer cylinder grows strongly with time and, eventually, dominates the flow in the third mode, of which axial velocity contours on the mid-height plane are shown in Fig. 15. This mode is rather simple and purely periodic, with its wavy structure aligned essentially radially, as opposed to the highly swept structures of the first and second modes. Note, on the other hand, that the azimuthal wave number is the same (m = 5). The oscillations affect nearly the entire axial extent of the flow and are very much strongest near the outer cylinder, as can be seen in Fig. 15 and in Fig. 16, which shows contours of axial velocity on unwrapped z-8 planes at various radial locations.

C. L. Streett, M. Y. Hussaini

Fig. 17. Fourth

/ Chaos in Ta_vlor-Couette flow

wavy mode axial velocity in r - 0 plane at mid-height;

61

r = 2.4, R = 2065

We believe that the regularly structured third mode is stable at some intermediate Reynolds number between 1790 and 2065, although we have not yet verified this. It is, however, a transient in the flow given the large Reynolds number jump taken in the simulation after the first mode was investigated. The next mode in the sequence at R = 2065 (the fourth mode) involves an increase in azimuthal structure complexity. Axial velocity contours on the mid-height plane for this mode (Fig. 17) show the appearance of a higher-frequency structure imposed on the radially oriented oscillations of the third mode. The complicated but regular oscillations span the entire axial extent, and are quite strong near the outer cylinder, as can be seen in Fig. 18. The time history of axial velocity shown in Fig. 19, was taken at a radial position of 0.97 gap. This history shows the progression from the essentially purely periodic third mode to the more complicated “quasiperiodic” fourth mode. Note that the magnitude of oscillation is almost two orders of magnitude stronger than for the first mode (Fig. 11). The temporal frequency spectrum for mode four, taken from the fourth-mode time history, is shown in Fig. 20, and indicates the two major frequency components present. The dominant frequency, approximately 0.2, was one of the low-energy components visible in the first-mode spectrum (Fig. 12). The quasiperiodic fourth mode is also probably stable at a somewhat lower Reynolds number; at R = 2065, however, the fourth mode is also a transient. After a period of time, the flow becomes disordered for a small interval, as can be seen in the time history shown in Fig. 21. After the aperiodic disruption, the flow again relaxes to the quasiperiodic fourth mode. Three

CL. Streett, M. Y. Hussnini

/ Chaos in Taylor-Couette

flow

.95 gap

.85 gap Fig. 18. Fourth wavy mode axial velocity in unwrapped

z - 8 plane; r = 2.4, R = 2065.

such ~sruptions were observed in this particular simulation, as shown in Fig. 22. The temporal frequency spectrum taken from the last period of disordered behavior (Fig. 23) shows the two major frequency components still visible, but against a broad-band noise background.

-6 k -10

0

11 .25

1

I*

30

11 -75

11 1.00

1

1.25

’

1

1.50

’

11

1.75

I”‘1

2.00

’

2.25

2.50

”

2.75

1

3.00

time Fig. 19. Third and fourth wavy mode time history of axial velocity at 0.98 gap; r = 2.4, R = 2065.

C. L. Street& M. Y. Hussaini

-5

I

1

-2.0

Fig. 20. Fourth

1

wavy mode frequency

1

1

-1.0

-1.5

/ Chaos in Taylor-Couette flow

1

1

1

-.5

spectrum

’

0

1

1

’

of axial velocity history

1

’

1.0

.5

63

J

’

1.5

2.0

at 0.98 gap; r = 2.4, R = 2065

After this third disruption, the flow relaxes into a more ordered quasiperiodic mode, with a completely different set of dominant frequency components. The time history of this transition is shown in Fig. 24. This lower-frequency quasiperiodic mode, the “sixth mode”, is an m = 2 azimuthal oscillation, as opposed to the m = 5 character of the previous mode, as can be seen in the mid-height contours of axial velocity shown in Fig. 25. The frequency content of this mode is fairly rich in discrete components, as shown in its history and corresponding frequency spectrum, Figs. 26 and 27. Visible in the spectrum are: the dominant frequency (wr - 0.25) and three higher harmonics; a very low frequency component (w,, - 0.02) and a harmonic; a component at a slightly lower frequency than dominant and one harmonic, which can be seen to be (w, - w,); and a rather broad component at about 0.8, again with a higher harmonic. The broad-band noise which characterized the previous mode has essentially disappeared, but broadening of the dominant component is apparent.

6

-6 -10

I

0

I

s

I

I

I

1I

1.5

1.0

I

2.0

2.5

5.0

3.5

I

4.0

I

I

4.5

I

I

5.0

time

Fig. 21. Appearance

of disorder

in fourth-mode

time history R = 2065.

of axial velocity

at 0.98 gap, early

time;

T = 2.4.

CL. Streett, M. Y. Hussaini / Chaos in Taylor-Couette flow

64

6

-6 -lot

I

0

I

I

I

I

I

2

1

I

3

I

I

4

1

I

I

5

6

I

I

I

I

7

I

8

I

I

I

9

10

I

I

1

11

I 12

time Fig. 22. Appearance

of disorder

in fourth-mode

time history of axial velocity at 0.98 gap; r = 2.4, R = 2065.

Animated graphic displays of cross-sectional streamlines show another qualitative difference between the m = 2 quasiperiodic sixth mode and the preceding modes. In the earlier modes, the azimuthal oscillations displaced the interface between the two vortices without significantly affecting the position of the centers of the vortices. In the sixth mode, however, the vortex centers assumed a “heaving” motion, moving together and apart a large distance. During the disordered fifth mode, the vortex centers alternated between being stationary and moving wildly about. The m = 2 quasiperiodic sixth mode is the final and stable mode for R = 2065; the simulation was carried for a very long period of time to verify this. This mode also remains stable on

-5t -3.0

Fig. 23. Frequency

spectrum

1 -2.5

I

I

-2.0

-1.5

,

1 -1.0

/ -.5

,

j 0

of axial velocity time history of disordered

,

I

.5

1.0

flow; r = 2.4, R = 2065.

C. L. Streett, M. Y. Hussaini

/ Chaos in Taylor-Cauette

flow

65

6

-6

-lOL-L~, .5

0

I

1.o

1s

,

2.5

2.0

I

I , I , I , J 3.5 4.0 4.5 5.0

3.0

time

Fig. 24. Time history of axial velocity,

transition

from disordered

flow to sixth mode;

Fig. 25. Sixth wavy mode axial velocity in r - 0 plane at mid-height;

0

I

I

1

2

I

I

I

3

4

1

I 5

)

I 6

1

I 7

I 8

I 9

r = 2.4, R = 2065.

r = 2.4, R = 2065.

I

I 10

1

I 11

time Fig. 26. Sixth wavy mode time history of axial velocity at 0.98 gap; r = 2.4, R = 2065

12

66

CL. Streett, M. Y. Hussaini / Chaos in Taylor-Couette

,;-3.o

-2.5

-2.0

-1.5

-.5

-1 .o

flow

0

.5

1.0

log We4 Fig. 21. Sixth wavy mode frequency

spectrum

of axial velocity history

at 0.98 gap; r = 2.4, R = 2065.

increase of Reynolds number up to at least R = 2250, with no qualitative change and slight quantitative frequency shift. On increase of Reynolds number to 2500 the aperiodic fifth mode, however, remains quite disordered, as can be seen from the long time history shown in Fig. 28; a small-scale portion of this history is shown in Fig. 29. Despite the large and aperiodic oscillations of this flow, the two-cell nature of the basic state remains intact (Fig. 30); the disordered nature of this flow is manifested in its temporal behavior, but is not strongly observed in its spatial character.

6

w 2

-2 -6 I

-lot 0

I 1

I

I 2

1

I

3

1

I

4

1

I

5

1

I

6

I

I

7

I

I

8

1

I

1

9

time Fig. 28. Chaotic wavy mode time history of axial velocity at 0.98 gap; r = 2.4, R = 2500.

I 10

C. L. Streett, M. Y. Hussaini

/ Chaos in Taylor-Couette

flow

67

6 w 2

-6

-10 6.0

6.2

6.4

6.6

6.8

7.0

7.2

7.4

7.6

7.8

8.0

time Fig. 29. Chaotic

wavy mode time history of axial velocity at 0.98 gap, detail;

r = 2.4, R = 2500.

The spectrum of the history from the R = 2500 flow shows no dominant frequency components; the broad-band noise drops off rapidly above a frequency of about 0.5 (Fig. 31). As pointed out by Brandstater and Swinney [g], the rate of decay of the broad-band component can be used to distinguish between chaotic but deterministic behavior and that of a stochastic process which can be considered fully turbulent. In the latter case the decay is algebraic, whereas the spectrum in the former case falls off faster than any finite algebraic power, which implies exponential decay. Plotted in semi-log form, the spectrum for this flow (Fig. 32) shows a clear linear behavior, indicating that the flow is chaotic in the dynamical systems theory sense.

Fig. 30. Chaotic

wavy mode cross-section

stream function;

T = 2.4, R = 2500.

68

L. Street&

-5t1

I

-3.0

-2.5

Fig. 31. Chaotic

I

Y. Hussaini

I1 -2.0

wavy mode frequency

I

Chaos in TaylorpCouette

II

I

-1.0

-1.5

spectrum

/

11l1lI -.5

0

of axial velocity history

.5

flow

J

111’

1.0

1.5

2.0

at 0.98 gap; r = 2.4, R = 2500.

0

-1

log (W) -2

-3

-4

I

-5

I

0

/

.5

1.0

11

1.5

I

2.0

I 2.5

II

II

I

3.0

3.5

4.0

I

’

1 4.5

’

1 5.0

freq Fig. 32. Chaotic

wavy mode frequency

spectrum

of axial velocity history

at 0.98 gap; T = 2.4, R = 2500.

4. Conclusions The spectral collocation techniques described here have proved to be accurate and efficient methods for solution of the incompressible Navier-Stokes equations. Algorithm developments during the process of this work have been carried over into a number of other simulations. The major feature of the time-split scheme described here is its simplicity, and the resulting low machine time per time step. The disadvantages are the requirement for very small time steps at high spatial resolution levels to keep splitting errors small, and errors due to the imposed approximate pressure boundary condition. A multi-domain implementation was also presented.

C.L. Streett, M. Y. Hussaini

/ Chaos in Taylor-Couette

flow

69

The staggered-mesh scheme requires no artificial pressure boundary conditions and is completely self-consistent, producing no spurious modes. Also, the coupled solution algorithm gives direct control over the level of mass conservation in the field. The multi-grid algorithm relieves the convergence difficulties of the conjugate gradient method used for solution of the Stokes operator in the time-dependent case. The time step for the coupled unsteady multi-grid algorithm is limited only by accuracy considerations of the explicit treatment of the nonlinear convection terms (up to O(1) Courant numbers), and may be orders of magnitude larger than the step size for which the splitting errors of the time-split scheme are acceptable. The algorithm is, however, very complicated and fairly expensive per time step. One out of our extensive investigations of axisymmetric bifuractions was described in Section 3. Reasonable agreement with experiment was shown for the very sensitive two-cell/one-cell bifurcation observed at low aspect ratios, including the change-over from supercritical to subcritical bifurcation as the aspect ratio was increased. A potential explanation of some of the systematic error between the experiment and computations was proposed. Three-dimensional unsteady simulations of the moderate Reynolds number wavy-vortex flow were carried out for a particular geometry. Seven distinct modes were observed in sequence as the Reynolds number was increased; many of these were transient modes, although we feel several would be stable at Reynolds numbers intermediate to those held for long periods in the simulations. Purely periodic and quasiperiodic states with varying degrees of spectral complexity appeared in the simulations. Both wave propagation speed and azimuthal mode changes were also observed. The final state simulated showed temporally aperiodic behavior; the temporal frequency spectrum of this mode was broadband, with exponential decay at high frequency. This feature indicates the flow to be low-order chaotic. We are in the process of analyzing these flow states. Animated color graphics provides necessary visualization of the simulation results. The dynamical systems character of these states, especially the chaotic one, is being studied using the now fairly standard tools of phase-space portraits and computations of Lyapunov exponents to estimate the dimension of the low-order attractor.

Acknowledgement

The authors would like to thank M.D. Salas, Head, Theoretical Aerodynamics Branch, NASA Langley Research Center, for his continued support and enthusiasm for this work, and the opportunity to present it here.

References [l] A. Aitta, G. Ahlers and D. Cannel, Tricritical phenomena in rotating Couette-Taylor Flow, Phyx Rev. Lett. 54 (1985) 673-676. [2] T. Alziary de Roquefort and G. Grillaud, Computation of Taylor vortex flow by a transient implicit method, Comput. &Fluids 6 (1978) 259-269. [3] C.D. Andereck, S.S. Liu and H.L. Swinney, Flow regimes in a circular Couette system with independently rotating cylinders, J. Fluid Mech. 164 (1986) 155-183.

70

[5] [6] [7] [8] [9] [lo] [ll] [12] [13] [14] [15]

[16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29]

C. L. Streett, M. Y. Hussaini

/ Chaos in Taylor-Couette

flow

Begue, Q. Dinh, B. Mantel, J. PCriaux, G. Terrasson, B. Cardot, F. El Dabaghi, F. Hecht, R. Muiioz, C. Pares, 0. Pironneau, M. Abdalas and R. Glowinski, Current progress on the numerical simulations of detached flow around airplanes, in: C. Taylor, W. Habashi and M. Hafez, eds., Numerical Methods in Laminar and Turbulent Flow (Pineridge, Swansea, 1987). T.B. Benjamin, Bifurcation phenomena in steady flow of a viscous fluid I: Theory, Proc. Roy. Sot. London A 359 (1978) l-26. T.B. Benjamin, Bifurcation phenomena in steady flow of a viscous fluid II: Experiment, Proc. Roy. Sot. London A 359 (1978) 27-43. C. Bernardi and Y. Maday, A staggered grid spectral method for the Stokes problem, in: Proceedings 6th International Symposium on Finite Element Methods in Flow Problems, Antibes, France (1986). A. Brandstater and H.L. Swinney, Strange attractors in weakly turbulent Couette-Taylor flow, Phys. Reu. A 35 (1987) 2207-2220. A. Brandt, Multi-level adaptive solutions to boundary value problems, Math. Comp. 31 (1977) 330-390. M.O. Bristeau, R. Glowinski and J. PCriaux, Numerical methods for the Navier-Stokes equations: Applications to the simulation of compressible and incompressible viscous flows, Comput. Phys. Report (1987). G.S. Bust, B.C. Domblaster and E.L. Koschmieder, Amplitudes and wavelengths of wavy Taylor vortices, Phys. Fluids 28 (1985) 1243-1247. C. Canuto, M.Y. Hussaini, A. Quarteroni and T.A. Zang, Spectral Methods in Fluid Dynamics (Springer, New York, 1988). K.A. Cliffe and T. Mullin, A numerical and experimental study of anomalous modes in the Taylor experiment, J. Fluid Mech. 153 (1985) 243-258. R.C. Di Prima, Transition in flow between rotating concentric cylinders, in: R.E. Meyer, ed., Transition and Turbulence (Academic Press, New York, 1981). R.C. Di Prima and H.L. Swinney, Instabilities and transitions in flow between concentric rotating cylinders, in: H.L. Swinney and J.P. Gollub, eds., Hydrodynamic Instabilities and the Transition to Turbulence (Springer, New York, 1981). P.R. Fenstermacher, H.L. Swinney and J.P. Gollub, Dynamical instabilities and the transition to chaotic Taylor vortex flow, J. Fluid Mech. 94 (1979) 103-128. M. Fortin, R. Peyret and R. Temam, Resolution numerique de equations de Navier-Stokes pour un fluide incompressible, J. Mecanique 10 (1971) 357-390. M. Gorman and H.L. Swinney, Spatial and temporal characteristics of modulated waves in the circular Couette system, J. Fluid Mech. 117 (1982) 123-142. F.H. Harlow and J.E. Welch, Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface, Phys. Fluids 8 (1965) 2182-2189. C.A. Jones, Nonlinear Taylor vortices and their stability, J. Fluid Mech. 102 (1981) 249-261. H.B. Keller, Numerical solutions of bifurcation and nonlinear eigenvalue problems, in: P.H. Rabinowitz, ed., Applications of Bifurcation Theory (Academic Press, New York, 1984) 359-384. G.P. King, Y. Li, W. Lee, H.L. Swinney and P.S. Marcus, Wave speeds in wavy Taylor-vortex flow, J. Fluid Mech. 141 (1984) 365-390. A. Lorenzen, G. Pfister and T. Mullin, End effects on the transition to time-dependent motion in the Taylor experiment, Phys. Fluids (1983) 10-13. M.G. Macaraeg and C.L. Streett, Improvements in spectral collocation discretization through a multiple domain technique, Appl. Numer. Math. 2 (1986) 95-108. P. Marcus, Simulation of Taylor-Couette flow: Numerical methods and comparison with experiment, J. Fluid Mech. 146 (1984) 45-64. R.D. Moser, P. Moin and A. Leonard, A spectral numerical method for the Navier-Stokes equations with applications to Taylor-Couette flow, J. Comput. Phys. 52 (1983) 524-544. T. Mullin and T.B. Benjamin, Transition to oscillatory motion in the Taylor experiment, Nature 288 (1980) 567-569. G.P. Neitzel, Numerical computation of time-dependent Taylor-vortex flows in finite length geometries, J. Fluid Mech. 141 (1984) 51-66. S.A. Orszag, M. Israeli and M.O. Deville, Boundary conditions for incompressible flows, J. Sci. Comput. 1 (1986) 75-111.

C.L. Streett, M. Y. Hussaini

/ Chaos in Taylor-Couette

flow

71

[30] G. Pfister, Deterministic chaos in rotational Taylor-Couette flow, in: Lecture Notes in Physics 235 (Springer, New York, 1985) 199-210. [31] M.P. Ross and A.K.M.F. Hussain, Effects of cylinder length on transition to doubly periodic Taylor-Couette flow, Phys. Fluids 30 (1987) 607-609. [32] C.L. Streett, Spectral methods and their implementation to solution of aerodynamic and fluid mechanic problems, Znternat. .Z. Numer. Meths. Fluids 7 (1987) 1159-1189. [33] CL. Streett, Two spectral collocation algorithms for the incompressible Navier-Stokes equations with two non-periodic directions, in: Proceedings First International Conference on Industrial and Apphed Mathematics, Paris, France (1987). [34] C.L. Streett and M.Y. Hussaini, A numerical simulation of finite-length Taylor-Couette flow, AIAA Paper 87-1444 (1987). [35] C.L. Streett and M.Y. Hussaini, A numerical simulation of finite-length Taylor-Couette flow II, in: G. deVah1 Davis and C. Fletcher, eds., Computational Fluid Dynamics (North-Holland, Amsterdam, 1988). [36] CL. Streett and M.Y. Hussaini, Numerical study of bifurcations in finite-length Taylor-Couette flow, Proc. Roy. Sot. London (submitted). [37] D. Walgraef, P. Borckmans and G. Dewel, Onset of wavy Taylor vortex flow in finite geometries, Phys. Rev. A 29 (1984) 1514-1519. [38] T.A. Zang and M.Y. Hussaini, On spectral multigrid methods for the time-dependent Navier-Stokes equations, Appl. Math. Comput. 19 (1986) 359-372.