Abstract

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

1. Introduction

There are many finite difference methods for the solution of the Navier-Stokes equations (NSEs) representing incompressible viscous flows. Some of these are schemes utilizing primitive variables (velocity-pressure) [14], vorticity-stream function [59], and stream function formulation [10, 11]. The practical estimation of any schemes may be different from the theoretical estimation because of the nonlinearity of the NSEs and the implicit characteristic of the continuity condition. There is no single method which is most suitable for all aspects. The solution of the full-nonlinear set of discretized fluid flow equations is usually obtained by solving a sequence of linear equations. The type of linearization which is used can significantly affect the rate of convergence of the iterations to the solution.

The primary difficulty in obtaining numerical solutions with primitive variable formulation is that there is no evolution equation for pressure variable. The pressure function can be removed from the equations by means of vorticity-stream function, but an explicit boundary condition for vorticity function is missing. The simplest way to avoid the case is to eliminate the vorticity and turn to a single nonlinear biharmonic equation for the stream function.

The condition of stability of an explicit scheme imposes very restrictive limitations on the time increment with the artificial time. For the fourth-order derivative, it is . That is why the problem of constructing an implicit scheme is of great significance. The straightforward implementation of an implicit scheme results in the very large linear system whose solution needs very large computer capacity. One of the most efficient approaches to reduce the computational time without compromising the stability of the scheme is the method of operator splitting [12]. The application of the splitting method to the NSEs in stream function formulation is not straightforward because of the fact that they are not a Cauchy-Kowalevska system and there is no time derivative of the stream function. What is present in the equation is the term which cannot serve as a basis for splitting [1].

The new scheme combines the computational simplicity of the implicit scheme in linear terms with semiexplicit approximation of nonlinear terms. The general idea when treating the nonlinear term is to represent it as implicit approximation and then to linearize it and to conduct internal iterations. After the inner iterations converge, one obtains, in fact, the solution for the new time stage. The explicit approximation of nonlinear terms accomplish severe requirement on time step. A single internal iteration on nonlinear terms induces sense of implicit approximation and reduce very severe band on time step. To improve stability properties of explicit approximation of nonlinear terms we require only 3 internal iterations and call such algorithm as a method with internal iteration.

As to the question of which method to be used, the answer is that it depends on the type of problems that to be solved. The objective of the present study is to validate the efficiency and accuracy of numerical implementations for the NSEs in term of stream function. The content of this paper is organized as follows. Section 2 contains the mathematical formulation of the problem. Section 3 deals with the discretization of the equations and detailed description of numerical algorithms. The methods are applied to the test problem of lid-driven cavity flow up to . Results and discussions of numerical solutions are presented in Section 4, where we make a detailed comparison with available numerical data.

2. Mathematical Formulation

Consider a closed 2D domain with a piecewise smooth boundary . The NSEs for a viscous incompressible flow in the terms of stream function are where .

Boundary and initial conditions are the following: where is a vector normal to domain boundary, and the Reynolds number is defined as , where is the characteristic velocity, is the characteristic length, and is the kinematic viscosity.

3. Difference Schemes

The lid-driven flow occupies the region (the cavity) No-slip boundary conditions take the following form: For the sake of simplicity, we assume a uniform grid with and spacing in and -direction, respectively, for and .

The mesh is staggered in direction on and in direction on with respect domain boundaries. The boundary conditions are approximated on two-point stencils with the second order of approximation as follows: where with .

3.1. Explicit Scheme in Nonlinear Terms

The simplest method is explicit in which all nonlinear terms are evaluated using known values at . The biharmonic operator is approximated by the Crank-Nicolson scheme. Thus, one has linear system to calculate the new value of the unknown at each node. It is clear that each term in (2.1) can be approximated using central-difference operators for all derivatives. The explicit finite difference scheme for (2.1) is where and . The difference formulas to be used are given in Table 1.

To combine equation (3.5) as a single linear system with band matrix, we introduce the new system of index as follows: Each node of grid associates with index . It is easy to see that Now we introduce a new grid function , which is defined on the composite grid. Substituting in lieu of in (3.5), we recast the algebraic system as follows: According to this idea, can be rewritten as follows: The general sequence of the algorithm is as follows: (i)Set the values of and choose an initial guess .(ii)Consider as known entities and calculate from the system of the algebraic equation (3.8) by the direct method to the solution-banded linear system. We used standard subroutines DGBSV and DGBSVX of LAPACK.(iii)If the following criterion is satisfied then the calculations are terminated. Otherwise the index of iterations is stepped up, and the computation is returned to step (ii). A disadvantage of explicit method is that convergence can be very slow if a small nonlinear residual is required.

3.2. An Operator Splitting Scheme

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

3.3. The Method with Internal Iteration in Nonlinear Terms

The general idea when treating the nonlinear term is to represent it as implicit approximation and then to linearize it and to conduct internal iterations. After the inner iterations converge, one obtains, in fact, the solution for the new time stage. The explicit approximation of nonlinear terms accomplished severe requirement on time step. A single internal iteration on nonlinear terms induced sense of implicit approximation and reduced very severe band on time step. An improvement in accuracy is achieved at small additional expense. The internal iteration scheme may be viewed as a modification of the explicit scheme (3.5). For the beginning, the internal iteration scheme used only 3 iterations. The reason behind our choices is that a single iteration is equivalent to one step of explicit scheme so, in order to take the economies of the explicit scheme, it is desirable not to perform to many iteration (we choose 3 internal iteration). The method with internal iteration of nonlinear terms can be recasted in the following form: Substituting in lieu of in (3.11), we recast the algebraic system as the following: The general sequence of the algorithm is as follows: (i)Set the values of and choose an initial guess . (ii)Let(a) then calculate from the system of the algebraic equation (3.12) by the direct method to solution-banded linear system. We used standard subroutines DGBSV and DGBSVX of LAPACK,(b) then calculate from the system of the algebraic equation (3.12),(c) then calculate from the system of the algebraic equation (3.12).(iii)If the following criterion is satisfied then the calculations are terminated. Otherwise the index of iterations is stepped up, or and the computation is returned to step (ii). Note that the linear system for the problem can be written as the following multidiagonal system for the composite grid function : where . The matrix of the algebraic system equation (3.14) is banded with lower and upper diagonals.

4. Results and Discussion

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

4.1. Convergence

To evaluate grid dependence on spatial variables the solution is obtained on a sequence of grids with , , and nodes. The finest grid is used as a reference solution. The rate of convergence is computed using two grids according to the formula In Tables 2 and 3   (appeared at the center of primary eddy) and (appeared at the center of bottom right-hand side eddy) are computed as minimum and maximum values of stream function, respectively. Parameters and correspond to the rate of convergence received from (4.1). Table 2 represents steady state. The solutions are qualified as steady when the absolute error between two time steps is less than . As shown in Tables 2 and 3 on particular choice of the parameters, the estimated rate is close to the theoretically predicted second rate of convergence. Figure 2 shows the absolute error versus the number of iterations when the largest value of time-step is used and It shows that the rate of convergence to solution of the method of internal decreased faster than that of explicit and the operator splitting method.

4.2. Time-Step Independence

For fixed values of , and , the number of iterations is needed for convergence to steady state depends on . For definiteness, we select for presentation a genuinely moderate and . The results for the number of iterations to steady state are shown in Tables 4 and 5 where the minimal number of iteration for each resolution is shown in bold face. In the case of , we minimal numbers of iterations on grid are 5276 iterations, 3581 iterations, and 1588 iterations for case explicit, operator splitting, and internal iteration schemes, respectively. We see that the minimal number of iterations of the method of internal iteration is less than other methods.

Figures 3 and 4 show the largest value of time-step for which it is possible to reach a steady solution which corresponds to other investigators. Figure 3 shows that the present schemes converge to correct solutions in all cases of 3000 considered. If the grid spacing is too coarse, then the operator splitting scheme will not compute a correct solution for . Moreover, the internal iteration method requires slightly larger time-step. An important feature of the internal iteration method is that it allows coarsening of the grid depending on the flow solution and this feature is extremely useful for accuracy in predicting flow fields in region with high Reynolds number and it is compared with the operator splitting scheme.

The following Figure 5 shows that the modification of explicit scheme can reduce the number of iterations for the higher Reynolds number. Figure 6 shows CPU time needed to find solution at instant of time with the largest value of parameter . From comparison, it is observed that the internal iteration method seems to be performing better than the other two methods. Computation was carried out on a personal computer with 2.4 GHz CPU and the CPU time per one time step (/one iteration) on the grid was about 2.22031, 2.20781, and 6.61875 sec for case explicit, operator splitting, and internal iteration scheme, respectively. The numerical results are demonstrating that the internal iteration is one method to compute the characteristic of flow for the high Reynolds number. We found that the largest time increment of 3 internal iterations gives convergence to solution roughly 3 times as explicit scheme.

4.3. Grid Independence

Figures 7 and 8 show the steady state at the secondary vortex at the bottom right corner of the cavity at . Figures 38 clearly illustrate the grid independence of the computational solution. The minimum and maximum of stream functions are listed in Tables 69 which includes a comparison to the high-accuracy solution of previous literatures. Tables 6 and 7 show the values of , and the space location of the primary vortex. As can be seen from Table 6, three finite-difference schemes give very similar quantities on the grid. The values of , , and their locations are in accordance with results observed in the literature [1, 4, 14, 15]. Table 7 shows the values of , and the space location of the primary eddy for . On the grid for , the method of operator splitting reaches steady solution, but this solution does not agree with the numerical solutions obtained by other authors ([1, 4, 14, 15]). At the same time, explicit and internal iteration schemes reach “correct” solution. For the primary eddy, results of operator splitting scheme with the grid agree within with those obtained in [4] with the grid and [1], but the two other schemes differ significantly up to .

Tables 8 and 9 report , and the space location of the bottom right secondary eddy. As can be seen from Table 8, the quantities of agree with those obtained by the other authors. For the bottom right secondary eddy, our results on the grid agree within 5% with those obtained by the other authors and in the case of grid they differ significantly up to 25% from those obtained in [15]. Table 7 shows the values of , and the space location of the bottom right secondary eddy for . For the bottom right secondary eddy, the quantities with the grid agree within 10% with those obtained in [4] with the grid and [1] and differ significantly up to 15% from those obtained [4] with the mesh.

At the high Reynolds number the flow becomes more complicated, and significantly finer grids are needed in the vicinity of the walls where the dynamics of the flow is dominated by viscosity. Coarse grid in general will not resolve the viscous layer near the boundary. For example, as can be seen on the 52 × 52 grid, operator splitting method cannot get correct solution. Moreover, the other two have value of differing significantly up to 67% from those obtained [4] with the mesh.

Figures 9, 10, and 11 show well-known stream function contours for , and on the grid. Whole three algorithms give the similar stream function that cannot be distinguished on the stream line patterns. All of these show secondary vortices at the bottom corners of the cavity as well as at the top left of the cavity for case and . The bottom right tertiary vortices become quite visible for . Figures 12 and 13 present comparisons of the -velocity profiles along a vertical line and the -velocity profiles along a horizontal line passing through the geometric center of the cavity for Reynolds numbers ranging from 100 to 5000. These profiles are in good agreement with that of [5]. While the data of [5] was obtained using that on the grid, our data is obtained using that on the grid.

5. Conclusions

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

Acknowledgments

This research is supported by the Centre of Excellence in Mathematics and the Commission on Higher Education, Thailand.