Turk. J. Math. Comput. Sci.

6(2017) 10–22
c MatDer
http://matder.dergipark.gov.tr/tjmcs/ MATDER

Numerical Solutions of 2-D Steady Incompressible Flow in a Driven Square

Cavity Using Stream Function-Vorticity Formulation

Vusala Ambethkara,∗ , Manoj Kumarb

a
Department of Mathematics, Faculty of Mathematical Sciences, University of Delhi, Delhi-110007, India.
b
Department of Mathematics, Faculty of Mathematical Sciences, University of Delhi, Delhi-110007, India.

Received: 15-02-2017 • Accepted: 24-06-2017

Abstract. In this study, a streamfunction-vorticity (ψ−ξ) method is suitably used to investigate the problem of 2-D
steady viscous incompressible flow in a driven square cavity with moving top and bottom wall. We used this method
to solve the governing equations along with no-slip and slip wall boundary conditions at low Reynolds number.
A general algorithm was used for this method in order to compute the numerical solutions for streamfunction
ψ, vorticityfunction ξ for low Reynolds numbers Re ≤ 100. We have executed this with the aid of a computer
programme developed and run in C++ compiler. We have also proposed the stability criterion of the numerical
scheme used. Streamline, vorticity and isobar contours have been depicted at different low Reynolds numbers. For
flows at Reynolds number Re=100, our numerical solutions are compared with established steady state results and
excellent agreement is obtained.

2010 AMS Classification: 35Q30, 76D05, 76D07, 76M20

Keywords: Stream function-vorticity formulation, numerical solutions, geometric center of square cavity, u-velocity,
v-velocity, pressure, streamline, isobars, low and moderate Reynolds number.

1. Introduction
The problem of 2-D steady incompressible flow in a square cavity plays an important role in a number of industrial
contexts such as coating, drying technologies, and melt spinning processes. It is also known that, the driven square
cavity flow problem is not only technically important but also of great scientific interest because it displays almost
all fluid mechanical phenomena in the simplest geometrical settings. A number of experimental and numerical studies
have been conducted to investigate the flow field of a driven square cavity flow in the last several decades.

Briefly discussing the related literature: Comparison of finite-difference computations of natural convection has
been investigated by Torrance [25]. Torrance and Rockett [26] have investigated a numerical study of natural convec-
tion in an enclosure with localized heating from below-creeping flow to the onset of laminar instability. Kopecky and
Torrance [14] have investigated about the initiation and structure of axisymmetric eddies in a rotating stream. Bozeman
and Dalton [3] have used the strongly implicit procedure to carried out the numerical solutions for the steady 2-D flow
of a viscous incompressible fluid in rectangular cavities. Gupta and Manohar [11] have discussed the boundary approx-
imations and accuracy in viscous flow computations. Ghia et al. [9] have used the vorticity-streamfunction formulation
for the 2-D incompressible Naiver-Stokes equations to study the effectiveness of the coupled strongly implicit multigrid

*Corresponding author
Email addresses: vambethkar@maths.du.ac.in, vambethkar@gmail.com (V. Ambethkar), manojransiwal@gmail.com (M. Kumar)
V. Ambethkar, M. Kumar, Turk. J. Math. Comput. Sci., 6(2017), 10–22 11

(CSI-MG) method in the determination of high-Re fine-mesh flow solutions. Schreiber and Keller [20] have discussed
on the spurious solutions in driven cavity calculations. Gustafson and Halsi [12] have investigated the unsteady vis-
cous incompressible Navier-Stokes flow in a driven cavity, with taking a particular attention to the formation and
evolution of vortices and eddies. Bruneau and Jouron [4] have obtained the numerical solutions for the steady incom-
pressible Navier-Stokes equations in a 2-D driven cavity in primitive variables by means of the multigrid method. Li
et al. [16] have proposed a compact fourth order finite difference scheme for the steady incompressible Naiver-Stokes
equations. Spotz [22] has investigated on the accuracy and performance of numerical wall boundary conditions for
steady 2-D incompressible Naiver-Stokes equations using streamfunction-vorticity formulation. Demir and Erturk [6]
have investigated numerically the wall driven flow of a viscoelastic fluid in rectangular cavity. Tian and Ge [23] have
used a fourth-order compact finite difference scheme on the nine-point 2-D stencil for solving the steady-state Navier-
Stokes/Boussinesq equations for 2-D incompressible fluid flow with heat transfer using the streamfunction-vorticity
formulation. Zhang [29] has used fourth-order compact finite difference schemes with multigrid techniques to simulate
the 2-D square driven cavity flow at low and high Reynolds numbers. Oztop and Dagtekin [17] have investigated the
steady state 2-D mixed convection problem in a vertical two sided lid-driven differentially heated square cavity. Erturk
et al. [7] have used the streamfunction-vorticity formulation to carry out the numerical computations of 2-D steady
incompressible flow in a driven cavity using a fine uniform grid mesh of 601 × 601. Numerical solutions of 2-D steady
incompressible flow in a driven skewed cavity was investigated by Erturk and Dursun [8]. Wahba [27] has investigated
the numerical simulations for incompressible flow in two sided and four sided lid driven cavities. Perumal and Dass [19]
have investigated the numerical solutions of 2-D unsteady viscous incompressible fluid flow in a two sided lid driven
square cavity with horizontal walls in motion using streamfunction-vorticity formulation. Tian and Yu [24] have used
an efficient compact finite difference approximation, named as five point constant coefficient second-order compact
(5PCC-SOC) scheme, was proposed for the streamfunction formulation of the steady incompressible Navier-Stokes
equations. Laminar mixed convection flow in the presence of magnetic field in a top sided lid-driven cavity heated by a
corner heater was investigated by Oztop et al. [18]. Chamkha and Nada [5] have discussed on the numerical modeling
of steady laminar mixed convection flow in single and double-lid square cavities filled with water- Al2 O3 nanofluid. Yu
and Tian [28] have developed an effective compact finite difference approximation which carries streamfunction and its
first derivatives (velocities) as the unknown variables for the streamfunction-velocity formulation of the steady two di-
mensional incompressible Navier-Stokes equations on non-uniform orthogonal cartesian grids. Numerical experiments
on 1-sided and 2-sided lid driven cavity with aspect ratio = 1 and with inclined wall were performed by Anthony et
al. [1]. Ismael et al. [13] have investigated the numerical solutions of the steady laminar mixed convection inside a
lid-driven square cavity filled with water.
The literature survey with regard to 2-D steady incompressible flow in a driven square cavity revealed that to obtain
high accurate and efficient numerical solutions, we need to depend on compact third and fourth order finite-difference
schemes, higher order Jenson formulation and multigrid techniques which have been developed during the last three
decades.
Although in most of the illustrated development vis-a-vis 2-D steady incompressible flow in a driven square cavity,
they have developed compact new finite difference schemes which are highly accurate in terms of the numerical results.
However, in most of the work cited above, they have not been able to show that they were successful in applying their
methods to test the benchmark problem of driven square cavity with two moving top and bottom walls. Another
lacking in the work cited above is that the behaviour of the velocity and pressure profiles and isobar contours along the
horizontal and vertical line through the geometric center of the square cavity has not been depicted and described.
Due to enormous scope of industrial applications mentioned in the beginning of this introduction and also due to
drawbacks and lacking mentioned above, we have been motivated to undertake the current study.
The main target of this work is to suitably use the streamfunction-vorticity method for solving the problem of 2-D
steady incompressible flow in a driven square cavity with moving top and bottom walls. We have used this method to
solve the governing equations along with no-slip and slip wall boundary conditions. The importance of these appli-
cations can be investigated only by determining numerical solutions of the unknown flow variables, streamfunction,
vorticity function for low Reynolds numbers up to 50.
The design of the current work is as follows: Section 2 provides the physical description, governing equations
of the 2-D steady incompressible flow in a driven cavity along with the boundary conditions for a square cavity,
determination of pressure for viscous flow. Section 3 describes numerical discretization of governing equations and
specification of boundary conditions. Section 4 provides proof of the stability of the numerical scheme used. Section 5
Numerical Solutions of 2-D Steady Incompressible Flow in a Driven Square Cavity Using Stream Function-Vorticity Formulation 12

provides numerical computations and general algorithm for computation of numerical solutions of the flow variables.
Section 6 discusses the numerical results. Section 8 lays out the conclusions of the study.

2. Problem Formulation
2.1. Physical Description. Figure 1 illustrates the geometry of the problem in this study along with the boundary
conditions. ABCD is a square cavity in which incompressible viscous flow is considered. Now, by sliding infinite
long plates lying on top and bottom of the cavity, vorticity along the walls of the cavity is generated. Suppose that all
variables are normalised so that the size of the cavity is 1×1 and the sliding velocities are 1 and -1 in the positive and
negative x-directions respectively.

A u = 1, v=0
D

u=0 u=0
v=0 v=0

B u = -1, v=0 C

Figure 1. Square cavity flow caused by moving plates.

The boundary conditions are defined as no slip on the stationary walls (AB and DC) and as slip on the moving walls
(AD and BC). We have assumed that, at all the four corner points of the square domain, the velocity components (u, v)
and pressure P vanish.

2.2. Governing Equations. In the present investigation, steady 2-D incompressible flow in a driven square cavity with
no-slip and slip wall boundary conditions has been considered. The governing equations of 2-D steady incompressible
viscous flow are the continuity equation, and the two components of the momentum equations. Using the Boussinesq
approximations, the dimensionless governing equations of this problem are expressed as follows:
∂u ∂v
continuity equation + = 0,
∂x ∂y
∂u ∂u ∂P ∂2 u ∂2 u
!
1
x-momentum equation u +v =− + + ,
∂x ∂y ∂x Re ∂x2 ∂y2
∂v ∂v ∂P ∂2 v ∂2 v
!
1
y-momentum equation u +v =− + + ,
∂x ∂y ∂y Re ∂x2 ∂y2
where u, v, P and Re are the velocity components along x and y axis, pressure and Reynolds number respectively. The
no-slip and slip wall boundary conditions are

u = 0, v = 0 on plate DC: u = 0, v = 0

on plate AB: 


on plate AD: u = 1, v = 0 on plate BC: u = −1, v = 0 


V. Ambethkar, M. Kumar, Turk. J. Math. Comput. Sci., 6(2017), 10–22 13

2.3. Determination of Pressure for Viscous Flow. In the streamfunction-vorticity method, to obtain pressure at each
grid point for viscous flow, it is necessary to solve an additional equation for pressure. This equation is derived by
differentiating x-momentum equation with respect to x and y-momentum equation with respect to y and adding them
together. The resulting equation is expressed as follows:

 ∂u ∂v 2
 !
∂u ∂v  ∂v ∂u ∂ ∂u ∂v ∂ ∂u ∂v
! ! ! ! ( !) ( !)
+ −2 +2 +u + +v +
∂x ∂y ∂x ∂y  ∂x ∂y ∂x ∂x ∂y ∂y ∂x ∂y


1 ∂2 ∂u ∂v ∂2 ∂u ∂v
" ! !#
= −O2 P + + + + (2.1)
Re ∂x2 ∂x ∂y ∂y2 ∂x ∂y
From the continuity equation for incompressible flow, we have
∂u ∂v
+ = 0,
∂x ∂y
We obtain from equation (2.1),
∂u ∂v ∂v ∂u
! " ! ! !#
O P=2 2
−
∂x ∂y ∂x ∂y
∂ ∂ψ ∂ ∂ψ ∂ ∂ψ ∂ ∂ψ
" ! ! ! !#
or O P = 2
2
− − −
∂x ∂y ∂y ∂x ∂x ∂x ∂y ∂y
or O2 P = S (2.2)
 2 
 ∂2 ψ ∂2 ψ ∂2 ψ 
! ! !
where S = 2  −
∂x2 ∂y2 ∂x∂y 

Equation (2.2) is known as Poisson equation for pressure. A suitable second-order difference representation for right
hand side of equation (2.2) is given as

ψi+1, j − 2ψi, j + ψi−1, j ψi, j+1 − 2ψi, j + ψi, j−1 ψi+1, j+1 − ψi+1, j−1 − ψi−1, j+1 + ψi−1, j−1 2 
" ! ! ! 
S i, j =2 −  .
∆x2 ∆y2 44x4y

3. Numerical Discretization
Discretization of the governing equations by finite difference method although a well-known technique we have
adopted this technique in the present study due to its compatibility with the regularly shaped geometry, flow in a
square cavity caused by moving plates. With the help of stream function-vorticity (ψ − ξ) formulation [10, p. 121],
the governing equations of unsteady 2-D Navier-Stokes equations and the equation of continuity can be expressed as
follows:
∂2 ψ ∂2 ψ
+ 2 = −ξ, (3.1)
∂x2 ∂y
∂ξ ∂ξ ∂ξ 1 ∂2 ξ ∂2 ξ
!
= −u − v + + , (3.2)
∂t ∂x ∂y Re ∂x2 ∂y2
where Re is the Reynolds number, x and y are the Cartesian coordinates. Essentially, the system is composed of the
Poisson equation for streamfunction (3.1) and the vorticity-transport equation (3.2). It is intended to obtain the steady
state solution from the discretized equations of (3.1) and (3.2) in a time marching fashion.

To obtain the numerical solutions, the coupled equations (3.1) and (3.2) need to be solved in an iterative manner.
Thus, we have used the method developed by Torrance [25] for solving natural convection (Torrance and Rockett [26])
and rotating flow (Kopecky and Torrance [14]) problems, to carry out the numerical computations of the unknown flow
variables: ψ, ξ, u, v, P for the present problem.

Consider a square numerical grid of size 1 × 1 having n horizontal interior grid lines and an equal number of vertical
grid lines as shown in Figure 2.
Numerical Solutions of 2-D Steady Incompressible Flow in a Driven Square Cavity Using Stream Function-Vorticity Formulation 14

u = 1, v=0
A D
j = n+1
(i,n+1)
(i-1,n) (i+1,n)
(i,n)

(i,n-1)

u=0 u=0
v=0 v=0

j=0
i=0 i = n+1
u = -1, v=0
B C

Figure 2. Square cavity flow caused by moving plates.

Using the three point central differences to both derivatives in equation (3.1), then the discretized Poisson equation
(3.1) for ψ can be written as
ψt+1
i+1, j − 2ψi, j + ψi−1, j
t+1 t+1
ψt+1
i, j+1 − 2ψi, j + ψi, j−1
t+1 t+1
+ = −ξi,t+1
j .
∆x2 ∆y2
We choose ∆x = ∆y = h so that the above discretized equation reduces to
ψt+1
i+1, j + ψi−1, j + ψi, j+1 + ψi, j−1 − 4ψi, j = −ξi, j h .
t+1 t+1 t+1 t+1 t+1 2

Now, to solve the vorticity-transport equation (3.2) the computationally stable upwind-differencing scheme is used
to approximate the first two terms on the right-hand side of this equation. We first define u f and ub as the average
x-directional velocities evaluated, respectively, at half a grid point forward and backward from the point (i, j) in x
direction, given as
1 1
u f = (ui+1, j + ui, j ), ub = (ui, j + ui−1, j )
2 2
and, similarly, for v
1 1
vf = (vi, j+1 + vi, j ), vb = (vi, j + vi, j−1 )
2 2
Further defining,
   
ξ1 = u f − u f ξi+1, j + u f + u f − ub + |ub | ξi, j − (ub + |ub |) ξi−1, j ,
   
ξ2 = v f − v f ξi, j+1 + v f + v f − vb + |vb | ξi, j − (vb + |vb |) ξi, j−1 ,
1
the upwind differencing form is preserved. The terms multiplied by Re are approximated by central-differencing
schemes. For them we let
ξ3 = ξi+1, j + ξi−1, j + ξi, j+1 + ξi, j−1 − 4ξi, j
Finally, a forward-differencing scheme is used to approximate the time derivative, so that
 0 
∂ξ
! ξi, j − ξi, j
=
∂t i, j ∆t
in which ∆t is the size of the time increment and a prime is used to denote the value of a variable evaluated at time
t + ∆t. Thus after rearranging terms, (3.2) becomes
0 ∆t  ξ3 
ξi, j = ξi, j + −ξ1 − ξ2 + 2 . (3.3)
2h Re h
V. Ambethkar, M. Kumar, Turk. J. Math. Comput. Sci., 6(2017), 10–22 15

3.1. Specification of Boundary Conditions.

2 ψ0, j −ψ1, j + ∂ψ
 
∂x |0, j ∆x

On the left plate AB: ξ0, j =


∆x 2 



2 ψn+1, j −ψn, j − ∂ψ
  
∂x |n+1, j ∆x

ξn+1, j =

On the right plate DC:


∆x2



2 ψi,n+1 −ψi,n − ∂ψ
 
∂y |i,n+1 ∆y

ξi,n+1 =


On the top plate AD: ∆y2





2 ψi,0 −ψi,1 + ∂ψ
  
∂y |i,0 ∆y


ξi,0 =

On the bottom plate BC:


∆y2

The homogeneous Neumann boundary conditions for pressure are given by

∂P
= 0,
∂n
where n is the normal direction. Elliptic equations with Neumann boundary conditions on all boundaries, such as the
pressure equation (2.2), present an indeterminate problem, as the coefficient matrix of the finite-difference representa-
tion of the equation has one zero eigenvalue [2, p. 269]. Consequently, the resulting system of equations are linearly
dependent and can not be solved uniquely. This can be alleviated by assigning a constant value to pressure at one refer-
ence point in the solution domain. We have assigned P=5 at one reference point in the solution domain. The resulting
equations will be linearly independent and will have a unique solution. Thus, the pressure solution will be off by this
constant, but the pressure gradient, which is the actual term that exists in the equation of motion, will be correctly
calculated.

4. Stability of the Numerical Scheme

In this section, we will prove the stability and convergence of the numerical scheme obtained under the section
numerical discretization based on the criteria suggested in [2, p. 239, [14,25,26]]. To prove the stability of the numerical
scheme given in equation (3.3), we first rewrite it in the form:
0
ξi, j = a1 ξi+1, j + a2 ξi−1, j + a3 ξi, j + a4 ξi, j+1 + a5 ξi, j−1 (4.1)
where
∆t  ∆t
"  #
a1 = − uf − uf + ,
2h Re h2
∆t ∆t
" #
a2 = (ub + |ub |) + ,
2h Re h2
" (  )#
1  4
a3 = 1 − ∆t u f + u f − ub + |ub | + v f + v f − vb + |vb | + ,
2h Re h2
∆t  ∆t
"  #
a4 = − vf − vf + ,
2h Re h2
∆t ∆t
" #
a5 = (vb + |vb |) + .
2h Re h2
Clearly, all the coefficients except for a3 are positive, no matter what the flow direction is. According to the quasilinear
analysis of Lax and Richtmyer [15] the scheme is stable if every coefficient in (4.1) is positive or, equivalently, if
a3 ≥ 0. Which shows that this requirement gives the stability criterion that
"  #−1
1  4
∆t ≤ u f + u f − ub + |ub | + v f + v f − vb + |vb | + .
2h Re h2
Hence, the numerical scheme we have used is conditionally stable. The scheme is consistent as the limiting value of
the local truncation is zero as ∆t, ∆x and ∆y −→ 0. So, by Lax’s equivalence theorem [21, p. 72], the given scheme is
convergent.
Numerical Solutions of 2-D Steady Incompressible Flow in a Driven Square Cavity Using Stream Function-Vorticity Formulation 16

5. Numerical Computations
We obtained the numerical computation of the unknown flow variables ψ, ξ, u, v, P for the present problem with
the aid of a computer programme developed and run in C++ compiler. The input data for the relevant parameters in
the governing equations like Reynolds number Re has been properly chosen incompatible with the present problem
considered.

5.1. The Algorithm for Obtaining Numerical Solutions by Streamfunction-Vorticity Formulation. The algorithm
for obtaining the numerical solutions by streamfunction-vorticity formulation consists of the following steps:
Step 1 Specify the initial values for ξ, ψ, u and v at time t = 0.
Step 2 Solve the vorticity transport equation (3.2) for ξ at each interior grid point at time t + ∆t.
Step 3 Iterate for new ψ values at all points by solving the Poisson equation (3.1) using new ξ values at interior points.
Step 4 Find the velocity components u = ∂ψ ∂ψ
∂y and v = − ∂x .
Step 5 Determine the values of ξ on the boundaries using ψ and ξ values at interior points.
Step 6 Return to step 2 if the solution is not converged.
Step 7 Solve the Poisson equation (2.2) for P using the calculated ψ values at each grid point.

6. Numerical Results and Discussion

We used the method developed by Torrance [25] for solving natural convection (Torrance and Rockett [26]) and
rotating flow (Kopecky and Torrance [14]) problems, to carry out the numerical computations of the unknown flow
variables: ψ, ξ, u, v, P for the present problem. We summarized an algorithm for streamfunction-vorticity formulation
as described under 5.1 above. We used a computer code developed and run in C++ compiler to execute this. To verify
our computer code, the numerical results obtained by the present method were compared with the bench mark results
reported in [9]. It is seen that the results obtained in the present work are in good agreement with those reported in [9]
at low Reynolds number Re=100. This indicates the validity of the numerical code that we developed.

Figure 3 illustrates the variation of u-velocity along the vertical line through the geometric center of the square
cavity for different Reynolds numbers Re=15, 30, 50 and 100. We observed that the absolute value of u-velocity first
decreases, then increase, and finally, decreases in the vicinity of the top wall and the same behavior is observed below
the geometric center for Reynolds number Re=15. We also observed that, the absolute value of u-velocity decreases
in the vicinity of the top wall and the same behaviour is observed below the geometric center of the square cavity for
Reynolds numbers Re=30, 50 and 100.

1.0

0.8

0.6 Re=15
Re=30
Re=50
0.4 Re=100

0.2

u
-0.5 0.0 0.5

Figure 3. u-velocity profiles along the vertical line through geometric center of the square cavity for
Re=15, 30, 50 and 100.
V. Ambethkar, M. Kumar, Turk. J. Math. Comput. Sci., 6(2017), 10–22 17

Figure 4 illustrates the variation of v-velocity along the horizontal line through the geometric center of the square
cavity at different Reynolds numbers Re=15, 30, 50 and 100. We observed that the absolute value of v-velocity first
increases uniformly, and finally, converges to the boundary of the right wall. The behaviour of v-velocity from the left
side wall towards the geometric center of the cavity is that it increases uniformly and converges at the geometric center
of the cavity. We also observed that the absolute value of v-velocity from the left side wall towards the geometric center
of the cavity first increases, then decreases as Reynolds number increases (from Re=15 to Re=100), and also the same
behaviour is observed to the right side of the geometric center.

v
0.4 Re=15
Re=30
Re=50
0.2
Re=100

0.0 x
0.2 0.4 0.6 0.8 1.0
-0.2

-0.4

Figure 4. v-velocity profiles along the horizontal line through geometric center of the square cavity
for Re=15, 30, 50 and 100.

Figure 5 illustrates the vorticity vector along the horizontal line through geometric center of the square cavity for
different Reynolds numbers Re=15, 30, 50 and 100. We observed that the absolute value of vorticity ξ first decreases,
then increases, in between the left wall boundary to the mid of the square domain and the same behaviour is observed
in between the midpoint to the right wall boundary of the square cavity for Reynolds numbers Re=15, 30 and 50. We
also observed that the absolute value of vorticity first decreases, then increases, and finally, decreases in the vicinity of
the right wall and the same behaviour is observed to left of the geometric center for Reynolds number Re=100.

Ξ
4
Re=15
Re=30
2 Re=50
Re=100

x
0.2 0.4 0.6 0.8 1.0

-2

-4

Figure 5. vorticity vector along horizontal line through geometric center of the square cavity for
different Reynolds numbers Re=15, 30, 50 and 100.

Figure 6 illustrates the variation of pressure along the horizontal line through geometric center of the square cavity
at different Reynolds numbers Re=15, 30, 50 and 100. We observed that, the pressure first increases then decreases till
the geometric center of the square cavity and the same behaviour is also observed to right of the geometric center for
Numerical Solutions of 2-D Steady Incompressible Flow in a Driven Square Cavity Using Stream Function-Vorticity Formulation 18

Re=15 and 30. The pressure value first increases, then decreases, and finally increases in between the left wall boundary
and the geometric center of the cavity and the same behaviour is observed in the vicinity of right wall boundary for
Re=50 and 100.

P
5.06 Re=15
Re=30
5.04 Re=50

Re=100
5.02

5.00

4.98

4.96
x
0.0 0.2 0.4 0.6 0.8 1.0

Figure 6. Pressure along the horizontal line through geometric center of the square cavity for differ-
ent Reynolds numbers Re=15, 30, 50 and 100.

Figure 7 illustrates the variation of pressure along the vertical line through geometric center of the square cavity at
different Reynolds numbers Re=15, 30, 50 and 100. We observed that, the pressure first increases, then decreases, again
it increases and finally decreases above the geometric center of the square cavity, and the same behaviour is observed
below the geometric center of the square cavity for Re=15, 30, 50 and 100.

y
1.0

0.8

0.6
Re=15
Re=30
0.4
Re=50
Re=100
0.2

P
4.96 4.98 5.00 5.02 5.04

Figure 7. Pressure along the vertical line through geometric center of the square cavity for different
Reynolds numbers Re=15, 30, 50 and 100.

The streamline contours for different Reynolds numbers Re=15, 30, 50 and 100 have been depicted in Figure 8.
We observed that two small streamline contours are generated, one above the geometric center and the other below the
geometric center in clockwise direction at Reynolds number Re=15 while a single streamline contour is formed at the
geometric center of the cavity for Reynolds numbers Re=30, 50 and 100.
V. Ambethkar, M. Kumar, Turk. J. Math. Comput. Sci., 6(2017), 10–22 19
1.0 1.0

0.8 0.8

0.6 0.6

y
0.4 0.4

0.2 0.2

0.0 0.0
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

x x

(a) Re = 15. (b) Re = 30.

1.0 1.0

0.8 0.8

0.6 0.6
y

y
0.4 0.4

0.2 0.2

0.0 0.0
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

x x

(c) Re = 50. (d) Re = 100.

Figure 8. Streamline contours at different Reynolds numbers :(a) Re=15; (b) Re=30; (c) Re=50; (d)
Re=100.

The vorticity contours for different Reynolds numbers Re=15, 30, 50 and 100 have been depicted in Figure 9. We
observed that, as the Reynolds number increases, the vorticity contours started rotating in the same direction as the
direction of the moving top and bottom walls and finally they rotate around the geometric center of the square cavity
in clockwise sense.

1.0 1.0

-6-10
-8 -8
-6 -6-10
-8 -8
-6
-4 -4
-2 0 -2 0
2 2
4 4
0.8 0.8

0.6 0.6
y

0.4 0.4

0.2 0.2
4 4
2 2
0 -2 0 -2
-4-6 -8-6 -4
-6 -8 -6
-8 -10 -8 -10
0.0 0.0
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

x x

(a) Re = 15. (b) Re = 30.

1.0 1.0
-6 -8 -8 -6 -8 -8
-6
-4 -6
-2 0 -2 -4
0
2
4 24
0.8 0.8

0.6 0.6
y

0.4 0.4

0.2 0.2
4 42
2
0
-2 -4 -2 0
-4
-6 -6
-8 -8 -6 -8 -8 -6
0.0 0.0
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

x x

(c) Re = 50. (d) Re = 100.

Figure 9. Vorticity contours at different Reynolds numbers :(a) Re=15; (b) Re=30; (c) Re=50; (d) Re=100.
Numerical Solutions of 2-D Steady Incompressible Flow in a Driven Square Cavity Using Stream Function-Vorticity Formulation 20

The Isobars for different Reynolds numbers Re=15, 30, 50 and 100 have been depicted in figure 10. Isobars are the
lines which connect points of equal pressure. Figure 10 shows a number of curves with a particular number on it, which
represent the pressure value on that curve. We observed that the maximum and minimum value of pressure are 5.1 and
4.9 respectively, for Re=50. In figure 10, a small blue colour contour shows the smallest value of pressure.
1.0 1.0
5.02 4.98 4.98 4.96 5.02 4.96 4.98
5.04 5.04 5.04 5.08
5.06 4.96 5.02 4.94
5.02 5 5.06
5.02 4.98 5 4.98
4.98
0.8 5 0.8 5.04
5.04 4.96 4.98
4.98
5 4.98
4.96 5
5.04
5 5.02
5 5.02 4.94
5.02 4.94 5
4.96 4.92
0.6 5.04 0.6 5.08 5.06
5 4.96 4.98
5.02 5.02 4.98 5.04 4.96
4.98 5 4.98 5
y

y
5.06
4.98 5.02 4.98
4.98 5.08
0.4 4.96 0.4 4.94 4.96 5.02 5.02
5.02 5.04
5.04 5.04
4.98
5.04 4.96 5.06
4.98 4.96 5.04 4.96
5.04 5.06
0.2 5 0.2 4.98 5
5.04
5 5.02 5
5.02 5.06 4.96
5.02
4.98 5.1
5.04 4.96 4.94 5.02
5.06 5 5.08 4.98
4.98 5 5.04 5.02
0.0 0.0
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

x x

(a) Re = 15. (b) Re = 30.

1.0 1.0
4.95 5.025 4.975 5.04 5.02 4.98
5.05 5.05 4.95 5.06 4.96 4.94 5.08 5.06
4.96
5.025
5.075
4.975 5.02 4.98 5.04 5
5 5.02 5 5.04
0.8 0.8
4.95 5.025 4.975 4.94 4.94
4.975 4.95 5.05 4.96 5.06 4.96 5.06
5.05 4.975 4.98 5.04 4.98 5.02 5
5.025 4.98
5 5.04 4.94
0.6 4.9 0.6 5.06 5.06
5.075 4.95 4.925 5.08 4.94 4.92
4.975 5.025
5.05 4.95 4.96 5 5.04 4.96
5.02
5 4.98

y
5
5 4.98 4.98 5.04
4.975 5.025 5.075
5.05 5.06 5.02 4.96
4.925 4.95 4.94 4.92 5.08 5.06
0.4 0.4 4.94
5.02
5.05 4.96 5.04
5
4.975 4.98 5
5.025 5 5.04 5.02
5.05 4.96
4.95 4.95 5.06
5.06
5.05 4.96 4.94 4.94
0.2 0.2 4.98
5 5.025 5.04 5
4.975 5
5.075 5.04 4.98 5.02
4.96
5.1 5.06
5 5.05 5.025 4.95 5.06 5.08 5.02 4.94 5.04
4.975 5.025 4.96 4.98
0.0 0.0
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

x x

(c) Re = 50. (d) Re = 100.

Figure 10. Isobars at different Reynolds numbers :(a) Re=15; (b) Re=30; (c) Re=50; (d) Re=100.

7. Code Validation
To check the validity of our present computer code used to obtain the numerical results of u-velocity and v-velocity,
we have compared our present results with those benchmark results given by Ghia et al. [9] and it has been found that
they are in good agreement.

1.1
0.2

Ghia et al. [7]

0.15
0.9
Present current code

0.1

0.7

0.05

0
0.5
-0.1 0.1 0.3 0.5 0.7 0.9 1.1

Ghia et al. [7] Present current code -0.05

0.3

-0.1

0.1 -0.15

-0.25 -0.05 0.15 0.35 0.55 0.75

-0.2

-0.1

-0.25

Figure 11. Comparison of the nu- Figure 12. Comparison of the nu-
merical results of u-velocity along merical results of v-velocity along
the vertical line through geomet- the horizontal line through geomet-
ric center of the square cavity for ric center of the square cavity for
Re=100. Re=100.
V. Ambethkar, M. Kumar, Turk. J. Math. Comput. Sci., 6(2017), 10–22 21

8. Conclusions
The main conclusions of this study are as follows:
(i) The absolute value of u-velocity in the vicinity of the top and bottom walls of the square cavity first decreases, then
increases.
(ii) The absolute value of v-velocity in the vicinity of left and right walls of the square cavity first increases, then
decreases.
(iii) The absolute value of vorticity in the vicinity of left and right walls of the square cavity first increases, then
decreases and finally, increases.
(iv) The absolute value of pressure along the horizontal and vertical line through geometric center of the square cavity
first decreases, then increases.
(v) We found that two small streamline contours are generated, one of them is above the geometric center and other is
below the geometric center in clockwise direction at Reynolds number Re=15.
(vi) As the Reynolds number increases, the vorticity contours started rotating in the same direction as the direction of
the moving top and bottom walls and finally they rotate around the geometric center of the square cavity in clockwise
sense.
(vii) The maximum and minimum value of pressure are 5.1 and 4.9 respectively, for Re=50.

Nomenclature
Re Reynolds number
x, y Cartesian Co-ordinates
∆x Grid spacing along x-axis
∆y Grid spacing along y-axis
∆t Time spacing
u x-Component of velocity
v y-Component of velocity
uf Average x-directional velocities evaluated at half a grid point
forward from the point (i, j)
ub Average x-directional velocities evaluated at half a grid point
backward from the point (i, j)
vf Average y-directional velocities evaluated at half a grid point
forward from the point (i, j)
vb Average y-directional velocities evaluated at half a grid point
backward from the point (i, j)
uf Absolute value of u f
|ub | Absolute value of ub
vf Absolute value of v f
|vb | Absolute value of vb
ψt+1
i, j Streamfunction at (i, j) node at time t + 1
ξi,t+1
j Vorticityfunction at (i, j) node at time t + 1
P Dimensionless pressure
t Time level
2
O Laplacin operator
Greek symbols
ψ Streamfunction
ξ Vorticityfunction
Numerical Solutions of 2-D Steady Incompressible Flow in a Driven Square Cavity Using Stream Function-Vorticity Formulation 22

Acknowledgement: The first and second authors acknowledge the financial support received from the University
of Delhi under Research and Development Scheme 2015-16 for providing research grant vide letter no.RC/2015/9677
dated 15th October 2015 as well as Non-Net Fellowship to M.phil. Scholars vide letter no.Ref.NO.Schs/Non-Net/139/EXT-
145/2015-16/254 Dated: 11th April,2015.

References
[1] Anthony, O.O., lyiola, O.O., Miracle, O.O., Numerical simulation of the lid driven cavity flow with inclined walls, International Journal of
Scientific and Engineering Research, 4(5)(2013). 1
[2] Biringen, S., Chow, C. Y., An Introduction to Computational Fluid Mechanics By Examples, John Wiley and Sons, Inc., Hoboken, New Jersey,
2011. 3.1, 4
[3] Bozeman, J.D., Dalton, C., Numerical study of viscous flow in a cavity, Journal of Computational Physics, 12(1973), 348–363. 1
[4] Bruneau, C.H., Jouron, C., An efficient scheme for solving steady incompressible Navier-Stokes equations, Journal of Computational Physics,
89(1990), 389–413. 1
[5] Chamkha, A.J., Nada, E.A., Mixed convection flow in single and double-lid driven square cavities filled with water-Al2 O3 nanofluid: Effect of
viscosity models, European Journal of Mechanics B/Fluids, 36(2012), 82–96. 1
[6] Demir, H., Erturk, V.S., A numerical study of wall driven flow of a viscoelastic fluid in Rectangular cavities, Indian Journal of Pure and applied
Mathematics, 32(10)(2001), 1581–1590. 1
[7] Erturk, E., Corke, T.C., Gokcol, C., Numerical solutions of 2-D steady incompressible driven cavity flow at high Reynolds numbers, Interna-
tional Journal for Numerical Methods in Fluids, 48(2005), 747–774. 1
[8] Erturk, E., Dursun, B., Numerical solutions of 2-D steady incompressible flow in a driven skewed cavity, Journal of Applied Mathematics and
Mechanics, 87(2007), 377–392. 1
[9] Ghia, U., Ghia, K.N., Shin, C.T., High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method, Journal
of Computational Physics, 48(1982), 387–411. 1, 6, 7
[10] Ghoshdastidar, P.S., Computer Simulation of Flow and Heat Transfer, Tata McGraw-Hill Publishing Company Limited, New Delhi, 1998. 3
[11] Gupta, M.M., Manohar, R.P., Boundary approximations and accuracy in viscous flow computations, Journal of Computational Physics,
31(1979), 265–288. 1
[12] Gustafson, K., Halasi, K., Vortex dynamics of cavity flows, Journal of Computational Physics, 4(1986), 279–319. 1
[13] Ismael, M.A., Pop, I., Chamkha, A.J., Mixed convection in a lid-driven square cavity with partial slip, International Journal of Thermal
Sciences, 82(2014), 47–61. 1
[14] Kopecky, R.M., Torrance, K.E., Initiation and structure of axisymmetric eddies in a rotating stream, Computers and Fluids, 1(1973), 289–300.
1, 3, 4, 6
[15] Lax, P.D., Richtmyer, R.D., Survey of the stability of linear finite difference equations, Communications on Pure Applied Mathematics, 9(1956),
267–293. 4
[16] Li, M., Tang, T., Fornberg, B., A compact fourth-order finite difference scheme for the steady incompressible Navier-Stokes equations, Interna-
tional Journal for Numerical Methods in Fluids, 20(1995), 1137–1151. 1
[17] Oztop, H.F., Dagtekin, I., Mixed convection in two sided lid-driven differentially heated square cavity, International Journal of Heat and Mass
Transfer, 47(2004), 1761–1769. 1
[18] Oztop, H.F., Salem, K.A., Pop, I., MHD mixed convection in a lid-driven cavity with corner heater, International Journal of Heat and Mass
Transfer, 54(2011), 3494–3504. 1
[19] Perumal, D.A., Dass, A.K., Simulation of incompressible flows in two-sided lid-driven square cavities, Part I - FDM, CFD Letters, 2(1)(2010).
1
[20] Schreiber, Keller, Spurious solutions in driven cavity calculations, Journal of Computational Physics, 49(1983), 165–172. 1
[21] Smith, G.D., Numerical Solution of Partial Differential Equations: Finite Difference Methods, Oxford University Press, New York, U.S.A.,
1985. 4
[22] Spotz, W.F., Accuracy and performance of numerical wall boundary conditions for steady 2-D incompressible streamfunction vorticity, Inter-
national Journal for Numerical Methods in Fluids, 28(1998), 737–757. 1
[23] Tian, Z., Ge, Y., A fourth-order compact finite difference scheme for the steady stream function-vorticity formulation of the Navier-
Stokes/Boussinesq equations, International Journal for Numerical Methods in Fluids, 41(2003), 495–518. 1
[24] Tian, Z.F., Yu, P.X., An efficient compact difference scheme for solving the streamfunction formulation of the incompressible Navier-Stokes
equations, Journal of Computational Physics, 230(2011), 6404–6419. 1
[25] Torrance, K.E., Comparison of finite-difference computations of natural convection, Journal of Research of the National Bureau of Standards,
72B (1968), 281–301. 1, 3, 4, 6
[26] Torrance, K.E., Rockett, J. A., Numerical study of natural convection in an enclosure with localized heating from below-creeping flow to the
onset of laminar instability, Journal of Fluid Mechanics, 36(1969), 33–54. 1, 3, 4, 6
[27] Wahba, E.M., Multiplicity of states for two sided and four sided lid driven cavity flows, Computers and Fluids, 38(2009), 247–253. 1
[28] Yu, P.X., Tian, Z.F., A compact streamfunction-velocity scheme on nonuniform grids for the 2-D steady incompressible Navier-Stokes equations,
Computers and Mathematics with Applications, 66(2013), 1192–1212. 1
[29] Zhang, J., Numerical simulation of 2-D square driven cavity using fourth-order compact finite difference schemes, Computers and Mathematics
with Applications, 45(2003), 43–52. 1