Research Article  Open Access
K. Poochinapan, "Numerical Implementations for 2D LidDriven Cavity Flow in Stream Function Formulation", International Scholarly Research Notices, vol. 2012, Article ID 871538, 17 pages, 2012. https://doi.org/10.5402/2012/871538
Numerical Implementations for 2D LidDriven Cavity Flow in Stream Function Formulation
Abstract
The aim of this paper is to study the properties of approximations to nonlinear terms of the 2D incompressible NavierStokes equations in the stream function formulation (timedependent 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 liddriven cavity is solved up to Reynolds number Re with secondorder 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 NavierStokes equations (NSEs) representing incompressible viscous flows. Some of these are schemes utilizing primitive variables (velocitypressure) [1–4], vorticitystream function [5–9], 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 fullnonlinear 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 vorticitystream 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 fourthorder 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 CauchyKowalevska 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 liddriven 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 liddriven flow occupies the region (the cavity) Noslip 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 twopoint 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 CrankNicolson 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 centraldifference 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 solutionbanded 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 solutionbanded 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 liddriven 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 timedependent 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 clockwiserotating primary, whose location occurs toward the geometric center of the square cavity, and several small eddies: the counterclockwise 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 righthand 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 timestep 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.


(a)
(b)
4.2. TimeStep 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 timestep 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 timestep. 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 3–8 clearly illustrate the grid independence of the computational solution. The minimum and maximum of stream functions are listed in Tables 6–9 which includes a comparison to the highaccuracy 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 finitedifference 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 .

(a)
(b)
(c)
(a)
(b)
(c)
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