# Numerical Implementations for 2D Lid-Driven Cavity Flow in Stream Function Formulation

Authors: K. Poochinapan
Journal: ISRN Applied Mathematics
 (2012)

---

## Abstract

The aim of this paper is to study the properties of approximations to nonlinear terms
of the 2D incompressible Navier-Stokes equations in the stream function formulation (time-dependent
biharmonic equation). The nonlinear convective terms are numerically solved by
using the method with internal iterations, compared to the ones which are solved by using
explicit and implicit schemes (operator splitting scheme Christov and Marinova; (2001)). Using schemes and algorithms,
the steady 2D incompressible flow in a lid-driven cavity is solved up to Reynolds number Re=5000 with second-order spatial accuracy. The schemes are thoroughly validated on
grids with different resolutions. The result of numerical experiments shows that the finite difference
scheme with internal iterations on nonlinearity is more efficient for the high Reynolds number.


---

## Body

## 1. Introduction

There are many finite difference methods for the solution of the Navier-Stokes equations (NSEs) representing incompressible viscous flows. Some of these are schemes utilizing primitive variables (velocity-pressure) [1–4], vorticity-stream function [5–9], and stream function formulation [10, 11]. The practical estimation of any schemes may be different from the theoretical estimation because of the nonlinearity of the NSEs and the implicit characteristic of the continuity condition. There is no single method which is most suitable for all aspects. The solution of the full-nonlinear set of discretized fluid flow equations is usually obtained by solving a sequence of linear equations. The type of linearization which is used can significantly affect the rate of convergence of the iterations to the solution.The primary difficulty in obtaining numerical solutions with primitive variable formulation is that there is no evolution equation for pressure variable. The pressure function can be removed from the equations by means of vorticity-stream function, but an explicit boundary condition for vorticity function is missing. The simplest way to avoid the case is to eliminate the vorticity and turn to a single nonlinear biharmonic equation for the stream function.The condition of stability of an explicit scheme imposes very restrictive limitations on the time increment with the artificial time. For the fourth-order derivative, it isτ≤0.24h4. That is why the problem of constructing an implicit scheme is of great significance. The straightforward implementation of an implicit scheme results in the very large linear system whose solution needs very large computer capacity. One of the most efficient approaches to reduce the computational time without compromising the stability of the scheme is the method of operator splitting [12]. The application of the splitting method to the NSEs in stream function formulation is not straightforward because of the fact that they are not a Cauchy-Kowalevska system and there is no time derivative of the stream function. What is present in the equation is the term (∇2ψ)t which cannot serve as a basis for splitting [1].The new scheme combines the computational simplicity of the implicit scheme in linear terms with semiexplicit approximation of nonlinear terms. The general idea when treating the nonlinear term is to represent it as implicit approximation and then to linearize it and to conduct internal iterations. After the inner iterations converge, one obtains, in fact, the solution for the new time stage. The explicit approximation of nonlinear terms accomplish severe requirement on time step. A single internal iteration on nonlinear terms induces sense of implicit approximation and reduce very severe band on time step. To improve stability properties of explicit approximation of nonlinear terms we require only 3 internal iterations and call such algorithm as a method with internal iteration.As to the question of which method to be used, the answer is that it depends on the type of problems that to be solved. The objective of the present study is to validate the efficiency and accuracy of numerical implementations for the NSEs in term of stream function. The content of this paper is organized as follows. Section2 contains the mathematical formulation of the problem. Section 3 deals with the discretization of the equations and detailed description of numerical algorithms. The methods are applied to the test problem of lid-driven cavity flow up to Re=5,000. Results and discussions of numerical solutions are presented in Section 4, where we make a detailed comparison with available numerical data.


## 2. Mathematical Formulation

Consider a closed 2D domainΩ with a piecewise smooth boundary ∂Ω. The NSEs for a viscous incompressible flow in the terms of stream function ψ are
(2.1)∂Δψ∂t+Re(∂ψ∂y∂Δψ∂x-∂ψ∂x∂Δψ∂y)-Δ2ψ=0,(x,y)∈Ω,t∈[0,T],
where Δ=(∂2/∂x2)+(∂2/∂y2).Boundary and initial conditions are the following:(2.2)ψ=0,∂ψ∂n=b(x,y),(x,y)∈∂Ω,ψ(x,y,0)=ψ0(x,y),(x,y)∈Ω,
where n is a vector normal to domain boundary, and the Reynolds number is defined as Re=UL/ν, where U is the characteristic velocity, L is the characteristic length, and ν is the kinematic viscosity.


## 3. Difference Schemes

The lid-driven flow occupies the region (the cavity)(3.1)Ω={0⩽x⩽1,0⩽y⩽1}.
No-slip boundary conditions take the following form:
(3.2)ψ=∂ψ∂x=0forx=0,x=1,ψ=∂ψ∂y=0fory=0,ψ=0,∂ψ∂y=1fory=1.
For the sake of simplicity, we assume a uniform grid with hx and hy spacing in x- and y-direction, respectively,
(3.3)hx=1M,hy=1N,(xi,yj)=((i-1.5)hx,(j-1.5)hy)
for i=1,…,M+2 and j=1,…,N+2.The mesh is staggered inx-direction on 0.5hx and in y-direction on 0.5hy with respect domain boundaries. The boundary conditions are approximated on two-point stencils with the second order of approximation as follows:
(3.4)ψ2,jn-ψ1,jn=0;ψ1,jn+ψ2,jn=0⇒ψ1,jn=ψ2,jn=0,ψM+2,jn-ψM+1,jn=0;ψM+2,jn+ψM+1,jn=0⇒ψM+2,jn=ψM+1,jn=0,ψi,2n-ψi,1n=0;ψi,1n+ψi,2n=0⇒ψi,1n=ψi,2n=0,ψi,N+2n-ψi,N+1n=hy;ψi,N+2n+ψi,N+1n=0⇒ψi,N+2n=-ψi,N+1n=hy2,
where ψi,jn=ψ(xi,yj,tn) with tn=(n-1)τ,i=1,…,M+2,j=1,…,N+2,n=1,2,….


### 3.1. Explicit Scheme in Nonlinear Terms

The simplest method is explicit in which all nonlinear terms are evaluated using known values attn. The biharmonic operator is approximated by the Crank-Nicolson scheme. Thus, one has linear system to calculate the new value of the unknown at each node. It is clear that each term in (2.1) can be approximated using central-difference operators for all derivatives. The explicit finite difference scheme for (2.1) is
(3.5)Δhψi,jn+1-Δhψi,jnτ+Re(δyψi,jn)(Δhψi+1,jn-Δhψi-1,jn2hx)-Re(δxψi,jn)(Δhψi,j+1n-Δhψi,j-1n2hy)-12(Δh2ψi,jn+1+Δh2ψi,jn)=0,
where Δh=δx2+δy2 and Δh2=δx4+2δx2y2+δy4. The difference formulas to be used are given in Table 1.Table 1
Second-order central difference approximations.
δ
x
ψ
i
,
j
(
ψ
i
+
1
,
j
-
ψ
i
-
1
,
j
)
/
2
h
x
δ
y
ψ
i
,
j
(
ψ
i
,
j
+
1
-
ψ
i
,
j
-
1
)
/
2
h
y
δ
x
2
ψ
i
,
j
(
ψ
i
+
1
,
j
-
2
ψ
i
,
j
+
ψ
i
-
1
,
j
)
/
h
x
2
δ
y
2
ψ
i
,
j
(
ψ
i
,
j
+
1
-
2
ψ
i
,
j
+
ψ
i
,
j
-
1
)
/
h
y
2
δ
x
3
ψ
i
,
j
(
ψ
i
+
2
,
j
-
2
ψ
i
+
1
,
j
+
2
ψ
i
-
1
,
j
-
ψ
i
-
2
,
j
)
/
2
h
x
3
δ
y
3
ψ
i
,
j
(
ψ
i
,
j
+
2
-
2
ψ
i
,
j
+
1
+
2
ψ
i
,
j
-
1
-
ψ
i
,
j
-
2
)
/
2
h
y
3
δ
x
2
y
ψ
i
,
j
(
(
ψ
i
+
1
,
j
+
1
-
ψ
i
+
1
,
j
-
1
)
-
2
(
ψ
i
,
j
+
1
-
ψ
i
,
j
-
1
)
+
(
ψ
i
-
1
,
j
+
1
-
ψ
i
-
1
,
j
-
1
)
)
/
2
h
y
h
x
2
δ
x
y
2
ψ
i
,
j
(
(
ψ
i
+
1
,
j
+
1
-
ψ
i
-
1
,
j
+
1
)
-
2
(
ψ
i
+
1
,
j
-
ψ
i
-
1
,
j
)
+
(
ψ
i
+
1
,
j
-
1
-
ψ
i
-
1
,
j
-
1
)
)
/
2
h
x
h
y
2
δ
x
4
ψ
i
,
j
(
ψ
i
+
2
,
j
-
4
ψ
i
+
1
,
j
+
6
ψ
i
,
j
-
4
ψ
i
-
1
,
j
+
ψ
i
-
2
,
j
)
/
h
x
4
δ
y
4
ψ
i
,
j
(
ψ
i
,
j
+
2
-
4
ψ
i
,
j
+
1
+
6
ψ
i
,
j
-
4
ψ
i
,
j
-
1
+
ψ
i
,
j
-
2
)
/
h
y
4
δ
x
2
y
2
ψ
i
,
j
(
(
ψ
i
+
1
,
j
+
1
-
2
ψ
i
+
1
,
j
+
ψ
i
+
1
,
j
-
1
)
-
2
(
ψ
i
,
j
+
1
-
2
ψ
i
,
j
+
ψ
i
,
j
-
1
)
+
(
ψ
i
-
1
,
j
+
1
-
2
ψ
i
-
1
,
j
+
ψ
i
-
1
,
j
-
1
)
)
/
h
x
2
h
y
2To combine equation (3.5) as a single linear system with band matrix, we introduce the new system of index as follows:
(3.6)m(i,j)=(i-1)(N+2)+j,i=1,…,M+2,j=1,…,N+2.
Each node (i,j) of grid Ωh associates with index m(i,j). It is easy to see that
(3.7)m(i+1,j)=m(i,j)+(N+2),m(i,j+1)=m(i,j)+1,m(i-1,j)=m(i,j)-(N+2),m(i,j-1)=m(i,j)-1.
Now we introduce a new grid function σm, which is defined on the composite grid. Substituting σm in lieu of ψi,j in (3.5), we recast the algebraic system as follows:
(3.8)Δhσmn+1-Δhσmnτ+Re(δyσmn)(Δhσm+(N+2)n-Δhσm-(N+2)n2hx)-Re(δxσmn)(Δhσm+1n-Δhσm-1n2hy)-12(Δh2σmn+1+Δh2σmn)=0.
According to this idea,  can be rewritten as follows:
(3.9)σmn=σm+(N+2)n=0,i=1,j=1,…,N+2,σmn=σm-(N+2)n=0,i=M+2,j=1,…,N+2,σmn=σm+1n=0,i=1,…,M+2,j=1,σmn=-σm-1n=hy2,i=1,…,M+2,j=N+2.
The general sequence of the algorithm is as follows:  (i)
Set the values ofRe,τ,ε,M,N and choose an initial guess σm0=0.(ii)
Considerσmn as known entities and calculate σmn+1 from the system of the algebraic equation (3.8) by the direct method to the solution-banded linear system. We used standard subroutines DGBSV and DGBSVX of LAPACK.(iii)
If the following criterion is satisfied(3.10)maxm|σmn+1-σmn|⩽ε,
then the calculations are terminated. Otherwise the index of iterations is stepped up, n:=n+1 and the computation is returned to step (ii). A disadvantage of explicit method is that convergence can be very slow if a small nonlinear residual is required.


### 3.2. An Operator Splitting Scheme

We use the operator splitting scheme developed in [12, 13]. In the paper [13], a new fully implicit unconditionally stable iterative scheme for the unsteady NSEs in terms of stream function has been developed.


### 3.3. The Method with Internal Iteration in Nonlinear Terms

The general idea when treating the nonlinear term is to represent it as implicit approximation and then to linearize it and to conduct internal iterations. After the inner iterations converge, one obtains, in fact, the solution for the new time stage. The explicit approximation of nonlinear terms accomplished severe requirement on time step. A single internal iteration on nonlinear terms  induced  sense of implicit approximation and reduced very severe band on time step. An improvement in accuracy is achieved at small additional expense. The internal iteration scheme may be viewed as a modification of the explicit scheme (3.5). For the beginning, the internal iteration scheme used only 3 iterations. The reason behind our choices is that a single iteration is equivalent to one step of explicit scheme so, in order to take the economies of the explicit scheme, it is desirable not to perform to many iteration (we choose 3 internal iteration). The method with internal iteration of nonlinear terms can be recasted in the following form:
(3.11)Δhψi,jn+1,k-Δhψi,jnτ+Re(δyψ~i,j)(Δhψ~i+1,j-Δhψ~i-1,j2hx)-Re(δxψ~i,j)(Δhψ~i,j+1-Δhψ~i,j-12hy)-12(Δh2ψi,jn+1,k+Δh2ψi,jn)=0.
Substituting σm in lieu of ψi,j in (3.11), we recast the algebraic system as the following:
(3.12)Δhσmn+1,k-Δhσmnτ+Re(δyσ~m)(Δhσ~m+(N+2)-Δhσ~m-(N+2)2hx)-Re(δxσ~m)(Δhσ~m+1-Δhσ~m-12hy)-12(Δh2σmn+1,k+Δh2σmn)=0.
The general sequence of the algorithm is as follows:  (i)
Set the values ofRe,τ,ε,M,N and choose an initial guess σm0=0.(ii)
Let
(a)
σ~m=σmn then calculate σmn+1,1 from the system of the algebraic equation (3.12) by the direct method to solution-banded linear system. We used standard subroutines DGBSV and DGBSVX of LAPACK,(b)
σ~m=(σmn+1,1+σmn)/2 then calculate σmn+1,2 from the system of the algebraic equation (3.12),(c)
σ~m=(σmn+1,2+σmn)/2 then calculate σmn+1,3 from the system of the algebraic equation  (3.12).(iii)
If the following criterion is satisfied(3.13)maxm|σmn+1,3-σmn|⩽ε,
then the calculations are terminated. Otherwise the index of iterations is stepped up, n:=n+1 or σmn=σmn+1,3 and the computation is returned to step (ii). Note that the linear system for the problem can be written as the following multidiagonal system for the composite grid function σ:
(3.14)Bm-2N-4σm-2N-4n+1,k+Bm-N-2σm-N-2n+1,k+Bm-2σm-2n+1,k+Bm-1σm-1n+1,k+Bmσmn+1,k+Bm+1σm+1n+1,k+Bm+2σm+2n+1,k+Bm+N+2σm+N+2n+1,k+Bm+2N+4σm+2N+4n+1,k=Fm,
where m=1,…,(M+2)(N+2). The matrix of the algebraic system equation (3.14) is banded with 2N+4 lower and upper diagonals.


## 3.1. Explicit Scheme in Nonlinear Terms

The simplest method is explicit in which all nonlinear terms are evaluated using known values attn. The biharmonic operator is approximated by the Crank-Nicolson scheme. Thus, one has linear system to calculate the new value of the unknown at each node. It is clear that each term in (2.1) can be approximated using central-difference operators for all derivatives. The explicit finite difference scheme for (2.1) is
(3.5)Δhψi,jn+1-Δhψi,jnτ+Re(δyψi,jn)(Δhψi+1,jn-Δhψi-1,jn2hx)-Re(δxψi,jn)(Δhψi,j+1n-Δhψi,j-1n2hy)-12(Δh2ψi,jn+1+Δh2ψi,jn)=0,
where Δh=δx2+δy2 and Δh2=δx4+2δx2y2+δy4. The difference formulas to be used are given in Table 1.Table 1
Second-order central difference approximations.
δ
x
ψ
i
,
j
(
ψ
i
+
1
,
j
-
ψ
i
-
1
,
j
)
/
2
h
x
δ
y
ψ
i
,
j
(
ψ
i
,
j
+
1
-
ψ
i
,
j
-
1
)
/
2
h
y
δ
x
2
ψ
i
,
j
(
ψ
i
+
1
,
j
-
2
ψ
i
,
j
+
ψ
i
-
1
,
j
)
/
h
x
2
δ
y
2
ψ
i
,
j
(
ψ
i
,
j
+
1
-
2
ψ
i
,
j
+
ψ
i
,
j
-
1
)
/
h
y
2
δ
x
3
ψ
i
,
j
(
ψ
i
+
2
,
j
-
2
ψ
i
+
1
,
j
+
2
ψ
i
-
1
,
j
-
ψ
i
-
2
,
j
)
/
2
h
x
3
δ
y
3
ψ
i
,
j
(
ψ
i
,
j
+
2
-
2
ψ
i
,
j
+
1
+
2
ψ
i
,
j
-
1
-
ψ
i
,
j
-
2
)
/
2
h
y
3
δ
x
2
y
ψ
i
,
j
(
(
ψ
i
+
1
,
j
+
1
-
ψ
i
+
1
,
j
-
1
)
-
2
(
ψ
i
,
j
+
1
-
ψ
i
,
j
-
1
)
+
(
ψ
i
-
1
,
j
+
1
-
ψ
i
-
1
,
j
-
1
)
)
/
2
h
y
h
x
2
δ
x
y
2
ψ
i
,
j
(
(
ψ
i
+
1
,
j
+
1
-
ψ
i
-
1
,
j
+
1
)
-
2
(
ψ
i
+
1
,
j
-
ψ
i
-
1
,
j
)
+
(
ψ
i
+
1
,
j
-
1
-
ψ
i
-
1
,
j
-
1
)
)
/
2
h
x
h
y
2
δ
x
4
ψ
i
,
j
(
ψ
i
+
2
,
j
-
4
ψ
i
+
1
,
j
+
6
ψ
i
,
j
-
4
ψ
i
-
1
,
j
+
ψ
i
-
2
,
j
)
/
h
x
4
δ
y
4
ψ
i
,
j
(
ψ
i
,
j
+
2
-
4
ψ
i
,
j
+
1
+
6
ψ
i
,
j
-
4
ψ
i
,
j
-
1
+
ψ
i
,
j
-
2
)
/
h
y
4
δ
x
2
y
2
ψ
i
,
j
(
(
ψ
i
+
1
,
j
+
1
-
2
ψ
i
+
1
,
j
+
ψ
i
+
1
,
j
-
1
)
-
2
(
ψ
i
,
j
+
1
-
2
ψ
i
,
j
+
ψ
i
,
j
-
1
)
+
(
ψ
i
-
1
,
j
+
1
-
2
ψ
i
-
1
,
j
+
ψ
i
-
1
,
j
-
1
)
)
/
h
x
2
h
y
2To combine equation (3.5) as a single linear system with band matrix, we introduce the new system of index as follows:
(3.6)m(i,j)=(i-1)(N+2)+j,i=1,…,M+2,j=1,…,N+2.
Each node (i,j) of grid Ωh associates with index m(i,j). It is easy to see that
(3.7)m(i+1,j)=m(i,j)+(N+2),m(i,j+1)=m(i,j)+1,m(i-1,j)=m(i,j)-(N+2),m(i,j-1)=m(i,j)-1.
Now we introduce a new grid function σm, which is defined on the composite grid. Substituting σm in lieu of ψi,j in (3.5), we recast the algebraic system as follows:
(3.8)Δhσmn+1-Δhσmnτ+Re(δyσmn)(Δhσm+(N+2)n-Δhσm-(N+2)n2hx)-Re(δxσmn)(Δhσm+1n-Δhσm-1n2hy)-12(Δh2σmn+1+Δh2σmn)=0.
According to this idea,  can be rewritten as follows:
(3.9)σmn=σm+(N+2)n=0,i=1,j=1,…,N+2,σmn=σm-(N+2)n=0,i=M+2,j=1,…,N+2,σmn=σm+1n=0,i=1,…,M+2,j=1,σmn=-σm-1n=hy2,i=1,…,M+2,j=N+2.
The general sequence of the algorithm is as follows:  (i)
Set the values ofRe,τ,ε,M,N and choose an initial guess σm0=0.(ii)
Considerσmn as known entities and calculate σmn+1 from the system of the algebraic equation (3.8) by the direct method to the solution-banded linear system. We used standard subroutines DGBSV and DGBSVX of LAPACK.(iii)
If the following criterion is satisfied(3.10)maxm|σmn+1-σmn|⩽ε,
then the calculations are terminated. Otherwise the index of iterations is stepped up, n:=n+1 and the computation is returned to step (ii). A disadvantage of explicit method is that convergence can be very slow if a small nonlinear residual is required.


## 3.2. An Operator Splitting Scheme

We use the operator splitting scheme developed in [12, 13]. In the paper [13], a new fully implicit unconditionally stable iterative scheme for the unsteady NSEs in terms of stream function has been developed.


## 3.3. The Method with Internal Iteration in Nonlinear Terms

The general idea when treating the nonlinear term is to represent it as implicit approximation and then to linearize it and to conduct internal iterations. After the inner iterations converge, one obtains, in fact, the solution for the new time stage. The explicit approximation of nonlinear terms accomplished severe requirement on time step. A single internal iteration on nonlinear terms  induced  sense of implicit approximation and reduced very severe band on time step. An improvement in accuracy is achieved at small additional expense. The internal iteration scheme may be viewed as a modification of the explicit scheme (3.5). For the beginning, the internal iteration scheme used only 3 iterations. The reason behind our choices is that a single iteration is equivalent to one step of explicit scheme so, in order to take the economies of the explicit scheme, it is desirable not to perform to many iteration (we choose 3 internal iteration). The method with internal iteration of nonlinear terms can be recasted in the following form:
(3.11)Δhψi,jn+1,k-Δhψi,jnτ+Re(δyψ~i,j)(Δhψ~i+1,j-Δhψ~i-1,j2hx)-Re(δxψ~i,j)(Δhψ~i,j+1-Δhψ~i,j-12hy)-12(Δh2ψi,jn+1,k+Δh2ψi,jn)=0.
Substituting σm in lieu of ψi,j in (3.11), we recast the algebraic system as the following:
(3.12)Δhσmn+1,k-Δhσmnτ+Re(δyσ~m)(Δhσ~m+(N+2)-Δhσ~m-(N+2)2hx)-Re(δxσ~m)(Δhσ~m+1-Δhσ~m-12hy)-12(Δh2σmn+1,k+Δh2σmn)=0.
The general sequence of the algorithm is as follows:  (i)
Set the values ofRe,τ,ε,M,N and choose an initial guess σm0=0.(ii)
Let
(a)
σ~m=σmn then calculate σmn+1,1 from the system of the algebraic equation (3.12) by the direct method to solution-banded linear system. We used standard subroutines DGBSV and DGBSVX of LAPACK,(b)
σ~m=(σmn+1,1+σmn)/2 then calculate σmn+1,2 from the system of the algebraic equation (3.12),(c)
σ~m=(σmn+1,2+σmn)/2 then calculate σmn+1,3 from the system of the algebraic equation  (3.12).(iii)
If the following criterion is satisfied(3.13)maxm|σmn+1,3-σmn|⩽ε,
then the calculations are terminated. Otherwise the index of iterations is stepped up, n:=n+1 or σmn=σmn+1,3 and the computation is returned to step (ii). Note that the linear system for the problem can be written as the following multidiagonal system for the composite grid function σ:
(3.14)Bm-2N-4σm-2N-4n+1,k+Bm-N-2σm-N-2n+1,k+Bm-2σm-2n+1,k+Bm-1σm-1n+1,k+Bmσmn+1,k+Bm+1σm+1n+1,k+Bm+2σm+2n+1,k+Bm+N+2σm+N+2n+1,k+Bm+2N+4σm+2N+4n+1,k=Fm,
where m=1,…,(M+2)(N+2). The matrix of the algebraic system equation (3.14) is banded with 2N+4 lower and upper diagonals.


## 4. Results and Discussion

The results from numerical simulations of the 2D lid-driven cavity flow are presented and compared with published observations. The standard benchmark problem for testing 2D plane NSEs is the driven cavity flow. The fluid contained inside a square cavity is set into motion by the upper wall which is sliding at constant velocity from left to right. The domain is the unit square cavity, and the viscous incompressible flow is governed by the 2D time-dependent incompressible NSEs and driven by the upper wall as seen in Figure1. An unexpected balance of viscous and pressure forces makes the fluid turn into the square cavity. The properties of these forces depend upon the Reynolds number, a hierarchy of eddies develops, the large clockwise-rotating primary, whose location occurs toward the geometric center of the square cavity, and several small eddies: the counter-clockwise rotating secondary eddies, the clockwise rotating tertiary eddies, whose locations occur at the three relevant corners of the square cavity: bottom left, bottom right, and top left.Figure 1
Sketch of the basic features in the 2D flow problem.


### 4.1. Convergence

To evaluate grid dependence on spatial variables the solution is obtained on a sequence of grids with32×32, 62×62, and 122×122 nodes. The finest grid is used as a reference solution. The rate of convergence is computed using two grids according to the formula
(4.1)rate=log2|ψ32×32-ψ122×122||ψ62×62-ψ122×122|.
In Tables 2 and 3   ψmin (appeared at the center of primary eddy) and ψmax (appeared at the center of bottom right-hand side eddy) are computed as minimum and maximum values of stream function, respectively. Parameters ratemin and ratemax correspond to the rate of convergence received from (4.1). Table 2 represents steady state. The solutions are qualified as steady when the absolute error between two time steps is less than 10-8. As shown in Tables 2 and 3 on particular choice  of the parameters, the estimated rate is close to the theoretically predicted second rate of convergence. Figure 2 shows the absolute error εn versus the number of iterations when the largest value of time-step  is used  and
(4.2)εn=maxi,j|ψi,jn-ψi,jn-1|.
It shows that the rate of convergence to solution of the method of internal decreased faster than  that of  explicit and the operator splitting method.Table 2
Convergence rate atRe=100 and τ=2×10-4.
Reference
Grid
ψ
min
ψ
max
rat
e
min
rat
e
max
32
×
32
-
0.1015785
3.787565
×
1
0
-
6
Explicit scheme
62
×
62
-
0.1030120
1.041364
×
1
0
-
5
2.284
2.265
122
×
122
-
0.1033825
1.215392
×
1
0
-
5
32
×
32
-
0.1024849
4.245849
×
1
0
-
6
Operator splitting
62
×
62
-
0.1032607
1.068302
×
1
0
-
5
2.384
2.372
122
×
122
-
0.1034445
1.222466
×
1
0
-
5
32
×
32
-
0.1015785
3.787567
×
1
0
-
6
Internal iteration
62
×
62
-
0.1030120
1.041365
×
1
0
-
5
2.285
2.265
122
×
122
-
0.1033822
1.215392
×
1
0
-
5Table 3
Convergence rate atRe=100, τ=h2 and t=0.1.
Reference
Grid
ψ
min
ψ
max
rat
e
min
rat
e
max
32
×
32
-
0.1014032
3.669264
×
1
0
-
6
Explicit scheme
62
×
62
-
0.1028318
1.041314
×
1
0
-
5
2.251
2.498
122
×
122
-
0.1032119
1.186333
×
1
0
-
5
32
×
32
-
0.1022921
4.128720
×
1
0
-
6
Operator splitting
62
×
62
-
0.1030805
1.014134
×
1
0
-
5
2.343
2.122
122
×
122
-
0.1032741
1.193415
×
1
0
-
5
32
×
32
-
0.1014053
3.673905
×
1
0
-
6
Internal iteration
62
×
62
-
0.1028323
1.014581
×
1
0
-
5
2.250
2.253
122
×
122
-
0.1032120
1.186404
×
1
0
-
5Convergence history: (a)Re=100 on 122 × 122 grid, (b) Re = 2000 on 102 × 102 grid.
(a)
(b)


### 4.2. Time-Step Independence

For fixed values ofhx,hy, and Re, the number of iterations is needed for convergence to steady state depends on τ. For definiteness, we select for presentation a genuinely moderate Re=100 and Re=1000. The results for the number of iterations to steady state are shown in Tables 4 and 5 where the minimal number of iteration for each resolution is shown in bold face. In the case of Re=1000, we minimal numbers of iterations on 52×52 grid are 5276 iterations, 3581 iterations, and 1588 iterations for case explicit, operator splitting, and internal iteration schemes, respectively. We see that the minimal number of iterations of the method of internal iteration is less than other methods.Table 4
The number of iterations to steady solution atRe=100 on  52×52 grid.
τ
Explicit scheme
Operator splitting
Internal iteration
0.0030
—
—
108
0.0010
243
—
243
0.0005
458
458
458
0.0001
1986
1985
1986Table 5
The number of iterations to steady solution atRe=1000 on  52×52 grid.
τ
Explicit scheme
Operator splitting
Internal iteration
0.00008
—
—
1588
0.00005
—
—
2404
0.00003
—
3581
3757
0.00002
5276
5082
5340
0.00001
9606
9174
9667
0.000005
17251
16365
17312Figures3 and 4 show the largest value of time-step for which it is possible to reach a steady solution which corresponds to other investigators. Figure 3 shows that the present schemes converge to correct solutions in all cases of Re≤ 3000 considered. If the grid spacing is too coarse, then the operator splitting scheme will not compute a correct solution for Re>3000. Moreover, the internal iteration method requires slightly larger time-step. An important feature of the internal iteration method is that it allows coarsening of the grid depending on the flow solution and this feature is extremely useful for accuracy in predicting flow fields in region with high Reynolds number and it is compared with the operator splitting scheme.Figure 3
Largest value of time-step on52×52 grid.Figure 4
Largest value of time-step on102×102 grid.The following Figure5 shows that the modification of explicit scheme can reduce the number of iterations for the higher Reynolds number. Figure 6 shows CPU time needed to find solution at instant of time t=0.42 with the largest value of parameter τ. From comparison, it is observed that the internal iteration method seems to be performing better than the other two methods. Computation was carried out on a personal computer with 2.4 GHz CPU and the CPU time per one time step (/one iteration) on the 52×52 grid was about 2.22031, 2.20781, and 6.61875 sec for case explicit, operator splitting, and internal iteration scheme, respectively. The numerical results are demonstrating that the internal iteration is one method to compute the characteristic of flow for the high Reynolds number. We found that the largest time increment of 3 internal iterations gives convergence to solution roughly 3 times as explicit scheme.Figure 5
Minimal number of iteration on52×52 grid at t=0.42.Figure 6
CPU time to reach instant of timet=0.42 on 52×52 grid.


### 4.3. Grid Independence

Figures7 and 8 show the steady state at the secondary vortex at the bottom right corner of the cavity at Re=1000. Figures 3–8 clearly illustrate the grid independence of the computational solution. The minimum and maximum of stream functions are listed in Tables 6–9 which includes a comparison to the high-accuracy solution of previous literatures. Tables 6 and 7 show the values of ψmin, ω and the space location of the primary vortex. As can be seen from Table 6, three finite-difference schemes give very similar quantities on the 102×102 grid. The values of ψmin, ω, and their locations are in accordance with results observed in the literature [1, 4, 14, 15]. Table 7 shows the values of ψmin, ω and the space location of the primary eddy for Re=5000. On the grid for 52×52, the method of operator splitting reaches steady solution, but this solution does not agree with the numerical solutions obtained by other authors ([1, 4, 14, 15]). At the same time, explicit and internal iteration schemes reach “correct” solution. For the primary eddy, results of operator splitting scheme with the 102×102 grid agree within 5% with those obtained in [4] with the 128×128 grid and [1], but the two other schemes differ significantly up to 10%.Table 6
Comparison of the three schemes on the primary vortex atRe=1000: minimum of the stream function and location.
Reference
Grid
ψ
min
ω
x
min
y
min
Explicit
52
×
52
-
0.1096
-
1.911
0.540
0.560
Operator
52
×
52
-
0.1163
-
2.077
0.520
0.580
Internal
52
×
52
-
0.1096
-
1.911
0.540
0.560
Explicit
102
×
102
-
0.1162
-
2.022
0.530
0.560
Operator
102
×
102
-
0.1186
-
2.074
0.530
0.560
Internal
102
×
102
-
0.1162
-
2.022
0.530
0.560
Ref. [14]
160
-
0.1189366
-
2.067753
0.5308
0.5652
Ref. [4]
128
×
128
-
0.11796
-
2.0508
0.53125
0.5625
Ref. [4]
1024
×
1024
-
0.11892
-
2.0674
0.53125
0.56543
Ref. [1]
512
×
512
-
0.116269
—
0.5316
0.5660
Ref. [15]
41
×
41
-
0.117240
-
2.05332
0.5250
0.5750Table 7
Comparison of the three schemes on the primary vortex atRe=5000: minimum of the stream function and location.
Reference
Grid
ψ
min
ω
x
min
y
min
Explicit
52
×
52
-
0.09299
-
1.422
0.540
0.520
Operator
52
×
52
—
—
—
—
Internal
52
×
52
-
0.09315
-
1.431
0.540
0.520
Explicit
102
×
102
-
0.10920
-
1.714
0.520
0.540
Operator
102
×
102
-
0.11490
-
1.848
0.520
0.550
Internal
102
×
102
-
0.10950
-
1.727
0.520
0.540
Ref. [4]
128
×
128
-
0.11731
-
1.8595
0.53125
0.5625
Ref. [4]
1024
×
1024
-
0.12193
-
1.9322
0.53125
0.56543
Ref. [1]
512
×
512
-
0.116120
—
0.5160
0.5357Table 8
Comparison of the three schemes on the bottom right secondary vortex atRe=1000: maximum of the stream function and location.
Reference
Grid
ψ
max
ω
x
max
y
max
Explicit
52
×
52
1.301
×
1
0
-
3
0.9891
0.880
0.120
Operator
52
×
52
1.597
×
1
0
-
3
0.9802
0.860
0.120
Internal
52
×
52
1.302
×
1
0
-
3
0.9891
0.880
0.120
Explicit
102
×
102
1.611
×
1
0
-
3
1.005
0.870
0.110
Operator
102
×
102
1.699
×
1
0
-
3
1.079
0.860
0.110
Internal
102
×
102
1.611
×
1
0
-
3
1.006
0.870
0.110
Ref. [14]
160
1.729717
×
1
0
-
3
1.109789
0.8640
0.1118
Ref. [4]
128
×
128
1.7322
×
1
0
-
3
1.1304
0.85937
0.10938
Ref. [4]
1024
×
1024
1.7297
×
1
0
-
3
1.1120
0.86426
0.11230
Ref. [1]
512
×
512
1.640
×
1
0
-
3
—
0.8651
0.1118
Ref. [15]
41
×
41
1.731
×
1
0
-
3
0.9847
0.8500
0.1250Table 9
Comparison of the three schemes on the bottom right secondary vortex atRe=5000: maximum of the stream function and location.
Reference
Grid
ψ
max
ω
x
max
y
max
Explicit
52
×
52
1.006
×
1
0
-
3
1.485
0.90
0.06
Operator
52
×
52
—
—
—
—
Internal
52
×
52
1.006
×
1
0
-
3
1.483
0.90
0.06
Explicit
102
×
102
2.622
×
1
0
-
3
2.359
0.83
0.07
Operator
102
×
102
2.576
×
1
0
-
3
2.207
0.80
0.08
Internal
102
×
102
2.624
×
1
0
-
3
2.357
0.83
0.07
Ref. [4]
128
×
128
2.9313
×
1
0
-
3
2.7718
0.80469
0.070313
Ref. [4]
1024
×
1024
3.0694
×
1
0
-
3
2.7245
0.80566
0.073242
Ref. [1]
512
×
512
2.890
×
1
0
-
3
—
0.8077
0.0736Secondary vortex at the bottom right corner of the cavity atRe=1000 predicted on 52×52 grid: (a) explicit scheme, (b) operator splitting, and (c) internal iteration.
(a)
(b)
(c)Secondary vortex at the bottom right corner of the cavity atRe=1000 predicted on 102×102 grid: (a) explicit scheme, (b) operator splitting, and (c) internal iteration.
(a)
(b)
(c)Tables8 and 9 report ψmax, ω and the space location of the bottom right secondary eddy. As can be seen from Table 8, the quantities of ψmax agree with those obtained by the other authors. For the bottom right secondary eddy, our results on the 102×102 grid agree within 5% with those obtained by the other authors and in the case of 52×52 grid they differ significantly up to 25%  from those obtained in [15]. Table 7 shows the values of ψmax, ω and the space location of the bottom right secondary eddy for Re=5000. For the bottom right secondary eddy, the quantities ψmax with the 102×102 grid agree within 10% with those obtained in [4] with the 128×128 grid and [1] and differ significantly up to 15%  from those obtained [4] with the 1024×1024 mesh.At the high Reynolds number the flow becomes more complicated, and significantly finer grids are needed in the vicinity of the walls where the dynamics of the flow is dominated by viscosity. Coarse grid in general will not resolve the viscous layer near the boundary. For example, as can be seen on the 52 × 52 grid, operator splitting method cannot get correct solution. Moreover, the other two have value ofψmax differing significantly up to 67% from those obtained [4] with the 1024×1024 mesh.Figures9, 10, and 11 show well-known stream function contours for Re=100,2000,  and 5000 on the 102×102 grid. Whole three algorithms give the similar stream function that cannot be distinguished on the stream line patterns. All of these show secondary vortices at the bottom corners of the cavity as well as at the top left of the cavity for case Re=2000 and 5000. The bottom right tertiary vortices become quite visible for Re=5000. Figures 12 and 13 present comparisons of the u-velocity profiles along a vertical line and the v-velocity profiles along a horizontal line passing through the geometric center of the cavity for Reynolds numbers ranging from 100 to 5000. These profiles are in good agreement with that of [5]. While the data of [5] was obtained using that on the 129×129 grid, our data is obtained using that on the 102×102 grid.Streamline contours forRe=100: (a) explicit scheme, (b) operator splitting, and (c) internal iteration.
(a)
(b)
(c)Streamline contours forRe=2000: (a) explicit scheme, (b) operator splitting, and (c) internal iteration.
(a)
(b)
(c)Streamline contours forRe=5000: (a) explicit scheme, (b) operator splitting, and (c) internal iteration.
(a)
(b)
(c)u-velocity along vertical line through the cavity: (a) explicit scheme, (b) operator splitting, and (c) internal iteration.
(a)
(b)
(c)v-velocity along horizontal line through the cavity: (a) explicit scheme, (b) operator splitting, and (c) internal iteration.
(a)
(b)
(c)


## 4.1. Convergence

To evaluate grid dependence on spatial variables the solution is obtained on a sequence of grids with32×32, 62×62, and 122×122 nodes. The finest grid is used as a reference solution. The rate of convergence is computed using two grids according to the formula
(4.1)rate=log2|ψ32×32-ψ122×122||ψ62×62-ψ122×122|.
In Tables 2 and 3   ψmin (appeared at the center of primary eddy) and ψmax (appeared at the center of bottom right-hand side eddy) are computed as minimum and maximum values of stream function, respectively. Parameters ratemin and ratemax correspond to the rate of convergence received from (4.1). Table 2 represents steady state. The solutions are qualified as steady when the absolute error between two time steps is less than 10-8. As shown in Tables 2 and 3 on particular choice  of the parameters, the estimated rate is close to the theoretically predicted second rate of convergence. Figure 2 shows the absolute error εn versus the number of iterations when the largest value of time-step  is used  and
(4.2)εn=maxi,j|ψi,jn-ψi,jn-1|.
It shows that the rate of convergence to solution of the method of internal decreased faster than  that of  explicit and the operator splitting method.Table 2
Convergence rate atRe=100 and τ=2×10-4.
Reference
Grid
ψ
min
ψ
max
rat
e
min
rat
e
max
32
×
32
-
0.1015785
3.787565
×
1
0
-
6
Explicit scheme
62
×
62
-
0.1030120
1.041364
×
1
0
-
5
2.284
2.265
122
×
122
-
0.1033825
1.215392
×
1
0
-
5
32
×
32
-
0.1024849
4.245849
×
1
0
-
6
Operator splitting
62
×
62
-
0.1032607
1.068302
×
1
0
-
5
2.384
2.372
122
×
122
-
0.1034445
1.222466
×
1
0
-
5
32
×
32
-
0.1015785
3.787567
×
1
0
-
6
Internal iteration
62
×
62
-
0.1030120
1.041365
×
1
0
-
5
2.285
2.265
122
×
122
-
0.1033822
1.215392
×
1
0
-
5Table 3
Convergence rate atRe=100, τ=h2 and t=0.1.
Reference
Grid
ψ
min
ψ
max
rat
e
min
rat
e
max
32
×
32
-
0.1014032
3.669264
×
1
0
-
6
Explicit scheme
62
×
62
-
0.1028318
1.041314
×
1
0
-
5
2.251
2.498
122
×
122
-
0.1032119
1.186333
×
1
0
-
5
32
×
32
-
0.1022921
4.128720
×
1
0
-
6
Operator splitting
62
×
62
-
0.1030805
1.014134
×
1
0
-
5
2.343
2.122
122
×
122
-
0.1032741
1.193415
×
1
0
-
5
32
×
32
-
0.1014053
3.673905
×
1
0
-
6
Internal iteration
62
×
62
-
0.1028323
1.014581
×
1
0
-
5
2.250
2.253
122
×
122
-
0.1032120
1.186404
×
1
0
-
5Convergence history: (a)Re=100 on 122 × 122 grid, (b) Re = 2000 on 102 × 102 grid.
(a)
(b)


## 4.2. Time-Step Independence

For fixed values ofhx,hy, and Re, the number of iterations is needed for convergence to steady state depends on τ. For definiteness, we select for presentation a genuinely moderate Re=100 and Re=1000. The results for the number of iterations to steady state are shown in Tables 4 and 5 where the minimal number of iteration for each resolution is shown in bold face. In the case of Re=1000, we minimal numbers of iterations on 52×52 grid are 5276 iterations, 3581 iterations, and 1588 iterations for case explicit, operator splitting, and internal iteration schemes, respectively. We see that the minimal number of iterations of the method of internal iteration is less than other methods.Table 4
The number of iterations to steady solution atRe=100 on  52×52 grid.
τ
Explicit scheme
Operator splitting
Internal iteration
0.0030
—
—
108
0.0010
243
—
243
0.0005
458
458
458
0.0001
1986
1985
1986Table 5
The number of iterations to steady solution atRe=1000 on  52×52 grid.
τ
Explicit scheme
Operator splitting
Internal iteration
0.00008
—
—
1588
0.00005
—
—
2404
0.00003
—
3581
3757
0.00002
5276
5082
5340
0.00001
9606
9174
9667
0.000005
17251
16365
17312Figures3 and 4 show the largest value of time-step for which it is possible to reach a steady solution which corresponds to other investigators. Figure 3 shows that the present schemes converge to correct solutions in all cases of Re≤ 3000 considered. If the grid spacing is too coarse, then the operator splitting scheme will not compute a correct solution for Re>3000. Moreover, the internal iteration method requires slightly larger time-step. An important feature of the internal iteration method is that it allows coarsening of the grid depending on the flow solution and this feature is extremely useful for accuracy in predicting flow fields in region with high Reynolds number and it is compared with the operator splitting scheme.Figure 3
Largest value of time-step on52×52 grid.Figure 4
Largest value of time-step on102×102 grid.The following Figure5 shows that the modification of explicit scheme can reduce the number of iterations for the higher Reynolds number. Figure 6 shows CPU time needed to find solution at instant of time t=0.42 with the largest value of parameter τ. From comparison, it is observed that the internal iteration method seems to be performing better than the other two methods. Computation was carried out on a personal computer with 2.4 GHz CPU and the CPU time per one time step (/one iteration) on the 52×52 grid was about 2.22031, 2.20781, and 6.61875 sec for case explicit, operator splitting, and internal iteration scheme, respectively. The numerical results are demonstrating that the internal iteration is one method to compute the characteristic of flow for the high Reynolds number. We found that the largest time increment of 3 internal iterations gives convergence to solution roughly 3 times as explicit scheme.Figure 5
Minimal number of iteration on52×52 grid at t=0.42.Figure 6
CPU time to reach instant of timet=0.42 on 52×52 grid.


## 4.3. Grid Independence

Figures7 and 8 show the steady state at the secondary vortex at the bottom right corner of the cavity at Re=1000. Figures 3–8 clearly illustrate the grid independence of the computational solution. The minimum and maximum of stream functions are listed in Tables 6–9 which includes a comparison to the high-accuracy solution of previous literatures. Tables 6 and 7 show the values of ψmin, ω and the space location of the primary vortex. As can be seen from Table 6, three finite-difference schemes give very similar quantities on the 102×102 grid. The values of ψmin, ω, and their locations are in accordance with results observed in the literature [1, 4, 14, 15]. Table 7 shows the values of ψmin, ω and the space location of the primary eddy for Re=5000. On the grid for 52×52, the method of operator splitting reaches steady solution, but this solution does not agree with the numerical solutions obtained by other authors ([1, 4, 14, 15]). At the same time, explicit and internal iteration schemes reach “correct” solution. For the primary eddy, results of operator splitting scheme with the 102×102 grid agree within 5% with those obtained in [4] with the 128×128 grid and [1], but the two other schemes differ significantly up to 10%.Table 6
Comparison of the three schemes on the primary vortex atRe=1000: minimum of the stream function and location.
Reference
Grid
ψ
min
ω
x
min
y
min
Explicit
52
×
52
-
0.1096
-
1.911
0.540
0.560
Operator
52
×
52
-
0.1163
-
2.077
0.520
0.580
Internal
52
×
52
-
0.1096
-
1.911
0.540
0.560
Explicit
102
×
102
-
0.1162
-
2.022
0.530
0.560
Operator
102
×
102
-
0.1186
-
2.074
0.530
0.560
Internal
102
×
102
-
0.1162
-
2.022
0.530
0.560
Ref. [14]
160
-
0.1189366
-
2.067753
0.5308
0.5652
Ref. [4]
128
×
128
-
0.11796
-
2.0508
0.53125
0.5625
Ref. [4]
1024
×
1024
-
0.11892
-
2.0674
0.53125
0.56543
Ref. [1]
512
×
512
-
0.116269
—
0.5316
0.5660
Ref. [15]
41
×
41
-
0.117240
-
2.05332
0.5250
0.5750Table 7
Comparison of the three schemes on the primary vortex atRe=5000: minimum of the stream function and location.
Reference
Grid
ψ
min
ω
x
min
y
min
Explicit
52
×
52
-
0.09299
-
1.422
0.540
0.520
Operator
52
×
52
—
—
—
—
Internal
52
×
52
-
0.09315
-
1.431
0.540
0.520
Explicit
102
×
102
-
0.10920
-
1.714
0.520
0.540
Operator
102
×
102
-
0.11490
-
1.848
0.520
0.550
Internal
102
×
102
-
0.10950
-
1.727
0.520
0.540
Ref. [4]
128
×
128
-
0.11731
-
1.8595
0.53125
0.5625
Ref. [4]
1024
×
1024
-
0.12193
-
1.9322
0.53125
0.56543
Ref. [1]
512
×
512
-
0.116120
—
0.5160
0.5357Table 8
Comparison of the three schemes on the bottom right secondary vortex atRe=1000: maximum of the stream function and location.
Reference
Grid
ψ
max
ω
x
max
y
max
Explicit
52
×
52
1.301
×
1
0
-
3
0.9891
0.880
0.120
Operator
52
×
52
1.597
×
1
0
-
3
0.9802
0.860
0.120
Internal
52
×
52
1.302
×
1
0
-
3
0.9891
0.880
0.120
Explicit
102
×
102
1.611
×
1
0
-
3
1.005
0.870
0.110
Operator
102
×
102
1.699
×
1
0
-
3
1.079
0.860
0.110
Internal
102
×
102
1.611
×
1
0
-
3
1.006
0.870
0.110
Ref. [14]
160
1.729717
×
1
0
-
3
1.109789
0.8640
0.1118
Ref. [4]
128
×
128
1.7322
×
1
0
-
3
1.1304
0.85937
0.10938
Ref. [4]
1024
×
1024
1.7297
×
1
0
-
3
1.1120
0.86426
0.11230
Ref. [1]
512
×
512
1.640
×
1
0
-
3
—
0.8651
0.1118
Ref. [15]
41
×
41
1.731
×
1
0
-
3
0.9847
0.8500
0.1250Table 9
Comparison of the three schemes on the bottom right secondary vortex atRe=5000: maximum of the stream function and location.
Reference
Grid
ψ
max
ω
x
max
y
max
Explicit
52
×
52
1.006
×
1
0
-
3
1.485
0.90
0.06
Operator
52
×
52
—
—
—
—
Internal
52
×
52
1.006
×
1
0
-
3
1.483
0.90
0.06
Explicit
102
×
102
2.622
×
1
0
-
3
2.359
0.83
0.07
Operator
102
×
102
2.576
×
1
0
-
3
2.207
0.80
0.08
Internal
102
×
102
2.624
×
1
0
-
3
2.357
0.83
0.07
Ref. [4]
128
×
128
2.9313
×
1
0
-
3
2.7718
0.80469
0.070313
Ref. [4]
1024
×
1024
3.0694
×
1
0
-
3
2.7245
0.80566
0.073242
Ref. [1]
512
×
512
2.890
×
1
0
-
3
—
0.8077
0.0736Secondary vortex at the bottom right corner of the cavity atRe=1000 predicted on 52×52 grid: (a) explicit scheme, (b) operator splitting, and (c) internal iteration.
(a)
(b)
(c)Secondary vortex at the bottom right corner of the cavity atRe=1000 predicted on 102×102 grid: (a) explicit scheme, (b) operator splitting, and (c) internal iteration.
(a)
(b)
(c)Tables8 and 9 report ψmax, ω and the space location of the bottom right secondary eddy. As can be seen from Table 8, the quantities of ψmax agree with those obtained by the other authors. For the bottom right secondary eddy, our results on the 102×102 grid agree within 5% with those obtained by the other authors and in the case of 52×52 grid they differ significantly up to 25%  from those obtained in [15]. Table 7 shows the values of ψmax, ω and the space location of the bottom right secondary eddy for Re=5000. For the bottom right secondary eddy, the quantities ψmax with the 102×102 grid agree within 10% with those obtained in [4] with the 128×128 grid and [1] and differ significantly up to 15%  from those obtained [4] with the 1024×1024 mesh.At the high Reynolds number the flow becomes more complicated, and significantly finer grids are needed in the vicinity of the walls where the dynamics of the flow is dominated by viscosity. Coarse grid in general will not resolve the viscous layer near the boundary. For example, as can be seen on the 52 × 52 grid, operator splitting method cannot get correct solution. Moreover, the other two have value ofψmax differing significantly up to 67% from those obtained [4] with the 1024×1024 mesh.Figures9, 10, and 11 show well-known stream function contours for Re=100,2000,  and 5000 on the 102×102 grid. Whole three algorithms give the similar stream function that cannot be distinguished on the stream line patterns. All of these show secondary vortices at the bottom corners of the cavity as well as at the top left of the cavity for case Re=2000 and 5000. The bottom right tertiary vortices become quite visible for Re=5000. Figures 12 and 13 present comparisons of the u-velocity profiles along a vertical line and the v-velocity profiles along a horizontal line passing through the geometric center of the cavity for Reynolds numbers ranging from 100 to 5000. These profiles are in good agreement with that of [5]. While the data of [5] was obtained using that on the 129×129 grid, our data is obtained using that on the 102×102 grid.Streamline contours forRe=100: (a) explicit scheme, (b) operator splitting, and (c) internal iteration.
(a)
(b)
(c)Streamline contours forRe=2000: (a) explicit scheme, (b) operator splitting, and (c) internal iteration.
(a)
(b)
(c)Streamline contours forRe=5000: (a) explicit scheme, (b) operator splitting, and (c) internal iteration.
(a)
(b)
(c)u-velocity along vertical line through the cavity: (a) explicit scheme, (b) operator splitting, and (c) internal iteration.
(a)
(b)
(c)v-velocity along horizontal line through the cavity: (a) explicit scheme, (b) operator splitting, and (c) internal iteration.
(a)
(b)
(c)


## 5.   Conclusions

In this paper, we introduce the finite difference scheme to treat nonlinear convective terms in the stream function formulation. After special renumbering of the grid points, the coupled system is formulated as a single system and solved by a LAPACK algorithm. The numerical model is applied to the flow in lid-driven cavity. This study shows that the internal iteration technique performs robustly and allows one to accurately follow the flow patterns. The results are quantitative in good agreement with [1, 4, 14, 15] in common ranges of the main parameters.



---

*Source: 871538-2012-10-31.xml*