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. If we discretize equations (1) and (2) using three point central differencing, the following finite difference equations provide second order (O (Dx^{2} ,Dy^{2})) accurate solutions
where subscripts denote second order (O (Dx^{2} ,Dy^{2})) three point central finite difference derivative approximations.
As explained in Erturk and Gokcol [7], with using second order (O (Dx^{2} ,Dy^{2})) three point central differencing, the solution of the following finite difference equations are fourth order (O (Dx^{4} ,Dx^{2} Dy^{2} ,Dy^{4})) accurate to the NavierStokes equations (1) and (2).
where
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 (Dx^{2} ,Dy^{2})) to the NavierStokes 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 (Dx^{4} ,Dx^{2} Dy^{2} ,Dy^{4})) to the NavierStokes 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 NavierStokes 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
where
As it is described briefly in Erturk and Gokcol [7], the numerical solutions of equations (8) and (9) are fourth order accurate (O (Dx^{4} ,Dx^{2} Dy^{2} ,Dy^{4})) to the NavierStokes equations (1) and (2), strictly provided that second order central discretizations (O (Dx^{2} ,Dy^{2})) 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 NavierStokes 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 NavierStokes equations (1) and (2). Since equations (8) and (9) and the NavierStokes equations (1) and (2) are in the same form, any numerical method suitable for the NavierStokes 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 METHODOLOGYWe will numerically solve both the NS equations (1) and (2) and the introduced compact fourth order formulation of the NavierStokes 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
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 MethodologyIn 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 tridiagonal system in xdirection
then we solve the following tridiagonal system in ydirection
After this we solve the vorticity equation and in order to do this we first solve the following tridiagonal system in xdirection
then we solve the following tridiagonal system in ydirection
In equations (13), (14), (15) and (16), the superscripts denote the iteration time level, d_{xx} and d_{yy} denote the second derivative three point central finite difference operators, and similarly d_{x} and d_{y} denote the first derivative three point central finite difference operators in x and ydirection respectively, for example
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 tridiagonal 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 NavierStokes equations.
3.1 Wide Fourth Order Solution MethodologyIn 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 pentadiagonal system in xdirection
then we solve the following pentadiagonal system in ydirection
Then, similarly, we solve the vorticity equation and in order to do this we first solve the following pentadiagonal system in xdirection
then we solve the following pentadiagonal system in ydirection
We note that, in equations (18), (19), (20) and (21) the second and first derivative operators d_{xx}, d_{yy}, d_{x} and d_{y} denote five point central finite difference operators in x and ydirection respectively, for example
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 pentadiagonal systems. Following the above numerical procedure defined in equations (18), (19), (20) and (21) we obtain spatially fourth order accurate solutions of the NavierStokes equations.
3.1 Boundary ConditionsThe 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 1a) as the following
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 1b.
For corner points, we again follow Stortkuhl et al. [22] and calculate the vorticity values at the corner grid points (diamond point in Figure 1a) as the following
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 1c. 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
where V_{w} 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 1a) as the following
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 1a) is calculated as the following
We note that, this approximation require the streamfunction values at grid points outside the computational domain (hexagon points in Figure 1a). In order to find the streamfunction values at these grid points we use the following fourth order relation
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 1a), we calculate the streamfunction and vorticity values, using GaussSeidel 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 1a.
4 RESULTS AND DISCUSSIONSFollowing the numerical procedure described in the previous section we obtain spatially fourth order solutions of the NavierStokes 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, tridiagonal matrices appear on the implicit LHS of the equations. In streamfunction equation the diagonal elements on the LHS matrices become 1+2Dt/Dh^{2}, and in vorticity equations the diagonal elements on the LHS matrices become 1+2Dt/(ReDh^{2}). 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=aDh^{2} for the streamfunction equation and similarly we use Dt=aReDh^{2} 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
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 Residual_{y} and Residual_{w} 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 NavierStokes equations (1) and (2).
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 (y_{ex},w_{ex}) given in equation (30) and the numerical solution (y_{i,j},w_{i,j}), defined as
where N being the total number of grid points, should be proportional to h^{m}where the power should be equal to m=4. Since these average absolute differences should have the following form
logarithm of both sides give
In Table 1 we have tabulated the average absolute differences (E_{y} and E_{w}) 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 loglog 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 (E_{y} and E_{w}) 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
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 (Dh^{4}) 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 ptype finite element scheme and presented highly accurate (Dh^{8} order) solutions. Schreiber and Keller [18], have presented highorder (Dh^{8} 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 highorder (Dh^{6} 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 tridiagonal systems, and these tridiagonal systems are solved efficiently using the Thomas algorithm. On the other hand, the wide formulation requires the solution of pentadiagonal systems. For these pentadiagonal systems we used the modified Thomas algorithm. The modified Thomas algorithm for pentadiagonal systems runs slower than the Thomas algorithm for tridiagonal 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 CONCLUSIONSIn this study we have numerically solved the NavierStokes 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 NS 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
File translated from T_{E}X by T_{T}H, version 3.63.


