This is a scientific web page about the two-dimensional steady incompressible flow in a driven cavity. The steady incompressible 2-D Navier-Stokes equations are solved numerically. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations. This is a scientific web page about the two-dimensional steady incompressible flow in a driven cavity. The steady incompressible 2-D Navier-Stokes equations are solved numerically. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations. Driven cavity flow, numerical methods, steady incompressible flow, finite difference, Navier Stokes equations.

Published in : International Journal for Numerical Methods in Fluids (2009)
Int. J. Numer. Meth. Fluids 2009, Vol 60, pp 992-1010

 

Comparison of Wide and Compact Fourth Order Formulations of the Navier-Stokes Equations

Ercan Erturk

Gebze Institute of Technology, Energy Systems Engineering Department,
Gebze, Kocaeli 41400, Turkey

 

Abstract

In this study the numerical performances of wide and compact fourth order formulation of the steady 2-D incompressible Navier-Stokes equations will be investigated and compared with each other. The benchmark driven cavity flow problem will be solved using both wide and compact fourth order formulations and the numerical performances of both formulations will be presented and also the advantages and disadvantages of both formulations will be discussed.

 

1  INTRODUCTION

In Computational Fluid Dynamics (CFD) field of study, most second order accurate finite difference approaches use three point discretizations. Depending on the flow problem, if a higher spatial accuracy is desired then fourth order accuracy can be chosen. In order to achieve fourth order spatial accuracy, standard five point discretization can be used. With five point discretization, wide fourth order formulation of the Navier-Stokes formulation has disadvantages near the boundaries such that when a wide fourth order formulation is used the points adjacent to the boundaries have to be treated specially.

Another way of achieving fourth order spatial accuracy is to use compact fourth order formulations. Compact fourth order formulations provide fourth order spatial accuracy in a 3×3 stencil. The main advantage of this type of formulation is that it could be easily used at the points adjacent to the boundary without a complexity.

In the literature it is possible to find many studies that uses "compact high order finite difference approximations". These studies can be categorized into two depending on the approach used for obtaining the compact high order finite difference approximations. In one category, the compact high order finite difference approximations are basically generalizations of the Pad scheme and the studies of Lele [13] and Visbal and Gaitonde [23] are examples for this type of "compact high order" studies. For this category of the compact high order finite difference approximations Zingg [26,27] compared the numerical efficiency of the compact versus non-compact finite difference methods.

In the second category, basically the Taylor series and the derivatives of the governing equations are used in obtaining the compact high order finite difference approximations. As examples of this type of "compact high order" studies, Spotz and Carey [19] and Li et al. [14], Zhang [25], Dennis and Hudson [5], MacKinnon and Johnson [15] and Gupta et al. [11] have demonstrated the efficiency of compact fourth order formulations of the streamfunction and vorticity equations for uniform grids, and also Ge and Zhang [9] and Spotz and Carey [20] have applied the fourth order compact formulation to nonuniform grids. This study is intended to contribute to this type of "compact high order" literature.

Recently Erturk and Gokcol [7] have presented a new fourth order compact formulation. The uniqueness of this formulation is that the presented compact fourth order formulation of the Navier-Stokes equations are in the same form with the Navier-Stokes equations with additional coefficients. In this formulation if these coefficients are chosen as zero the Navier-Stokes equations are obtained. Therefore if a numerical code is written for the solution of the introduced compact fourth order formulation, the obtained numerical solution is spatially fourth order (O (Dx4)) accurate to the Navier-Stokes equations. In this code if the coefficients are chosen as zero then the solution of the introduced compact fourth order formulation is spatially second order (O (Dx2)) accurate to the Navier-Stokes equations. Therefore with this formulation, using a single equation one can either obtain a second accurate or a fourth order accurate numerical solution of the Navier-Stokes equations just by setting up some coefficients. Moreover, using the introduced formulation one can easily change an existing second order accurate numerical code to provide fourth order accuracy and to do this the only thing have to be done is to add some coefficients into the existing second order accurate numerical code. When these coefficients are added into a second order code to obtain fourth accurate solutions, the code will run slower because of the extra CPU work of evaluating these inserted coefficients. Erturk [8] have shown that this extra CPU work is slightly dependent on the numerical method used. Erturk [8] also showed that for the driven cavity flow when ADI method is used in order to achieve numerical solutions with the same level of convergence (defined as convergence to the same maximum percent change in flow variables), a fourth order accurate numerical code requires 1.8 times more CPU time than a second order accurate numerical code.

Using the Taylor series expansion, one can easily prove that both compact fourth order formulations and five point wide fourth order formulations indeed provide fourth order spatial accurate numerical solutions. In the literature, to the best of authors' knowledge, there is no study that compares the solutions of compact fourth order formulation and the solutions of wide fourth order formulation in terms of accuracy and numerical performance.

The aim of this study is then to solve the steady 2-D incompressible Navier-Stokes equations with fourth order accuracy using both compact fourth order formulation and wide fourth order formulation. First in order to compare the spatial accuracy of both formulations we will consider an analytical test case. Then in order to compare both formulations in terms of numerical performances, as a test case we will consider the benchmark driven cavity flow for different Reynolds numbers with different grid mesh and we will compare the CPU time and number of iterations needed for convergence for both formulations. We will also document the maximum allowable time step that could be used in both formulations as a sign of stability characteristics of the formulations. Finally we will document the advantages and disadvantages of each formulation.

 

2  FOURTH ORDER COMPACT FORMULATION

We consider the steady 2-D incompressible Navier-Stokes equations in streamfunction (y) and vorticity (w) formulation. In non-dimensional form they are given as


2 y

x2
+ 2 y

y2
= - w
(1)

1

Re
2 w

x2
+ 1

Re
2 w

y2
= y

y
w

x
- y

x
w

y
(2)

where x and y are the Cartesian coordinates and Re is the Reynolds number.

If we discretize equations (1) and (2) using three point central differencing, the following finite difference equations provide second order (O (Dx2 ,Dy2)) accurate solutions

yxx + yyy = - w
(3)

1

Re
wxx + 1

Re
wyy = yy wx - yx wy
(4)
where subscripts denote second order (O (Dx2 ,Dy2)) three point central finite difference derivative approximations.

As explained in Erturk and Gokcol [7], with using second order (O (Dx2 ,Dy2)) three point central differencing, the solution of the following finite difference equations are fourth order (O (Dx4 ,Dx2 Dy2 ,Dy4)) accurate to the Navier-Stokes equations (1) and (2).

yxx +yyy = - w+ A
(5)

1

Re
(1+B)wxx + 1

Re
(1+C)wyy = (yy +D)wx -(yx +E)wy +F
(6)
where

A = - Dx2

12
wxx - Dy2

12
wyy -
Dx2

12
+ Dy2

12

yxxyy

B = - Re Dx2

6
yxy + Re2 Dx2

12
yy yy

C = Re Dy2

6
yxy + Re2 Dy2

12
yx yx

D =
Dx2

12
+ Dy2

12

yxxy - Re Dx2

12
yy yxy + Re Dy2

12
yx yyy

E =
Dx2

12
+ Dy2

12

yxyy - Re Dx2

12
yy yxx + Re Dy2

12
yx yxy

F =
Dx2

12
+ Dy2

12

yy wxyy -
Dx2

12
+ Dy2

12

yx wxxy - Dx2

6
yxx wxy

+ Dy2

6
yyy wxy + Re
Dx2

12
+ Dy2

12

yx yy wxy -
Dx2

12
- Dy2

12

wx wy

- 1

Re

Dx2

12
+ Dy2

12

wxxyy
(7)
We note that equations (5) and (6) are in the same form with equations (3) and (4) except with additional coefficients A, B, C, D, E and F. In equations (5) and (6) if the coefficients A, B, C, D, E and F are chosen as zero, we identically obtain equations (3) and (4). Therefore the solution of equations (5) and (6)become second order accurate (O (Dx2 ,Dy2)) to the Navier-Stokes equations when the coefficients are chosen as zero. On the other hand, if the coefficients A, B, C, D, E and F in equations (5) and (6) are calculated as they are defined in equation (7), then the solution of these equations ((5) and (6)) are fourth order accurate (O (Dx4 ,Dx2 Dy2 ,Dy4)) to the Navier-Stokes equations. As demonstrated in Erturk [8], equations (5) and (6) provide easy way of obtaining either second order or fourth order spatial accurate solutions of the Navier-Stokes simply just by using the appropriate coefficients.

In the limit when Dx 0 and Dy 0, the finite difference derivative approximations in equations (5), (6) and (7) could be written as partial derivatives as the following

2 y

x2
+ 2 y

y2
= - w+A
(8)

1

Re
(1+B) 2 w

x2
+ 1

Re
(1+C) 2 w

y2
=
y

y
+ D
w

x
-
y

x
+ E
w

y
+ F
(9)
where

A = - Dx2

12
2 w

x2
- Dy2

12
2 w

y2
-
Dx2

12
+ Dy2

12

4 y

x2 y2

B = - Re Dx2

6
2 y

xy
+ Re2 Dx2

12
y

y
y

y

C = Re Dy2

6
2 y

xy
+ Re2 Dy2

12
y

x
y

x

D =
Dx2

12
+ Dy2

12

3 y

x2 y
- Re Dx2

12
y

y
2 y

xy
+ Re Dy2

12
y

x
2 y

y2

E =
Dx2

12
+ Dy2

12

3 y

xy2
- Re Dx2

12
y

y
2 y

x2
+ Re Dy2

12
y

x
2 y

xy

F =
Dx2

12
+ Dy2

12

y

y
3 w

xy2
-
Dx2

12
+ Dy2

12

y

x
3 w

x2 y
- Dx2

6
2 y

x2
2 w

xy

+ Dy2

6
2 y

y2
2 w

xy
+ Re
Dx2

12
+ Dy2

12

y

x
y

y
2 w

xy
-
Dx2

12
- Dy2

12

w

x
w

y

- 1

Re

Dx2

12
+ Dy2

12

4 w

x2 y2
(10)
As it is described briefly in Erturk and Gokcol [7], the numerical solutions of equations (8) and (9) are fourth order accurate (O (Dx4 ,Dx2 Dy2 ,Dy4)) to the Navier-Stokes equations (1) and (2), strictly provided that second order central discretizations (O (Dx2 ,Dy2)) are used and also strictly provided that a uniform grid mesh with Dx and Dy is used.

One thing is important to note here that, the Navier-Stokes equations (1) and (2) are a subset of equations (8) and (9) such that, when the coefficients A, B, C, D, E and F in equations (8) and (9) are chosen as zero, we obtain the Navier-Stokes equations (1) and (2). Since equations (8) and (9) and the Navier-Stokes equations (1) and (2) are in the same form, any numerical method suitable for the Navier-Stokes equations (1) and (2) can easily be applied to solve equations (8) and (9).

Moreover, if we have an existing code that solve the NS equations (1) and (2) with second order accuracy, we can easily alter the existing code by adding the coefficients A, B, C, D, E and F defined in equation (10) into the code, then the existing second order accurate code will turn into a fourth order accurate code. Therefore a single numerical code for the solution of equations (8) and (9) can provide both second order and fourth order spatial accuracy just by setting some coefficients.

 

3  NUMERICAL SOLUTION METHODOLOGY


We will numerically solve both the NS equations (1) and (2) and the introduced compact fourth order formulation of the Navier-Stokes equations (8) and (9). Both of these equation sets are nonlinear and therefore they need to be solved in an iterative manner. In order to have an iterative numerical algorithm we assign pseudo time derivatives to these equations, such that as an example for equations (8) and (9) we obtain the following

y

t
= 2 y

x2
+ 2 y

y2
+ w-A
(11)

w

t
= 1

Re
(1+B) 2 w

x2
+ 1

Re
(1+C) 2 w

y2
-
y

y
+ D
w

x
+
y

x
+ E
w

y
-F
(12)
We solve these equations (11) and (12) in the pseudo time domain until the solutions converge to steady state.

For the solution method in the pseudo time, we decided to use the Alternating Direction Implicit (ADI) method. ADI method is a very widely used numerical method and in this method a two dimensional problem is solved in two sweeps while solving the equation implicitly in one dimension in each sweep. The reader is referred to [16], [1] and [4] for details.

 

3.1  Compact Fourth Order Solution Methodology

In order to obtain solutions of the fourth order compact formulation, when we apply the ADI method to solve equations (8) and (9), first we solve the streamfunction equation and in order to do this we first solve the following tri-diagonal system in x-direction


1- Dt

2
dxx
yn+1/2 = yn + Dt

2
dyy yn + Dt

2
w- Dt

2
A
(13)
then we solve the following tri-diagonal system in y-direction


1- Dt

2
dyy
yn+1 = yn+1/2 + Dt

2
dxx yn+1/2 + Dt

2
w- Dt

2
A
(14)
After this we solve the vorticity equation and in order to do this we first solve the following tri-diagonal system in x-direction


1- Dt

2
1

Re
(1+B)dxx + Dt

2
(yy +D)dx
wn+1/2 = wn

+ Dt

2
(1+C) 1

Re
dyy wn + Dt

2
(yx +E)dy wn - Dt

2
F
(15)
then we solve the following tri-diagonal system in y-direction


1- Dt

2
1

Re
(1+C)dyy - Dt

2
(yx +E)dy
wn+1 = wn+1/2

+ Dt

2
1

Re
(1+B)dxx wn+1/2 - Dt

2
(yy +D)dxx wn+1/2 - Dt

2
F
(16)
In equations (13), (14), (15) and (16), the superscripts denote the iteration time level, dxx and dyy denote the second derivative three point central finite difference operators, and similarly dx and dy denote the first derivative three point central finite difference operators in x- and y-direction respectively, for example

dx q = qi+1,j -qi-1,j

2Dx
+ O ( Dx2 )

dxx q = qi+1,j - 2 qi,j +qi-1,j

Dx2
+O ( Dx2 )
(17)
where i and j are the grid index and q denote any differentiable quantity.

We note that we use the Thomas algorithm for the solution of these tri-diagonal systems . Using the above numerical procedure defined in equations ((13), (14), (15) and (16)) we solve equations (8) and (9) and obtain spatially fourth order accurate solutions of the Navier-Stokes equations.

 

3.1  Wide Fourth Order Solution Methodology

In order to obtain solutions of the fourth order wide formulation, when we apply the ADI method to solve the NS equations (1) and (2), similarly first we solve the streamfunction equation and in order to do this we first solve the following penta-diagonal system in x-direction


1- Dt

2
dxx
yn+1/2 = yn + Dt

2
dyy yn + Dt

2
w
(18)
then we solve the following penta-diagonal system in y-direction


1- Dt

2
dyy
yn+1 = yn+1/2 + Dt

2
dxx yn+1/2 + Dt

2
w
(19)
Then, similarly, we solve the vorticity equation and in order to do this we first solve the following penta-diagonal system in x-direction


1- Dt

2
1

Re
dxx + Dt

2
yy dx
wn+1/2 = wn + Dt

2
1

Re
dyy wn + Dt

2
yx dy wn
(20)
then we solve the following penta-diagonal system in y-direction


1- Dt

2
1

Re
dyy - Dt

2
yx dy
wn+1 = wn+1/2 + Dt

2
1

Re
dxx wn+1/2 - Dt

2
yy dxx wn+1/2
(21)
We note that, in equations (18), (19), (20) and (21) the second and first derivative operators dxx, dyy, dx and dy denote five point central finite difference operators in x- and y-direction respectively, for example

dx q = -qi+2,j +8qi+1,j -8qi-1,j +qi-2,j

12Dx
+O ( Dx4 )

dxx q = -qi+2,j +16qi+1,j -30qi,j +16qi-1,j -qi-2,j

12Dx2
+O ( Dx4 )
(22)
where i and j are the grid index and q denote any differentiable quantity.

We note that, in numerical solution of the wide formulation we use the modified Thomas algorithm for the solution of the penta-diagonal systems. Following the above numerical procedure defined in equations (18), (19), (20) and (21) we obtain spatially fourth order accurate solutions of the Navier-Stokes equations.

 

3.1  Boundary Conditions

The compact formulations require a 3×3 stencil and therefore this formulation can be easily used at the first set of grid points near the boundaries. However, wide formulations can not be applied directly to the first set of grid points near the boundaries. We are going to compare the numerical solutions of wide and compact formulations, therefore in order to isolate the effect of the boundary conditions on both formulations we would like to use the same boundary conditions in both cases, so that the effect of the boundary conditions will be the same in both cases.

Stortkuhl et al. [22] have presented an analytical asymptotic solution near the corners of the cavity and using finite element bilinear shape functions they also have presented a singularity removed boundary condition for vorticity at the corner grid points as well as at the wall grid points. In this study, we follow Stortkuhl et al. [22] and calculate the vorticity values at the wall grid points (circle points in Figure 1-a) as the following

wb = - 9 V

2 Dh
- 3

2 Dh2
(yd +ye +yf )- 1

8
(2 wa + 2 wc +wd + 4 we +wf )
(23)
where V is the speed of the wall which is equal to 1 for the moving top wall and equal to 0 for the three stationary walls. The grid points a, b, c, d, e and f are shown in Figure 1-b.

For corner points, we again follow Stortkuhl et al. [22] and calculate the vorticity values at the corner grid points (diamond point in Figure 1-a) as the following

wb = - 9 V

2 Dh
- 3

Dh2
yf - 1

4
(2 wc + wf +2 we )
(24)
where again V is equal to 1 for the upper two corners and it is equal to 0 for the bottom two corners. The grid points b, c, e and f are also shown in Figure 1-c. The reader is referred to Stortkuhl et al. [22] for details on the boundary conditions

          

In this study, for the first set of grid points adjacent to the wall we decided to use the Computational Boundary Method (CBM). For details on the CBM, the reader is referred to Huang [12], Yang [24], Gresho [10] and Spotz [21]. Using a fourth order one sided approximation for the velocity on the wall we obtain the following expression

Vw = y

n



0 
= -25 y0 +48 y1 - 36 y2 + 16 y3 - 3 y4

12h
+O (h4 )
(25)
where Vw denotes the wall velocity, n is the normal wall direction, h is the grid spacing, 0 denotes the grid points on the wall, 1 denotes the first set of grid points adjacent to the wall and similarly 2, 3 and 4 denotes the second, third and fourth set of grid points adjacent to the wall respectively. From this relation we can calculate the streamfunction at the grid points near the boundary (rectangle points in Figure 1-a) as the following

y1 = 12 h Vw + 25 y0 + 36 y2 - 16 y3 + 3 y4

48
(26)
In order to calculate the vorticity values at these grid points, we used the streamfunction equation (1). Using a five point wide fourth order formulation for the streamfunction equation, the vorticity at the first set of grid points adjacent to the boundary (rectangle points in Figure 1-a) is calculated as the following

wi,j = -yi+2,j + 16 yi+1,j - 30 yi,j + 16 yi-1,j - yi-2,j

12 Dx2

+ -yi,j+2 + 16 yi,j+1 - 30 yi,j + 16 yi,j-1 - yi,j-2

12 Dy2
+ O ( Dx4 ,Dy4 )
(27)
We note that, this approximation require the streamfunction values at grid points outside the computational domain (hexagon points in Figure 1-a). In order to find the streamfunction values at these grid points we use the following fourth order relation

Vw = y

n



0 
= -3 y-1 - 10 y0 + 18 y1 - 6 y2 + y3

12h
+ O (h4 )
(28)
where -1 denotes the grid point outside the computational domain. From this relation we calculate the value of streamfunction at the exterior grid points (y-1) and use it in equation (27) to calculate the vorticity at the grid points adjacent to the boundary.

At the first diagonal grid points near the corners (triangle point in Figure 1-a), we calculate the streamfunction and vorticity values, using Gauss-Seidel method on the five point fourth order wide formulation of the streamfunction equation (1) and vorticity equations (2) respectively.

We note that in both wide and compact formulations, we used the above described boundary conditions for the grid points on the wall and for the first set of grid points adjacent to the wall. Then we used the wide and compact formulations to obtain the numerical solutions at interior grid points shown in Figure 1-a.

 

4  RESULTS AND DISCUSSIONS

Following the numerical procedure described in the previous section we obtain spatially fourth order solutions of the Navier-Stokes equations obtained using both compact and wide formulations.

For the choice of time steps in solving the governing equations, we decided to use different time steps, Dt, for streamfunction and vorticity equations, as it was also done in Erturk [8]. If we solve the streamfunction and vorticity equations using three point second order accurate discretization using the ADI method, tri-diagonal matrices appear on the implicit LHS of the equations. In streamfunction equation the diagonal elements on the LHS matrices become 1+2Dt/Dh2, and in vorticity equations the diagonal elements on the LHS matrices become 1+2Dt/(ReDh2). We decided to choose different time steps for streamfunction and vorticity equations such that these different time steps would make the diagonal elements the same. Therefore we use Dt=aDh2 for the streamfunction equation and similarly we use Dt=aReDh2 for the vorticity equation, where a is a coefficient we can choose. We solve both the wide and compact formulations using the same time steps.

In this study we would also like to document the maximum allowable time step Dt that we can use in both formulation for convergence as a sign of numerical stability characteristics of both formulations. In order to this, we first solve the wide and compact formulations using the same time steps. Then we increase a, i.e. increase the time step Dt, and solve the wide and compact formulations again. We continue doing this until at some Dt the solution does not converge. Therefore we would document the maximum allowable Dt for convergence for wide and compact formulation for a given Reynolds number and grid mesh.

In this study, in all of the cases considered, we start the iterations from a homogenous initial guess and continue until a certain condition of convergence is satisfied. As the measure of the convergence to the same level, one option is to use the residual of the equations as it was also used by Erturk et al. [6]. However, we are actually solving two different equations, the NS equations (1) and (2) and also the fourth order equations (8) and (9), and we are trying to compare the numerical performance of these two different equations. Therefore the residual of these two different equations may not indicate the convergence to the same level. Alternatively, we can use the difference of the streamfunction and vorticity variables between two time steps as the measure of convergence for the two equations. However, the solutions of the two different equations are slightly different. Since the solutions are different, the difference of the streamfunction and vorticity variables between two time steps may not show the same convergence level for these equations also. Therefore in this study, as convergence criteria we decided to use the difference of the streamfunction and vorticity variables between two time steps normalized by the previous value of the corresponding variable, such that

Residualy = max


yi,jn+1 - yi,jn

yi,jn



Residualw = max


wi,jn+1 -wi,jn

wi,jn


(29)
These residuals provide an indication of the maximum percent change in y and w variables in each iteration step and in this study we let the iterations converge until both Residualy and Residualw are less than 10-8.

Using the Taylor series expansion, mathematically it is straight forward to show that the first leading truncation error term in both fourth order compact formulations and five point fourth order wide formulations is indeed fourth order. In order to document this numerically, we decided to compare the numerical results with a known analytical solution using a model problem introduced by Richards and Crane [17]. Inside the domain (x,y)=([0,1],[0,1]) the following analytical solutions satisfy the Navier-Stokes equations (1) and (2).

y = y-x

Re
- e(x+y)

w = 2e(x+y)
(30)
For this model problem, as the boundary conditions we decided to use the analytical solutions defined in equation (30) at the grid points on the boundaries and at the first set of grid points adjacent to the boundaries. This way we would be able to avoid any effect of numerical boundary condition approximations on the numerical solution and concentrate on the accuracy of the solution of both formulations in the interior domain for a given analytical values at the boundaries.

We note that, in equation (30) the vorticity is independent of Re and the solution of the streamfunction at different Re numbers looks almost the same in a contour plot, therefore for this model problem we only considered the case where Re=1000. We solve this model problem with wide and compact fourth order formulations using different grid mesh with Dh=1/16  ,  1/32  ,  1/64  ,  1/128 where Dx=Dy=Dh. In these fourth order solutions the average absolute difference between the exact solution (yex,wex) given in equation (30) and the numerical solution (yi,j,wi,j), defined as




Ey =



i,j 
| yex (xi ,yj )-yi,j |

N




Ew =



i,j 
| wex (xi ,yj )-wi,j |

N

 

(31)

where N being the total number of grid points, should be proportional to hmwhere the power should be equal to m=4. Since these average absolute differences should have the following form

E=CDhm
(32)
logarithm of both sides give

logE=logC+mlogDh
(33)
In Table 1 we have tabulated the average absolute differences (Ey and Ew) for wide and compact fourth order solutions for different grid mesh and also in Figure 2, these average absolute differences are plotted with respect to the grid spacing in a log-log scale.

                     

In Figure 2 the average absolute differences of the wide and compact formulations are shown by square and circle symbols respectively. In order to have an idea on how these average absolute differences change with respect to the grid spacing, two linear lines with blue and red colors with slopes of 4 are drawn in the same figure. From Figure 2 we can clearly see that both wide and compact fourth order formulations (shown by square and circle in the figure), indeed provide fourth order accurate solutions, such that for both formulations the slope between logE and logDh is very close to m 4.

Using the average absolute differences (Ey and Ew) tabulated in Table 1, since the grid size is decreased by a factor of 2, we can calculate the convergence rate m using the following formula

E(Dh = Dh1 )

E
Dh = Dh1

2

= 2m
(34)
The convergence rate m for both wide and fourth order formulations is also tabulated in Table 1. From this table, we can again see that when the grid spacing is decreased progressively by half , the convergence rates of both formulations is very close to m 4.

Next, we considered the benchmark driven cavity flow problem. For the driven cavity flow we obtain fourth order numerical solutions using both wide and compact fourth order formulations using the boundary conditions described in the previous section and also using both 128×128 and 256×256 grid mesh.

Figures 3, 4 and 5 show the streamfunction contours of the compact and wide formulation solutions for a variety of Reynolds numbers (Re=100, 1000, 2500) obtained with using 256×256 grid points.

                                

From these figures we can see that both compact and wide fourth order formulation solutions are very close to each other. The difference between them is most visible in the region of the center of the main circulation. In the literature, for the benchmark driven cavity flow, the value of the streamfunction and the vorticity at the center of the main circulation and also the location of this center is widely accepted as the flow parameters to compare the accuracy of the numerical solutions. In Table 2 we tabulated the value of the streamfunction and the vorticity at the center of the main circulation and the location of the center for the Reynolds numbers considered.

          

In Table 2 we also have tabulated solutions found in the literature that we believe to be very accurate. Erturk and Gokcol [7] have presented compact fourth order (Dh4) solutions obtained using a very fine grid mesh (600×600). Botella and Peyret [3] have used a Chebyshev collocation method and obtained highly accurate spectral solutions for the driven cavity flow. Barragy and Carey [2] have used a p-type finite element scheme and presented highly accurate (Dh8 order) solutions. Schreiber and Keller [18], have presented high-order (Dh8 order in theory) solutions obtained using repeated Richardson extrapolation using the solutions obtained on different grid mesh sizes. Erturk et al. [6] have also used repeated Richardson extrapolation and presented high-order (Dh6 order in theory) solutions. From Table 2 we can see that the compact and wide formulation solutions are very close to each other and also all of the solutions agree well with the highly accurate solutions found in the literature.

In order to test the numerical performances of the wide and compact formulations, we then solve the benchmark driven cavity flow problem several times with different time steps Dt by changing a. In order to see the effect of the number of grids on the numerical performances of both formulations, we have tabulated the numerical performances of wide and compact formulations to achieve same level of convergence. Table 3 shows the CPU time and the number of iterations necessary for convergence of the wide and compact fourth order formulations for a 128×128 grid mesh and also Table 4 shows the same for a 256×256 grid mesh.

                     

From Tables 3 and 4 we can see that, for the same Reynolds number and for the same time step (a) both wide fourth order formulation and compact fourth order formulation converge to the same level in about the same number of iteration such that the ratio of the number of iteration necessary for convergence for wide formulation to the number of iteration necessary for convergence for compact formulation is approximately 1. For wide and compact fourth order formulations, the convergence rate in the pseudo time is approximately the same.

From Tables 3 and 4 we also can see that, for the same Reynolds number and for the same time step (a), the CPU time necessary for convergence of compact fourth order formulation is higher than the CPU time necessary for convergence of wide fourth order formulation. For the solution of the compact formulation, in each time step the coefficients A, B, C, D, E and F in equations (5) and (6) have to be calculated as they are defined in equation (7) and this would require a CPU time. The compact formulation requires the solution of tri-diagonal systems, and these tri-diagonal systems are solved efficiently using the Thomas algorithm. On the other hand, the wide formulation requires the solution of penta-diagonal systems. For these penta-diagonal systems we used the modified Thomas algorithm. The modified Thomas algorithm for penta-diagonal systems runs slower than the Thomas algorithm for tri-diagonal systems since it requires more mathematical instructions. Even though the modified Thomas algorithm runs slower, from Tables 3 and 4 we see that, for convergence the compact formulation requires more CPU time than that of the wide formulation. This shows that the extra CPU time necessary to calculate the coefficients A, B, C, D, E and F is higher than the extra CPU time necessary for the modified Thomas algorithm. For the same Reynolds number and for the same time step (a), in terms of the CPU time necessary for convergence, the wide fourth order formulation is more advantageous than the compact forth order formulation such that the compact fourth order code requires almost 1.5 more CPU time than the wide fourth order code.

For a numerical formulation, the largest time step that the numerical code would not diverge is an indicative of the numerical stability of the numerical formulation. In order to find the largest time step for the formulations, we progressively increased the time step (a) with 0.01 and run the code several times for a given Reynolds number until the solution is no longer convergent such that we could not obtain a solution, i.e. the solution is divergent. In Tables 3 and 4 we can see that, in all cases, we can obtain numerical solution using larger time steps in compact formulation compared to that of wide formulation. This would indicate that the compact formulation is numerically more stable than the wide formulation and allows us to use larger time steps.

When running a code, one would like to use the maximum possible time step (Dt) for a faster convergence. In Table 3 and 4, comparing the CPU time of the compact formulation when maximum time step is considered (a = 1.3) with the CPU time of the wide formulation when maximum time step is considered (a = 1.1), we can see that when maximum possible time steps are used, the compact formulation runs approximately 1.3 times slower than the wide formulation in terms of the CPU time.

 

5  CONCLUSIONS

In this study we have numerically solved the Navier-Stokes equations using both fourth order compact formulation and fourth order wide formulation and compared the numerical performances of both fourth order formulations.

Solving a model problem which has an exact analytical solution to N-S equations with several different grid mesh showed that the solution of both wide and compact formulations, indeed converge with fourth degree with respect to the grid spacing. Also solving the benchmark driven cavity flow showed that the numerical solution of both formulations are very close to each other and produce spatially very accurate results.

We see that both wide and compact formulations converge to a specified convergence level in almost the same number of iterations for a given Reynolds number and time step. In terms of number of iterations necessary for convergence, both formulations have the similar convergence rate.

For a given Reynolds number and time step, wide formulation requires less CPU time than the compact formulation. In terms of computing time, wide formulation runs faster and is advantageous over the compact formulation.

In compact formulation we were able to obtain convergence using larger time steps than the time steps we use in wide formulation. This would show that the compact formulation is numerically more stable than the wide formulation, i.e. allows us to use larger time steps. In terms of numerical stability compact formulation is advantageous over wide formulation.

We would like to point out that in order to use the wide formulation, the points adjacent to boundaries have to be treated specially and doing this adds a great complexity in coding stage. On the other hand, even though compact formulation does not have this complexity and can be easily used for the points adjacent to boundaries, in order to use a compact formulation one has to deal with the extra coefficients in the equations and this can be counted as the complexity of the compact formulation in coding stage. Therefore we can fairly say that both wide and compact formulation requires almost the same level of coding effort.

 

References

[1]
Anderson DA, Tannehill JC, Pletcher RH, 1984. Computational Fluid Mechanics and Heat Transfer. McGraw-Hill, New York.

[2]
Barragy E, Carey GF, 1997. Stream function-vorticity driven cavity solutions using p finite elements. Computers and Fluids 26, 453- 468.

[3]
Botella O, Peyret R, 1998. Benchmark spectral results on the lid-driven cavity flow. Computers and Fluids 27, 421- 433.

[4]
Chung TJ, 2002. Computational Fluid Dynamics. Cambridge University Press, Cambridge.

[5]
Dennis SC, Hudson JD, 1989. Compact h4 Finite Difference Approximations to Operators of Navier-Stokes Type, Journal of Computational Physics 85, 390-416.

[6]
Erturk E, Corke TC, Gokcol C, 2005. Numerical Solutions of 2-D Steady Incompressible Driven Cavity Flow at High Reynolds Numbers, International Journal for Numerical Methods in Fluids 48, 747-774.

[7]
Erturk E, Gokcol C, 2006. Fourth Order Compact Formulation of Navier-Stokes Equations and Driven Cavity Flow at High Reynolds Numbers, International Journal for Numerical Methods in Fluids 50, 421-436.

[8]
Erturk E, 2008. Numerical Performance of Compact Fourth-Order Formulation of the Navier-Stokes Equations, Communications in Numerical Methods in Engineering, In Press, DOI: 10.1002/cnm.1090.

[9]
Ge L, Zhang J, 2001. High Accuracy Iterative Solution of Convection Diffusion Equation with Boundary Layers on Nonuniform Grids, Journal of Computational Physics 171, 560-578.

[10]
Gresho PM, 1991. Incompressible Fluid Dynamics: Some Fundamental Formulation Issues, Annual Review of Fluid Mechanics 23, 413-453

[11]
Gupta MM, Manohar RP, Stephenson JW, 1984. A Single Cell High Order Scheme for the Convection-Diffusion Equation with Variable Coefficients, International Journal for Numerical Methods in Fluids 4, 641-651.

[12]
Huang H, Yang H, 1990. The Computational Boundary Method for Solving the Navier-Stokes Equations, Institute of Applied Mathematics and Statistics Technical Report. No: Iam90-05, http://www.iam.ubc.ca/pub/techreports/.

[13]
Lele SK, 1992. Compact finite difference schemes with spectral-like resolution, Journal of Computational Physics 103, 16-42.

[14]
Li M, Tang T, Fornberg B, 1995. A Compact Forth-Order Finite Difference Scheme for the Steady Incompressible Navier-Stokes Equations, International Journal for Numerical Methods in Fluids 20, 1137-1151.

[15]
MacKinnon RJ, Johnson RW, 1991. Differential-Equation-Based Representation of Truncation Errors for Accurate Numerical Simulation, International Journal for Numerical Methods in Fluids 13, 739-757.

[16]
Peaceman DW, Rachford Jr. HH, 1955. The numerical solution of parabolic and elliptic differential equations, Journal of the Society for Industrial and Applied Mathematics 3, 28-41.

[17]
Richards CW, Crane CM, 1979. The accuracy of finite difference schemes for the numerical solution of the Navier-Stokes equations, Applied Mathematical Modelling 3, 205-211.

[18]
Schreiber R, Keller HB, 1983. Driven cavity flows by efficient numerical techniques. Journal of Computational Physics 49, 310 -333.

[19]
Spotz WF, Carey GF, 1995. High-Order Compact Scheme for the Steady Streamfunction Vorticity Equations, International Journal for Numerical Methods in Engineering 38, 3497-3512.

[20]
Spotz WF, Carey GF, 1998. Formulation and Experiments with High-Order Compact Schemes for Nonuniform Grids, International Journal of Numerical Methods for Heat and Fluid Flow 8, 288-303.

[21]
Spotz WF, 1998. Accuracy and Performance of Numerical Wall Boundary Conditions For Steady, 2D, Incompressible Streamfunction Vorticity, International Journal for Numerical Methods in Fluids 28, 737-757.

[22]
Stortkuhl T, Zenger C, Zimmer S, 1994. An Asymptotic Solution for the Singularity at the Angular Point of the Lid Driven Cavity, International Journal of Numerical Methods for Heat and Fluid Flow 4, 47-59.

[23]
Visbal MR, Gaitonde DV, 1999. High-Order-Accurate Methods for Complex Unsteady Subsonic Flows, AIAA Journal 37, 1231-1239.

[24]
Yang HH, 1993. The Computational Boundary Method for Solving the Incompressible Flows, Applied Mathematics Letters 6, 3-7.

[25]
Zhang J, 1998. An Explicit Fourth-Order Compact Finite Difference Scheme for Three-Dimensional Convection-Diffusion Equation, Communications in Numerical Methods in Engineering 14, 209-218.

[26]
Zingg DW, 1996. On finite-difference methods for linear acoustic wave propagation, AIAA Paper, Paper No:1996-280.

[27]
Zingg DW, 2000. Comparison of High-Accuracy Finite-Difference Methods for Linear Wave Propagation, SIAM Journal on Scientific Computing 22, 476-502.



File translated from TEX by TTH, version 3.63.

 

www.cavityflow.com