JOURNAL

OF COMPUTATIONAL

PHYSICS48, 387-411 (1982)

High-Re Solutions for Incompressible
Using the Navier-Stokes
Equations
Multigrid Method*

Flow
and a

U. GHIA, K. N. GHIA, AND C. T. SHIN
Universityof Cincinnati,

Cincinnati,

Received January

Ohio 45221

15, 1982

The vorticity-stream
function formulation of the two-dimensional incompressible NavierStokes equations is used to study the effectiveness of the coupled strongly implicit multigrid(CSI-MG)method in the determination of high-Re fine-mesh flow solutions. The driven flow
in a square cavity is used as the model problem. Solutions are obtained for configurations
with Reynolds number as high as 10.000 and meshes consisting of as many as 257 x 257
points. For Re = 1000, the (129 x 129) grid solution required 1.5 minutes of CPU time on the
AMDAHL470 V/6 computer. Because of the appearance of one or more secondary vortices
in the flow field, uniform mesh refinement was preferred to the use of one-dimensional gridclustering coordinate transformations.

1. INTRODUCTION
The past decade has witnessed a great deal of progress in the area ofcomputational fluid dynamics. Developments in computer technology hardware aswell as in advanced numerical algorithms have enabled attempts to be made towards
analysis and numerical solution of highly complex flow problems. For some of theseapplications, the use of simple iterative techniques to solve the Navier-Stokes
equations leads to a rather slow convergence rate for the solutions. The solution
convergence rate can be seriously affected if the coupling among the various
governing differential equations is not properly honored either in the interior of thesolution domain or at its boundaries. The rate of convergence is also generally
strongly dependent on such problem parameters as the Reynolds number, the meshsize, and the total number of computational points. This has led several researchers to
examine carefully the recently emerging multigrid (MG) technique as a useful meansfor enhancing the convergence rate of iterative numerical methods for solving
discretized equations at a number of computational grid points so large as to beconsidered impractical previously.

* This research was supported
Technical Monitor.

in part by AFOSR

Grant 80-0160, with Dr. James D. Wilson as

387
0021.9991/82/120387-25$02,00/O
Copyright
C 1982 by Academic
Press, Inc.
All rights of reproduction
in any form reserved.

388

GHIA, GHIA,

AND SHEN

The theoretical potential of the multigrid method has been adequately exposed forthe system of discretized equations arising from a single differential equation (e.g.,[ 15, 221). In fact, for 1-D problems, Merriam [ 151 has shown the likeness of themultigrid method to the direct solution procedure of cyclic reduction. This potential
has been realized and demonstrated in actual solutions of sample problems [3, 5,
8, 131. Application of the multigrid technique to the solution of a system of coupled
nonlinear differential equations still poses several questions, however, that are
currently being studied by various investigators [ 7, 21, 221.
The present study represents an effort to employ the multigrid method in thesolution of the Navier-Stokes equations for a model flow problem with a goal of
obtaining solutions for Reynolds numbers and mesh refinements as high as possible.
The fundamental principle of the multigrid procedure is first described briefly, then its
application to the governing equations is discussed in detail. Finally, the results
obtained for the shear-driven flow in a square cavity at Reynolds number as high as
5000 and 10,000 are presented, together with the particular details that needed to beobserved in obtaining these solutions.

2. BASIC PRINCIPLE OF MULTIGRID TECHNIQUE
Following Brandt and Dinar [7], the continuous differential problem considered isa system of I partial differential equations represented symbolically as
Lj CT(X)= Fj(~),

j=

1,2 ,..., 1, YED,

(2.1)

with the m boundary conditions
Bi o(2) = G,(Z),i = 1, 2,..., m,

X E aD,(2.2)

where 0 = (U, , U, ,..., U,) are the unknown variables, X = (x,, x2,..., xd) are the dindependent variables of the d-dimensional problem, Fj and Gi are known functions
on domain D and its boundary c?D, respectively, and Lj and Bi are general differential
operators.
A finite-difference solution to the problem described by Eqs. (2.1) and (2.2) isdesired in a computational domain with grid spacing h. With a superscript h todenote the finite-difference approximation, the linear system of algebraic equations
resulting from a selected difference scheme can be represented as(2.3)
A conventional iterative technique for solving Eq. (2.3) consists of repeated sweeps ofsome relaxation scheme, the simplest being the Gauss-Seidel scheme, untilconvergence is achieved. It is often experienced that the convergence of the method isfast only for the first few iterations. This phenomenon can be explained if one
considers a Fourier analysis of the error. Brandt (51 has thus estimated theHIGH RE INCOMPRESSIBLE

FLOW

389

magnitude of the smoothing rate p defined as the factor by which each error
component is decreased during one relaxation sweep of the Gauss-Seidel procedure.
It is observed that Gauss-Seidel relaxation produces a good smoothing rate for thoseerror components whose wave length is comparable to the size of the mesh; thesmoothing rate of more slowly varying Fourier components of the error is relatively
poor. The multigrid method is based primarily on this feature. It recognizes that awavelength which is long relative to a tine mesh is shorter relative to a coarser mesh.
Hence, after the first two or three iterations on a given fine mesh, the multigrid
method switches to a coarser mesh with step size 2h, where the error components
with wavelength comparable to 2h can be rapidly annihilated. The fine-grid solution
determined in the first step then needs to be corrected to reflect appropriately theremoval of the Z&wavelength content from the error. Repeated application of thisprocess over a sequence of grids constitutes the basic idea of the multigrid method.
Accordingly, the multigrid method makes use of a hierarchy of computational
grids Dk with the corresponding grid functions ok, k = 1, 2,..., M. The step size on Dk
is h, and hk+, = fhk, so that as k decreases, Dk becomes coarser. On the kth grid,
Eq. (2.1) has the discretized approximate form
LjkUk = Fj”.

(2.4)

The operations of transfer of functions from tine to coarse grids or from coarse tofine grids, have been termed “interpolations”
by Brandt 151. This terminology issomewhat unconventional when referring to transfer from fine to coarse grids. The
alternative terminology of restriction and prolongation, used by Hackbush [ 131 andWesseling [22], for example, is preferred here. The restriction operator Rip’ transfersa fine-grid function fik to a coarse-grid function ok-‘. On the other hand, theprolongation operator, denoted as Pi- I, transfers a coarse-grid function ok-’ to atine-grid function fl”.
For the restriction operator, the simplest possible form is “injection,” whereby thevalues of a function in the coarse grid are taken to be exactly the values at thecorresponding points of the next fine grid, i.e.,(Rt-lUk)i+l,j+l=“‘?i+I,2j+l*(2.5)

Being computationally efficient, injection has been used very frequently, particularlyin the initial phases of development of a multigrid program. In general, however, therestriction operator R :- ’ may be formulated as one of many possible weighted
averages of neighboring fine-grid function values. Two such operators are theoptimal-weighted averaging and the full-weighted averaging operators defined by
Brandt [6]. It is significant to note that these two are equivalent for 1-D problems.
For 2-D problems, optimal-weighted averaging involves fewer points than fullweighted averaging, which as the name indicates, involves all eight points (i f V,
j f- a), v, u = 0, 1, adjacent to a given point (i, j). Hence, optimal weighting may becomputationallymore efficient than full weighting, but the latter provides better

390

GHIA,

GHIA,

AND

SHIN

stability and convergence properties to a multigrid technique, particularlyfor
problems with rapidly varying coefficients. Full weighting is also preferred by
Wesseling [22] who termed it “9-point restriction” because of the number of points it
employs, i.e.,(Ri-‘Uk)i+

l,j+

1 =

auk+

l,*j+

I

+ QLU:i+*,*j+l + Ul;i+l,*j+2 + u:i,*j+l+ uti+l,*jl+ 3kl”ti+2,*j+2+ u:i,*,jl*+ Uti,*j+Z+ u:i+*.*j(2.6)

Wesseling also tested a 7-point modification of the above 9-point restriction operator
and found it to be almost equally suitable. The optimal-weighted averaging operatorof Brandt [6] is a 5-point restriction operator derivable from Eq. (2.6) by eliminating
the influence of the four corner-point values and doubling the center-point influence.
For the prolongation operator, the simplest form is derived using linear interpolation. This has been indicated by Brandt to be suitable for second-orderdifferential equations. Prolongation by linear interpolation introduces no ambiguity
when the interpolated value is desired at the midpoints of the boundaries of a meshcell. Two options are possible, however, for obtaining the interpolated value at thecenter of a cell. The choice of Rip ’ as defined by Eq. (2.6) suggests that the
prolongation operator Pi_, also involve nine points so that the value at the cell
center is obtained as the arithmetic mean of the four corner points. This leads to the
9-point prolongation operator defined by Wesseling as(Pk1Uk-‘)*i+2,2j+l=+[“ftt,j+l+ ufT:,j+l17(2.7)

1+ uf<:,j+Zl,(P~-,Uk-‘)*i+l,zj+*=~[U:,:,ji(Pf-lUk-1)2i+2,2j+2=

fIUfYll,j+l+

ufY:,j+l+

ufT:,,j+2+

uFT:,j+*l*

The operators defined in Eqs. (2.6) and (2.7) have the important
pfi-I is the adjoint of Rt-‘, i.e.,(R;-‘uk,

vk-’ )k-1 = (Uk,P;-,vk-‘)k,property that

(2.8)

where ( ) denotes an inner product defined as(2.9)
This properly maintains the coarse-grid equation to be a “homogenization”
of thefine-grid equation [2]. This is particularly important for nonlinear problems. In liniteelement discretization, this property leads to coarse-grid operators defined as
Lkp1=Rk-LLkpk
k

k-,

=

(P;-,)*Lk&.

(2.10)

HIGH RE INCOMPRESSIBLE FLOW

391

Brandt and Dinar [7] indicate that the choice of the restriction operators is guided
by rather definite rules, so that the only flexibility in a multigrid procedure is in theselection of the smoothing technique, i.e., the relaxation technique. While this may betrue to an extent, even the limited experience of the present authors with the multigrid
method indicates that, within the prescribed guidelines, some modifications in the
restriction and the prolongation operators do influence the efficiency of the overall
algorithm. Also, the definition of convergence in the finer grids appears to influence
the final solution obtained.
3. APPLICATION TO NAVIER-STOKESSHEAR-DRIVEN

EQUATIONS FORCAVITY FLOW

The laminar incompressible flow in a square cavity whose top wall moves with auniform velocity in its own plane has served over and over again as a model problem
for testing and evaluating numerical techniques, in spite of the singularities at two of
its corners. For moderately high values of the Reynolds number Re, published results
are available for this flow problem from a number of sources (e.g., [9, 17, 19]), using
a variety of solution procedures, including an attempt to extract analytically thecorner singularities from the dependent variables of the problem [IO]. Some results
are also available for high Re [ 161, but the accuracy of most of these high-Re
solutions has generally been viewed with some skepticism because of the size of thecomputational mesh employed and the difftculties experienced with convergence of
conventional iterative numerical methods for these cases.Possible exceptions to these
may be the results obtained by Benjamin and Denny (41 for Re = 10,000 using anonuniform 151 X 151 grid such that Ax = Ay ‘v l/400 near the walls and those of
Agarwal [ 1] for Re = 7500 using a uniform 121 x 121 grid together with a higherorder accurate upwind scheme. The computational time required in these studies,
however, is of the order of one hour or more for these high-Re solutions. The present
study aims to achieve these solutions in a computational time that is considerably
smaller, thereby rendering tine-mesh high-Re solutions more practical to obtain.
Governing D@erential Equations and Boundary Conditions

With the nomenclature shown in Fig. 1, the two-dimensional flow in the cavity canbe represented mathematically in terms of the stream function and the vorticity as
follows, with the advective terms expressed in conservation form:
Stream Function Equation:

yfxx + yyy + u = 0.

(3.1)

wx.x+ UYY- ReKvpL - (v,~>,l = Rew,.

(3.2)

Vorticity Transport Equation:

Boundary conditions.
The zero-slip condition at the nonporous walls yields that w
and its normal derivatives vanish at all the boundaries. As is well known, this392

GHIA, GHIA, AND SHIN

FIG. 1. Cavity flow configuration, coordinates, nomenclature, and boundary conditions.

provides no direct condition for o at the walls. Theoretically, this should not pose adifficulty if the equations for u and v are solved simultaneously and if all boundary
conditions are imposed implicitly. In practice, however, the boundary conditions foro are derived from the physical boundary conditions together with the definition of w
as given by Eq. (3.1).
Thus, at the moving wall y = 1,j = J:
(3.3)

vJ=o,

*J = -‘(/yy = -(vJ+

1-

~VJ

+

WJ-

,)/Ay2,(3.4)

where wJ+, is evaluated from a third-order accurate finite-difference expression forw,,, which is a known quantity at the boundary, i.e.,WY,=c2y/,+

1 +3~J-661//~-l+

VJ-I)/(~AY).

(3.5)

HIGH

RE INCOMPRESSIBLE

FLOW

393

The resulting expression (3.4) for wJ is then second-order accurate (see 1121).
Expressions for w at other boundaries are obtained in an analogous manner.
Discretization
The discretization is performed on a uniform mesh. In fact, with a multigrid
solution technique, nonuniform mesh or grid-clustering coordinate transformations
are not essential since local mesh refinement may be achieved by defining
progressively liner grids in designated subdomains of the computational region.
Second-order accurate central finite-diferrence approximations are employed for all
second-order derivatives in Eqs. (3.1) and (3.2). The convective terms in Eq. (3.2) arerepresented via a first-order accurate upwind difference scheme including its secondorder accurate term as a deferred correction, as formally suggestedby Khosla andRubin [ 14). This ensures diagonal dominance for the resulting algebraic equations,thus lending the necessary stability property to the evolving solutions while restoring
second-order accuracy at convergence.
Relaxation Scheme (Smoothing Operator)
In the multigrid method, the role of the iterative relaxation schemeis not so muchto reduce the error as to smooth it, i.e., to eliminate the high-frequency error
components. Due to the coupling between governing equations (3.1) and (3.2) as wellas through the vorticity boundary conditions (Eq. (3.4)), sequential relaxation of theindividual equations (3.1) and (3.2) would have poor smoothing rate. For example,
smooth errors in w could produce high-frequency error components in the vorticity
solution via the boundary condition for w. On the other hand, a convergent solutionof each equation at each step would constitute a very inefficient procedure. An
appropriate approach consists of relaxing the coupled governing equations (3.1) and(3.2) simultaneously and incorporating the vorticity boundary conditions implicitly.
A coupled Gauss-Seidel procedure or a coupled alternating-directing implicit scheme
may be used for this purpose. Rapid convergence, however, of the coarsest-grid
solution as required in the full multigrid algorithm [6] can be safely assured by theuse of these methods only when the coarsest mesh is not too line. On the other hand,for the driven-cavity flow at high Re, too coarse a grid does not retain enough of thesolution features and cannot, therefore, provide an appropriate initial approximation
to the fine-grid solution for high-Re flow. Hence, the present work employs thecoupled strongly implicit (CSI) procedure of Rubin and Khosla [ 181.This scheme isa two-equation extension of the strongly implicit procedure developed by Stone 1201
for a scalar elliptic equation and may be viewed as a generalization of the Thomasalgorithm to two-dimensional implicit solutions. Its effectiveness has been
demonstrated by Rubin and Khosla [ 181 via application to a number of flow
problems. The present authors have also found it to be useful in conjunction with themultigrid technique [ 111.The procedure may be approximately likened to incomplete
lower-upper (ILU) decomposition which is considerably more efficient, manifesting
lower values of the smoothing factor p than the simple Gauss-Seidel relaxation
procedure.

394

GHIA,

Prolongation

GHIA,

AND

SHIN

and Restriction Operators

The prolongation operator was always chosen to be the 9-point operator given in
Eq. (2.7) except for converged coarse-grid solutions where cubic interpolations wereused as suggestedby Brandt [6]. For the restriction operator, simple injection as wellas the 5-point operator (optima1 weighting) and the 9-point operator (Eq. (2.6)) were
employed. While the first two generally provided convergent solutions, 9-point
restriction led to improved convergence for the very high Re cases computed.
Multigrid

Procedure

For the present nonlinear problem, the full approximation scheme (FAS) was
employed, rather than a correction scheme.Also, the full multigrid (FMG) algorithm
was preferred over the cycling algorithm since a converged coarse-grid solution is
generally obtainable by the CSI procedure used for relaxation. It is possible that, forthe higher-Re cases computed, the cycling algorithm could also be used with the first
approximation of the finest-grid solution being provided by the solution of apreceding calculation with a lower value of Re. Finally, the “accommodative” versionof the multigrid procedure was used so that convergence as well as convergence rate
were monitored during the process of relaxation on a given grid in order to control
the computational procedure, particularly with respect to switching from one grid toanother. The accommodative FAS-FMG procedure used here follows that detailed
by Brandt [6]. This procedure is briefly described below.
The solution on grid Dk is denoted as uk. This is prolongated to the next finer grid
Dkt ’ using the prolongation operator to provide an estimate for ukt’ asU k+’est -pt+‘uk*
-

(3.6)

This estimate is used as the initial guessfor the solution on grid Dkt I, i.e., for solving
the equation
Lk+lUktl
=Fktl
(3.7)
Convergence is defined to occur when the norm ekt, of the dynamic residuals of
Eq. (3.7) is below a specified tolerance, sk+, , i.e.,ek+l<‘ktl?

(3-g)

with ek+, taken to be 10m4 in the present work. It must be recognized that thisassigned value for ek+, remains in effect only until the first multigrid cycle is
executed, during which Ed+I is redefined as described later in Eqs. (3.13) and (3.14).
During the relaxation process, the convergence rate for Eq. (3.7) is also monitored
and compared with the theoretical smoothing rate of the relaxation procedure used.
Accordingly, if the ratio of the residuals eI+ I and et:: for two consecutive iterations
n and (n + 1) is smaller than the smoothing rate of the scheme,the iteration process
is continued. Following convergence, k is incremented by unity and the steps
indicated by Eqs. (3.6) and (3.7) are repeated. This is continued until k = M, i.e.,
convergenceon the desired finest grid is achieved, yielding the required final solution.

HIGH RE INCOMPRESSIBLE FLOW

395

If, at any stage of the relaxation process for k # 1, i.e., for all but the coarsest grid,
the convergence rate is not satisfactory, i.e., if
flt1
ektl

>

let+,

(3‘9)

v,

where q = ,LJ,the scheme smoothing rate, then a multigrid cycle is interjected in theprocedure. The multigrid cycle consists of computing a coarse-grid correction u:+ , tothe evolving line-grid solution u$,’ by solving the equation
LkUkk+lwhere(3.10)-P,-

fk & Lk(Rkkt,~;&‘)+R;+l(fk+l-Lk+‘u;&‘).

(3.11)If (k + 1) is currently the finest level, then f””= Fkt ’ . This correction is used to
improve the old line-grid solution according to the relationu k+’
new= z&l+ P;+l(u;t1 -R;,,u::,‘).

(3.12)

Convergence of Eq. (3.10) is defined to occur when the residual norm ek for thisequation is smaller than the residual norm ek+l for the finer grid, i.e., whenek < ck= de,+,(3.13)

where 6 < 1; the value used was 6 = 0.2.
Following the correction according to Eq. (3.12), the solution of Eq. (3.7)
proceeds as before.
If the solution of Eq. (3.10) itself does not exhibit a satisfactory convergence rate,
defined in a manner analogous to Eq. (3.9), then a second multigrid cycle may be
performed by going to a yet coarser grid Dk-’ to enhance the convergence rate ofEq. (3.10). Thus, a sequence of multigrid cycles may be nested, one inside another, tosolve the current finest-grid equation effkiently. On the currently finest grid Dkt ’ inthis nest, convergence should be attained to within the estimated truncation error rk
so that, corresponding to Eq. (3.8) the convergence criterion used isektl

cEktl=aqrkT(3.14)

where Tk is the norm of (Fk -fk),
a = (hk+ /hk)‘, and q = 1 for second-orderaccurate discretization. As will be discussed in the next section, modification of
Eq. (3.14) to include aq, with q > 1, appears to influence the converged numericalvalues of the solution.

4. FINE-GRID

AND HIGH-RE

RESULTS FOR DRIVEN CAVITY

The correctness of the analysis, the solution procedure, and the computer program
were assessed by first obtaining fine-mesh solutions for the case with Re = 100 for
which ample reliable results are available in the literature. This case was also
intended for experimentation with some of the parameters associated with the

a

“9.c.
-0.4

-0.2

0.0

0.0

0.2

0.4

0.6

0.8

1.0

1.0
0.8
I

PROFiLES THROUGHGEOMETRIC CENTER

0.6
0.4
0.2

0.0

- 0.8

Y

. 0.6
PRIMARY-VORTEX(Locations
Listed-0.6
FIG. 2a.
center.

b-0.4Comparison

-0.2

i1::

0.0

of u-velocity

0.0

0.2

CENTER (p.v.c)

inof p.v.c.
Table
5)

0.4

0.6

- 0.4

0.8

1.0

-

along vertical lines through geometric center and primary vortex

0.5 ,PRIMARYVORTEX

0.0

0.0

-0.1

-0.1

-0.3

0.0

0.2

0.4

0.6

0.8

I
0.0

0.2

x

FIG. 2b. Comparison
primary vortex center.

of profiles of tr-velocity

along horizontal

0.4

0.8
1
0.6
x
lines through geometric center andRE INCOMPRESSIBLE

HIGH

397

FLOW

multigrid procedure, namely, selection of q in Eq. (3.9), values qr and q, for q in
Eq. (3.14), the coarsest mesh width, the finest mesh width h,, and the prolongation
and restriction operators. It was observed that h, is the most important parameter,
especially for high Re. Also, as Re increased, very coarse grids could not be includedin the procedure. The smoothing factor ,u for the CSI procedure used is expected tobe smaller than that for the Gauss-Seidel scheme. Nevertheless, for Re = 100 and400, r = 0.5 was used. The value of q had to be increased gradually with Re; for
Re = 10,000, v = 0.7 was needed. Similarly, a time step of infinity could be used in
the vorticity equation for Re up to 3200 but had to be reduced rather rapidly as Re
increased. For Re = 10,000, At = 0.1 was required. The corresponding values At,
used by Benjamin and Denny [4] in conjunction with the AD1 solution procedure forthe case with Re = 10,000 were smaller than this by several orders of magnitude; thismay also partly explain the reduction in computational time achieved by the present
solution technique. Initially, q = 1 was used in Eq. (3.14). For Re = 100, this proved
adequate in the respect that the results agreed well with available solutions. But for
Re > 1000, the values obtained, for instance, for 1vmin/ at the center of the primaryvortex were somewhat lower than the published solutions. Better comparison resulted
from the use of q > 1. This is because increasing q enforces continuation of theiteration process and leads to some further reduction of the dynamic residuals ek+, ,while also modifying the actual solution simultaneously. The results published bymost previous investigators have been obtained subject to the convergence criterion
that the relative change in two successive iterates of the solution at eachcomputational point be below a prescribed small value. Frequently, the choice of thisvalue is not related to the truncation error in the finite-difference approximation. On
the other hand, in the present computational procedure, convergenceon the finest gridis defined in terms of the truncation error. Hence, the results presentedhere employedq, = 4 and q,,,= 5, except for the cases with Re = 7500 and 10,000 which used
q* = 4.

Figures 2a and b show the velocity profiles for u along vertical lines and c’ along

LEGEND

Present
Rubin
Khosla

100

400

1000

3200

5000

7500

10000

a

b

c

d

e

f

9

-.-

-

-

and
[I9771A

A

Nallasamy&
Prasad[I97710aAgarwal119811

j0

0

0

-

-

__

0

0

0

398

GHIA, GHIA,

AND SHIN

horizontal lines passing through the geometric center of the cavity and through thecenter of the primary vortex in the flow. The origins for these graphs for various
values of Re have been displaced for clarity of the profiles. The thinning of the wallboundary layers with increase in Re is evident from these profiles, although the rateof this thinning is very slow for Re > 5000. The near-linearity of these velocity
profiles in the central core of the cavity is indicative of the uniform vorticity region
that develops here for large Re. The high-Re profiles of u exhibit a kink near y = 1,
while a similar behavior is observed for the u profiles near x = 1. Such behavior has
been reported by some previous investigators, and is seen to persist in the present
line-grid solutions. This would imply that the velocity distributions near these walls
are not extremely sensitive to mesh size. The values of I// and w at the vortex centerare observed to be considerably more sensitive to mesh size than these velocity
profiles.
Included in Fig. 2 are the available results of [ 1, 16, 171. For Re = 100, all results
agree well with one another as well as with the present solutions, indicating that forthis value of Re, the coarser grids employed by the previous investigators were quite
adequate. As Re increases, however, the inadequacy of coarse meshes gradually
becomesapparent. This is particularly evident in the solutions reported by Nallasamy
TABLE
Results for u-velocity

I

along Vertical Line through Geometric Center of Cavity
Re

129.
grid
pt. no.

y129

1.ooooo

1.ooooo

1.ooooo

1.ooooo

1.ooooo

1.ooooo

1.ooooo

I .ooooo

126

0.9766

0.841230.758370.659280.532360.48223

0.472440.472211250.96880.7887 1

0.684390.57492

0.48296

0.461200.470480.47783124

0.96090.737220.61756

0.511170.465470.45992

0.473230.480701230.953 10.687170.55892

0.46604

0.461010.460360.47167

0.47804110

0.85160.231510.290930.33304

0.34682

0.33556

0.34228

0.34635950.73440.00332

0.16256

0.187190.19791

0.20087

0.2059 I

0.2067380

0.6172

AI136410.02135

0.057020.07156

0.08 183

0.083420.0834465
590.50000.453 1-0.2058 1-0.21090-0.06080-0.10648-0.04272-0.86636

a.03039-0.07404

0.03 111-0.07540-0.27805

a.24427-0.22855

XJ.03800
AI.07503
AI.23 176

4.38289-0.29730-0.33050-0.32393-0.32709-0.20196
AI.18109-0.35344-0.32407-0.40435-0.43643
a.42901-0.41165
0.00000

AI.38324-0.43025
a.43590
AI.43154
0.00000-0.38000

a.22220

a.34323a.41933-0.37827

100

400

31

0.2813AI.15662-0.11477-0.17119-0.3272623
140.1719Al.10150-0.24299

0.1016a.06434-0.146121098
10.07030.06250.0547-0.04775
a.04192-0.03717
0.00000-0.10338
AI.09266-0.08 186
0.00000

0.00001000

0.00000

3200

5000

7500

10,000-0.23 186

-0.41657-0.42537-0.42735
0.00000

HIGH RE INCOMPRESSIBLE

399

FLOW

and Prasad [ 161. Nevertheless, the fourth-order accurate spline method of Rubin and
Khosla [ 171 remains satisfactory with a 17 x 17 mesh at Re = 1000. Also, the thirdorder accurate scheme of Agarwal [ 1] performs well with a 121 X 121 grid at
Re = 7500, but the corresponding computer time is quite large. Unfortunately,
Benjamin and Denny (41 did not present any velocity data, although their solutions
are considered to be very accurate for high Re.
In view of the above remarks, the present line-mesh results should be very useful.
Consequently, Tables I and II list the numerical values corresponding to the velocity
profiles shown in Fig. 2 for lines passing through the geometric center of the cavity.
Only typical points, rather than the entire large set of computational points, along
these profiles have been listed. Care has been taken to include the points of local
maxima and minima for all values of Re; these points are underscored.
The streamline contours for the cavity flow configurations with Re increasing from100 to 10,000 are shown in Fig. 3. A magnified view of the various secondary
vortices is also included. The values of I,V along the contours shown are listed inTable III. For Re = 400, the results from a 129 X 129 grid as well as a 257 x 257
grid are presented in order to demonstrate that the 129 x 129 grid is adequate formoderate values of Re. Although a comparison is not shown in this figure, the extent
TABLE
Results for u-Velocity

along Horizontal

II

Line through Geometric Center of Cavity

129m
grid
pt. no.

Re
x

100

400

1000

3200

5000

129

1.oooo0.00000

0.00000

0.00000

0.00000

0.00000

0.00000

0.00000125

0.9688-0.05906

a.12146

a.21388

xJ.39017-0.49774

a.53858-0.54302

124

0.9609-0.0739 1-0.15663-0.27669

a.47425a.55069-0.55216-0.52987

123

0.953 1-0.08864-0.19254

a.33714

a.52357-0.55408-0.52347-0.49099

122

0.9453-0.10313-0.22847-0.39188-0.54053

XI.52876

a.48590-0.45863

117

0.9063-0.16914-0.23827

a.5 1550

a.44307a.41442

Jx41050-0.4 1496

111

0.8594-0.22445

a.44993

a.42665-0.37401-0.36214

a.36213-0.36737

104

0.8047-0.24533-0.38598-0.3 1966-0.31184

X).30018

a.30448-0.307 19650.50000.054540.05 188

0.02526

0.009990.009450.008240.0083 1310.23440.175270.30174

0.32235

0.28 188

0.27280

0.273480.27224300.22660.302030.330750.290300.28066

0.28 1170.2800321

0.15630.175070.16077

0.28124

0.35368

0.09380.078 10.07030.06250.123170.108900.100910.229650.20920
0.197130.37095
0.32627

0.37119

13
11
10
9

0.303530.290120.42768
0.4 1906
0.409 170.4295 1
0.436480.433290.35060
0.418240.435640.440300.350700.41487
0.43 124
0.43733

0.09233

0.18360

0.27485

0.39560

0.424470.43979

0.4398310.00000.00000

0.00000

0.00000

0.00000

0.00000

0.00000

0.00000

7500

10,000400

GHIA, GHIA, AND SHIN

RE-100.

RE=400,UNIFORMGRID

UNIFORM

GRID(129x129)(129x125

RE=400,RE = 1000,

Eddy

UNIFORMGRID

UNIFORM

GRID(129~129)

EL,

.d

x

FIG. 3.

0.2

0.4

Streamline pattern for primary,

0.6*

secondary, and additional

0.8

1.0

corner vortices.

(257x257)

HIGH

Eddy

KE

RE = 3200,

TL,

EddiesBL,.

INCOMPRESSIBLE

FLOW

UNIFORM

GRID

(129x

129)

UNIFORM

GRID (257x 257)

BL2

Y

0.6

Y

0.2

0.0

Eddy

0.0

TL,

x

RE = 5000,

0.2

1.0

Y

0.8

0.6

FIGURE

3 (conhued)

401

GHIA, GHIA, AND SHIN

402

Eddy

TL,

RE = 7500,

UNIFORM

GRID

(257x

257)

RE = 10000,

UNIFORM

GRID

(257x

257)

0.6

Y

Y

FIGURE 3 (concluded)

403

HIGH RE INCOMPRESSIBLE FLOWTABLE III
Values for Streamline and Vorticity Contours in Figs. 3 and 4
Stream function

Vorticity

ContourletterValue of (I/numberValue of v

kIm

-1.0 x lo- IL7-1.0 x lo-’-1.0 x loms
-1.0 x 1om4
-0.0 100-0.0300-0.0500-0.0700-0.0900-0.1000-0.1100-0.1150-0.1175

0
1
2
3
4
5
6
7
8
9
10

1.0 x 10-n
1.0 x 10-7
1.0 x lo-”
1.0 x 10. 5
5.0x lo-5
1.0 x 1om4
2.5 x 10mJ
5.0 x 10-J
1.0 x 10-j
1.5 x lo-’
3.0 x 10-j

ContourContournumberValue of w

0
fl
+2
+3
f4
5
6

0.0
*to.5
*1.0
zt2.0
f3.0
4.0
5.0

of the various secondary vortices is in excellent agreement with that reported by
Benjamin and Denny [4]. The present results, however, are computationally more
efficient.
In Fig. 4 we show the vorticity contours corresponding to the streamline patterns
presented in Fig. 3. Again, the values of w along these contours are listed inTable III. As Re increases, several regions of high vorticity gradients, indicated byconcentration of the vorticity contours, appear within the cavity. It is seen from
Fig. 4 that these regions are not aligned with the geometric boundaries of the cavity.
It is for these reasons that uniform mesh refinement was used in the present study.
Possible suitable alternatives appear to be the use of a basically modified nonCartesian coordinate system and of a solution-adaptive local mesh refinement. An
often-compared quantity for cavity flows is the vorticity at the midpoint of themoving wall or the minimum value of o at this boundary. Hence, the values of w atseveral selected points along this boundary are listed in Table IV, with the minimum
value indicated by the underscore. These values of mrnin agree very well with theresults tabulated in [4].
As seen from Figs. 3 and 4, fine-mesh solutions exhibit additional counter-rotating
vortices in or near the cavity corners as Re increases. The effect of Re on the locationof the centers of these vortices is shown in Fig. 5. In terms of the notation shown in
Fig. 1, the letters T, B, L, and R denote top, bottom, left, and right, respectively; thesubscript numeral denotes the hierarchy of these secondary vortices. Thus, BR, refers
to the second in the sequence of secondary vortices that occur in the bottom rightcorner of the cavity. As is well known, the center of the primary vortex is offset

GHIA, GHIA, AND SHIN
RE = 100,

UNIFORM

GRID(129x129)

RE = 400,

UNIFORM

GRID(129x129)

RE = 400,

UNIFORM

GRID(257x257)

RE = 1000,

UNIFORM

GRID(129x129)

FIG. 4.

Vorticity

contours for flow in driven cavity.

towards the top right corner at Re = 100. It moves towards the geometric center ofthe cavity with increase in Re. Its location becomes virtually invariant for Re > 5000.
All the secondary vortices appear initially very near the corners (or near the wall, inthe case of the vortex TL,) and their centers also move, though very slowly, towards
the cavity center with increase in Re. At the larger values of Re considered, theconvection of these secondary eddies is evidenced by the direction of movement of
the centers of these vortices.
The computational advantage gained by use of the MG procedure is best illustrated
in terms of the behavior of the root-mean square (RMS) value of the dynamic
residuals of the discretized governing equations in the finest grid. In Fig. 6 we show
the finest-grid RMS residuals for w and o obtained during a single-grid computation
with h = & (solid curve) as well as a multigrid calculation with h, = & and M = 6

HIGH RE INCOMPRESSIBLE

405

FLOW

RE = 3200,

UNIFORM

GRID(129x129)

RE = 5000,

UNIFORM

GRID

RE = 7500,

UNIFORM

GRID(257x257)

RE= 10000,

UNIFORlVl GRID

(257x257)(257 x 257)

FIGURE 4 (continued)(solid and dashed lines). Flow configurations with Re = 100 and Re = 1000 havebeen examined. In both cases, even the single-grid calculations exhibit a rapid initial
decay of the RMS residuals for y as well as w during the first 4-6 iterations (workunits). Thereafter, the solid curves show a marked decrease in their slope. Employing
the multigrid process after these first 4-6 work units tends to retain the initial decayrate for the errors during the overall computation.
It is important to mention two points with respect to the MG curves in Fig. 6.
First, the solid portions of the MG-curves correspond to the relaxation step(smoothing) on the finest grid while the dashed portions correspond to the coarse-gridcorrection due to the MG cycle. Second, although convergence was defined on thebasis of the arithmetic average of the RMS residual in u and u/, the convergence rate
was examined in terms of the RMS residual in o alone. It is perhaps for this reason

TABLE
Results for VorticityIVw along Moving Boundary
Re

100

X

10,0001000

3200

75.5980
51.0557
40.5437
32.2953
25.4341
20.2666
16.8350
14.8901
14.0928
14.1374
14.8061
16.0458
18.312023.870742.1124-

126.670
89.3391
75.640161.786447.144335.8795
28.941325.3889
24.1457
24.463925.8572
27.95 1430.4779
34.232749.96640.00004O.OllO
22.5378
16.2862
12.7844
10.41998.69628
7.43218
6.5745 1
6.13973
6.18946
6.82674
8.2211010.741415.6591
30.7923

0.06250.12500.18750.25000.31250.37500.43750.50000.56250.62500.68750.75000.81250.87500.9375
1.oooo53.6863
34.635 1
26.5825
21.0985
16.8900
13.7040
11.4537
10.05459.388899.345999.88979
11.2018
13.9068
19.685935.0773

Vortex146.702103.436
91.5682
77.950960.0065
45.862237.360933.0115
3 1.3793
31.5791
33.048635.3504
38.0436
41.339456.7091-

209.452
145.073
127.928116.275
90.023 1
67.1400
53.5905
46.827144.3287
44.630346.8672
50.3792
54.3725
57.7756
66.0352

180.927125.131111.115
98.2364
75.633456.934545.9128
40.3982
38.3834
38.695 140.612343.564146.890150.0769
61.4046

TL,

0.92
Y

Primary

Vortex

0.91

0.70

0.90
ow891

Re=:200

,

Y

1

0.65
0.05

0.06

0.07

0.08

x
0.60

.50

0.55

0.60

0.65

x

Vortex

BL,

0.20

.
0

Vortex
Vortex

0.151

1

I

0.75

0.80

0.85

BR,
BR
2
I

0.90

I

1

0.95

1.00

x

FIG. 5.

Effect of Reynolds number on location of vortex centers.

HIGH

RE INCOMPRESSIBLE

FIG.
6. Convergence of single grid and multigrid
Re= 100: (O)e,, (W)e,; Re= 1000: (A)e,, (*)e,.
De ,; Re= 1000: (A)e,, (O)e,.

FLOW

407

computational procedures. Single grid (h = &)Multigrid (h,, = &, A4 = 6) Re = 100: (0) e,.

that w exhibits a much more desirable convergence behavior than v because theconvergence rate is indeed the parameter that comprises the basis for interjecting an
MG cycle in the solution procedure. Some further improvement in the overall
convergence process may be possible by also including the convergence rate of v/ inthe criterion controlling switching to the coarse-grid correction step.
Finally, a comprehensive survey of the properties of the primary and secondary
vortices in the driven-cavity flow is provided in Table V. Some of these are directlycomparable with the numerical data listed in [ 1, 41. In particular, attention is drawn
to the values of vmin and cc),,,, for the primary vortex. The present calculations for
Re = 7500 with a 257 X 257 grid exhibit a stronger secondary vortex BR, thanreported by Agarwal [ 11. Consequently, the present primary vortex is somewhat
weakened. Nevertheless, the approach of o,.~, to the infinite-Re value of 1.886 isclear, although this value is approached “from below” for the present solutions.
581/48/3-7

408

GHIA,

GHIA,

AND

SHIN

TABL’E
Properties of Primar:yNumberPrimary

Property
Vmin

w V.E.
Location, x, y
First
T

Wmax

100
-0.103423
3.16646
0.6172, 0.7344

-

400

1000

-0.113909
2.29469
0.5547, 0.6055-0. I 17929
2.04968
0.5313, 0.5625-

W V.C.

--

1.74877 x lo-6
-1.55509 x lo-’0.03 13, 0.039 1
0.078 10.078 11.41951 x lomJ
-5.69697 x lo-’0.0508, 0.0469
0.12730.10812.31129 x 10 -’-0.36175
0.0859, 0.078 I
0.21880.16801.25374 x lo-’
-3.30749 x 10-l
0.9453, 0.0625
0.13280.14846.42352 x IO-’
-4.33519 x 10-l
0.8906, 0.1250
0.26170.32031.75 102 x 10 mi-1.15465
0.8594, 0: 1094
0.30340.3536‘Ymin
wv.,.
Location, x, y

-

4.
VI

-

-7.67738 x IO-“’
9.18377 x IO-’
0.0039, 0.0039
0.00390.0039Location, x, y
4.
VL
BL

vmax
W Y.C.

Location, x, y
HL
VI.
BR

Wmax
W V.C.

Location, x, y4
VL
SecondBLBRvlnli”W V.E.

Location, x, y
HI.
VL
Third
BR

Wmax
Location, x, y
HLVIWork units
CPU seconds
Mesh points--1.86595 x lo-@
4.38726 x 10-j
0.9922, 0.0078
0.01560.0156-9.31929 x IO-’
8.52782 x 10-j
0.9922, 0.0078
0.00780.0078

18.84
53.59
129

18.08
215.05
257

31.56
92.27
129

409

HIGH RE INCOMPRESSIBLE FLOWV
and Secondary Vortices

3200

5000

7500

10,000-0.120377
1.98860
0.5 165, 0.5469-0.118966
1.86016
0.5117, 0.5352-0.119976
1.87987
0.5117, 0.5322-0.119731
1.88082
0.5117, 0.5333

7.27682 x lo-”
-1.71161
0.0547, 0.8984
0.08590.20571.45641 x 10-l
-2.08843
0.0625, 0.9 102
0.12110.26932.04620 x lo-’
-2.15507
0.0664, 0.9141
0.14450.29932.42103 x 10 ’
-2.18276
0.0703, 0.9 14 1
0.15890.32039.7823 x 10m4
‘1.06301
0.0859, 0.1094
0.28440.23051.36119 x lo-1.53055
0.0703, 0.1367
0.3 1840.2643

1.46709 x 10-l
-1.78511
0.0645, 0.1504
0.33390.2793

1.51829 x 10 ’
-2.08560
0.0586, 0.1641
0.34380.2891

3.13955 x 10-l
-2.27365
0.8125, 0.0859
0.3406
0.4 1023.08358 x lo-’
-2.66354
0.8086, 0.0742
0.35650.41803.28484 x 10-j
-3.49312
0.7813, 0.0625
0.37790.43753.41831 x lo-’
-4.05 3 1
0.7656. 0.0586
0.39060.4492-6.33001 x lo-’1.44550 x 10-2
0.0078, 0.0078
0.00780.0078-7.08860 x lOm8
1.88395 x lo-*0.0117, 0.0078
0.01560.0163-1.83167 x 10-l
1.72980 x 10. 2
0.0117, 0.0117
0.02340.0254-7.75652 x lo-’
2.75450 x 10 -’0.0156, 0.0195
0.03520.044 1-2.51648 x lo9.74230 x IO-’
0.9844, 0.0078
0.02540.0234-1.43226 x 10m6
3.19311 x lo-*0.9805, 0.0195
0.05280.0417-3.28148 x 10.’1.41058 x 10-l
0.9492, 0.0430
0.12700.0938-1.31321 x 10 -’3.12583 x 10-l
0.9336, 0.0625
0.17060.1367-

1.58111 x lO-9
0.9961, 0.0039
0.00390.00395.66830 x 10 my0.9961, 0.0039
0.00390.003970.8125734.49257

68.50
705.62257

99.5
986.65257

78.25
207.26129

410

GHIA, GHIA, AND SHIN

SUMMARY

Fine-mesh solutions have been obtained very efficiently for high-Re flow using thecoupled strongly implicit and multigrid methods. The various operators andparameters in the multigrid procedure were examined, especially for high-Re flow.
The use of 9-point restriction, or full-weighting, was found to be superior to 5-point
restriction, or optimal weighting. The finest mesh size employed in the grid sequence
continues to be a very significant parameter. The smoothing factor of the iteration
scheme was seen to be influenced by the physical problem parameters, namely, Re.
The definition used for convergence on current tine grids was also observed toinfluence the final solutions.
The robustness and the efficiency of the overall solution technique has been
demonstrated using the model problem of flow in a driven square cavity. Detailed
accurate results have been presented for this problem. Up to 257 x 257, i.e., 66049
computational points and Re as high as 10,000 have been considered, with CPU timeof 16 to 20 minutes on the AMDAHL 470 V/6 computer. The present results agree
well with published fine-grid solutions but are about four times as efficient.
Future effort includes consideration of primitive-variable formulation; true 3-D
solutions may be possible in the foreseeable future, with practical CPU timerequirements, by use of multigrid techniques.
REFERENCES
1. R. K. AGARWAL,“A Third-Order-Accurate Upwind Scheme for Navier-Stokes Solutions at HighReynolds Numbers,” AIAA Paper No. 8 l-01 12, 1981.
2. 1. BABLJSKA,in “Numerical Solutions of Partial Differential Equations III” (B. Hubbard, Ed.),
Academic Press, New York, 1975.
3. N. S. BAKHVALOV,USSR Comput. Math. Phys. 6(5) (1966), 101.
4. A. S. BENJAMINAND V. E. DENNY, J. Comput. Phys. 33 (1979), 340.
5. A. BRANDT,Math. Comput. 31 (1977), 333.
6. A. BRANDT,“Multi-Level Adaptive Computations in Fluid Dynamics,” AIAA Paper No. 79-1455,
1979.
7. A. BRANDT AND N. DINAR, in “Numerical Methods for Partial Differential Equations” (S. Parter,Ed.), Academic Press, New York, 1979.
8. R. P. FEDORENKO,USSR Comput. Math. Phys. 1 (1962), 1092.
9. K. N. GHIA, W. L. HANKEY, AND J. K. HODGE,“Study of Incompressible Navier-Stokes Equationsin Primitive Variables Using Implicit Numerical Technique,” AIAA Paper No. 77-648, 1977; AIAA
J. 17(3) (1979), 298.

10. K. N. GHIA, C. T. SHIN, AND U. GHIA, “Use of Spline Approximations for Higher-Order Accurate
Solutions of Navier-Stokes Equations in Primitive Variables,” AIAA Paper No. 79-1467, 1979.
1I. K. N. GHIA, U. GHIA, C. T. SHIN, AND D. R. REDDY, in “Computers in Flow Predictions and Fluid
Dynamics Experiments” (K. Ghia et al., Eds.) ASME Publication, New York, 1981.
12. U. GHIA, K. N. GHIA, S. G. RUBIN AND P. K. KHOSLA, Comput. Fluids 9 (198 l), 123.
13. W. HACKBUSCH,Computing 20 (1978), 291.
14. P. K. KHOSLA AND S. G. RUBIN, Comput. Fluids 2 (!974), 207.
15. M. L. MERRIAM, “Formal Analysis of Multi-Grid Techniques Applied to Poisson’s Equations inThree Dimensions,” AIAA Paper No. 81-1028, 1981.

HIGH RE INCOMPRESSIBLE

FLOW

411

16. M. NALLASAMYAND K. K. PRASAD,J. Fluid Mech. 79(2) (1977), 391.
17. S. G. RUBIN AND P. K. KHOSLA,J. Comput. Phys. 24(3) (1977) 217.
18. S. G. RUBIN AND P. K. KHOSLA, Comput. Fluids 9 (1981), 163.
19. R. E. SMITH AND A. KIDD, “Comparative Study of Two Numerical Techniques for the Solution of
Viscous Flow in a Driven Cavity,” pp. 61-82, NASA SP-378, 1975.
20. H. L. STONE,SIAM J. Numer. Anal. 5(3) (1968), 530.
21. T. THUNELL AND L. FUCHS,in “Proceedings of Symposium on Numerical Methods in Laminar andTurbulent Flow” (C. Taylor, K. Morgan, and B. A. Schrefler, Eds.), pp. 141-152, Venice, Italy,
1981.
22. P. WESSELING,“Theoretical and Practical Aspects of a Multi-Grid Method,” Report NA-37, Delft
University of Technology, The Netherlands, 1980.