This is a scientific web page about the twodimensional steady incompressible flow in a driven cavity. The steady incompressible 2D NavierStokes 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 twodimensional steady incompressible flow in a driven cavity. The steady incompressible 2D NavierStokes 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.
where x and y are the Cartesian coordinates and Re is the Reynolds number. Erturk and Gokcol [11] have introduced the Fourth Order NavierStokes (FONS) equations. The introduced FONS equations are the following
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 NavierStokes equations and , the following finite difference equations provide second order (O (Dx^{2} )) accuracy
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 (Dx^{4} )) accurate to the NavierStokes.
where
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 (Dx^{2} ,Dy^{2} )) 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 (Dx^{4} ,Dx^{2} Dy^{2} ,Dy^{4} )) 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 semiimplicit 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
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 tridiagonal system in xdirection
then we solve the following tridiagonal system in ydirection
Similarly when we apply the ADI method to solve equation , we first solve the following tridiagonal system in xdirection
then we solve the following tridiagonal system in ydirection
where d_{xx} and d_{yy} denote the second order finite difference operators, and similarly d_{x} and d_{y} denote the first order finite difference operators in x and ydirection respectively, for example
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
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
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 Dt^{2} 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 Dt^{2} terms due to the factorization, Erturk et al. [10] have added the same amount of Dt^{2} 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
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 xdirection in the following tridiagonal system
When the solution for f is obtained, the streamfunction variable is advanced into the next time level by solving the following tridiagonal system in ydirection
Similarly, for the solution of equation , we first solve for the introduced variable g in xdirection in the following tridiagonal system
When the solution for g is obtained, the vorticity variable is advanced into the next time level by solving the following tridiagonal system in ydirection
Störtkuhl 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 Störtkuhl et al. [13] and use the following expression for calculating vorticity values at the wall
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 Störtkuhl et al. [13] and use the following expression for calculating the vorticity values
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
Similarly, for the corner points also shown in Figure 1, the vorticity is calculated as the following.
The reader is referred to Störtkuhl 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, tridiagonal 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/Dh^{2}, also in vorticity equations the diagonal elements on the LHS matrices become 1+2Dt/(ReDh^{2}). 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=aDh^{2} and for vorticity equation we use Dt=aReDh^{2}, 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
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 Residual_{y} and Residual_{w} 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 (Dx^{4} )). 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 NavierStokes equations can be easily applied to the FONS equations in order to obtain fourth order accurate solutions (O (Dx^{4} )). One of the features of the introduced FONS equations is that any existing code that solve the NavierStokes equations with second order accuracy (O (Dx^{2} )) can be easily altered to provide fourth order accuracy (O (Dx^{4} )) 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
File translated from T_{E}X by T_{T}H, version 3.63.


