Published in : ZAMM  Journal of Applied Mathematics and Mechanics (2007)
ZAMM  Z. Angew. Math. Mech. 2007, Vol 87, pp 377392
Numerical Solutions of 2D Steady Incompressible Flow in a Driven Skewed Cavity
Ercan Erturk and Bahtiyar Dursun
Energy Systems Engineering Department, Gebze Institute of Technology,
Gebze, Kocaeli 41400, Turkey
Abstract
The benchmark test case for nonorthogonal grid mesh, the "driven
skewed cavity flow", first introduced by Demirdzic et al.
(1992, IJNMF, 15, 329) for skew angles of a=30^{°} and
a=45^{°}, is reintroduced with a more variety of skew
angles. The benchmark problem has nonorthogonal, skewed grid mesh
with skew angle (a). The governing 2D steady incompressible
NavierStokes equations in general curvilinear coordinates are
solved for the solution of driven skewed cavity flow with
nonorthogonal grid mesh using a numerical method which is efficient
and stable even at extreme skew angles. Highly accurate numerical
solutions of the driven skewed cavity flow, solved using a fine grid
(512×512) mesh, are presented for Reynolds number of 100 and
1000 for skew angles ranging between
15^{°} £ a £ 165^{°}.
1 INTRODUCTION
In the literature, it is possible to find many numerical methods
proposed for the solution of the steady incompressible NS
equations. These numerical methods are often tested on several
benchmark test cases in terms of their stability, accuracy as well
as efficiency. Among several benchmark test cases for steady
incompressible flow solvers, the driven cavity flow is a very well
known and commonly used benchmark problem. The reason why the driven
cavity flow is so popular may be the simplicity of the geometry. In
this flow problem, when the flow variables are nondimensionalized
with the cavity length and the velocity of the lid, Reynolds number
appears in the equations as an important flow parameter. Even though
the geometry is simple and easy to apply in programming point of
view, the cavity flow has all essential flow physics with counter
rotating recirculating regions at the corners of the cavity. Among
numerous papers found in the literature, Erturk et al.
[ 6], Botella and Peyret [ 4], Schreiber and
Keller [ 21], Li et al. [ 12], Wright and
Gaskel [ 30], Erturk and Gokcol [ 7], Benjamin and
Denny [ 2] and Nishida and Satofuka [ 16] are
examples of numerical studies on the driven cavity flow.
Due to its simple geometry, the cavity flow is best solved in
Cartesian coordinates with Cartesian grid mesh. Most of the
benchmark test cases found in the literature have orthogonal
geometries therefore they are best solved with orthogonal grid mesh.
However often times the real life flow problems have much more
complex geometries than that of the driven cavity flow. In most
cases, researchers have to deal with nonorthogonal geometries with
nonorthogonal grid mesh. In a nonorthogonal grid mesh, when the
governing equations are formulated in general curvilinear
coordinates, cross derivative terms appear in the equations.
Depending on the skewness of the grid mesh, these cross derivative
terms can be very significant and can affect the numerical stability
as well as the accuracy of the numerical method used for the
solution. Even though, the driven cavity flow benchmark problem
serve for comparison between numerical methods, the flow is far from
simulating the real life fluid problems with complex geometries with
nonorthogonal grid mesh. The numerical performances of numerical
methods on orthogonal grids may or may not be the same on
nonorthogonal grids.
Unfortunately, there are not much benchmark problems with
nonorthogonal grids for numerical methods to compare solutions with
each other. Demirdzic et al. [ 5] have introduced
the driven skewed cavity flow as a test case for nonorthogonal
grids. The test case is similar to driven cavity flow but the
geometry is a parallelogram rather than a square. In this test case,
the skewness of the geometry can be easily changed by changing the
skew angle ( a). The skewed cavity problem is a perfect test
case for body fitted nonorthogonal grids and yet it is as simple as
the cavity flow in terms of programming point of view. Later
Oosterlee et al. [ 17], Louaked et al.
[ 13], Roychowdhury et al. [ 20], Xu and Zhang
[ 31], Wang and Komori [ 28], Xu and Zhang
[ 32], Tucker and Pan [ 27], Brakkee et al.
[ 3], Pacheco and Peck [ 18], Teigland and
Eliassen [ 25], Lai and Yan [ 11] and Shklyar and
Arbel [ 22] have solved the same benchmark problem. In all
these studies, the solution of the driven skewed cavity flow is
presented for Reynolds numbers of 100 and 1000 for only two
different skew angles which are a=30 ^{°} and
a=45 ^{°}, and also the maximum number of grids used in
these studies is 320×320.
Peric [ 19] considered the 2D flow in a skewed cavity and
he stated that the governing equations fail to converge for
a < 30 ^{°}. The main motivation of this study is then to
reintroduce the skewed cavity flow problem with a wide range of skew
angle (15 ^{°} £ a £ 165 ^{°}) and present
detailed tabulated results obtained using a fine grid mesh with
512×512 points for future references.
Erturk et al. [ 6] have introduced an efficient,
fast and stable numerical formulation for the steady incompressible
NavierStokes equations. Their method solve the streamfunction and
vorticity equations separately, and the numerical solution of each
equation requires the solution of two tridiagonal systems. Solving
tridiagonal systems are computationally efficient and therefore they
were able to use very fine grid mesh in their solution. Using this
numerical formulation, they have solved the very well known
benchmark problem, the steady flow in a square driven cavity, up to
Reynolds number of 21000 using a 601×601 fine grid mesh.
Their formulation proved to be stable and effective at very high
Reynolds numbers ([ 6], [ 7], [ 8]).
In this study, the numerical formulation introduced by Erturk
et al. [ 6] will be applied to NavierStokes equations in
general curvilinear coordinates and the numerical solutions of the
driven skewed cavity flow problem with body fitted nonorthogonal
skewed grid mesh will be presented. By considering a wide range of
skew angles, the efficiency of the numerical method will be tested
for grid skewness especially at extreme skew angles. The numerical
solutions of the flow in a skewed cavity will be presented for
Reynolds number of 100 and 1000 for a wide variety of skew angles
ranging between a=15 ^{°} and a=165 ^{°}
with Da=15 ^{°} increments.
2 NUMERICAL FORMULATION
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
y_{y} w_{x}  y_{x} w_{y} = 
1
Re

( w_{xx} +w_{yy} ) 
 (2) 
where, Re is the Reynolds number, and x and y are
the Cartesian coordinates. We consider the governing NavierStokes
equations in general curvilinear coordinates as the following
(x_{x}^{2} + x_{y}^{2}) y_{xx} + (h_{x}^{2} + h_{y}^{2})y_{hh} + (x_{xx} + x_{yy}) y_{x} + (h_{xx} +h_{yy}) y_{h} + 2 (x_{x} h_{x} + x_{y} h_{y})y_{xh} = w 
 (3) 
(x_{x} h_{y}) y_{h} w_{x}  (x_{x} h_{y}) y_{x}w_{h} = 
1
Re


æ è

(x_{x}^{2} + x_{y}^{2}) w_{xx}+ (h_{x}^{2} + h_{y}^{2}) w_{hh} 

+ (x_{xx} + x_{yy}) w_{x} + (h_{xx} + h_{yy})w_{h} + 2 (x_{x} h_{x} + x_{y} h_{y}) w_{xh} 
ö ø


 (4) 
Following Erturk et al. [ 6], first pseudo time
derivatives are assigned to streamfunction and vorticity equations
and using an implicit Euler time step for these pseudo time
derivatives, the finite difference formulations in operator notation
become the following
æ è

1  Dt (x_{x}^{2} + x_{y}^{2}) d_{xx}  Dt(h_{x}^{2} + h_{y}^{2}) d_{hh}  Dt (x_{xx} +x_{yy}) d_{x}  Dt (h_{xx} + h_{yy}) d_{h} 
ö ø

y^{n+1} 

= y^{n} + 2 Dt (x_{x} h_{x} + x_{y} h_{y})y^{n}_{xh} + Dt w^{n} 

 (5) 
æ è

1  
Dt
Re

(x_{x}^{2} + x_{y}^{2}) d_{xx}  
Dt
Re

(h_{x}^{2} + h_{y}^{2}) d_{hh}  
Dt
Re

(x_{xx} + x_{yy}) d_{x}  
Dt
Re

(h_{xx} + h_{yy}) d_{h} 


+ Dt (x_{x} h_{y}) y_{h}^{n} d_{x}  Dt (x_{x}h_{y}) y_{x}^{n} d_{h} 
ö ø

w^{n+1} 

= w^{n}+ 
2 Dt
Re

(x_{x} h_{x} + x_{y} h_{y})w^{n}_{xh} 

 (6) 
where d_{xx} and d_{hh} denote the
second order finite difference operators, and similarly
d_{x} and d_{h} denote the first order finite
difference operators in x and hdirection respectively.
The equations above are in implicit form and require the solution of
a large matrix at every pseudo time iteration which is
computationally inefficient. Instead these equations are spatially
factorized such that
æ è

1  Dt (x_{x}^{2} + x_{y}^{2}) d_{xx}  Dt(x_{xx} + x_{yy}) d_{x} 
ö ø


æ è

1  Dt (h_{x}^{2}+ h_{y}^{2}) d_{hh}  Dt (h_{xx} + h_{yy})d_{h} 
ö ø

y^{n+1} 

= y^{n} + 2 Dt (x_{x} h_{x} + x_{y} h_{y})y^{n}_{xh} + Dt w^{n} 

 (7) 
æ è

1  
Dt
Re

(x_{x}^{2} + x_{y}^{2}) d_{xx}  
Dt
Re

(x_{xx} + x_{yy}) d_{x} + Dt(x_{x} h_{y}) y_{h}^{n} d_{x} 
ö ø



× 
æ è

1  
Dt
Re

(h_{x}^{2} + h_{y}^{2})d_{hh}  
Dt
Re

(h_{xx} + h_{yy})d_{h}  Dt (x_{x} h_{y}) y_{x}^{n} d_{h} 
ö ø

w^{n+1} 

= w^{n}+ 
2 Dt
Re

(x_{x} h_{x} + x_{y} h_{y})w^{n}_{xh} 

 (8) 
The advantage of these equations are that each equation require the
solution of a tridiagonal systems that can be solved very
efficiently using the Thomas algorithm. It can be shown that
approximate factorization introduces additional second order terms
( O( Dt ^{2})) in these equations. In order for the
equations to have the correct physical representation, to cancel out
the second order terms due to factorization the same terms are added
to the right hand side of the equations. The reader is referred to
Erturk et al. [ 6] for more details of the numerical
method. The final form of the equations take the following form
æ è

1  Dt (x_{x}^{2} + x_{y}^{2}) d_{xx}  Dt(x_{xx} + x_{yy}) d_{x} 
ö ø


æ è

1  Dt (h_{x}^{2}+ h_{y}^{2}) d_{hh}  Dt (h_{xx} + h_{yy})d_{h} 
ö ø

y^{n+1} 


= y^{n} + 2 Dt (x_{x} h_{x} + x_{y} h_{y})y^{n}_{xh} + Dt w^{n} 

æ è

Dt (x_{x}^{2} + x_{y}^{2}) d_{xx} + Dt(x_{xx} + x_{yy}) d_{x} 
ö ø


æ è

Dt (h_{x}^{2} +h_{y}^{2}) d_{hh} + Dt (h_{xx} + h_{yy})d_{h} 
ö ø

y^{n} 

 (9) 
æ è

1  
Dt
Re

(x_{x}^{2} + x_{y}^{2}) d_{xx}  
Dt
Re

(x_{xx} + x_{yy}) d_{x} + Dt(x_{x} h_{y}) y_{h}^{n} d_{x} 
ö ø



× 
æ è

1  
Dt
Re

(h_{x}^{2} + h_{y}^{2})d_{hh}  
Dt
Re

(h_{xx} + h_{yy})d_{h}  Dt (x_{x} h_{y}) y_{x}^{n} d_{h} 
ö ø

w^{n+1} 


= w^{n}+ 
2 Dt
Re

(x_{x} h_{x} + x_{y} h_{y})w^{n}_{xh} 


æ è


Dt
Re

(x_{x}^{2} + x_{y}^{2}) d_{xx} + 
Dt
Re

(x_{xx} + x_{yy}) d_{x}  Dt(x_{x} h_{y}) y_{h}^{n} d_{x} 
ö ø


× 
æ è


Dt
Re

(h_{x}^{2} + h_{y}^{2})d_{hh} + 
Dt
Re

(h_{xx} + h_{yy})d_{h} + Dt (x_{x} h_{y}) y_{x}^{n} d_{h} 
ö ø

w^{n} 

 (10) 
3 DRIVEN SKEWED CAVITY FLOW
Figure 1 illustrates the schematic view of the benchmark problem,
the driven skewed cavity flow. We will consider the most general
case where the skew angle can be a > 90^{°} or
a < 90^{°}.
In order to calculate the metrics, the grids in the physical domain
are mapped onto orthogonal grids in the computational grids as shown
in Figure 2.
The inverse transformation metrics are calculated using
central differences, as an example [(¶x)/(¶x)] = [(x_{i+1,j}x_{i1,j})/(2 Dx)] = [(2[1/N])/2] = [1/N]. Similarly, inverse transformation metrics are
calculated as the following
x_{x} = 
1
N

, x_{h} = 
cosa
N

, y_{x} = 0 , y_{h} = 
sina
N


 (11) 
where N is the number of grid points. We consider a
(N×N) grid mesh. The determinant of the Jacobian matrix
is found as
 J  = x_{x} y_{h}  x_{h} y_{x} = 
sina
N^{2}


 (12) 
The transformation metrics are defined as
x_{x} = 
1
 J 

y_{h} , x_{y} = 
1
 J 

x_{h} , h_{x} = 
1
 J 

y_{x} , h_{y} = 
1
J 

x_{x} 
 (13) 
Substituting Equations (11) and (12) into (13), the
transformation metrics are obtained as the following
x_{x} = N , x_{y} = 
N cosa
sina

, h_{x} = 0 , h_{y} = 
N
sina


 (14) 
Note that since we use equal grid spacing, the second
order transformation metrics will be all equal to zero such that for
example
x_{xx} = 
¶(x_{x})
¶x

= x_{x} 
¶(x_{x})
¶x

+ h_{x} 
¶(x_{x})
¶h

= 0 
 (15) 
Hence
x_{xx} = x_{yy} = h_{xx} = h_{yy} = 0 
 (16) 
These calculated metrics are substituted into Equations
(9) and (10) and the final form of the numerical equations become as
the following
æ è

1  Dt 
æ è


N^{2}
sin^{2} a


ö ø

d_{xx} 
ö ø


æ è

1  Dt 
æ è


N^{2}
sin^{2} a


ö ø

d_{hh} 
ö ø

y^{n+1} 

= y^{n} + 2 Dt 
æ è


N^{2} cosa
sin^{2}a


ö ø

y^{n}_{xh} + Dt w^{n} + 
æ è

Dt 
æ è


N^{2}
sin^{2} a


ö ø

d_{xx} 
ö ø


æ è

Dt 
æ è


N^{2}
sin^{2} a


ö ø

d_{hh} 
ö ø

y^{n} 

 (17) 
æ è

1  
Dt
Re


æ è


N^{2}
sin^{2} a


ö ø

d_{xx} + Dt 
æ è


N^{2}
sina


ö ø

y_{h}^{n} d_{x} 
ö ø


æ è

1  
Dt
Re


æ è


N^{2}
sin^{2} a


ö ø

d_{hh}  Dt 
æ è


N^{2}
sina


ö ø

y_{x}^{n} d_{h} 
ö ø

w^{n+1} 

= w^{n}+ 
2 Dt
Re


æ è


N^{2} cosa
sin^{2} a


ö ø

w^{n}_{xh} 

+ 
æ è


Dt
Re


æ è


N^{2}
sin^{2} a


ö ø

d_{xx}  Dt 
æ è


N^{2}
sina


ö ø

y_{h}^{n} d_{x} 
ö ø


æ è


Dt
Re


æ è


N^{2}
sin^{2} a


ö ø

d_{hh} + Dt 
æ è


N^{2}
sina


ö ø

y_{x}^{n} d_{h} 
ö ø

w^{n} 
 (18) 
The solution methodology of each of the above two equations,
Equations (17) and (18), involves a twostage timelevel updating.
First the streamfunction equation (17) is solved, and for this, the
variable f is introduced such that
æ è

1  Dt 
æ è


N^{2}
sin^{2} a


ö ø

d_{hh} 
ö ø

y^{n+1} = f 
 (19) 
where

æ è

1  Dt 
æ è


N^{2}
sin^{2} a


ö ø

d_{xx} 
ö ø

f = y^{n} + 2 Dt 
æ è


N^{2}cosa
sin^{2} a


ö ø

y^{n}_{xh} + Dtw^{n} 


+ 
æ è

Dt 
æ è


N^{2}
sin^{2} a


ö ø

d_{xx} 
ö ø


æ è

Dt 
æ è


N^{2}
sin^{2}a


ö ø

d_{hh} 
ö ø

y^{n} 
 (20) 
In Equation (20) f is the only unknown variable. First,
this Equation (20) is solved for f at each grid point. Following
this, the streamfunction (y) variable is advanced into the new
time level using Equation (19).
Then the vorticity equation (18) is solved, and in a similar
fashion, the variable g is introduced such that

æ è

1  
Dt
Re


æ è


N^{2}
sin^{2} a


ö ø

d_{hh}  Dt 
æ è


N^{2}
sina


ö ø

y_{x}^{n} d_{h} 
ö ø

w^{n+1} = g 
 (21) 
where
æ è

1  
Dt
Re


æ è


N^{2}
sin^{2} a


ö ø

d_{xx} + Dt 
æ è


N^{2}
sina


ö ø

y_{h}^{n} d_{x} 
ö ø

g = w^{n}+ 
2 Dt
Re


æ è


N^{2} cosa
sin^{2} a


ö ø

w^{n}_{xh} 

+ 
æ è


Dt
Re


æ è


N^{2}
sin^{2} a


ö ø

d_{xx}  Dt 
æ è


N^{2}
sina


ö ø

y_{h}^{n} d_{x} 
ö ø


æ è


Dt
Re


æ è


N^{2}
sin^{2} a


ö ø

d_{hh} + Dt 
æ è


N^{2}
sina


ö ø

y_{x}^{n} d_{h} 
ö ø

w^{n} 
 (22) 
As with f, first the variable g is determined at every
grid point using Equation (22), then vorticity (w) variable
is advanced into the next time level using Equation (21).
3.1 Boundary Conditions
In the computational domain the velocity components are defined as
u = y_{y} = x_{y} y_{x} + h_{y} y_{h} = 
N cosa
sina

y_{x}+ 
N
sina

y_{h} 
 (23) 
v =  y_{x} =  x_{x} y_{x}  h_{x} y_{h} =  N y_{x} 
 (24) 
On the left wall boundary we have
y_{0,j}=0 , y_{h} 
ê ê

0,j

=0 , y_{hh} 
ê ê

0,j

=0 
 (25) 
where the subscripts 0 and j are the grid indexes.
Also on the left wall, the velocity is zero (u=0 and v=0).
Using Equations (23) and (24) we obtain
and also
y_{xh}= 
¶(y_{x})
¶h


ê ê

0,j

= 0 
 (27) 
Therefore, substituting these into the streamfunction Equation (3)
and using Thom's formula [ 26], on the left wall boundary the
vorticity is calculated as the following
w_{0,j} =  
2 N^{2} y_{1,j}
sin^{2} a


 (28) 
Similarly the vorticity on the right wall (w_{N,j}) and the
vorticity on the bottom wall (w_{i,0}) are defined as the
following
w_{N,j} =  
2 N^{2} y_{N1,j}
sin^{2} a

, w_{i,0} =  
2 N^{2}y_{i,1}
sin^{2} a


 (29) 
On the top wall the uvelocity is equal to u=1. Following the
same procedure, the vorticity on the top wall is found as
w_{i,N} =  
2 N^{2} y_{i,N1}
sin^{2} a

 
2 N
sina


 (30) 
We note that, it is well understood ([ 10],
[ 15], [ 23], [ 29]) that, even though
Thom's method is locally first order accurate, the global solution
obtained using Thom's method preserves second order accuracy.
Therefore in this study, since three point second order central
difference is used inside the skewed cavity and Thom's method is
used at the wall boundary conditions, the presented solutions are
second order accurate.
In the skewed driven cavity flow, the corner points are singular
points for vorticity. We note that due to the skew angle, the
governing equations have cross derivative terms and because of these
cross derivative terms the computational stencil includes 3×3
grid points. Therefore, the solution at the first diagonal grid
points near the corners of the cavity require the vorticity values
at the corner points. For square driven cavity flow Gupta et
al. [ 9] have introduced an explicit asymptotic solution in
the neighborhood of sharp corners. Similarly, Stortkuhl et
al. [ 24] have presented an analytical asymptotic
solutions near the corners of cavity and using finite element
bilinear shape functions they also have presented a singularity
removed boundary condition for vorticity at the corner points as
well as at the wall points. In this study we follow Stortkuhl
et al. [ 24] and use the following expression for
calculating vorticity values at the corners of the skewed cavity
N^{2}
3 sin^{2} a


é
ê
ê
ê
ê
ê
ê
ê
ë 


ù
ú
ú
ú
ú
ú
ú
ú
û 
y 
+ 
1
9


é
ê
ê
ê
ê
ê
ê
ê
ë 


ù
ú
ú
ú
ú
ú
ú
ú
û 
w =  
V N
2 sina


(31) 
where V is the speed of the wall which is equal to 1 for
the upper two corners and it is equal to 0 for the bottom two
corners. The reader is referred to Stortkuhl et al.
[ 24] for details.
4 RESULTS
The steady incompressible flow in a driven skewed cavity is
numerically solved using the described numerical formulation and
boundary conditions. We have considered two Reynolds numbers,
Re=100 and Re=1000. For these two Reynolds numbers we have
varied the skewangle (a) from a=15^{°} to
a=165^{°} with Da=15^{°} increments.
We have solved the introduced problem with a 512×512 grid
mesh, for the two Reynolds number and for all the skew angles
considered.
During the iterations as a measure of the convergence to the steady
state solution, we monitored three error parameters. The first error
parameter, ERR1, is defined as the maximum absolute residual of the
finite difference equations of the steady streamfunction and
vorticity equations in general curvilinear coordinates, Equations
(3) and (4). These are respectively given as
ERR1_{y} = max 
æ è

abs 
æ è


ê ê


N^{2}
sin^{2} a

y_{xx} + 
N^{2}
sin^{2}a

y_{hh}  2 
N^{2} cosa
sin^{2}a

y_{xh} + w 
ê ê

n i,j


ö ø


ö ø

 (32) 
ERR1_{w} = max 
æ è

abs 
æ è


ê ê


1
Re


N^{2}
sin^{2} a

w_{xx} + 
1
Re


N^{2}
sin^{2} a

w_{hh} 

 
N^{2}
sina

y_{h} w_{x} + 
N^{2}
sina

y_{x} w_{h}  
2
Re


N^{2} cosa
sin^{2} a

w_{xh} 
ê ê

n i,j


ö ø


ö ø


 (33) 
The magnitude of ERR1 is an indication of the degree to
which the solution has converged to steady state. In the limit ERR1
would be zero.
The second error parameter, ERR2, is defined as the maximum absolute
difference between an iteration time step in the streamfunction and
vorticity variables. These are respectively given as

ERR2_{y} = max 
æ è

abs 
æ è

y_{i,j}^{n+1}  y_{i,j}^{n} 
ö ø


ö ø



ERR2_{w} = max 
æ è

abs 
æ è

w_{i,j}^{n+1}  w_{i,j}^{n} 
ö ø


ö ø

 (34) 
ERR2 gives an indication of the significant digit of the
streamfunction and vorticity variables are changing between two time
levels.
The third error parameter, ERR3, is similar to ERR2, except that it
is normalized by the representative value at the previous time step.
This then provides an indication of the maximum percent change in
y and w at each iteration step. ERR3 is defined as

ERR3_{y} = max 
æ è

abs 
æ è


y_{i,j}^{n+1}  y_{i,j}^{n}
y_{i,j}^{n}


ö ø


ö ø



ERR3_{w} = max 
æ è

abs 
æ è


w_{i,j}^{n+1}  w_{i,j}^{n}
w_{i,j}^{n}


ö ø


ö ø

 (35) 
In our computations, for every Reynolds numbers and for every skew
angles, we considered that convergence was achieved when both
ERR1_{y} and ERR1_{w} were 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,
ERR2_{y}, was in the order of 10^{17} and for
vorticity, ERR2_{w}, it was in the order of
10^{15}. And also at these convergence levels, between two time
steps the maximum absolute normalized difference in streamfunction,
ERR3_{y}, and in vorticity, ERR3_{w}, was
in the order of 10^{14}, and 10^{13} respectively.
We note that at extreme skew angles, convergence to such low
residuals is necessary. For example, at skew angle
a=15 ^{°} at the bottom left corner, and at skew angle
a=165 ^{°} at the bottom right corner, there appears
progressively smaller counter rotating recirculating regions in
accordance with Moffatt [ 14]. In these recirculating
regions confined in the sharp corner, the value of streamfunction
variable is getting extremely smaller as the size of the
recirculating region gets smaller towards the corner. Therefore, it
is crucial to have convergence to such low residuals especially at
extreme skew angles.
Before solving the skewed driven cavity flow at different skew
angles first we have solved the square driven cavity flow to test
the accuracy of the solution. The square cavity is actually a
special case for skewed cavity and obtained when the skew angle is
chosen as a=90^{°}. For the square driven cavity flow,
the streamfunction and the vorticity values at the center of the
primary vortex and the location of this center are tabulated in
Table 1 for Reynolds numbers of Re=1000, together with results
found in the literature.
The present results are almost exactly the
same with that of Erturk et al. [ 6]. This was
expected since in both studies the same number of grid points were
used and also the spatial accuracy of both the boundary condition
approximations and the solutions were the same. Furthermore the
presented results are in very good agreement with that of highly
accurate spectral solutions of Botella and Peyret [ 4] and
extrapolated solutions of Schreiber and Keller [ 21] and
also fourth order solutions of Erturk and Gokcol [ 7] with
approximately less than 0.18% and 0.14% difference in
streamfunction and vorticity variables respectively. For all the
skew angles considered in this study
(15 ^{°} £ a £ 165 ^{°}) we expect to have
the same level of accuracy we achieved for a=90 ^{°}.
With Li et al. [ 12], Wright and Gaskel [ 30],
Benjamin and Denny [ 2] and Nishida and Satofuka
[ 16] again our solutions compare good.
After validating our solution for a=90 ^{°}, we decided
to validate our solutions at different skew angles. In order to do
this we compare our results with the results found in the
literature. At this point, we would like to note that in the
literature among the studies that have solved the skewed cavity flow
([ 5], [ 17], [ 13], [ 20],
[ 31], [ 28], [ 32], [ 27],
[ 3], [ 18], [ 25], [ 11] and
[ 22]), only Demirdzic et al. [ 5],
Oosterlee et al. [ 17], Shklyar and Arbel
[ 22] and Louaked et al. [ 13] have
presented tabulated results therefore we will mainly compare our
results with those studies.
As mentioned earlier, Demirdzic et al. [ 5] have
presented solutions for skewed cavity for Reynolds number of 100 and
1000 for skewed angles of a=45 ^{°} and
a=30 ^{°}. Figure 3 compares our results of uvelocity
along line AB and vvelocity along line CD with that of Demirdzic et al. [ 5] for Re=100 and 1000 for
a=45 ^{°}, and also Figure 4 compares the same for
a=30 ^{°}. Our results agree excellent with results of
Demirdzic et al. [ 5].
Table 2 compares our results of the minimum and also maximum
streamfunction value and also their location for Reynolds numbers of
100 and 1000 for skew angles of a=30 ^{°} and
a=45 ^{°} with results of Demirdzic et al.
[ 5], Oosterlee et al. [ 17], Louaked
et al. [ 13] and Shklyar and Arbel [ 22]. The
results of this study and the results of Demirdzic et al.
[ 5] and also those of Oosterlee et al.
[ 17], Shklyar and Arbel [ 22] and Louaked
et al. [ 13] agree well with each other, although we
believe that our results are more accurate since in this study a
very fine grid mesh is used.
Figure 5 to Figure 8 show the streamline and also vorticity contours
for Re=100 and Re=1000 for skew angles from
a=15^{°} to a=165^{°} with Da=15^{°} increments. As it is seen from these contour
figures of streamfunction and vorticity, the solutions obtained are
very smooth without any wiggles in the contours even at extreme skew
angles.
We have solved the incompressible flow in a skewed driven cavity
numerically and compared our numerical solution with the solutions
found in the literature for a=90 ^{°}, 30 ^{°} and
45 ^{°}, and good agreement is found. We, then, have presented
solutions for 15 ^{°} £ a £ 165 ^{°}. Since
we could not find solutions in the literature to compare with our
presented solutions other than a=90 ^{°}, 30 ^{°} and
45 ^{°}, in order to demonstrate the accuracy of the numerical
solutions we presented, a good mathematical check would be to check
the continuity of the fluid, as suggested by Aydin and Fenner
[ 1]. We have integrated the uvelocity and vvelocity
profiles along line AB and line CD, passing through the geometric
center of the cavity shown by the red dotted line in Figure 1, in
order to obtain the net volumetric flow rate through these sections.
Through section AB, the volumetric flow rate is Q _{AB}=  òu d y + òv d x  and through section CD it
is Q _{CD}=  òv d x . Since the flow is
incompressible, the net volumetric flow rate passing through these
sections should be equal to zero, Q=0. Using Simpson's rule for
the integration, the volumetric flow rates Q _{AB} and Q _{CD} are
calculated for every skew angle ( a) and every Reynolds number
considered. In order to help quantify the errors, the obtained
volumetric flow rate values are normalized by the absolute total
flow rate through the corresponding section at the considered Re
and a. Hence Q _{AB} is normalized by òu  d y + òv  d x and similarly Q _{CD} is
normalized by òv  d x. Table 3 tabulates the
normalized volumetric flow rates through the considered sections.
We note that in an integration process the numerical errors will add
up, nevertheless, the normalized volumetric flow rate values
tabulated in Table 3 are close to zero, such that they can be
considered as Q_{AB} » Q_{CD} » 0. This
mathematical check on the conservation of the continuity shows that
our numerical solution is indeed very accurate at the considered
skew angles and Reynolds numbers.
We note that, to the authors best knowledge, in the literature there
is not a study that considered the skewed cavity flow at the skew
angles used in the present study other than a=30^{°}
and a=45^{°}. The solutions presented in this study are
unique therefore, for future references, in Table 4 we have
tabulated the minimum and also maximum streamfunction values and
their locations and also the vorticity value at these points for
Reynolds number of 100 and 1000 for all the skew angles considered,
from a=15^{°} to a=165^{°} with Da=15^{°} increments.
In this table the interesting point
is, at Reynolds number of 1000, the strength of vorticity (absolute
value of the vorticity) at the center of the primary vortex decrease
as the skew angle increase from a=15^{°} to
a=90^{°}, having the minimum value at
a=90^{°}. As the skew angle increase further from
a=90^{°} to a=165^{°}, the strength of
vorticity at the center of the primary vortex increase. At this
Reynolds number, Re=1000, the streamfunction value at the center
of the primary vortex also show the same type of behavior, where the
value of the streamfunction start to increase as the skew angle
increase until a=90^{°}, then start to decrease as the
skew angle increase further. However at Reynolds number of 100, the
minimum value of the strength of vorticity and also the maximum
value of the streamfunction at the center of the primary vortex
occur at a=105^{°}. In order to explain this behavior,
we decided to look at the location of the center of the primary
vortex at different skew angles. The strength of vorticity at the
center of the primary vortex is proportional with the vertical
distance between the center of the primary vortex and the top moving
lid. As the center of the primary vortex move away from the top
moving wall, the strength of vorticity at the center should
decrease. Figure 9 shows the vertical distance between the center of
the primary vortex and the top moving lid and as a function of the
skew angle for both Re=100 and 1000.
As this figure show at
Re=100 the eye of the primary vortex is at its farthest position
from the top moving lid at a=105^{°} where as at
Re=1000 the same occurs at a=90^{°}, explaining the
minimum vorticity strength we obtain for the primary vortex at
a=105^{°} for Re=100 and also at
a=90^{°} for Re=1000. For future references, in Table
5 and 6 we have tabulated the uvelocity profiles along line AB
(shown in Figure 1) for Reynolds number of 100 and 1000
respectively, and similarly, in Table 7 and 8 we have tabulated the
vvelocity profiles along line CD (also shown in Figure 1) for
Reynolds number of 100 and 1000 respectively.
5 CONCLUSIONS
In this study the benchmark test case for nonorthogonal grid mesh,
the skewed cavity flow introduced by Demirdzic et al.
[ 5] for skew angles of a=30 ^{°} and
a=45 ^{°}, is reintroduced with a variety of skew
angles. The skewed cavity flow is considered for skew angles ranging
between 15 ^{°} £ a £ 165 ^{°} with Da=15 ^{°} increments, for Re=100 and Re=1000. The
governing NavierStokes equations are considered in most general
form, in general curvilinear coordinates. The nonorthogonal grids
are mapped onto a computational domain. Using the numerical
formulation introduced by Erturk et al. [ 6], fine
grid solutions of streamfunction and vorticity equations are
obtained with very low residuals. The numerical formulation of
Erturk et al. [ 6] have proved to be very effective
on nonorthogonal problems with nonorthogonal grid mesh even at
extreme skew angles. The driven skewed cavity flow problem is a
challenging problem and it can be a perfect benchmark test case for
numerical methods to test performances on nonorthogonal grid
meshes. For future references detailed results are tabulated.
References
 [1]

M. Aydin and R. T. Fenner, Boundary Element Analysis of Driven
Cavity Flow for Low and Moderate Reynolds Numbers, International
Journal for Numerical Methods in Fluids, 37 (2001),
4564.
 [2]

A. S. Benjamin and V. E. Denny, On the Convergence of
Numerical Solutions for 2D Flows in a Cavity at Large Re, Journal
of Computational Physics 33 (1979), 340358.
 [3]

E. Brakkee, P. Wesseling and C. G. M. Kassels, Schwarz Domain
Decomposition for the Incompressible Navier–Stokes Equations in
General Coordinates, International Journal for Numerical Methods
in Fluids 32, (2000) 141173.
 [4]

O. Botella and R. Peyret, Benchmark Spectral Results on the
LidDriven Cavity Flow, Computers and Fluids 27, (1998)
421433.
 [5]

I. Demirdzic, Z. Lilek and M. Peric, Fluid Flow and
Heat Transfer Test Problems for Nonorthogonal Grids: Benchmark
Solutions, International Journal for Numerical Methods in Fluids
15, (1992) 329354.
 [6]

E. Erturk, T. C. Corke and C. Gokcol, Numerical Solutions of
2D Steady Incompressible Driven Cavity Flow at High Reynolds
Numbers, International Journal for Numerical Methods in Fluids
48, (2005) 747774.
 [7]

E. Erturk and C. Gokcol, Fourth Order Compact Formulation of
NavierStokes Equations and Driven Cavity Flow at High Reynolds
Numbers, International Journal for Numerical Methods in Fluids
50, (2006) 421436.
 [8]

E. Erturk, O. M. Haddad and T. C. Corke, Numerical Solutions
of Laminar Incompressible Flow Past Parabolic Bodies at Angles of
Attack, AIAA Journal 42, (2004) 22542265.
 [9]

M. M. Gupta, R. P. Manohar and B. Noble, Nature of Viscous
Flows Near Sharp Corners, Computers and Fluids 9, (1981)
379388.
 [10]

H. Huang and B. R. Wetton, Discrete Compatibility in Finite
Difference Methods for Viscous Incompressible Fluid Flow, Journal
of Computational Physics 126, (1996) 468478.
 [11]

H. Lai and Y. Y. Yan, The Effect of Choosing Dependent
Variables and Cellface Velocities on Convergence of the SIMPLE
Algorithm Using NonOrthogonal Grids, International Journal of
Numerical Methods for Heat & Fluid Flow 11, (2001)
524546.
 [12]

M. Li, T. Tang and B. Fornberg, A Compact ForthOrder Finite
Difference Scheme for the Steady Incompressible NavierStokes
Equations International Journal for Numerical Methods in Fluids
20, (1995) 11371151.
 [13]

M. Louaked, L. Hanich and K. D. Nguyen, An Efficient Finite
Difference Technique For Computing Incompressible Viscous Flows,
International Journal for Numerical Methods in Fluids 25,
(1997) 10571082.
 [14]

H. K. Moffatt, Viscous and resistive eddies near a sharp
corner, Journal of Fluid Mechanics 18, (1963) 118.
 [15]

M. Napolitano, G. Pascazio and L. Quartapelle, A Review of
Vorticity Conditions in the Numerical Solution of the zy Equations, Computers and Fluids 28, (1999) 139185.
 [16]

H. Nishida and N. Satofuka, HigherOrder Solutions of Square
Driven Cavity Flow Using a VariableOrder MultiGrid Method,
International Journal for Numerical Methods in Fluids 34,
(1992) 637653.
 [17]

C. W. Oosterlee, P. Wesseling, A. Segal and E. Brakkee,
Benchmark Solutions for the Incompressible NavierStokes
Equations in General Coordinates on Staggered Grids, International
Journal for Numerical Methods in Fluids 17, (1993)
301321.
 [18]

J. R. Pacheco and R. E. Peck, Nonstaggered BoundaryFitted
Coordinate Method For Free Surface Flows, Numerical Heat Transfer
Part B 37, (2000) 267291.
 [19]

M. Peric, Analysis of PressureVelocity Coupling on
Nonorthogonal Grids, Numerical Heat Transfer Part B 17,
(1990) 6382.
 [20]

D. G. Roychowdhury, S. K. Das and T. Sundararajan, An
Efficient Solution Method for Incompressible NS Equations Using
NonOrthogonal Collocated Grid, International Journal for Numerical
Methods in Engineering 45, (1999) 741763.
 [21]

R. Schreiber and H. B. Keller, Driven Cavity Flows by
Efficient Numerical Techniques, Journal of Computational Physics
49, (1983) 310333.
 [22]

A. Shklyar and A. Arbel, Numerical Method for Calculation of
the Incompressible Flow in General Curvilinear Coordinates With
Double Staggered Grid, International Journal for Numerical Methods
in Fluids 41, (2003) 12731294.
 [23]

W. F. Spotz, Accuracy and Performance of Numerical Wall
Boundary Conditions for Steady 2D Incompressible Streamfunction
Vorticity, International Journal for Numerical Methods in Fluids
28, (1998) 737757.
 [24]

T. Stortkuhl, C. Zenger and S. Zimmer, An Asymptotic Solution
for the Singularity at the Angular Point of the Lid Driven Cavity,
International Journal of Numerical Methods for Heat & Fluid Flow
4, (1994) 4759.
 [25]

R. Teigland and I. K. Eliassen, A Multiblock/Multilevel Mesh
Refinement Procedure for CFD Computations, International Journal
for Numerical Methods in Fluids 36, (2001) 519538.
 [26]

A. Thom, The Flow Past Circular Cylinders at Low Speed,
Proceedings of the Royal Society of London Series A 141,
(1933) 651669.
 [27]

P. G. Tucker and Z. Pan, A Cartesian Cut Cell Method for
Incompressible Viscous Flow, Applied Mathematical Modelling
24, (2000) 591606.
 [28]

Y. Wang and S. Komori, On the Improvement of the SIMPLELike
method for Flows with Complex Geometry, Heat and Mass Transfer
36, (2000) 7178.
 [29]

E. Weinan and L. JianGuo, Vorticity Boundary Condition and
Related Issues for Finite Difference Schemes, Journal of
Computational Physics 124, (1996) 368382.
 [30]

N. G. Wright and P. H. Gaskell, An Efficient Multigrid
Approach to Solving Highly Recirculating Flows, Computers and
Fluids 24, (1995) 6379.
 [31]

H. Xu and C. Zhang, Study Of The Effect Of The
NonOrthogonality For NonStaggered Grids—The Results,
International Journal for Numerical Methods in Fluids 29,
(1999) 625644.
 [32]

H. Xu and C. Zhang, Numerical Calculation of Laminar Flows
Using Contravariant Velocity Fluxes, Computers and Fluids
29, (2000) 149177.
File translated from
T_{E}X
by
T_{T}H,
version 3.63.
