Numerical calculations of the 2D steady incompressible driven cavity flow are presented. The NavierStokes equations in streamfunction and vorticity formulation are solved numerically using a fine uniform grid mesh of 601×601. The steady driven cavity solutions are computed for Re £ 21,000 with a maximum absolute residuals of the governing equations that were less than 10^{10}. A new quaternary vortex at the bottom left corner and a new tertiary vortex at the top left corner of the cavity are observed in the flow field as the Reynolds number increases. Detailed results are presented and comparisons are made with benchmark solutions found in the literature.
Numerical methods for 2D steady incompressible NavierStokes (NS) equations are often tested for code validation, on a very well known benchmark problem; the liddriven cavity flow. Due to the simplicity of the cavity geometry, applying a numerical method on this flow problem in terms of coding is quite easy and straight forward. Despite its simple geometry, the driven cavity flow retains a rich fluid flow physics manifested by multiple counter rotating recirculating regions on the corners of the cavity depending on the Reynolds number. In the literature, it is possible to find different numerical approaches which have been applied to the driven cavity flow problem [134]. Though this flow problem has been numerically studied extensively, still there are some points which are not agreed upon. For example; 1) an interesting point among many studies is that different numerical solutions of cavity flow yield about the same results for Re £ 1,000 however, start to deviate from each other for larger Re, 2) also another interesting point is that while some studies predict a periodic flow at a high Reynolds number, some others present steady solutions for even a higher Reynolds number. The objective of this work is then to investigate these two points and present accurate very fine grid numerical solutions of steady 2D driven cavity flow.
Among the numerous studies found in the literature, we will give a brief survey about some of the significant studies that uses different types of numerical methods on the driven cavity flow. In doing so the emphasis will be given on three points; on the numerical method used, on the spatial order of the numerical solution and on the largest Reynolds number achieved.
Recently, Barragy & Carey [2] have used a ptype finite element scheme on a 257 ×257 strongly graded and refined element mesh. They have obtained a highly accurate (Dh^{8} order) solutions for steady cavity flow solutions up to Reynolds numbers of Re=12,500. Although Barragy & Carey [2] have presented qualitative solution for Re=16,000 they concluded that their solution for Re=16,000 was underresolved and needed a greater mesh size.
Botella & Peyret [6] have used a Chebyshev collocation method for the solution of the liddriven cavity flow. They have used a subtraction method of the leading terms of the asymptotic expansion of the solution of the NS equations in the vicinity of the corners, where velocity is discontinuous, and obtained a highly accurate spectral solutions for the cavity flow with a maximum of grid mesh of N=160 (polynomial degree) for Reynolds numbers Re £ 9,000. They stated that their numerical solutions exhibit a periodic behavior beyond this Re.
Schreiber & Keller [28] have introduced an efficient numerical technique for steady viscous incompressible flows. The nonlinear differential equations are solved by a sequence of Newton and chord iterations. The linear systems associated with the Newton iteration is solved by LUfactorization with partial pivoting. Applying repeated Richardson extrapolation using the solutions obtained on different grid mesh sizes (maximum being 180 ×180), they have presented highorder accurate (Dh^{8} order in theory) solutions for Reynolds numbers Re £ 10,000.
Benjamin & Denny [5] have used a method of relaxing the algebraic equations by means of the ADI method, with a nonuniform iteration parameter. They have solved the cavity flow for Re £ 10,000 with three different grid mesh sizes (maximum being 101 ×101) and used a Dh^{n} extrapolation to obtain the values when Dh ® 0.
Wright & Gaskell [36] have applied the Block Implicit Multigrid Method (BIMM) to the SMART and QUICK discretizations. They have presented cavity flow results obtained on a 1024 ×1024 grid mesh for Re £ 1,000.
Nishida & Satofuka [24] have presented a new higher order method for simulation of the driven cavity flow. They have discretized the spatial derivatives of the NS equations using a modified differential quadrature (MDQ) method. They have integrated the resulting system of ODEs in time with fourth order RungeKuttaGill (RKG) scheme. With this they have presented spatially Dh^{10} order accurate solutions with grid size of 129 ×129 for Re £ 3,200.
Liao [20], and Liao & Zhu [21], have used a higher order streamfunctionvorticity boundary element method (BEM) formulation for the solution of NS equations. With this they have presented solutions up to Re=10,000 with a grid mesh of 257 ×257.
Hou et. al. [17] have used Lattice Boltzmann Method for simulation of the cavity flow. They have used 256 ×256 grid points and presented solutions up to Reynolds numbers of Re=7,500.
Goyon [13] have solved the streamfunction and vorticity equations using Incremental Unknowns. They have presented steady solutions for Re £ 7,500 on a maximum grid size of 256×256.
Rubin & Khosla [27] have used the strongly implicit numerical method with 2×2 coupled streamfunctionvorticity form of the NS equations. They have obtained solutions for Reynolds numbers in the range of Re £ 3,000 with a grid mesh of 17 ×17.
Ghia et. al. [12] later have applied a multigrid strategy to the coupled strongly implicit method developed by Rubin & Khosla [27]. They have presented solutions for Reynolds numbers as high as Re=10,000 with meshes consisting of as many as 257 ×257 grid points.
Gupta [15] has used a fourth (Dh^{4}) order compact scheme for the numerical solution of the driven cavity flow. He has used a 9 point (3×3) stencil in which the streamfunction and vorticity equations are approximated to fourth order accuracy. He has used pointSOR type of iteration and presented steady cavity flow solutions for Re £ 2,000 with a maximum of 41 ×41 grid mesh.
Li et. al. [19] have used a fourth (Dh^{4}) order compact scheme which had a faster convergence than that of Gupta [15]. They have solved the cavity flow with a grid size of 129 ×129 for Re £ 7,500.
Bruneau & Jouron, [7] have solved the NS equations in primitive variables using a full multigridfull approximation storage (FMGFAS) method. With a grid size 256 ×256, they have obtained steady solutions for Re £ 5,000.
Grigoriev & Dargush [14] have presented a BEM solution with improved penalty function technique using hexagonal subregions and they have discretized the integral equation for each subregion as in FEM. They have used a nonuniform mesh of 5040 quadrilateral cells. With this they were able to solve driven cavity flow up to Re=5,000.
Aydin & Fenner [1] have used BEM formulation with using central and upwind finite difference scheme for the convective terms. They have stated that their formulation lost its reliability for Reynolds numbers greater than 1,000.
To the authors best knowledge, in all these studies among with other numerous papers found in the literature, the maximum Reynolds number achieved for the 2D steady incompressible flow in a liddriven cavity is Re=12,500 and is reported by Barragy & Carey [2]. Although Barragy & Carey [2] have presented solutions for Re=16,000, their solution for this Reynolds number display oscillations related to boundary layer resolution issues due to a coarse mesh. Nallasamy & Prasad [22] have presented steady solutions for Reynolds numbers up to Re £ 50,000, however, their solutions are believed to be inaccurate as a result of excessive numerical dissipation caused by their first order upwind difference scheme.
A very brief discussion on computational as well as experimental studies on the liddriven cavity flow can be found in Shankar & Deshpande [29].
Many factors affect the accuracy of a numerical solution, such as, the number of grids in the computational mesh (Dh), and the spatial discretization order of the finite difference equations, and also the boundary conditions used in the solution.
It is an obvious statement that as the number of grids in a mesh is increased (smaller Dh) a numerical solution gets more accurate. In this study the effect of number of grid points in a mesh on the accuracy of the numerical solution of driven cavity flow is investigated, especially as the Reynolds number increases. For this, the governing NS equations are solved on progressively increasing number of grid points (from 128 ×128 to 600×600) and the solutions are compared with the highly accurate benchmark solutions found in the literature.
Order of spatial discretization is another factor on the accuracy of a numerical solution. Finite element and spectral methods provide high spatial accuracy. This study also investigates the effect of spacial accuracy on the numerical solution of driven cavity flow. For this, following [28], repeated Richardson extrapolation is used and then the extrapolated results are compared with the benchmark solutions.
Another factor on the accuracy of a numerical solution is the type and the order of the numerical boundary conditions used in computation. For driven cavity flow, this subject was briefly discussed by Weinan & JianGuo [34], by Spotz [30], and also by Gupta & Manohar [16].
Until Barragy & Carey [2], many studies have presented steady solutions of driven cavity flow for Re=10,000 ([28], [5], [21] and [12]). They ([2]) have presented solutions for Re=12,500. In all these studies that have presented high Reynolds number numerical solutions, the largest grid size used was 256×256. In this study, the grid size will be increased up to 600 ×600 and the effect of the grid size on the largest computable Reynolds number solution of the driven cavity flow will be investigated.
This paper presents accurate, very fine grid (600 ×600) numerical solutions of 2D steady incompressible flow in a liddriven cavity for Reynolds numbers up to Re=21,000. A detailed comparisons of our results with mainly the studies mentioned above as well as with other studies not mentioned here, will be done.
For twodimensional and axisymmetric flows it is convenient to use the streamfunction (y) and vorticity (w) formulation of the NavierStokes equations. In nondimensional form, they are given as
¶^{2}y
¶x^{2}
+
¶^{2}y
¶y^{2}
= w
(1)
1
Re
æ
è
¶^{2} w
¶x^{2}
+
¶^{2} w
¶y^{2}
ö
ø
=
¶y
¶y
¶w
¶x

¶y
¶x
¶w
¶y
(2)
where, Re is the Reynolds number, and x and y are the Cartesian coordinates. The NS equations are nonlinear, and therefore need to be solved in an iterative manner. In order to have an iterative numerical algorithm, pseudo time derivatives are assigned to these equations (1 and 2) such as
¶y
¶t
=
¶^{2}y
¶x^{2}
+
¶^{2}y
¶y^{2}
+ w
(3)
¶w
¶t
=
1
Re
¶^{2}w
¶x^{2}
+
1
Re
¶^{2}w
¶y^{2}
+
¶y
¶x
¶w
¶y

¶y
¶y
¶w
¶x
(4)
Since we are seeking a steady solution, the order of the pseudo time derivatives will not affect the final solution. Therefore, an implicit Euler time step is used for these pseudo time derivatives, which is first order (Dt) accurate. Also for the nonlinear terms in the vorticity equation, a first order (Dt) accurate approximation is used and the governing equations (3 and 4) are written as the following
y^{n+1}  Dt
¶^{2}y
¶x^{2}
n+1
 Dt
¶^{2}y
¶y^{2}
n+1
= y^{n} + Dt w^{n}
(5)
w^{n+1}  Dt
1
Re
¶^{2}w
¶x^{2}
n+1
 Dt
1
Re
¶^{2}w
¶y^{2}
n+1
 Dt
¶y
¶x
n
¶w
¶y
n+1
+ Dt
¶y
¶y
n
¶w
¶x
n+1
= w^{n}
(6)
In operator notation the above equations look like the following.
æ
è
1  Dt
¶^{2}
¶x^{2}
 Dt
¶^{2}
¶y^{2}
ö
ø
y^{n+1} = y^{n}+ Dt w^{n}
(7)
æ
è
1  Dt
1
Re
¶^{2}
¶x^{2}
 Dt
1
Re
¶^{2}
¶y^{2}
 Dt
æ
è
¶y
¶x
ö
ø
n
¶
¶y
+ Dt
æ
è
¶y
¶y
ö
ø
n
¶
¶x
ö
ø
w^{n+1} = w^{n}
(8)
The above equations are in fully implicit form where each equation requires the solution of a large banded matrix which is not computationally efficient. Instead of solving the equations (7 and 8) in fully implicit form, more efficiently, we spatially factorize the equations, such as
æ
è
1  Dt
¶^{2}
¶x^{2}
ö
ø
æ
è
1  Dt
¶^{2}
¶y^{2}
ö
ø
y^{n+1} = y^{n} + Dt w^{n}
(9)
æ
è
1  Dt
1
Re
¶^{2}
¶x^{2}
+ Dt
æ
è
¶y
¶y
ö
ø
n
¶
¶x
ö
ø
æ
è
1  Dt
1
Re
¶^{2}
¶y^{2}
 Dt
æ
è
¶y
¶x
ö
ø
n
¶
¶y
ö
ø
w^{n+1} = w^{n}
(10)
Here, we split the single LHS operators into two operators, where each only contains derivatives in one direction. The only difference between equations (7 and 8) and (9 and 10) are the Dt^{2} terms produced by the factorization. The advantage of this process is that each equation now require the solution of a tridiagonal system, which is numerically more efficient.
As the solution converges to the steady state solution, the Dt^{2} terms due to the factorization will not become zero and will remain in the solution. To illustrate this, equation (9) can be written in explicit form as
y^{n+1}  Dt
¶^{2} y^{n+1}
¶x^{2}
 Dt
¶^{2} y^{n+1}
¶y^{2}
+ Dt^{2}
¶^{4} y^{n+1}
¶x^{2} ¶y^{2}
= y^{n} + Dt w^{n}
(11)
As the solution converges to a steady state, y^{n+1} becomes equal to y^{n} and they cancel each other from the above equation and at convergence, equation (11) therefore appears as the following
¶^{2} y
¶x^{2}
+
¶^{2} y
¶y^{2}
+ w =  Dt
¶^{4} y
¶x^{2} ¶y^{2}
(12)
In equation (12), there is no guarantee that the RHS, which is the consequence of the factorization, will be small. Only a very small time step (Dt) will ensure that the RHS of equation (12) will be close to zero. This, however, will slow down the convergence. In addition, the product may not be small since [(¶^{4} y)/(¶x^{2} ¶y^{2})] can be O ( [1/(Dt)] ) or higher in many flow situations.
In order to overcome this problem, we add a similar term at the previous time level to the RHS of equation (11) to give the following finite difference form
y^{n+1}  Dt
¶^{2} y^{n+1}
¶x^{2}
 Dt
¶^{2} y^{n+1}
¶y^{2}
+ Dt^{2}
¶^{4} y^{n+1}
¶x^{2} ¶y^{2}
= y^{n} + Dt w^{n} + Dt^{2}
¶^{4} y^{n}
¶x^{2} ¶y^{2}
(13)
As a result, at steady state, y^{n+1} converges to y^{n} and [(¶^{4} y^{n+1})/(¶x^{2}¶y^{2})] converges to [(¶^{4}y^{n})/(¶x^{2} ¶y^{2})], so that they cancel each out, and the final solution converges to the correct physical form
¶^{2} y
¶x^{2}
+
¶^{2} y
¶y^{2}
+ w = 0
(14)
We would like to emphasize that the added Dt^{2} term on the RHS is not an artificial numerical diffusion term. Such terms are often used to stabilize the numerical schemes so that convergence can be achieved. Here, the added Dt^{2} term on the RHS is not meant to stabilize the numerical scheme. The only reason for this approach is to cancel out the terms introduced by factorization so that the equations are the correct physical representation at the steady state.
One of the most popular spatially factorized schemes is undoubtedly the Beam & Warming [4] method. In that method, the equations are formulated in delta (D) form. The main advantage of a delta formulation in a steady problem is that there will be no second order (Dt^{2}) terms due to the factorization in the solution at the steady state. However in our formulation the second order (Dt^{2}) terms due to factorization will remain even at the steady state. Adding a second order (Dt^{2}) term to the RHS of the equations to cancel out the terms produced by factorization, easily overcomes this problem. One of the main differences between the presented formulation and a delta formulation is in the numerical treatment of the vorticity equation. In a delta formulation of the vorticity equation, there appears Dt order convection and dissipation terms on the RHS (explicit side) of the numerical equations. In the presented formulation, on the RHS of the equation there are no Dt order terms. Instead cross derivative terms ([(¶^{4})/(¶x^{2}¶y^{2})], [(¶^{3})/(¶x^{2}¶y)], [(¶^{3})/(¶x¶y^{2})], and [(¶^{2})/(¶x¶y)]) appear on the RHS. Note that these terms are Dt^{2} order. Since Dt is small, the contribution of these Dt^{2} terms are even smaller. Our extensive numerical tests showed that this approach has better numerical stability characteristics and is more effective especially at high Reynolds numbers compared to a delta formulation of the steady streamfunction and vorticity equations.
Applying the approach of removing second order terms due to factorization on the RHS, to both equations (9 and 10), the final form of the numerical formulation we used for the solution of streamfunction and vorticity equation becomes
æ
è
1  Dt
¶^{2}
¶x^{2}
ö
ø
æ
è
1  Dt
¶^{2}
¶y^{2}
ö
ø
y^{n+1} = y^{n} + Dt w^{n}
+
æ
è
Dt
¶^{2}
¶x^{2}
ö
ø
æ
è
Dt
¶^{2}
¶y^{2}
ö
ø
y^{n}
(15)
æ
è
1  Dt
1
Re
¶^{2}
¶x^{2}
+ Dt
æ
è
¶y
¶y
ö
ø
n
¶
¶x
ö
ø
æ
è
1  Dt
1
Re
¶^{2}
¶y^{2}
 Dt
æ
è
¶y
¶x
ö
ø
n
¶
¶y
ö
ø
w^{n+1} = w^{n}
+
æ
è
Dt
1
Re
¶^{2}
¶x^{2}
 Dt
æ
è
¶y
¶y
ö
ø
n
¶
¶x
ö
ø
æ
è
Dt
1
Re
¶^{2}
¶y^{2}
+ Dt
æ
è
¶y
¶x
ö
ø
n
¶
¶y
ö
ø
w^{n}
(16)
The solution methodology for the two equations involves a twolevel updating. For the streamfunction equation, the variable f is introduced such that
æ
è
1  Dt
¶^{2}
¶y^{2}
ö
ø
y^{n+1} = f
(17)
and
æ
è
1  Dt
¶^{2}
¶x^{2}
ö
ø
f = y^{n} + Dt w^{n} +
æ
è
Dt
¶^{2}
¶x^{2}
ö
ø
æ
è
Dt
¶^{2}
¶y^{2}
ö
ø
y^{n}
(18)
In equation (18) f is the only unknown. First, this equation is solved for f at each grid point. Following this, the streamfunction (y) is advanced into the new time level using equation (17).
In a similar fashion for the vorticity equation, the variable g is introduced such that
æ
è
1  Dt
1
Re
¶^{2}
¶y^{2}
 Dt
æ
è
¶y
¶x
ö
ø
n
¶
¶y
ö
ø
w^{n+1} = g
(19)
and
æ
è
1  Dt
1
Re
¶^{2}
¶x^{2}
+ Dt
æ
è
¶y
¶y
ö
ø
n
¶
¶x
ö
ø
g = w^{n}
+
æ
è
Dt
1
Re
¶^{2}
¶x^{2}
 Dt
æ
è
¶y
¶y
ö
ø
n
¶
¶x
ö
ø
æ
è
Dt
1
Re
¶^{2}
¶y^{2}
+ Dt
æ
è
¶y
¶x
ö
ø
n
¶
¶y
ö
ø
w^{n}
(20)
As with f, g is determined at every grid point using equation (20), then vorticity (w) is advanced into the next time level using equation (19).
We note that, in this numerical approach the streamfunction and vorticity equations are solved separately. Each equation is advanced into a next time level by solving two tridiagonal systems, which allows the use of very large grid sizes easily. The method proved to be very effective on flow problems that require high accuracy on very fine grid meshes (Erturk & Corke [9] Erturk, Haddad & Corke [10]).
The boundary conditions and a schematics of the vortices generated in a driven cavity flow are shown in Figure 1.
In this figure, the abbreviations BR, BL and TL refer to bottom right, bottom left and top left corners of the cavity, respectively. The number following these abbreviations refer to the vortices that appear in the flow, which are numbered according to size.
We have used the wellknown Thom's formula [32] for the wall boundary condition, such that
y_{0} = 0
w_{0} =
2 y_{1}
Dh^{2}

2 U
Dh
(21)
where subscript 0 refers to points on the wall and 1 refers to points adjacent to the wall, Dh refers to grid spacing and U refers to the velocity of the wall with being equal to 1 on the moving wall and 0 on the stationary walls.
We note that, it is well understood ([34], [30], [18], [23]) that, even though Thom's method is locally first order accurate, the global solution obtained using Thom's method preserve second order accuracy. Therefore in this study, since three point second order central difference is used inside the cavity and Thom's method is used at the wall boundary conditions, the presented solutions are second order accurate O (Dh^{2}).
During our computations we monitored the residual of the steady streamfunction and vorticity equations (1 and 2) as a measure of the convergence to the steady state solution, where the residual of each equation is given as
R_{y} =
y_{i1,j}^{n+1}2y_{i,j}^{n+1}+y_{i+1,j}^{n+1}
Dx^{2}
+
y_{i,j1}^{n+1}2y_{i,j}^{n+1}+y_{i,j+1}^{n+1}
Dy^{2}
+w_{i,j}^{n+1}
(22)
R_{w} =
1
Re
w_{i1,j}^{n+1}2w_{i,j}^{n+1} +w_{i+1,j}^{n+1}
Dx^{2}
+
1
Re
w_{i,j1}^{n+1}2 w_{i,j}^{n+1}+w_{i,j+1}^{n+1}
Dy^{2}

y_{i,j+1}^{n+1}y_{i,j1}^{n+1}
2Dy
w_{i+1,j}^{n+1}w_{i1,j}^{n+1}
2Dx
+
y_{i+1,j}^{n+1}y_{i1,j}^{n+1}
2 Dx
w_{i,j+1}^{n+1}w_{i,j1}^{n+1}
2 Dy
(23)
The magnitude of these residuals are an indication of the degree to which the solution has converged to steady state. In the limit these residuals would be zero. In our computations, for all Reynolds numbers, we considered that convergence was achieved when for each equation (22 and 23) the maximum of the absolute residual in the computational domain (max(R_{y}) and max(R_{w})) was less than 10^{10}. Such a low value was chosen to ensure the accuracy of the solution. At these residual levels, the maximum absolute difference in streamfunction value between two time steps, (max(y^{n+1}y^{n})), was in the order of 10^{16} and for vorticity, (max(w^{n+1}w^{n})), it was in the order of 10^{14}. And also at these convergence levels, between two time steps the maximum absolute normalized difference in streamfunction, (max([(y^{n+1}y^{n})/(y^{n})])), and in vorticity, (max([(w^{n+1}w^{n})/(w^{n})])), were in the order of 10^{13}, and 10^{12} respectively.
We have started our computations by solving the driven cavity flow from R=1,000 to Re=5,000 on a 129×129 grid mesh. With using this many number of grids, we could not get a steady solution for Re=7,500. The first natural conclusion was that the cavity flow may not be steady at Re=7,500, and therefore steady solution may not be computable. Even though we have used a pseudotime iteration, the time evolution of the flow parameters, either vorticity or streamfunction values at certain locations, were periodic suggesting that the flow indeed may be periodic, backing up this conclusion. Besides, in the literature many studies claim that the driven cavity flow becomes unstable between a Reynolds number of 7,000 and 8,000, for example, Peng et. al. [25], Fortin et. al. [11] and Poliashenko & Aidun [26] predict unstable flow beyond Re of 7402, 7998,5 and 7763 respectively. However, on the other hand there were many studies ([2], [5], [12], [21], and [28]) that have presented solutions for a higher Reynolds number, Re=10,000. This was contradictory. We then have tried to solve the same case, Re=7,500, with a larger grid size of 257×257. This time we were able to obtain a steady solution. In fact, with 257×257 number of grids, we carry our steady computations up to Re=12,500. Barragy & Carey [2] have also presented steady solutions for Re=12,500. Interesting point was, when we tried to increase Re furthermore with the same grid mesh, we have again obtained a periodic behavior. Once again the natural conclusion was that the steady solutions beyond this Re may not be computable. Having concluded this, we could not help but ask the question what happens if we increase the number of grids furthermore. In their study Barragy & Carey [2] have presented a steady solution for Re=16,000. Their solution for this Reynolds number displayed oscillations. They concluded that this was due to a coarse mesh. This conclusion encouraged us to use larger grid size than 257×257 in order to obtain solutions for Reynolds numbers greater than 12,500. At this point, since the focus was on the number of grids, we realized that among the studies that presented steady solutions at very high Reynolds numbers ([2], [5], [12], [21], and [28]), the maximum number of grids used was 257×257. Larger grid sizes haven't been used at these high Reynolds numbers. This was also encouraging us to use more grid points. We, then, decided to increase the number of grids and tried to solve for Re=15,000 with 401×401 number of grids. This time we were able to obtain a steady solution. This fact suggests that in order to obtain a numerical solution for 2D steady incompressible driven cavity flow for Re ³ 12,500, a grid mesh with more than 257×257 grid points is needed. With using 401×401 grid points, we continued our steady computations up to Re=21,000. Beyond this Re our computations again displayed a periodic behavior. This time the first natural thing came to mind was to increase the number of grids. We have tried to use 513×513 grids, however again we could not get a steady solution beyond Re=21,000. Thinking that the increase in number of grids may not be enough, we then again increased the number of grids and have used 601×601 grids. The situation was the same, and the maximum Reynolds number that we can obtain a steady solution with using 601×601 grids was Re=21,000. We did not try to increase the number of grids furthermore since the computations became time consuming. Whether or not steady computations of driven cavity flow are possible beyond Re=21,000 with grid numbers larger than 601×601, is still an open question.
One of the reasons, why the steady solutions of the driven cavity flow at very high Reynolds numbers become computable when finer grids are used, may be the fact that as the number of grids used increases, Dh gets smaller, then the cell Reynolds number or so called Peclet number defined as Re_{c}=[(uDh)/(n)] decreases. This improves the numerical stability characteristics of the numerical scheme used ([31], [34]), and allows high cavity Reynolds numbered solutions computable. Another reason may be the fact that finer grids would resolve the corner vortices better. This would, then, help decrease any numerical oscillations that might occur at the corners of the cavity during iterations.
We note that all the figures and tables present the solution of the finest grid size of 601×601 unless otherwise stated.
The accuracy of a finite difference solution is set by the mesh size, and by the spatial order of the finite difference equations and the boundary approximations. At low Reynolds numbers (Re=1,000), different numerical method solutions found in the literature agree with each other. However, the solutions at higher Reynolds numbers (Re > 1,000) have noticeable discrepancies. We believe it is mainly due to different spatial orders and different grid mesh sizes and different boundary conditions used in different studies.
Table 1 tabulates the minimum streamfunction value, the vorticity value at the center of the primary vortex and also the center location of the primary vortex for Re=1,000 along with similar results found in the literature with the most significant ones are underlined.
In Table 1 among the most significant (underlined) results, Schreiber & Keller [28] have used Richardson extrapolation in order to achieve high spatial accuracy. For Re=1,000 their extrapolated solutions are Dh^{6} order accurate and these solutions are obtained by using repeated Richardson extrapolation on three different mesh size solutions.
Following the same approach, we have solved the cavity flow on three different grid mesh (401×401, 513×513, and 601×601) separately for every Reynolds number (Re=1,000 to Re=21,000). Since the numerical solution is second order accurate O (Dh^{2}), the computed values of streamfunction and vorticity have an asymptotic error expansion of the form
We can use Richardson extrapolation in order to obtain highorder accurate approximations (Schreiber & Keller [28]). We have used repeated Richardson extrapolations using the three different mesh size solutions. The extrapolated values of the streamfunction and the vorticity at the center of the primary vortex are given in Table 2 and Table 3 respectively. The extrapolated results tabulated in Table 2 and Table 3 are, in theory, Dh^{6} order accurate.
Looking back to Table 1, for Re=1,000, our extrapolated results are in very good agrement with the extrapolated results of Schreiber & Keller [28]. In fact our extrapolated streamfunction value at the primary vortex differs only 0.013 % and the vorticity value differs only 0.095 % from the extrapolated results of Schreiber & Keller [28].
With spectral methods, very high spatial accuracy can be obtained with a relatively smaller number of grid points. For this reason, the tabulated results from Botella & Peyret [6] are believed to be very accurate. Our results are also in very good agreement such that our extrapolated results differs 0.016 % in the streamfunction value and 0.098 % in the vorticity value from the results of Botella & Peyret [6].
Barragy & Carey [2] have only tabulated the location of the vortex centers and the streamfunction values at these locations yet they have not presented any vorticity data. Nevertheless, their solutions based on ptype finite element scheme also have very high spatial accuracy (Dh^{8} order accurate). They have presented the maximum Reynolds number solutions (Re=12,500) for the cavity flow, to the best of our knowledge, found in the literature. We will show that the agreement between our results and results of Barragy & Carey [2] is excellent through the whole Reynolds number range they considered (1,000 £ Re £ 12,500). Comparing solutions for Re=1,000 in Table 1, the difference is 0.022 % in the streamfunction value.
Wright & Gaskell [36] have used a very fine grid mesh with 1024×1024 grid points in their solutions. With this many grids their solution should be quite accurate. Their streamfunction and vorticity value are very close to our solutions with a difference of 0.113 % and 0.114 % respectively.
Nishida & Satofuka [24] have presented higher order solutions in which for Re=1,000 their solutions are Dh^{8} order accurate. Comparing with our results we find a very good agreement such that our extrapolated results differs 0.040 % in streamfunction value and 0.136 % in vorticity value.
Grigoriev & Dargush [14] have solved the governing equations with a boundary element method (BEM). They have only tabulated the streamfunction values and the location of the centers of the primary and secondary vortices for driven cavity flow for Reynolds numbers up to 5,000. Their solution for Re=1,000 is quite accurate.
Benjamin & Denny [5] have used a Dh^{n} extrapolation to obtain the values when Dh ® 0 using the three different grid mesh size solutions. Their extrapolated results are close to other underlined data in Table 1. It is quite remarkable that even a simple extrapolation produces such good results.
Hou et. al. [17] have simulated the cavity flow by the Lattice Boltzmann Method. They have presented solutions for Re £ 7,500 obtained by using 256×256 lattice. However their results have significant oscillations in the two upper corners of the cavity. Especially their vorticity contours have wiggles aligned with lattice directions, for Re ³ 1,000.
Li et. al. [19] have used a compact scheme and presented fourth order accurate (Dh^{4}) results for a 129×129 grid mesh. Ghia et. al. [12] have used a second order method (Dh^{2}) with a mesh size of 129×129. These results, together with that of Goyon [13] and Liao & Zhu [21], might be considered as somewhat underresolved. Our calculations showed that for a second order (Dh^{2}) spatial accuracy, even a 401×401 grid mesh solutions can be considered as underresolved for Re=1,000.
The results of Bruneau & Jouron, [7] and Vanka [33] appear overdiffusive due to upwind differencing they have used.
From all these comparisons we can conclude that even for Re=1,000 higher order approximations together with the use of fine grids are necessary for accuracy.
Figures 2 through 6 show streamfunction contours of the finest (601×601) grid solutions, for various Reynolds numbers. These figures exhibit the formation of the counterrotating secondary vortices which appear as the Reynolds number increases.
Table 4 compares our computed streamfunction and vorticity values at the center of the primary vortex with the results found in the literature for various Reynolds numbers.
The agreement between our Richardson extrapolated results (Dh^{6} order accurate) with results of Barragy & Carey [2] based on ptype finite element scheme (Dh^{8} order accurate) is quite remarkable. One interesting result is that the streamfunction value at the center of the primary vortex decreases at first as the Reynolds number increases until Re=10,000, however then starts to increase as the Reynolds number increases further. This behavior was also documented by Barragy & Carey [2].
In terms of quantitative results, Gupta & Manohar [16], stated that for the numerical solutions of driven cavity flow, the value of the streamfunction (y) at the primary vortex center appears to be a reliable indicator of accuracy, where as the location of the vortex center (x and y) is limited by the mesh size, and the value of the vorticity (w) at the vortex center is sensitive to the accuracy of the wall boundary conditions. Nevertheless, the locations of the primary and secondary vortices, as well as the values of the streamfunction (y) and vorticity (w) at these locations are presented in Table 5.
These are comparable with the solutions of Barragy & Carey [2] and Ghia et. al. [12]. It is evident that as the Reynolds number increases, the center of the primary vortex moves towards the geometric center of the cavity. Its location is almost invariant for Re ³ 17,500. Our computations indicate the appearance of a quaternary vortex at the bottom left corner (BL3) at Reynolds number of 10,000. We would like to note that, this quaternary vortex (BL3) did not appear in our solutions for Re=10,000 with grid sizes less than 513×513. This would suggest that, especially at high Reynolds numbers, in order to resolve the small vortices appear at the corners properly, very fine grids have to be used. Moreover, our solutions indicate another tertiary vortex at the top left side of the cavity (TL2) at Re=12,500. These vortices (BL3 and TL2) were also observed by Barragy & Carey [2].
Figures 7 and 8 present the uvelocity profiles along a vertical line and the vvelocity profiles along a horizontal line passing through the geometric center of the cavity respectively and also Tables 6 and 7 tabulates u and vvelocity values at certain locations respectively, for future references. These profiles are in good agrement with that of Ghia et. al. [12] shown by symbols in Figures 7 and 8.
We have solved the incompressible flow in a driven cavity numerically. A very good mathematical check on the accuracy of the numerical solution would be to check the continuity of the fluid, as suggested by Aydin & Fenner [1], although this is done very rarely in driven cavity flow papers. We have integrated the uvelocity and vvelocity profiles, considered in Figures 7 and 8, along a vertical line and horizontal line passing through the geometric center of the cavity, in order to obtain the net volumetric flow rate, Q, through these sections. Since the flow is incompressible, the net volumetric flow rate passing through these sections should be equal to zero, Q=0. The uvelocity and vvelocity profiles at every Reynolds number considered are integrated using Simpson's rule, then the obtained volumetric flow rate values are divided by a characteristic flow rate, Q_{c}, which is the horizontal rate that would occur in the absence of the side walls (Plane Couette flow), to help quantify the errors, as also suggested by Aydin & Fenner [1]. In an integration process the numerical errors will add up. Nevertheless, the volumetric flow rate values (Q_{1}=[(ò^{1}_{0} u dy)/(Q_{c})] and Q_{2}=[(ò^{1}_{0} v dx)/(Q_{c})]) tabulated in Table 8 are close to zero such that, even the largest values of Q_{1} = 0.000000705 and Q_{2} = 0.000002424 in Table 8 can be considered as Q_{1} » Q_{2} » 0. This mathematical check on the conservation of the continuity shows that our numerical solution is indeed very accurate.
Distinctively than any other paper, Botella & Peyret [6] have tabulated highly accurate Chebyshev collocation spectral vorticity data from inside the cavity, along a vertical line and along a horizontal line passing through the geometric center of the cavity, for Re=1,000. Comparing our solutions with theirs ([6]) in Figure 9 and 10, we find that the agreement with Botella & Peyret [6] is remarkable, with the maximum difference between two solutions being 0.18 %.
Between the whole range of Re=1,000 and the maximum Re found in the literature, Re=12,500, our computed results agreed well with the published results. Since there is no previously published results (to the best of our knowledge) of steady driven cavity flow beyond Re > 12,500 in the literature, we decided to use Batchelor's [3] model in order to demonstrate the accuracy of our numerical solutions at higher Reynolds numbers. As seen in Figures 7 and 8, the u and vvelocity profiles change almost linearly in the core of the primary vortex as Reynolds number increases. This would indicate that in this region the vorticity is uniform. As Re increases thin boundary layers are developed along the solid walls and the core fluid moves as a solid body with a uniform vorticity, in the manner suggested by Batchelor [3]. At high Reynolds numbers the vorticity at the core of the eddy is approximately constant and the flow in the core is governed by
¶^{2}y
¶x^{2}
+
¶^{2}y
¶y^{2}
= C
(25)
where C is the constant vorticity value. With streamfunction value being y = 0 on the boundaries, inside the domain (x,y)=([0.5,0.5],[0.5,0.5]) the following expression is a solution to equation (25)
y = C
æ
è
y^{2}
2

1
8
+
¥
å
n=1
A_{n} cosh(2n1) px cos(2n1) py
ö
ø
(26)
where
A_{n}=
4(1)^{n+1}
(2n1)^{3} p^{3}
(cosh(n
1
2
)p)^{1}
(27)
The mean square law proposed by Batchelor [3] relates the core vorticity value with the boundary velocities. From using the relation
ó
(ç)
õ
[(y_{x})^{2}+(y_{y})^{2}] d s =
ó
(ç)
õ
U^{2}d s = 1
(28)
where U is the boundary velocity, we can have
1
C^{2}
=
1
6

16
p^{4}
¥
å
n=1
1
(2n1)^{4}(cosh(n
1
2
)p)^{2}
+ 2
ó
õ
[1/2]
[1/2]
æ
è
1
2
+p
¥
å
n=1
(1)^{n} (2n1) A_{n} cosh(2n1) px
ö
ø
2
d x
(29)
From the above expression, the numerical value of C (vorticity value at the core) can be found approximately as 1.8859, in accordance with the theoretical core vorticity value of 1.886 at infinite Re, calculated analytically by Burggraf [8] in his application of Batchelor's model [3], which consists of an inviscid core with uniform vorticity, coupled to boundary layer flows at the solid surface.
In Figure 11 we have plotted the variation of the vorticity value at the center of the primary vortex, tabulated in Table 2, as a function of Reynolds number. In the same figure, the theoretical value of 1.8859 is shown with the dotted line. This figure clearly shows that our computations, in fact, in strong agrement with the mathematical theory (Batchelor [3], Wood [35], and Burggraf [8]) such that computed values of the vorticity at the primary vortex asymptote to the theoretical infinite Re vorticity value as the Reynolds number increases.
Fine grid numerical solutions of the steady driven cavity flow for Reynolds numbers up to Re=21,000 have been presented. The steady 2D incompressible NavierStokes equations in streamfunction and vorticity formulation are solved computationally using the numerical method described. The streamfunction and vorticity equations are solved separately. For each equation, the numerical formulation requires the solution of two tridiagonal systems, which allows the use of large grid meshes easily. With this we have used a fine grid mesh of 601×601. At all Reynolds numbers, the numerical solutions converged to maximum absolute residuals of the governing equations that were less than 10^{10}.
Our computations indicate that fine grid mesh is necessary in order to obtain a steady solution and also resolve the vortices appear at the corners of the cavity, as the Reynolds number increases. Several unique features of the driven cavity flow have been presented such that the appearance of a quaternary vortex at the bottom left corner (BL3) at Reynolds number of 10,000 on fine grid mesh, and a tertiary vortex at the top left corner (TL2) at Reynolds number of 12,500. Detailed results were presented and they were compared with results found in the literature. The presented results were found to agree well with the published numerical solutions and with the theory as well.
Acknowledgement
This study was funded by Gebze Institute of Technology with project no BAP2003A22.
E. Ertürk is grateful for this financial support.
M. Aydin, R.T. Fenner. Boundary Element Analysis of Driven Cavity Flow for Low and Moderate Reynolds Numbers. Int. J. Numer. Meth. Fluids 2001; 37:4564.
E. Erturk, O.M. Haddad, T.C. Corke. Numerical Solutions of Laminar Incompressible Flow Past Parabolic Bodies at Angles of Attack. AIAA J. 2004; 42:22542265.
A. Fortin, M. Jardak, J.J. Gervais, R. Pierre. Localization of Hopf Bifurcations in Fluid Flow Problems. Int. J. Numer. Methods Fluids 1997; 24:11851210.
U. Ghia, K.N. Ghia, C.T. Shin. HighRe Solutions for Incompressible Flow Using the NavierStokes Equations and a Multigrid Method. J. Comp. Physics 1982; 48:387411.
O. Goyon. HighReynolds Number Solutions of NavierStokes Equations Using Incremental Unknowns. Computer Methods in Applied Mechanics and Engineering 1996; 130:319335.
M. Li, T. Tang, B. Fornberg. A Compact ForthOrder Finite Difference Scheme for the Steady Incompressible NavierStokes Equations. Int. J. Numer. Methods Fluids 1995; 20:11371151.
S.J. Liao, J.M. Zhu. A Short Note on HigherOrder StremfunctionVorticity Formulation of 2D Steady State NavierStokes Equations. Int. J. Numer. Methods Fluids 1996; 22:19.
M. Napolitano, G. Pascazio, L. Quartapelle. A Review of Vorticity Conditions in the Numerical Solution of the z  y Equations. Computers and Fluids 1999; 28:139–185.
H. Nishida, N. Satofuka. HigherOrder Solutions of Square Driven Cavity Flow Using a VariableOrder MultiGrid Method. Int. J. Numer. Methods Fluids 1992; 34:637653.