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 : Communications in Numerical Methods in Engineering (2008)
Commun. Numer. Meth. Engng. 2008, Vol 24, pp 2003-2019

 

Numerical Performance of
Compact Fourth Order Formulation
of the Navier-Stokes Equations

Ercan Erturk

Energy Systems Engineering Department, Gebze Institute of Technology,
Gebze, Kocaeli 41400, Turkey
and
Energy Systems Engineering Department, Bahcesehir University,
Besiktas, Istanbul 41400, Turkey

 

Abstract

In this study the numerical performance of the fourth order compact formulation of the steady 2-D incompressible Navier-Stokes equations introduced by Erturk et al. (Int. J. Numer. Methods Fluids, 50, 421-436) will be presented. The benchmark driven cavity flow problem will be solved using the introduced compact fourth order formulation of the Navier-Stokes equations with two different line iterative semi-implicit methods for both second and fourth order spatial accuracy. The extra CPU work needed for increasing the spatial accuracy from second order (O (Dx2 )) to fourth order (O (Dx4 )) formulation will be presented.

 

1  INTRODUCTION

In Computational Fluid Dynamics (CFD) field of study High Order Compact formulations are becoming more popular. Compact formulations provide more accurate solutions in a compact stencil.

In finite difference, in order to achieve fourth order spatial accuracy, standard five point discretization can be used. When a five point discretization is used, the points near the boundaries have to be treated specially. Another way to achieve fourth order spatial accuracy is to use High Order Compact schemes. High Order Compact schemes provide fourth order spatial accuracy in a 3×3 stencil, hence this type of formulations can be used near the boundaries without a complexity.

In the literature, Zhang [14], Dennis and Hudson [3], MacKinnon and Johnson [6], Gupta et al. [5], Spotz and Carey [7] and Li et al. [9] have demonstrated the efficiency of high order compact schemes on the streamfunction and vorticity formulation of 2-D steady incompressible Navier-Stokes equations for uniform grids. Also in the literature the studies of Ge and Zhang [4] and Spotz and Carey [8] are example studies on the application of the high order compact scheme to nonuniform grids. The advantage of the high order compact schemes is that for a given flow problem and for a chosen grid mesh, high order compact formulation provides more accurate solutions (O (Dx4 )) compared to standard second order formulation (O (Dx2 )). Also for a given flow problem, the same level of accuracy of the solution obtained by standard second order formulation (O (Dx2 )) using a certain grid mesh, can be obtained with a smaller grid mesh when high order compact (O (Dx4 )) formulation is used.

Recently Erturk and Gokcol [11] 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" (FONS) equations are in the same form with the Navier-Stokes (NS) equations with additional coefficients. In fact the NS equations are a subset of the FONS equations and obtained when the additional coefficients are chosen as zero. Therefore any numerical method that solves the NS equations can be easily applied to the introduced FONS equations to obtain fourth order spatial accurate solutions. Moreover, the most important feature of the FONS equations is that any existing code with second order accuracy (O (Dx2 )) can easily be changed to provide fourth order accuracy (O (Dx4 )) by just adding some coefficients into the existing code at the expense of extra CPU work of evaluating these coefficients. This is an important feature such that if one already has a second order accurate code and wants to increase the accuracy of it to fourth order, instead of writing a new code, one can use the FONS equations and just by inserting some coefficient into the existing second order code, the existing second order code can turn into a fourth order code. Therefore the FONS equations introduced by Erturk and Gokcol [11] provide a very easy way to convert an existing second order accurate code into a fourth order accurate code and this is as simple as inserting some coefficients into the existing code. Of course, when this is done, the code will run a little slower because of the extra CPU work of evaluating the inserted coefficients. It will be good to estimate the CPU time needed for convergence of a converted fourth order accurate code compared to the CPU time needed for a second order accurate code.

Zhang [15] have studied the convergence and performance of iterative methods with fourth order compact discretization schemes. To the best of the author's knowledge in the literature there is not a study that documents the numerical performance of high order compact formulation the Navier-Stokes equations compared to regular second order formulation of the Navier-Stokes equations in terms of numerical stability and convergence for a chosen iterative method. In this study using the FONS equations introduced by Erturk and Gokcol [11], we will numerically solve the Navier-Stokes equations for both fourth order (O (Dx4 )) and second order (O (Dx2 )) spatial accuracy. This way we will be able to compare the convergence and stability characteristics of both formulations. In this study we will also document the extra CPU work that is needed for convergence when a second order accurate code is converted into a fourth order accurate code using the introduced FONS equations by Erturk and Gokcol [11]. The stability and convergence characteristics of both formulations and also the extra CPU work can show variation depending on the iterative numerical method used for the solution therefore in this study we will use two different line iterative semi-implicit numerical methods. Using these two numerical methods we will solve the benchmark driven cavity flow problem. First we will solve the cavity flow with second order (O (Dx2 )) spatial accuracy then we will solve the same flow with fourth order (O (Dx4 )) spatial accuracy. We will document the stability characteristics, such as the maximum allowable time increment (Dt), convergence characteristics, such as the number of iterations and the CPU time necessary for a chosen convergence criteria and also the extra CPU work that is needed to increase the spatial accuracy of the numerical solution from second order to fourth order using the FONS equations.

 

2  FOURTH ORDER COMPACT FORMULATION

In non-dimensional form, steady 2-D incompressible Navier-Stokes equations in streamfunction (y) and vorticity (w) formulation 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.

Erturk and Gokcol [11] have introduced the Fourth Order Navier-Stokes (FONS) equations. The introduced FONS equations are the following

2 y

x2
+ 2 y

y2
=-w+A
(3)

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
(4)

A
=
- Dx2

12
2w

x2
- Dy2

12
2w

y2
-
Dx2

12
+ Dy2

12

4y

x2y2
B
=
- Re Dx2

6
2y

x y
+ Re2 Dx2

12
y

y
y

y
C
=
Re Dy2

6
2y

xy
+ Re2 Dy2

12
y

x
y

x
D
=

Dx2

12
+ Dy2

12

3y

x2 y
- Re Dx2

12
y

y
2y

x y
+ Re Dy2

12
y

x
2y

y2
E
=

Dx2

12
+ Dy2

12

3y

x y2
- Re Dx2

12
y

y
2 y

x2
+ Re Dy2

12
y

x
2y

x y
F
=

Dx2

12
+ Dy2

12

y

y
3w

x y2
-
Dx2

12
+ Dy2

12

y

x
3w

x2 y
- Dx2

6
2 y

x2
2w

x y
+ Dy2

6
2y

y2
2w

x y
+ Re
Dx2

12
+ Dy2

12

y

x
y

y
2w

x y
-
Dx2

12
- Dy2

12

w

x
w

y
- 1

Re

Dx2

12
+ Dy2

12

4w

x2 y2
(5)

As it is described briefly in Erturk and Gokcol [11], the numerical solutions of FONS equations and are fourth order accurate to NS equations and , strictly provided that second order central discretizations shown in Table 1 are used and also strictly provided that a uniform grid mesh with Dx and Dy is used.

       

We note that NS equations are a subset of FONS equations and obtained when the coefficientsA, B, C, D, E and F in FONS equations are chosen as zero. The FONS equations are in the same form with the NS equations therefore any iterative numerical method applied to streamfunction and vorticity equation and can be easily applied to fourth order streamfunction and vorticity equation and . Moreover if there is an existing code that solves the streamfunction and vorticity equation and with second order spatial accuracy, by just adding some coefficients A, B, C, D, E and F into the existing code using the FONS equations, the same existing code can provide fourth order spatial accuracy. A single numerical code for the solution of the FONS equations can provide both second order and fourth order spatial accuracy by just setting some coefficients. Of course there is a pay off switching from second order to fourth order spatial accuracy, that is the extra cost of CPU work of calculating the coefficients A, B, C, D, E and F as defined in equation .

 

3  FINITE DIFFERENCE EQUATIONS

For numerical solutions of the Navier-Stokes equations and , the following finite difference equations provide second order (O (Dx2 )) accuracy

yxx +yyy = -w
(6)

1

Re
wxx + 1

Re
wyy = yy wx -yx wy
(7)

where subscripts denote derivatives as defined in Table 1.

As explained in Erturk and Gokcol [11] the solution of the following finite difference equations are fourth order (O (Dx4 )) accurate to the Navier-Stokes.

yxx +yyy = -w+A
(8)

1

Re
(1+B)wxx + 1

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

where

A
=
- Dx2

12
wxx - Dy2

12
wyy -
Dx2

12
+ Dy2

12

yxxyy
B
=
- Re Dx2

6
yxy + Re2 Dx2

12
yyyy
C
=
Re Dy2

6
yxy + Re2 Dy2

12
yxyx
D
=

Dx2

12
+ Dy2

12

yxxy - Re Dx2

12
yyyxy + Re Dy2

12
yxyyy
E
=

Dx2

12
+ Dy2

12

yxyy - Re Dx2

12
yyyxx + Re Dy2

12
yxyxy
F
=

Dx2

12
+ Dy2

12

yywxyy -
Dx2

12
+ Dy2

12

yxwxxy - Dx2

6
yxxwxy
+ Dy2

6
yyywxy + Re
Dx2

12
+ Dy2

12

yxyywxy -
Dx2

12
- Dy2

12

wxwy
- 1

Re

Dx2

12
+ Dy2

12

wxxyy
(10)

Note that equations and are in the same form with equations and except with additional coefficients A, B, C, D, E and F. In equations and if the coefficients A, B, C, D, E and F are chosen to be zero then the solution of these equations are second order accurate (O (Dx2 ,Dy2 )) to NS equations, since when the coefficients are zero, equations and are identical with equations and . However, in equations and if the coefficients A, B, C, D, E and F are calculated as they are defined in equation , then the solution of these equations and are fourth order accurate (O (Dx4 ,Dx2 Dy2 ,Dy4 )) to NS equations. When FONS equations are used, one can easily switch from second order to fourth order spatial accuracy using a single equation by just using the appropriate coefficients. Computationally, calculating the coefficients defined in equation when fourth order accuracy is desired will require an extra CPU work compared to second order accuracy. In order to quantify the extra CPU work to switch from equations and to equations and , i.e. switch from second order accuracy to fourth order accuracy, we solve the above equations with two different line iterative semi-implicit numerical method and document the CPU time for comparison.

We note that the equations and are nonlinear equations therefore they need to be solved in an iterative manner. In order to have an iterative numerical algorithm we assign pseudo time derivatives to equations and , thus we have

y

t
=yxx +yyy +w-A
(11)

w

t
= 1

Re
(1+B)wxx + 1

Re
(1+C)wyy -(yy +D)wx +(yx +E)wy -F
(12)

We solve these equations and in the pseudo time domain until the solution converges to steady state.

One of the numerical methods we will use to solve equations and is 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 [12], [1] and [2] for details. When we apply the ADI method to solve equation , first we solve the following tri-diagonal system in x-direction


1- Dt

2
dxx
yn+1/2 = yn + Dt

2
yyyn + 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
yxxn+1/2 + Dt

2
w- Dt

2
A
(14)
Similarly when we apply the ADI method to solve equation , 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
wyyn + Dt

2
(yx +E)wyn - 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)wxxn+1/2 - Dt

2
(yy +D)wxn+1/2 - Dt

2
F
(16)

where dxx and dyy denote the second order finite difference operators, and similarly dx and dy denote the first order 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 -2qi,j +qi-1,j

Dx2
+O (Dx2 )
(17)

where i and j are the grid index and q denote any differentiable quantity.

The second numerical method we will use is the efficient numerical method proposed by Erturk et al. [10]. Following Erturk et al. [10], first using an implicit Euler approximation for the pseudo time derivatives in equations and , we obtain the following finite difference formulations


1- Dt
dxx -Dtdyy
yn+1 = yn +Dtwn -DtAn
(18)


1-Dt 1

Re
(1+B)dxx -Dt 1

Re
(1+C)dyy
+ Dt
(yy +D)dx -Dt(yx +E)dy
wn+1 = wn -DtFn
(19)

We note that equations and are in fully implicit form and each equation requires the solution of a large banded matrix which is not computationally efficient. Instead, we spatially factorize these equations and thus we obtain the following finite difference equations


1- Dt
dxx

1- Dt
dyy
yn+1 = yn +Dtwn -DtAn
(20)


1-Dt 1

Re
(1+B)dxx +Dt(yy +D)dx
×

1-Dt 1

Re
(1+C)dyy -Dt(yx +E)dy
wn+1 = wn -DtFn
(21)

The advantage of these equations and is that each equation requires the solution of tridiagonal systems which is computationally very efficient using the Thomas algorithm. However, spatial factorization introduces Dt2 terms into the left hand side of equations and , also these terms remain in the solution even at the steady state. To cancel out these Dt2 terms due to the factorization, Erturk et al. [10] have added the same amount of Dt2 terms to the right hand side of the equations so that the equations recover the correct physical representation at the steady state. The final form of the finite difference equations take the following form


1- Dt
dxx

1- Dt
dyy
yn+1 = yn +Dtwn -DtAn +
Dt
dxx

Dt
dyy
yn
(22)


1-Dt 1

Re
(1+B)dxx +Dt(yy +D)dx
×

1-Dt 1

Re
(1+C)dyy -Dt(yx +E)dy
wn+1 = wn -DtFn

Dt 1

Re
(1+B)dxx -Dt(yy +D)dx
×

Dt 1

Re
(1+C)dyy +Dt(yx +E)dy
wn
(23)

The reader is referred to Erturk et al. [10] for details. The solution methodology for equations and involves a two stage time level updating. For example, for the solution of equation , we first solve for the introduced variable f in x-direction in the following tri-diagonal system


1- Dt
dxx
f=yn +Dtwn -DtAn +
Dt
dxx

Dt
dyy
yn
(24)

When the solution for f is obtained, the streamfunction variable is advanced into the next time level by solving the following tri-diagonal system in y-direction


1- Dt
dyy
yn+1 = f
(25)

Similarly, for the solution of equation , we first solve for the introduced variable g in x-direction in the following tri-diagonal system


1-Dt 1

Re
(1+B)dxx +Dt(yy +D)dx
g=wn -DtFn

Dt 1

Re
(1+B)dxx -Dt(yy +D)dx
×

Dt 1

Re
(1+C)dyy +Dt(yx +E)dy
wn
(26)

When the solution for g is obtained, the vorticity variable is advanced into the next time level by solving the following tri-diagonal system in y-direction


1-Dt 1

Re
(1+C)dyy -Dt(yx +E)dy
wn+1 = g
(27)

Strtkuhl et al. [13] have presented an analytical asymptotic solution near the corners of cavity and using finite element bilinear shape functions they also have presented a singularity removed boundary condition for vorticity at the corner points as well as at the wall points. For the boundary conditions, in both of the numerical methods described above we follow Strtkuhl et al. [13] and use the following expression for calculating vorticity values at the wall

1

3 Dh2

é
ê
ê
ê
ê
ê
ê
ê
ë

1

2,
-4
1

2
1
1
1
ù
ú
ú
ú
ú
ú
ú
ú
û
y + 1

9

é
ê
ê
ê
ê
ê
ê
ê
ë

1

2
2
1

2
1

4
1
1

4
ù
ú
ú
ú
ú
ú
ú
ú
û
w = - V

h
(28)

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.

For corner points, we again follow Strtkuhl et al. [13] and use the following expression for calculating the vorticity values

1

3 Dh2

é
ê
ê
ê
ê
ê
ê
ê
ë

-2
1

2
1

2
1
ù
ú
ú
ú
ú
ú
ú
ú
û
y + 1

9

é
ê
ê
ê
ê
ê
ê
ê
ë

1
1

2
1

2
1

4
ù
ú
ú
ú
ú
ú
ú
ú
û
w = - V

2h
(29)

where again V is equal to 1 for the upper two corners and it is equal to 0 for the bottom two corners.

In explicit notation, for the wall points shown in Figure 1, the vorticity is calculated as the following

wb = - 9V

2Dh
- 3

2Dh2
(yd +ye +yf )- 1

8
(2wa +2wc +wd +4we +wf )
(30)

Similarly, for the corner points also shown in Figure 1, the vorticity is calculated as the following.

wb = - 9V

2Dh
- 3

Dh2
yf - 1

4
(2wc +wf +2we )
(31)

       

The reader is referred to Strtkuhl et al. [13] for details on the boundary conditions

 

4  RESULTS

In order to quantify the extra CPU work needed when a second order accuracy code is converted into a fourth order accuracy code using the FONS equations introduced by Erturk and Gokcol [11], we have solved both the second order and fourth order accurate equations , , and for the solution of the driven cavity flow. We consider the driven cavity flow for Reynolds numbers of Re=100, 1000 and 3200, with using a grid mesh of 128×128 (Dx=Dy=Dh). We note that since we are mainly interested in finding the ratio of CPU time needed for convergence of a fourth order accuracy code to CPU time needed for convergence of a second order accuracy code, the choice of grid mesh size is not important, such that the ratio will be the same whether a coarse or fine grid mesh is used.

In solving the equations we decided to use the Alternating Direction Implicit (ADI) method and the numerical method proposed by Erturk et al. [10]. By using two different numerical methods we would be able to see if the extra CPU work is dependent on the numerical method used. While doing this, as a corollary, we would also be able to compare the ADI method and the method proposed by Erturk et al. [10] in terms of numerical performance. In both of the numerical method we use, for both second order and fourth order accuracy, the two equations, i.e. the streamfunction and the vorticity equations, are solved separately. In order to document the extra CPU work when a fourth order accuracy is desired what we do is, we first solve for second order accuracy and solve equations and . Then keeping the number of grids, the time step Dt and boundary conditions the same, we solve for fourth order accuracy thus we calculate and insert the coefficientsA, B, C, D, E and F into the equations and solve for equations and . While we solve the same flow problem, i.e. the driven cavity flow, for both second and fourth order accuracy we document the necessary number of iterations and the CPU time needed for a certain defined convergence criteria. This way we would be able to compare the convergence characteristics of both second and fourth order formulations in terms of the number of iterations and the CPU time, and also we would be able to document the extra CPU time needed if a second order code is converted into a fourth order code using the FONS equations.

For the choice of the time steps in solving the governing equations, we decided to use different time steps, Dt, for streamfunction and vorticity equations. In both of the numerical methods we use, while solving both the streamfunction and vorticity equations, tri-diagonal matrices appear on the implicit LHS of the equations. When second order accuracy is considered, in streamfunction equation the diagonal elements on the LHS matrices become 1+2Dt/Dh2, also in vorticity equations the diagonal elements on the LHS matrices become 1+2Dt/(ReDh2). We choose different time steps for streamfunction and vorticity equations that would make the diagonal elements the same in both equations. Therefore for streamfunction equation we use Dt=aDh2 and for vorticity equation we use Dt=aReDh2, where a is a coefficient we can choose. For fourth order accuracy we use the same time steps we use in second order accuracy. By using the same time steps in second and fourth order accuracy we would be able to compare the numerical stability of the FONS equations and the NS equations. In order to do this first we both solve the NS and the FONS equations using the same time steps. Then we increasea, i.e. increase the time stepDt, and solve the NS and the FONS equations 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 both the NS and the FONS equations for a given Reynolds number and grid mesh. This maximum allowable Dt for convergence is an indicative of the numerical stability. For example, using either of the numerical method, i.e. the ADI method or the Erturk method, we solve the same flow problem using both second and fourth order formulations. Therefore for a chosen numerical method, the maximum allowable time step for second and fourth order formulations will be indicative of the numerical stability characteristics of the second order formulation compared to that of the fourth order formulation. Also, using either of the formulations, i.e. second or fourth order formulations, we solve the same flow problem using both the ADI method and the Erturk method. Therefore, for a chosen formulation, the maximum allowable time step for the ADI and the Erturk methods will be indicative of the numerical stability characteristics of the ADI method compared to that of the Erturk method.

Our extensive numerical studies show that the increase in the extra CPU work is dependent on the computer and the compiler used. In this study, we run the codes on a 64 bit HP ES45 machine with EV68 AlphaChip 1.25 GHz processors with HP Tru64 UNIX operating system. We run the codes with both compiled normally and also compiled using maximum compiler optimization (-fast -O5).

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, the residual of the equations can be used as it was also used in Erturk et al. [10]. However we are solving two different equations, the NS and the FONS equations, and try to compare the CPU time of convergence for each equation to the same level. Therefore the residual of these equations may not show the same convergence level. Alternatively, we can use the difference of the streamfunction and vorticity variables between two time steps as the measure of convergence. However the solutions of the two different equations are slightly different since one is spatially second order and the other is fourth order accurate. Since the solutions are different, the difference of the streamfunction and vorticity variables between two time steps may not also show the same convergence level for those equations. Therefore as the 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


(32)

These residuals provide an indication of the maximum percent change in y and w variables in each iteration step. In all of the data presented in this study, for both the solutions of the NS and the FONS equations, obtained using both of the numerical methods, we let the iterations converge until both Residualy and Residualw are less than 10-8 . At this convergence level, this would indicate that the variables y and w are changing less than 0.000001% of their value between two iterations at every grid point in the mesh.

Figure 2, 3 and 4 show the streamline and vorticity contours of the driven cavity flow for Re=100, 1000 and 3200 respectively, obtained using the method proposed by Erturk et al. [10] applied to FONS equations (O (Dx4 )). We note that both second order and fourth order accurate solutions of ADI method and the Erturk method, agree well with the solutions found in the literature especially with Erturk et al. [10],[11].

                    

Using both of the numerical methods, we solve the driven cavity flow using different coefficients for time (a), i.e. using different time steps, and document the CPU time and iteration number needed to the desired convergence level explained above. Table 2 shows the CPU time and iteration numbers for ADI method for different Reynolds numbers using various a values. Table 3 shows the same for the method proposed by Erturk et al. [10]. Looking at Table 2 and 3, for both of the numerical methods the number of iterations for convergence is almost the same for second order and fourth order accuracy. However for both of the numerical methods the CPU time for fourth order accuracy is greater than the CPU time for second order accuracy as expected since the coefficients A, B, C, D, E and F have to be calculated at each iteration in fourth order accuracy which will result in an increase in the CPU time. The ratio of the CPU times of fourth order accuracy to second order accuracy show the increase in CPU time when we switch from second order accuracy to fourth order accuracy. From Tables 2 and 3 we see that this ratio seems to increase slightly when Reynolds number increases.

             

It seems that for both of the numerical methods, at a given Reynolds number, the ratios of CPU time and iteration number for second and fourth order accuracy is constant and it is independent of the time step, i.e. a.

For the ADI method, in Table 2, when the order of accuracy is increased to fourth order from second order, the CPU time increases almost 1.76 times for Re=100. This number increases as the Reynolds number increases and the increase in CPU time becomes 1.90 times for Re=3200.

For the method proposed by Erturk et al. [10], in Table 3, when the order of accuracy is increased to fourth order from second order, the CPU time increases almost 1.60 times for Re=100 and it is almost 1.68 times for Re=3200.

We note that, when the order of accuracy is increased from second order to fourth order, the 1.6 and 1.68 times increase in CPU time for Re=100 and 3200 respectively for the method proposed by Erturk et al. [10] are less than the equivalent 1.76 and 1.90 times increase in CPU time for the same Reynolds numbers for the ADI method. This shows that, the extra CPU time needed for fourth order accuracy when FONS equations are used, is dependent on the numerical method used and the extra CPU time for the method proposed by Erturk et al. [10] is much lower than the extra CPU time for the ADI method.

In Table 2 and 3, comparing the two methods, for the same Reynolds numbers and for the same a values (Dt), the iteration numbers for convergence are almost the same for both ADI and the method proposed by Erturk et al. [10], however the CPU time for ADI is less than that of the method proposed by Erturk et al. [10]. The reason for this is that in the method proposed by Erturk et al. [10] on the RHS of the finite difference equations more terms have to be calculated at each iteration and this increases the CPU time compared to the ADI method.

For faster convergence one can use larger time steps, if the numerical method used has a higher numerical stability limit. Therefore, for a numerical method, the maximum allowable time step (Dt) for convergence gives an indication of the numerical stability limit of the numerical method. Since in this study we have used two different numerical methods for the same flow problem, we decided to compare the numerical stability limit of the two methods applied to both the NS and the FONS equations, by finding the maximum allowable time step for convergence for both numerical methods. In order to find the maximum allowable time step for convergence, for a given Reynolds number we solve the second and fourth order equations using both of the numerical methods several times while increasing a with 0.01 increments each time, until the solution no longer converges.

For the ADI method the maximum allowable a for convergence is 0.79 for second order accuracy and it is 0.78 for fourth order accuracy for Re=1000. For the method proposed by Erturk et al. [10] the maximum a values was 1.89 for second order accuracy and it was 1.75 for fourth order accuracy. This would indicate that one can use much larger time steps in Erturk method compared to the ADI method, for example, for Re=1000 the method proposed by Erturk et al. [10] allows to use 2.4 times larger time step for the NS equations and 2.2 times larger time step for the FONS equations than the ADI method. From this we can conclude that the Erturk method has better numerical stability characteristics compared to the ADI method. When the maximum allowable time steps are used, the required CPU time for the method proposed by Erturk et al. [10] is almost 0.53 of the required CPU time for the ADI method. This means that the method proposed by Erturk et al. [10] converge almost twice faster than the ADI method when the maximum allowable time steps are used.

Comparing the numerical stability of the NS and the FONS equations, we see that for a chosen numerical method the FONS equations have slightly less stability limit than that of the NS equations. For example, for the ADI method and for Re=1000 the value of 0.79 for maximum allowable a for convergence for second order accuracy drop down to 0.78 when fourth order accuracy is used. Also, for the Erturk method and for Re=1000 the maximum allowable a value of 1.89 for second order accuracy drop down to 1.75 if we switch to fourth order accuracy. This would indicate that, for fourth order formulations the maximum allowable time step for convergence is lower than the maximum allowable time step for convergence for second order formulations.

We then decided to run the same codes compiled with using the maximum compiler optimization (-fast -O5). Table 4 and 5 document the CPU time and iteration numbers when compiler optimization is used for the same conditions documented in Table 2 and 3 respectively. Comparing the numbers in Table 4 and 2 for ADI method, compiler optimization decreases the necessary CPU time for convergence about 0.71 times for second order accuracy and about 0.51 times for fourth order accuracy. Also comparing the same in Table 5 and 3 for the method proposed by Erturk et al. [10], compiler optimization decreases the necessary CPU time for convergence about 0.54 and 0.45 for second order and fourth order accuracy respectively. These numbers show that compiler optimization decreases the CPU time significantly and for the numerical method proposed by Erturk et al. [10] the codes almost run twice faster in terms of CPU time when compiler optimization is used.

             

 

5  CONCLUSIONS

The FONS equations introduced by Erturk and Gokcol [11] are in the same form of the NS equations, therefore any numerical method that solve the Navier-Stokes equations can be easily applied to the FONS equations in order to obtain fourth order accurate solutions (O (Dx4 )). One of the features of the introduced FONS equations is that any existing code that solve the Navier-Stokes equations with second order accuracy (O (Dx2 )) can be easily altered to provide fourth order accuracy (O (Dx4 )) just by adding some coefficients into the existing code using the FONS equations. This way, the accuracy of any second order code can be easily increased to fourth order, however there is a pay off for this increased accuracy, that is the extra CPU time for calculating the added coefficients.

In this study we have solved the NS equations and the FONS equations to document the extra CPU time necessary for convergence when an existing second order accurate code is altered to provide fourth order accurate solutions. For this we have used two different numerical methods. We find that the extra CPU time is slightly dependent on the numerical method used. For the ADI method to obtain fourth order accurate solutions of driven cavity flow, the CPU time increases 1.8 times compared to second order accurate solutions for Re=1000. Also for the numerical method proposed by Erturk et al. [10], with the cost of 1.64 times the CPU time necessary for second order accuracy, a fourth order accurate solutions can be obtained for Re=1000, using the FONS equations.

The FONS equations introduced by Erturk and Gokcol [11] provide a very easy way of obtaining fourth order accurate solutions by just adding some coefficients into an existing second order accurate code, at the expense of a minor increase in the CPU time.

 

References

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

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

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

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

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

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

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

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

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

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

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

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

[13]
Strtkuhl T, Zenger C, Zimmer S. 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 1994; 4:47-59.

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

[15]
Zhang J. On Convergence and Performance of Iterative Methods with Fourth-Order Compact Schemes, Numerical Methods for Partial Differential Equations 1998; 14:263-280.



File translated from TEX by TTH, version 3.63.

 

www.cavityflow.com