Linear stability of lid-driven cavity flow

Cite as: Physics of Fluids 6, 2690 (1994); https://doi.org/10.1063/1.868158

Submitted: 15 October 1993 • Accepted: 29 April 1994 • Published Online: 04 June 1998

Natarajan Ramanan and George M. Homsy

ARTICLES YOU MAY BE INTERESTED IN

Reynolds number and end-wall effects on a lid-driven cavity flow

Physics of Fluids A: Fluid Dynamics 1, 208 (1989); https://doi.org/10.1063/1.857491

Global stability of a lid-driven cavity with throughflow: Flow visualization studies

Physics of Fluids A: Fluid Dynamics 3, 2081 (1991); https://doi.org/10.1063/1.857891

Simulation of lid-driven cavity with top and bottom moving boundary conditions using
implicit finite difference method and staggered grid
AIP Conference Proceedings 2021, 020002 (2018); https://doi.org/10.1063/1.5062719

Physics of Fluids 6, 2690 (1994); https://doi.org/10.1063/1.868158 6, 2690

© 1994 American Institute of Physics.

Linear stability of lid-driven cavity flow
Natarajan Ramanan a ) and George M. Homsy
Department of Chemical Engineering, Stanford University, Stanford, California 94305
(Received 15 October 1993; accepted 29 April 1994)
Previous experimental studies indicate that the steady two-dimensional flow in a lid-driven cavity
becomes unstable and goes through a sequence of transitions before becoming turbulent. In this
study, an analysis of this instability is undertaken. The two-dimensional base flow is computed
numerically over a range of Reynolds numbers and is perturbed with three-dimensional
disturbances. The partial differential equations governing the evolution of these perturbations are
then obtained using linear stability analysis and normal mode analysis. Using a finite difference
discretization, a generalized eigenvalue problem is formulated from these equations whose solution
gives the dispersion relation between complex growth rate and wave number. An eigenvalue solver
using simultaneous iteration is employed to identify the dominant eigenvalue which is indicative of
the growth rate of these perturbations and the associated eigenfunction which characterizes the
secondary state. This paper presents stability curves to identify the critical Reynolds number and the
critical wavelength of the neutral mode and discusses the mechanism of instability through energy
calculations. This paper finds that the loss of stability of the base flow is due to a long wavelength
mode at a critical Reynolds number (Re) of 594. The mechanism is analyzed through a novel
application of the Reynolds-Orr equations and shown to be due to a Goertler type instability. The
stability curves are relatively flat indicating that this state will be challenged by many shorter
wavelength modes at a slightly higher Reynolds number. In fact, a second competing mode with a
wavelength close to the cavity width was found to be unstable at Re=730. The present results of the
reconstructed flow based on these eigenfunctions at the neutral state, show striking similarities to the
experimental observations.

I. INTRODUCTION The lid-driven cavity problem (see Fig. 1) was studied

The structure of the fluid flow in cavity driven by a by BurggraF who was interested in the applicability of the
moving lid is a classical problem that has been the subject of Prandtl-Batchelor theorem. The theorem states that a flow
wide-spread research for more than three decades. Such with closed streamlines, will, at infinite Reynolds number,
prominent interest can be attriouted to three reasons. First, contain a core of constant vorticity. Burggraf set out to verify
from a practical point of view, this problem has direct rel- this theorem numerically using finite differences. He ob-
evance to many of the manufacturing processes in the coat- tained converged solutions up to a Reynolds number of 400,
ing industry. Second, it is a prototype problem in which beyond which his numerical scheme failed to converge. The
many of the flow features (a shear flow, boundary layers, main characteristics of the lid-driven cavity flow in this Rey-
eddies, and a core vortex) present in flows with closed nolds number range are a core vortex, and vortices in the
streamlines may be studied. Finally, being simple in geom- corners (BR and BL vortex). The BR vortex is also referred
etry while exhibiting complex flow behavior, it is an ideal to as the downstream eddy (DSE). To infer the asymptotic
setting for testing the accuracy and efficacy of new numeri- behavior of flows in square and rectangular cavities, Pan and
cal methods to the solution of the Navier-Stokes equations. Acrivos4 studied the flow experimentally. They conducted
Our interest in this problem stems from the fact that the experiments in the Reynolds number range between 20 and
flow structure of high Marangoni number thermocapillary 4000. Their observations about the flow in a cubical box can
flow in a differentially heated cavity exhibits global flow be summarized as follows. (1) As the Reynolds number in-
structures akin to lid-driven cavity flow. Earlier numerical creases, the center of the core vortex shifts to the geometric
studies of Carpenter and Homsyl clearly illustrate this claim. center of the cavity and occupies the bulk of the cavity. (2)
The experiments of Gillon and Homsy2 discovered that be- The size of the downstream eddy increases in size up to a
yond a critical Marangoni number, the steady two- Reynolds number of 500 and soon after it diminishes in size.
dimensional thermocapillary flows in a box with a square The chief assumption of their study was that the flows ob-
cross section and a spanwise ratio of 3:1 becomes unstable served were two dimensional.
and three dimensional. The resulting three-dimensional A number of studies since then have focused on steady
structures also bear strong resemblance to those observed in two-dimensional cavity flows. Most notable among them are
lid-driven cavity flows. To better understand the stability of the accurate numerical simulations of Ghia, Ghia, and ShinS
thermocapillary flows, we felt compelled to first analyze the and Schreiber and KelIer6 which provide the following infor-
stability of lid-driven cavity flows. mation about steady two-dimensional flows. Up to a Rey-
nolds number of 500, the size of the downstream eddy com-
a)Currently, Research Scientist, Fluid Dynamics lnternational, Evanston, pares favorably with the experimental results of Pan and
Illinois 60201. Acrivos. 4 As the Reynolds number is increased further, the

2690 Phys. Fluids 6 (8), August 1994 1070-6631/94/6(8)/2690/12/$6.00 © 1994 American Institute of Physics
 v axial wave number and do not indicate if these vortices are
- - - - . . . . Moving Lid
steady or unsteady. In marked contrast, the calculations of
Ku et al. 13 are devoid of TGL vortices in a cubic cavity with
no-slip boundary conditions on the end walls. They used
pseudospectral methods and reported solutions up to Rey-
nolds numbers of 1000. It is probable that the TGL vortices
core vortex were suppressed by the stabilizing effects of sidewalls at
these low Reynolds numbers or improper initial disturbances
failed to stimulate these vortices.
Recently, Aidun et al. 14 conducted a careful experimen-
BR Vortex tal study with a spanwise aspect ratio of 3:1. A distinct fea-
BL Vortex
o 0 (DO)
ture of their experimental setup is that they allow for an
~-~------------~--~
Y,v
inflow near the bottom left portion of the cavity to replace
~x,u the fluid lost with the moving belt. The inflow removes one
of the lower eddies, however, the global flow characteristics
z
are assumed to be unaltered. Based on the visualization, they
draw the following conclusions: (1) Below Re=500, the flow
FIG. 1. Schematic of a lid-driven cavity flow.
is two dimensional. (2) Between Re=825 and 925, the flow
undergoes a transition to a three-dimensional, small-
DSE monotonically grows in size in contrast with the obser- amplitude, time periodic state. From their figures, it appears
vations of Pan and Acrivos. 4 Asymptotically, as Re increases, that three pairs of vortices appear in the span wise plane near
the center of the primary/core vortex shifts to the geometric the downstream eddy. In a later study, using a power-
center of the cavity. To understand the discrepancy between spectrum analysis, Benson et aI., 15 estimate the dimension-
the numerical and experimental results in regard to the size less frequency of this mode to be 0.1112. (3) Between Re
of the downstream eddy, Koseff and Streee and Freitas = 1000 and 1300, traveling waves appear, and the dynamics
et al. 8 conducted experimental studies. Working with cavities are more complex and do not appear to have a single fre-
with a spanwise aspect ratio of 3:1, they visualized flows quency. (4) At Re>1900, irregularly spaced, oscillatory
between Reynolds numbers of 1000 and 10 000. They con- mushroom-like structures appear. (5) They also noticed that
cluded that at the lower Reynolds numbers, the flow was multiple steady states are possible if the lid speed is suddenly
initially two dimensional and became observably three di- reduced from a higher value.
mensional beyond Re=3200. At that stage they also noticed From the experimental and numerical studies, we make
the presence of counter-rotating vortices with an axis parallel the following observations. It appears that the flow is most
to the moving lid, and located near the downstream eddy. definitely two dimensional below 500. At some critical Rey-
The vortices, called the Taylor-GoertIer-like (TGL) vortices, nolds number below 1000, there is a transition to a secondary
were present in the visualization plane normal to the base state. The only reliable experiments in this range of Re are
flow plane. They contended that the TGL vortices are a con- those due to Aidun et al. 14 The numerical studies conducted
sequence of the instability of the flow near the separating for this range of Re were limited to a wavelength range in the
streamline of the primary vortex and the downstream eddy. spanwise direction (i.e., less than or equal to the cavity
The flow near the separation streamline is similar to the flow width). The experiments conducted by Aidun et al. 14 with a
down a concave wall, which undergoes a centrifugal Goertler spanwise length of 3 allowed slightly larger wavelength
instability.9 The TGL vortices were nonstationary. At even modes to be present but if the flow is unstable to modes of
higher Reynolds numbers (>5000), they noticed that the even longer wavelength (>3), they would be suppressed. It
flow was turbulent. Prasad et al. 10 and Perng and Streetll is also not clear what causes the two-dimensional flow to
performed numerical simulations using finite differences in become unstable. The arguments put forth by Koseff and
cubic cavities at Re=3200 to show the presence of TGL Streee are based on visualization and currently, there is no
vortices observed in these earlier experiments. The vortices theoretical explanation of the mechanism of the instability.
were nonstationary and had a wavelength close to half the The objectives of this work are to analyze the linear stability
width of the cavity. of steady, two-dimensional cavity flow, to establish the criti-
The simulations of Kim and Moin 12 are noteworthy in cal value of the Reynolds number above, which the flow is
that they observed the presence of TGL vortices in a square unstable, and analyze the mechanism of the instability. The
cavity with periodic boundary conditions in the spanwise only previous relevant analysis (as distinct from the three-
direction. Their calculations showed that the TGL vortices do dimensional direct numerical simulations discussed above) is
not require the presence of end walls to initiate them. Work- a recent paper by Poliashenko and Aidun. 16 They analyze the
ing with a cubic cavity they first computed the two- stability of the two-dimensional steady cavity flow to two-
dimensional base flow, to which a random perturbation ve- dimensional oscillatory modes. For aspect ratios close to
locity was added. Assuming periodicity in the transverse unity, they find a subcritical Hopf bifurcation at Reynolds
direction, they found upon time integration, the presence of numbers on the order of 7000. We discuss below the rela-
two pairs of TGL vortices. The vortices appear only when tionship of this work to ours. Our analysis allows for three-
the Re>900. Unfortunately, they report results only for one dimensional disturbances, which all the available numerical

Phys. Fluids, Vol. 6, No.8, August 1994 N. Ramanan and G. M. Homsy 2691
and experimental work suggest are the first modes to become Substituting into the Navier-Stokes equations (1) and (2),
unstable. The main assumption of our work is that the cavity subtracting the base flow equations (3) and (4) and lineariz-
is of infinite axial extent, which allows for the eigenvalue ing, we obtain the following equations for the perturbation
problem to be decomposed into normal modes in the axial velocity subject to no-slip conditions on all boundaries:
direction. The problem is formulated in Sec. II, and Sec. III
V·V=o, (7)
explains and validates our numerical procedures. Section IV
contains our main results and the analysis of the instability v,[+v·VU+ U·Vv= ~ Vp' + (1/Re)V 2 v. (8)
mechanism, and Sec. V concludes the paper with a short
discussion. The perturbation vorticity transport equation is found in the
usual way taking the curl of (8) to obtain
II. THEORY (.rJ,t- VX(Uxw) - VX(vxO) = (lIRe)V 1 w. (9)
A. Formulation Here, the perturbation vorticity is represented by w, and the
base flow vorticity, D by Ok. By choosing a set of potentials
We consider lid-driven cavity flow in a channel that is
infinitely long in the transverse (Le., z) direction (see Fig. 1).
1/1 and ~
Let us assume that the fluid is incompressible, isothermal, V= VX(iI/l) + VX(j~) (10)
Newtonian, and that the motion in the fluid is driven by the
to force the velocity field to be solenoidal, we automatically
top boundaries moving at a constant speed of V. The cavity
satisfy the continuity equation (Kessler I7 ). In terms of the
is square in cross section with a width of L. If the nondimen-
potentials, the vorticity can be expressed as
sional scales for length, velocity and pressure are L, V and
pV2 , the dimensionless form of the Navier-Stokes equations w= (c/J,Xy - "'SY - "',Il' - c/J,zz - c/J,xx + "',yx, c/J,xz + c/J,yz)'
is . (11)

V·u=o, (1) The perturbation vorticity transport equation (9) can now be
expressed in terms of these potentials as
2
1
, e V u,
ut+u,Vu=-Vp+-R (2) Re[wl,t- (Ulwl- U1WI),y- U 1 w3,z+D~,zzJ= V1WI,
(12)
where V is the gradient operator (V=ia,x +ja,y+ ka,z), u the
velocity, p the pressure, t the time, and a comma indicates Re[ wz. t + (U 1 Wl - Uzwd.x - U z lJ)3,z-DI/I,zz] = Vzwz, (13)
differentiation. The Reynolds number, Re(=pVL/ p.,) is the where (WI ,W2) are the x and y components of secondary
only dimensionless parameter of the problem. No-slip condi- vorticity. Only Eqs. (12) and (13) are needed to determine
tions are imposed on all boundaries. the potentials as the z component equation for secondary
vorticity can be shown to be a linear combination of (12) and
B. Base flow (13). If the vorticity components (II) are substituted into the
vorticity transport equations we have a coupled set of fourth-
As the cavity is infinitely long in z direction, the base
order partial differential equations for the potentials. The
flow whose stability is being examined is two dimensional
boundary conditions on the potentials are obtained from (10)
and steady, with the result that the equations for this flow
as
simplify to
I/Iz=~z=o, (14a)
V·U=o, (3)
~x= I/Iy. (14b)
(4) For our problem of a cavity of infinite extent in the axial
direction, we represent the potential in terms of normal
where V * (V * =ia,x +ja,j) is the two-dimensional gradient modes in the usual way:
operator, U the base flow velocity, and P the base flow pres-
sure. These equations are subject to the boundary conditions 1/1= 'IJI'(x,y )e iKZ + + complex
lTt conjugates, (15)
that
~=<P(x,y)eiKz+<Tt+complex conjugates, (16)
U=(1,O), at y= 1,
where K is the axial wave number and 0' the complex growth
(5)
U=(O,O), at x=O,x= 1, and y=O. rate. The assumed form of the eigenfunction is completely
general and allows for both steady and oscillatory modes,
C. Perturbation equations depending on whether the eigenvalue 0' is real or complex,
respectively. If 0' is real, the disturbances either grow or de-
To study the stability of the base flow to disturbances,
cay monotonically, the critical Reynolds number is that for
we need the equations that govern the evolution of these
which 0'=0. If 0' is complex, then the neutral condition is
perturbations. To this end, we perturb the base flow by a
O'r=O, and the onset of instability is oscillatory with dimen-
disturbance velocity v[ =v(x,y,z,t)] and the pressure by p'.
sionless frequency O'i' This normal mode form also includes
The total velocity u is then written as
the time-dependent two-dimensional instability of the steady
u=U+v and p=P+p'. (6) flow, (corresponding to K=O), examined previously by Po-

2692 Phys. Fluids, Vol. 6, No.8, August 1994 N. Ramanan and G. M. Homsy
liashenko and Aidun. 16 Substituting (15) and (16) into (12)
and (13), presents an eigenvalue problem with the growth
rate being the eigenvalue
Re[<TA 1- (U 1A 2 - UzAd,y+ K2 U 1A 3 - K
2
n<I>]
=(a,xx+a.yy-K2)Al' (17)

Re[<TA 2+ (U1A z - U 2A 1),x+ K2UzA3+ K2nqr]

=(a,xx+a,yy-K2)Az, (18)
where
_ ·2
A 1 - <I> ,xy - 'I' ,yy + L '1',
..1.2 = qr,xy - <I> ,xx + P<I>, (19)

..1.3 = ike <I>,x + <t> ,y),

subject to the following boundary conditions on all bound-
aries; FIG. 2. 2-D Steady flow at Re=1000. (Streamfunctlon contours in the BR

<t> ='1' = ° and <I>x = "¥y' (20)

vortex are not at the same intervals as the core vortex.)

The equations can be written compactly as

in either direction. Details of the solution methodology on a
12 related problem has been reported elsewhere. 19 The solutions
b 22
<t> = 0,
b ][qr] (21)
are compared with those of Ghia et al. 5 and excellent agree-
ment is observed. For instance, at Re= 1000 the following
where al1, a12, a21, an, b ll , b 12 , b 21 , and b 22 are
comparison is made from our base flow solutions with those
functions of the Reynolds number, the associated base flow
field and the wave number of the normal mode. By choosing of Ghia et al. 5 Figure 2 shows the flow structure at this Rey-
nolds number which, as is well known, consists of a primary
the two parameters, Re and K, and computing the base flow
vortex, secondary vortices at the bottom corners, and an in-
corresponding to the Re, the coefficients to these linear equa-
cipient vortex in the upper left-hand corner. In addition to
tions are determined.
these global features, Table I demonstrates excellent agree-
ment of local quantities such as maximum streamfunction
III. NUMERICAL METHODS and the location of centers of primary and BR vortices with
In the last two sections, the formulation of the problem those of Ghia et al. 5
and the equations that govern the evolution of the distur- The base flow was computed as discrete values of Rey-
bances have been described. In this section, the numerical nolds numbers of 10, 50, 100, 150, 200, 300, 400, 500; 600,
methods used to compute the base flow and solve the eigen- and 800. At the lower Reynolds numbers comparable or bet-
value problem for the identification of the critical state are ter accuracy was obtained. All the derivatives required to
outlined. compute the coefficients of the linear equations (21) are cal-
culated using central differences.
A. Base flow computations
The evaluation of the two-dimensional base flow is per- B. Formulation of the eigenvalue problem
formed using the streamfunction-vorticity form of Eqs. (3) Equations (21) governing the evolution of the perturba-
and (4). These equations are discretized using finite differ- tion potentials !/J and <p, pose a generalized eigenvalue prob-
ences on a uniform grid. All terms in the equations are rep- lem for the growth rate <T. The formulation of this eigenvalue
resented using central differences except for the convective problem as a function of Reynolds number and wave number
terms, which are approximated using Khosla-Rubin (KR) involves two steps; discretization of the differential equa-
differencing,18 a variant of upwind and central differences. tions, and evaluation of the coefficients for the matrices.
This difference scheme is second order accurate in space
[O(~x2)] with the property that it reduces to central differ- 1. Difference equations for the partial differential
ences when a converged solution is reached. A multigrid equations
method was used to accelerate convergence and help attain The fourth-order partial differential equations for the po-
solutions in very fine grids with as many as 257 grid points tentials !/J and <p, are discretized on a uniform grid using

TABLE I. Comparison of current study and Ghia et at. (in brackets) at Re=1000.

.pmax (primary vortex) .pm•• (BR vortex) Center (primary vortex) Center (BR vortex)

0.119 (0.118) 0.5313,0.5634 (0.5313,0.5625) 0.8593,0.1111 (0.8594,0.1094)

Phys. Fluids, Vol. 6, No.8, August 1994 N. Ramanan and G. M. Homsy 2693
finite differences. Central differences are used for all the C. Solution of the eigenvalue problem
terms. The asymmetric part of the equations involves a 25 The eigenvalue pencil of the real un symmetric eigen-
point stencil and the symmetric part a nine-point stencil. If value problem [Eq. (23)] contains real values and complex
(i,j) is a grid location, the algebraic equations are written as
conjugate pairs. Since the growth rate of a mode is of the
i+2 j+2 opposite sign as the eigenvalue, eigenvalues with positive
2: 2: [a k,''IJI' k,t + b k,t<I> k,il real parts indicate decaying modes and vice versa. The neu-
k=i-2 '=j-2 tral modes are those with a zero real part.
HI j+l To detect the onset of instability, one needs to identify
those modes that possess eigenvalues closest to the imagi-
+ 0' 2: 2: [Ck,,'IJI'k,t+ dk,,<I>k,l] =0, (22a)
nary axis. The leading or dominant eigenvalues that we seek
k=i-I '=j-l
then are those with the smallest real value. When only one or
i+2 j+2 a few eigenvalues are desired, the commonly used methods
2: [ek ,,'IJI'kl+
, . .
fkl<I>kl] are inverse iteration, subspace iteration and iterative tech-
k=i-2 '=j-2 niques based on Lanczos methods. These methods find the
i+l j+l eigenvalues with the smallest norm and the eigenvectors as-
sociated with them. To identify complex eigenvalues with the
+0' 2: 2: . [gk,,'IJI'k,l+hk,l<I>k,l]=O, (22b)
smallest real value, either a bilinear or an exponential trans-
k=i-l t=j-l
formation is necessary which while retaining the same eigen-
we used MACSYMA20 to evaluate the coefficients of these vectors modifies the system so that the leading modes of the
equations. Mter defining a derivative using central differ- transformed problem are the desired ones. Such a method
ences symbolically, we were able to carry out the difference maps the eigenvalues closest to the imaginary axes to those
equations for all the derivatives symbolically. MACSYMA then with the smallest magnitude of the new system. There is a
combined the algebraic equations, collected like terms and direct relationship between the eigenvalues of the trans-
generated the FORTRAN code for each of the coefficients. We formed matrix system and the original one. A short descrip-
used one-sided differences that were second-order accurate tion of these methods is given below.
at the boundaries. We found symbolic computation to be so
powerful that we defined the Navier-Stokes equations in its
1. Inverse iteration
vector form in the MACSYMA script. A change of base flow or
the form of the divergence free potential resulted in modifi- If the leading eigenvalue is real, then it is also the small-
cation of only the vector form of the definitions and MAC- est in absolute magnitude. Inverse iteration is then an effi-
SYMA carried out all the symbolic computations and manipu- cient technique to obtain both the eigenvalue and the eigen-
lations. The resulting difference equations can be written in vector. If i is the iteration count, Xi the current guess of the
the matrix form eigenvector, we solve repeatedly the set of linear equations,
Ax=ABx, (23) (24)
where x is a vector that contains the unknown values of the normalize x i+ 1 and evaluate the eigenvalue A(i + 1) ,
potentials 'IJI' and <1>, and A is the decay rate (= - if, the
Ai+l = (xJ+ lAxi+ d/(xf+ I Bx i + 1)' (25)
growth rate). It can also be seen from the equations that the
matrix A is unsymmetric and real, and matrix B is symmetric If the relative change in the eigenvalues at the i and i + 1 th
and real. iteration are smaller than a prescribed tolerance, we assume
convergence. The convergence rate of inverse iteration is
proportional to O(IAII/IA2J), so that convergence is rapid
only if the first two eigenvalues are well separated. The so-
2. Computation of the coefficients
lution to the linear system at each iteration can be achieved
The matrix elements depend on the parameters of the by a direct solver such as LU decomposition with partial
problem (Re and K) and the derivatives of the base flow pivoting or an iterative solver using conjugate gradients. The
quantities. These quantities are computed at each and every use of an iterative solver is computationally efficient when
grid point location in the base flow grid using central differ- solving a large set of algebraic equations. So, we attempted
ences of the base flow solution. The eigenvalue problem was to implement an iterative solver based on the biconjugate
discretized in a uniform grid much like the one used for the gradient solver, LSQR,zl But, we found that the convergence
base flow. However, for reasons of computational efficiency, rates are critically dependent on the choice of the precondi-
we resorted to coarser grids than the base flow. Since grids tioner and were unable to find a suitable preconditioner that
for both the base flow and eigenvalue problems are uniform could render efficient convergence. Consequently, we re-
and the number of mesh intervals are multiples of two, all sorted to the use of direct methods which we found to be
grid points in any given coarse grid have an equivalent grid quite efficient even with the use of LV decomposition per-
point (at the same coordinate location) in the finer grid of formed out of core. However, the disk space necessary to
base flow computations. This enables us to compute all the solve the problem out of core became a limiting factor in the
flow variables and the associated derivatives at all grid points size of the grids that we could utilize for the eigenvalue
of the eigenvalue problem. problem.

2694 Phys. Fluids, Vol. 6, No.8, August 1994 N. Ramanan and G. M. Homsy
2. Bilinear and exponential transformation 3.5 -'f-Ra= l750

--.."'~ A
To determine the eigenvalue with the smallest real part, --Ra=1800
bilinear or exponential transformation can be used. Storer22
3.0

\
/;
2.5 - - - Ra= 2000 /;
utilized bilinear transformation while attempting to capture
the leading modes of a magnetohydrodynamic problem. In S 2.0
'"
....
.
?:'i
this method a transformed eigenvalue problem of the type OJ)

~
.S [.5

(B- r{)A)x= a(B+ r(JA)x (26) S 1.0

0'"
is solved where rand () are parameters. The eigenvectors of
this transformed system are the same as those of the original
0.5
~~~/
0.0
system. The true eigenvalue A of the original system is re-
lated to the eigenval~e a as .{l.5
~-/'
A= (1- a)!r{)(1 + a). (27) 2.0 2.5 3.0 3.5 4.0

Wavenumber
We combined this with inverse iteration to find the smallest
eigenvalue of the transformed system. We found that the
leading eigenvalues are very sensitive to the choice of the FIG. 3. Damping rate as function of wave number for normal modes of
parameter r. In fact, different choices of r can lead to differ- steady periodic rolls of Rayleigh-Benard convection.
ent predictions of leading eigenvalues. Storer22 and
Christodoulou and Scriven23 also report similar experiences
found that a simultaneous iteration with 10-12 eigenvalues
in their studies.
is significantly more efficient than inverse iteration in most
With exponential transformation, one chooses a transfor-
cases.
mation matrix of type
(28) D. Validation study
whence the eigenvalues of the transformed system are etA. Since we are interested here in the three-dimensional
Polynomial approximation of the exponent are typically used instability of two-dimensional flows dependent upon a pa-
and we followed those recommended by Erikkson and rameter, we chose a related thermal convection problem to
Rizzi.24 To be precise, we solve the eigenvalue problem validate the code. To this end, we examined the stability of
[B- A( 1- (J)r]kx= a(B+ A{)r)kx , (29) finite amplitude two-dimensional Rayleigh-Benard convec-
tion. It is well known that a fluid heated from below is qui-
where A=(I-a)h[{)a+(I-{))], and e and r are parameters escent only below a critical Rayleigh number, Ra. Beyond
that control the iteration. this critical Rayleigh number a recirculating flow develops in
We found exponential iteration to be more reliable than the form of steady periodic rolls, with a flow strength and
bilinear iteration, but it is extremely expensive since it in- structure dependent on the value of Ra. As the Rayleigh
volves solution of the linear system k times for every itera- number is raised further the flow undergoes a sequence of
tion. We used it to confirm the fact that the leading eigenval- transitional states. The mechanism of the instabilities of
ues are in the subspace of the simultaneous iteration and these transitional states has been studied extensively by
used simultaneous iteration to obtain solutions at finer grids. Busse.26 We chose the stability of finite amplitude periodic
rolls between no-slip boundaries as our test case, since it is a
nontrivial two-dimensional flow. Clever and Busse27 show
3. Simultaneous iteration
that at low Rayleigh numbers (Le., close to the critical Ra
This is the desired approach to obtain a limited number = 1708) these steady flows undergo a cross-roll instability.
of eigenvalues and eigenvectors. The method as described by Since for infinite plane layers, the axes of rotation for the
Stewart and Jennings25 captures the prescribed number of rolls is arbitrary, the cross-roll instability replaces the rolls by
eigenvalues (say, m) with the smallest magnitude. The con- another set of rolls with a wavelength close to the original
vergence rate of the method is D(IAII/IAml). The algorithm is one but with an axis of rotation 90° to the original state. For
briefly described as follows. For details, the reader is referred our tests, we assumed that the Prandtl number of the fluid is
to Stewart and Jennings.25 Given an initial guess or random 7.0 and the wavelength of the periodic steady rolls is 71'. The
guess for the eigenvectors, do the following: (1) premultipli- base flow was computed as discrete values of Rayleigh num-
cation, Z=(A -lB)X, executed k times; (2) do matrix multi- ber using the multigrid method detailed above. We solved the
plications of the sort, G=ZTZ and H=XTZ; (3) solve the eigenvalue problem using inverse iteration for a range of
system GB=H; (4) obtain the full eigensolution of B, sort its wave numbers of the perturbation flow. The leading eigen-
eigenvalues in ascending order; (5) perform matrix multipli- values were always real indicating that the new state will
cation W=ZP, where P is a sorting matrix that locks the also be steady. The calculations were performed with two
vectors that have already converged; (6) normalize eigenvec- uniform grids of 17x17 and 33X33. The computed damping
tors and test for convergence; (7) return to 1. rate based on Richardson extrapolation is shown in Fig. 3 as
The convergence rate of the leading eigenvalue is clearly a function of wave number for various Rayleigh numbers. It
more favorable than that of inverse iteration. In fact, we is clear from the curves that neutral modes have a wave-

Phys. Fluids, Vol. 6, No.8, August 1994 N. Ramanan and G. M. Homsy 2695
 1.4
6.00000

1.2 ~-Re=IOO

--Re=200
~Re=300
4.00000

CampI
0
s-g 1.0
----Re=400
----Re=500
O.B

---
0 0 0 ~
0 ....
'"
2.00000 ~
, - --
0 0
0
a u
.S>
Co
0.6
a 0 0 a 0
0 0 0
0
0
0 §
0 0 o 0 a a 0 0.4
0 o 0 ,J),Q
0
0
0 0 0 0
.00000 o0 0 0 0 0 0
0 0 o 0 0 0
0 0 0.2
0 0 0 0
0 0 0 0
0 0
a a a
-2.00000 ,- n 0.0'
0 2 4 6 8 10
a a 0 (a)
I Wavenumber

-4.00000 f- r.-----!-
--/--- -
0
- -- 08
~-Re=IOO

07 --Re = 200
.00000 .20000 .40000 .60000 .80000 1. 00000
~-A-Re=300
Real value (X10+ 1)
:>6 --'-Re = 400
--Re = 500
'':')

FIG. 4. Distribution of eigenvalues for normal mode, K=3 at Re=50. " ",5
'"
-,-
L4
,...
,~

length close to 'TT. To determine the critical parameters, we 0.3

extrapolated the data near Ra=lSOO and K=3 using bicubic
interpolation. 28 To this end, we fitted a bicubic polynomial U2

for the damping rate in terms of the two parameters, Ra and 0,1
K. The coefficients of the polynomial are computed, once we
specify the function value (i.e., damping rate), the first de- 3 4 5 6 7 9

rivatives of the function with respect to the two parameters (b) Wave number
and the cross-derivative with respect to both parameters at
four points in the parameter space. We estimated the deriva- FIG. 5. Ca) Real part of damping rate versus wave number for various Re.
tives of the function using central differences. Based on the (b) Imaginary part of damping rate (wave speed) versus wave number for
extrapolation, the critical parameters are Ra = 1725 and various Re.
K=3.2S, which are in agreement with Clever and Busse's27
results for cross-roll instability.
grid was adequate for obtaining grid invariant results and for
IV. RESULTS FOR LID-DRIVEN CAVITY FLOW
the higher Reynolds numbers the finer 65X65 was necessary.
We now present the results for the driven cavity flow. Unlike the validation study, the leading eigenvalues in
The base flow was computed at Re=lO, 50, 100, 150, 200, this problem can either be real or complex. To determine
300,400,500,600 and SOO on grids as fine as 257X257. As how the eigenvalues are distributed in the complex plane we
mentioned earlier, the base flow derivatives were computed first used simultaneous iteration to capture all the eigenval-
on this grid and used for stability calculations. For comput- ues on a small grid (9X9) at a low Reynolds number of 50
ing the eigenvalues/damping rates, we performed grid con- and wave number of 3.0. Figure 4 shows the distribution of
vergence tests on grids from 9X9 to 65X65. We found that at the eigenvalues in the complex plane. It can be seen for these
the lower Reynolds numbers (Le., below Re=200) a 33X33 parameters, that the leading eigenvalue is a Hopf mode fol-

TABLE II. Leading eigenvalue (A) versus wave number (K) for various Re, ( ... ) indicates unconverged solutions.

K Re=100 200 300 400 500 600 800

0 0.73
0.50 0.34 0.28 0.18:!::0.26i
2 0.48 0.23 0.11 0.10 0.03 -0.05
3 0.55XO.25i 0.29:!::0.11i O.19:!::O.14i 0.13:!::O.15i 0.1O:!::O.16i
4 O.59:!::o.26i 0.30:!::0.28i O.18:!::0.29i O.13:!::O.31i 0.11:!::0.33i
5 O.68:!::O.36i 0.34:!::0.43i O.24:!::O.45i O.17:!::O.46i O.16:!::O.43i
6 0.77:!::0.48i 0.39:!::O.58i 0.21:!::0.63i 0.13:!::O.62i 0.08:!::0.62i 0.04:!::O.59i -0.01:!::0.57i
7 0.90:!::0.58i 0.41:!::O.67i 0.21:!::0.67i O.14±O.64i O.09±O.63i
8 1.05:!::0.63i O.45:!::0.72i 0.27:!::0.70i O.21:!::O.63i 0.15:!::0.57i
9 l.23±O.71i O.54±O.76i O.39±O.75i O.23:!::O.58i 0.15:!::0.33i

2696 Phys. Fluids, Vol. 6, No.8, August 1994 N. Ramanan and G. M. Homsy
lowed by some real values and complex conjugate pairs and 0.8
that all modes are damped. Since the magnitude of the lead-
0.7
ing eigenvalues are very similar, we conclude that at least at
low Reynolds numbers, a simultaneous iteration with 10-15 0.6
eigenvalues would capture the most significant ones (i.e., 0.5
those with the smallest real values). We confirmed this for B
t:: 0.4
the higher Reynolds numbers by employing exponential it- OJ)

eration. 'i'" 03
It was observed that to obtain convergence at the lower CI
0.2
wave numbers close to zero a grid finer than 65 X65 was
OJ
necessary. However, for the rest of the wave-number spec-
trum, a 65X65 grid was adequate for obtaining eigenvalues 0.0
with an accuracy of 10- 4. In fact similar behavior was also -0.1
observed at the lower Reynolds numbers. Grid convergence 0 @ D B ~ D 0 ~ B ~

required fewer grid refinements at the intermediate and larger Reynolds number
wave numbers. Since critical modes are in the intermediate
range (based on our observations at lower Re), we did not FIG. 6. Damping rates (real part) versus Re at K=2.0 and 6.0 for extrapo-
feel compelled to increase grid refinement any more than lation to critical state.
65X65. The eigenvalues were computed using Richardson
extrapolation to extrapolate the results to zero grid size. All
the computations were obtained using both simultaneous it- real and imaginary parts of the eigenvalue, respectively. Over
eration and inverse iteration. While simultaneous iteration the range of Re, the damping rates in Fig. 5(a) show a mini-
converged in all cases with less than 30 matrix multiplica- mum at K=2.0. It is seen that the flow is stable over the
tions for six leading eigenvalues with tolerance less than Reynolds number covered, but there is a clear trend toward
10-4 and better than 5X 10- 4 for the next six eigenvalues. instability as Re increases. Furthermore, the dispersion rela-
Inverse iteration, on the other hand, took more than 30 ma- tion ["-r(K)] becomes flatter as Re increases, indicating a ten-
trix multiplications and failed to converge at low- and high- dency for several modes to become unstable simultaneously.
wave-number extremes at the higher Reynolds numbers. The imaginary part of the eigenvalue is shown in Fig. 5(b). It
However, when inverse iteration did converge, the results is essentially zero for K~2.0 for all Reynolds number studied
compared favorably with simultaneous iteration. indicating stationary modes. For K~3.0, the eigenvalues are
Figures 5(a) and 5(b) show our principal results for the complex, indicating (damped) Hopf modes, with the dimen-

,~
~
~
§
~
~
~
~ ?Jr
1f
! ~
" " "/
\ I
\
, ,1
~ ~ ~
~ ~ - ::. ; -~ ~
~ ~ ~
~
~ ~ ;; :::::
%::
'Z -- -"'
~ ~ -~ ~

~
?
r =- ~
~
~ ~ ~
~ ~ : I
1
.
- - - -- - --
~

2
I \
, ...
-
-::
I

::-
I
~
;
.;:::
I
;

-
(a)

.. ......................
~ --
... , ; " " .. .
........... -._-------, ................. ..
~

(b)

FIG. 7. Reconstructed total flow field at the near critical state at Re=600, K=2.0. (a) View from the side plane of a section along midplane of base flow. (b)
View of base flow plane at spanwise locations of 'Tr/4 and 'Tr/2.

Phys. Fluids, Vol. 6, No.8, August 1994 N. Ramanan and G. M. Homsy 2697
 FIG. 9. Base flow at Re=600(streamfunction contours).

(a) To find the critical Reynolds number for any given wave
number, we perform an interpolation by plotting Reynolds
number against the growth rate (see Fig. 6) for various wave
numbers. We do this since at the critical Reynolds number,
the eigenvalue problem is singular and the estimation of the
true eigenvalues and eigenvectors is corrupted by spurious
ones. We observed that the appearance of spurious eigenval-
ues increased as we got closer to critical states. Based on the
interpolations of Fig. 6, we estimate the critical Reynolds
number as 600 at a critical wave number of 2.0. This neutral
mode is stationary. Using bicubic interpolation as described
in the validation study, the refined estimates of the critical
parameters are Rec =594 and K=2.12-2.1S. The shorter
wavelength mode (K=6) crosses the axis at a Re close to
700. This mode has a wavelength close to unity and a dimen-
sionless frequency (= 'A)27il of 0.1. As observed from Fig.
5(a), the stability curves are quite flat indicating that at
Re> Rec the secondary flow will be challenged by the Hopf
modes nearby. We comment here that all the three-
dimensional modes are significantly more unstable than the
two-dimensional modes (K=O) computed by Poliashenko
(b) and Aidun. 16 Critical Reynolds numbers are 0(600) for
three-dimensional modes while they are 0(7000) for two-
FIG. 8. Production (a) and dissipation (b) energies of the near-critical mode
dimensional modes in the same range of aspect ratios near
at Re=600, K=2.0. Production contours (A = - 0.64,B= -0.45, ... ,1 unity. This is in accord with the numerical and experimental
= 1.07) are plotted at intervals of 0.19 and dissipation contours (A =-0.88, evidence that three-dimensional modes are the first to appear.
B= -0.793, ... ,J=-O.047) at intervals of 0.093. It is instructive to examine the flow field that corre-
sponds to the critical mode. We reconstructed the total three-
dimensional flow field by combining the two-dimensional
sionless frequency increasing with wave number. Table II base flow with an arbitrary amount of the disturbance eigen-
summarizes the data in Figs. 5(a) and 5(b). vector. The reconstructed flow is then fully three dimen-
At the higher Reynolds numbers, (i.e., beyond Re=300) sional, which we show for the case of Re=600, K=2.0.
we found that another local minimum appears in the disper- While weakly damped, we consider it a good approximation
sion relation. This minimum occurs at a wave number of 6.0 to the steady three-dimensional flow. Figure 7(a) shows the
and is only slightly more damped than the mode at K=2.0. velocity vector field in the (y-z) plane along the axial mid-
The K= 2 mode crosses the imaginary axis first and is imme- plane of the cavity. As can be seen, it is periodic in the axial
diately followed by the K=6 mode at a slightly higher Re. It direction, and consists of a cellular structure similar to that
is also noted that the later mode is a Hopf mode with a observed in flow visualizations of Aidun et al. 14 and Koseff
dimensionless wave speed of 0.59. and Street7 and in the direct numerical simulations of Kim

2698 Phys. Fluids, Vol. 6, No.8, August 1994 N. Ramanan and G. M. Homsy
 (a)

L, L
(b)

FIG. 10. Reconstructed total flow field at the near critical state at Re=800, K=6.0. Ca) View from the side plane of a section along midplane of base flow. (b)
View of base flow plane at spanwise locations of 0.25 and 0.5.

and MoinY Figure 7(b) shows the flow structure in the (x- y) where
plane at two axial positions corresponding to the quarter
wavelengths of the secondary flow. These patterns show the
distortion of the base circulation due to out-of-plane motion.
These reconstructions illustrate how the vortex center can be
moved closer or further from the moving lid by significant In this equation, the first term is the production of energy
vertical velocities of either sign, and the corresponding effect While the second term represents the dissipation. While the
on the strength of the corner eddies. first term may be positive or negative, the dissipation is al-
To understand the mechanism of instability, we compute ways negative. Equation (31) may be further manipulated as
the production and dissipation energies of the eigenmode. To follows. Since all perturbed quantities may be written in a
do so we identify the net energy associated with the second- form similar to Eq. (15), and since both the Reynolds stress
ary flow by taking a scalar product of the perturbation veloc- and the dissipation are quadratic functionals of the distur-
ity with momentum equation (8). If we integrate this over the bances, the integrands may be broken into a z-independent
domain, V, corresponding to the parallelpiped containing one part, and a part proportional to e 2iKZ • Since this latter part
axial wavelength of the perturbed flow and define the net integrates to zero over V, we can reduce the right-hand side
energy E as of Eq. (31) to area integrals over Q:sox:sol, O:SOy:so1. We can
therefore study the spatial distribution of the integrands giv-
~
E=fv 2 u·u·
L L
dV (30) ing the mean production and mean dissipation rates, which
we do by computing them by second-order accurate differ-
we can obtain the classical Reynolds-Orr energy equation ences at each grid point.
(in indicial form) as We plot in Figs. Sea) and 8(b) the spatial variation of
these production and dissipation integrands, respectively.
~~ =- Iv llilljDij dV- ~e Iv dijd ij dV, (31)
The production terms show very clearly that the mechanism
of instability is related to regions of the flow near the down-

Phys. Fluids, Vol. 6, No.8, August 1994 N. Ramanan and G. M. Homsy 2699
stream eddy. It is seen by superimposing these contours with
the base flow shown in Fig. 9 that the production of distur-
bance energy is highly localized near the dividing streamline
between the core vortex and the downstream eddy. This is as
it should be if the mechanism is associated with a Goertler
type instability. In this interpretation, there is an unstable
stratification of centrifugal force associated with the free-
shear layer near the dividing streamline, which approxi-
mately mimics the flow on a concave wall. On the other
hand, the dissipation contours indicate that the energy is dis-
sipated near the walls and the shear layer at the top. This
local field information provides convincing support to the
speculation of Koseff and Street 7 that the instability is asso-
ciated with a Goertler mechanism.
As we mentioned before, there is a shorter wavelength
mode that corresponds to a wave number of 6.0 (wavelength
-1) that becomes critical at a Reynolds number of 730.
Since this secondary mode is very close to the first unstable
mode, we also investigated its properties. On observing the
velocity vector plots of the reconstructed flow based on this
eigenmode (Fig. 10), it is clear that this mode significantly (a)
distorts the flow near the upstream eddy. The flow near this
eddy (namely the BL vortex) is quite similar to that near the
DSE and is susceptible to the same mechanisms of centrifu-
gal instability. Convincing evidence of this appears from the
contours of the production and dissipation (Fig. 11) for this
mode. As can be seen the energy production is highest along
the dividing streamline of the upstream eddy and the shear
layer near the moving lid. The dissipation is primarily in this
region and in the viscous boundary layer on the left wall.
Why this mode is oscillatory while the primary mode is sta-
tionary is an open question.

V. DISCUSSION AND CONCLUSIONS

In this study, we examine the stability of the classical
lid-driven cavity flow. The steady, two-dimensional base
flow is computed using a multigrid method and the equations
governing perturbations to the base flow pose a generalized
eigenvalue problem for the growth rate of the normal mode.
We solve the discrete version of this eigenvalue problem
through inverse iteration and simultaneous iteration to iden-
tify the leading modes. Our results indicate that the flow (b)
becomes unstable at Re-594 to a long-wave stationary
mode with a wave number of -2.12. A second minimum in FIG. 11. Production (a) and dissipation (b) energies of the near-critical
the dispersion relation is observed at K=6. This dimension- mode at Re=800, K=6.0. Production contours (A=-0.12,B=0.0, ... ,J
less secondary mode K=6 is of shorter wavelength and is "'0.994) are plotted at intervals of 0.12 and dissipation contours (A
oscillatory with a frequency of 0.1. The new flow, visualized = -1.13,B = -1.01, ... ,J= -0.06) at intervals of 0.12.

by reconstructing the flow by addition of the base flow and

an arbitrary amount of the leading eigenmode shows striking
resemblance to the flow observed in the experiments of Ko- latory state at Re-825. The oscillatory state had a character-
seff and Street7 and Aidun et al. 14 and the numerical simula- istic dimensionless frequency at the onset of 0.11 and a
tions of Kim and Moin.l0 A novel application of energy wavelength estimated from Fig. 4(b) of Aidun et at. 14 of ap-
analysis shows very clearly that both modes are related to a proximately 1.0. Transition to a nonperiodic state is also ob-
Goertler type instability of the flow near the streamline served as Re increases (Aid un et at. 14). Interestingly, steady
which separates the core vortex from the corner eddies. axially periodic states are achieved in the experiments by a
Comparison with theory is ambiguous. The recent ex- sudden decrease in Re. Figure 9 of Aidun et at. 14 shows flow
periments of Aidun and co_workers 14,15 show the following visualizations of such states which are strikingly similar to
characteristics. Working with an apparatus with a spanwise our reconstructed 3-D flow fields. Aidun et al. 14 interpret
aspect ratio of 3:1, they first observe a transition to an oscil- their experiments in terms of perturbed bifurcations of mul-

2700 Phys. Fluids, Vol. 6, No.8, August 1994 N. Ramanan and G. M. Homsy
ticelled steady solution branches, some of which exhibit a three dimensional flow in a cavity," Proceedings of the 4th International
Hopf bifurcations, the imperfection being the weak axial Conference of Numerical Methods on Laminar and Turbulent Flow, 1985,
p.503.
flow driven by the no-slip end walls.
~P. G. Drazin and W. H. Reid, Hydrodynamic Instability (Cambridge Uni-
Our calculations which represent an unperturbed bifur- versity Press, New York, 1985).
cation, show transition to a steady 3-D flow at Re~594 with 1UA. K. Prasad, C. Y. Perng, and J. R. Koseff, "Some observations on the
an axial dimensionless wave number at approximately 2.1, influence of longitudinal vortices in a lid-driven cavity flow," AIAA Paper
No. 88-3654, 1988.
immediately followed by a Hopf mode becoming critical at II C. Y. Perng and R. L. Street, "Three dimensional unsteady flow simula-
Re~730, wave number ~6, and dimensionless frequency of tions: alternative strategies for a volume averaged calculation," Int. J.
0.1. Since the wavelength of the primary stationary mode is Num. Meth. Fluids 9, 341 (1989).
equal or slightly larger than the box size employed in the 12J. Kim and P. Moin, "Application of a fractional step method to incom-
pressible Navier-Stokes equations," J. Comput. Phys. 59, 308 (1985).
experiments, it is entirely likely to have been stabilized by I3H.C. Ku, R. S. Hirsch, and T. D. Taylor, "A pseudo-spectral method for
the no-slip sidewalls. We speculate that the experiments first the solution of the three-dimensional incompressible Navier-Stokes equa-
accessed the Hopf mode at Re-825 with a wavelength of tions," J. Comput. Phys. 70, 439 (1987).
1.0, since the transitional parameters are so similar, and that 14c.K. Aidun. N. G. Triantafillopoulos, and J. D. Benson, "Global stability
of lid-driven cavity with throughflow: Flow visualization studies," Phys.
the SUdden decrease in Re was able to access the super- Fluids A 3.2081 (1991).
critical stationary mode. These speCUlations can only be re- 15J. D. Benson and C. K. Aidun, "Transition to unsteady nonperiodic state
solved by further experimentation and a detailed bifurcation in a throughflow driven cavity," Phys. Fluids A 4,2316 (1992).
analysis, however. 16M. Poliashenko and C. K. Aidun, "A direct method for computation of
simple bifurcations," submitted to J. Comput. Phys. (1994).
17R. Kessler, "Non-linear transition in three dimensional convection," J.
ACKNOWLEDGMENTS Fluid Mech. 174,357 (1987).
18R. Khosla and S. G. Rubin, "A diagonally dominant second order accurate
This work was begun during the postdoctoral study of implicit scheme," Comput. Fluids 2, 207 (1974).
NR and was supported by NASA Grant No. NAG-8-715-51. ION. Ramanan and S. A. Korpela, "Multi-grid solutions of natural convec-
NR acknowledges Fluid Dynamics International for provid- tion in a tall vertical cavity," Num. Heat Trans. A 15, 323 (1989).
ing computational facilities and GMH acknowledges contin- 2oMACSYMA, "A symbolic computation package," Macsyma Inc., Arlington,
Massachusetts, [993.
ued support of NASA Grant No. NAG-3-1475. 21 C. Paige and M. S. Saunders, "LSQR: Sparse linear equations of least
squares problems" ACM Trans. Math. Software 8, 195 (1982).
I B.Carpenter and G. M. Homsy, "High Marangoni number convection in a 22R. G. Storer, "Numerical studies of toroidal resistive magneto-
square cavity, Part II," Phys. Fluids A 2, 137 (1990). hydrodynamic instabilities," J. Comput. Phys. 66, 294 (1986).
Ip. Gillon and G. M. Homsy, "Combined buoyant-thermocapillary Flow: ~3K. N. Christodoulou and L. E. Scriven, "Finding leading modes of a
An experimental study" (private communication, 1989). viscous free surface flow: An unsymmetric generalized eigenvalue prob-
~O. R. Burggraf, "Analytical and numerical studies of the structure of lem," J. Sci. Com put. 3, 355 (1988).
steady separated flows," J. Fluid Mech. 24, 113 (1966). 24L. E. Erikkson and A. Rizzi, "Computer aided analysis of the convergence
4F. Pan and A. Acrivos, "Steady flows in rectangular cavities," J. Fluid to steady state discrete approximation to the Euler equations," J. Comput.
Mcch. 28, 643 (1967). Phys. 57, 90 (1985).
5U. Ghia, K. N. Ghia, and C. T. Shin, "High Re solutions for incompress- 25W. J. Stewart and A. Jennings, "A simultaneous iteration algorithm for
ible flows using Navier-Stokes equations and a multi-grid method," J. real matrices," ACM Trans. Math. Software 7, 184 (1981).
Comput. Phys. 48, 387 (1982). 26F. H. Busse, "Non-linear properties of thermal convection," Rep. Prog.
"R. Schreiber and H. B. Keller, "Driven cavity flows by effective numerical Phys. 41, 1929 (19715).
techniques," J. Comput. Phys. 49, 310 (1983). 27R. M. Gever and F. H. Busse, "Transition to time dependent convection,"
7J. R. Koseff and R. L. Street, "The lid-driven cavity flow: a synthesis of J. Fluid Mech. 65, 625 (1974).
qualitative and quantitative observations," J. Fluids Eng. 106,390 (1984). 2RW. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Nu-
8C. J. Freitas, A. N. Findikakis, and R. L. Street, "Numerical simulations of merical Recipes (Cambridge University Press, New York, 1986), p. 95.

Phys. Fluids, Vol. 6, No.8, August 1994 N. Ramanan and G. M. Homsy 2701