## Services on Demand

## Journal

## Article

## Indicators

- Cited by SciELO
- Access statistics

## Related links

- Cited by Google
- Similars in SciELO
- Similars in Google

## Share

## Brazilian Journal of Chemical Engineering

##
*Print version* ISSN 0104-6632

### Braz. J. Chem. Eng. vol.30 no.4 São Paulo Oct./Dec. 2013

#### http://dx.doi.org/10.1590/S0104-66322013000400023

**PROCESS SYSTEMS ENGINEERING**

**Finite volume simulation of 2-D steady square lid driven cavity flow at high reynolds numbers**

**K. Yapici ^{I, *}; Y. Uludag^{II}**

^{I}Cumhuriyet University, Department of Chemical Engineering, Phone: + 90.346.219 1010/1338, Fax: + 90.346.219 1179, 58140 Sivas, Turkey. E-mail: kyapici@cumhuriyet.edu.tr

^{II}Middle East Technical University, Department of Chemical Engineering, 06531 Ankara, Turkey

**ABSTRACT**

In this work, computer simulation results of steady incompressible flow in a 2-D square lid-driven cavity up to Reynolds number (*Re*) 65000 are presented and compared with those of earlier studies. The governing flow equations are solved by using the finite volume approach. Quadratic upstream interpolation for convective kinematics (QUICK) is used for the approximation of the convective terms in the flow equations. In the implementation of QUICK, the deferred correction technique is adopted. A non-uniform staggered grid arrangement of 768x768 is employed to discretize the flow geometry. Algebraic forms of the coupled flow equations are then solved through the iterative SIMPLE (Semi-Implicit Method for Pressure-Linked Equation) algorithm. The outlined computational methodology allows one to meet the main objective of this work, which is to address the computational convergence and wiggled flow problems encountered at high Reynolds and Peclet (*Pe*) numbers. Furthermore, after *Re* __>__ 25000 additional vortexes appear at the bottom left and right corners that have not been observed in earlier studies.

**Keywords:** Finite volume method; QUICK; Driven cavity flow; High Reynolds number.

**INTRODUCTION**

Lid-driven cavity flow of Newtonian fluids is one of the most well-known problems in the Computational Fluid Dynamics (CFD) literature due to its peculiar challenges in the form of singularities in spite of its simple geometry (Figure 1) (Botella and Peyret, 1998). In addition, the availability of both analytical solutions and experimental results for the lid-driven cavity flow field has enabled researchers to test and improve their computational methods through this benchmark geometry.

Recently, physical, mathematical and numerical aspects of the steady incompressible lid-driven cavity flow were discussed in detail by Erturk (2009). The findings of the study were also compared with those of earlier publications.

The main outcome of the study can be listed as follows. Due to the Taylor-Görtler-Like (TGL) vortices appearing even at moderate Reynolds (*Re*) around 1000, it is misleading to assume that flow in a lid-driven cavity is steady and 2-D. At even higher *Re* 2-D flow becomes completely fictitious (Shankar, 2000). On the other hand, 2-D incompressible steady driven cavity flow can be computable at high *Re* (Erturk, 2009) with its peculiar computational challenges. Therefore, steady lid-driven cavity flow, especially at high *Re* (*Re* __>__ 10000), has been exploited by many researchers in order to test and improve robustness and stability of their computational methods (Erturk, 2009; Erturk *et al.*, 2005; Ramšak and Škerget, 2004; Erturk and Gökçel, 2006; Sahin and Owens, 2003; Barragy and Carey, 1997; Schreiber and Keller, 1983; Ghia *et al.*, 1982).

In their report, Erturk *et al.* (2005) used stream function-vorticity formulation for the solution of 2-D steady incompressible flow in a lid-driven cavity. With a uniform grid size of 601x601 they obtained a second-order accurate steady solution up to *Re *of 21000. Rausak and Skerget (2004) introduced a new formulation of the integral boundary element method (BEM) using subdomain technique. They used a non- uniform grid of 100x60 sub-domains. With this technique, they were able to obtain a steady solution of driven cavity flow up to *Re* = 50000 by transient computation. However, their results suffer from oscillations at this high *Re* because of the coarse mesh structure. To the best knowledge of the authors, this value of Re is the highest one in the literature when steady solution of lid-driven cavity flow is concerned. The main objective of the present study is to demonstrate that the numerical solution of 2-D steady incompressible flow in a lid-driven cavity can be obtained at even higher *Re* (*Re* __<__ 65000) by using high-order linear schemes such as the quadratic upstream interpolation for convective kinematics (QUICK) proposed by Leonard (1979).

It has been well established that accuracy of a numerical solution is improved by using a smaller mesh in the regions of high gradients than the mesh size of bulk flow (Hartmann *et al*., 1990; Wang *et al*., 2005; Lilek and Perić, 1995). More specifically, in the case of flow in the lid-driven cavity, adopting a finer grid structure near the lid than in the bulk region enables one to resolve high gradients (Erturk and Gökçel, 2006) and to obtain an oscillation-free solution at high *Re* numbers (Erturk *et al.*, 2005). Use of a non-uniform grid structure entails modification of the numerical schemes. In the current study, a non-uniform version of the QUICK scheme (Yuguo and Baldacchino, 1995; Arampatzis *et al.*, 1994; Freitas *et al.*, 1985) is used for the approximation of convective terms. Generally the implementation of the high-order schemes is carried out by using the deferred correction method that was proposed by Khosla and Rubin (1974). A well-documented study that applied this approach for the formulation of the uniform QUICK scheme is given by Hayese *et al.* (1992). They tested and compared their formulation of QUICK that satisfies the five sets of rules proposed by Patankar (1980), with the other formulations on 2-D lid-driven square cavity flow. They concluded that their formulation is both numerically robust and faster than the others. Moreover, they pointed out that the low-order boundary condition approximation reduces the overall accuracy of the solution. For that reason, they recommended that at least second-order accurate schemes should be employed at the boundaries. Recently, Nacer *et al.* (2007) proposed a new formulation for the uniform QUICK scheme. The layout treatment of the lower-order term in the deferred correction equation appears to be the main difference between the formulations of Hayese *et al.* (1992) and Nacer *et al.* (2007). In their formulation, Nacer *et al. *(2007) employed the central difference scheme (CDS) instead of the upwind difference scheme (UDS) provided that the magnitude of Peclet (Pe) number |*Pe*| __<__ 2. They reported that this modification of QUICK leads to decreased computational time due to the smaller source term in CDS compared to the UDS.

In the numerical study on 3-D lid-driven cavity flow, Lilek *et al.* (1997) reported that the use of CDS is not necessarily restricted to a Peclet number of 2. They showed that, with their computational methodology, lid-driven cavity flows having Peclet number up to 100 can be handled without any significant convergence problem or wiggles. Similarly, Shyy *et al.* (1992) also used CDS and second-order upwind schemes for convective transport terms in the flow equations for 2-D lid-driven cavity flow up to *Re* __<__ 3200. They emphasized that both CDS and second-order upwind schemes result in oscillation-free solutions in the studied range of *Re* number.

In the light of the reported findings, if the mesh structure is sufficiently refined near the boundaries, it seems possible to use CDS instead of UDS to obtain deferred correction coefficients with two important consequences. The first one is that oscillations and convergence problems appearing at high *Re* are eliminated. No Pe number restriction can be considered as the second benefit of CDS as opposed to the stringent criterion of |*Pe*| __<__ 2 for UDS. The authors then believe that exploring the promising CDS to materialize these potential benefits would be a major contribution to the CFD literature. Therefore, in the current computational study on 2-D incompressible steady lid-driven cavity flow the promising approach, i.e., the use of CDS along with near boundary mesh refinement, is adopted for the formulation of a non-uniform version of QUICK scheme at high Reynolds numbers (*Re* __<__ 65000). It should also be emphasized that this methodology also enables one to treat the boundaries through second order accuracy.

The rest of the paper is organized as follows. In the following sections governing equations and new formulation of the QUICK scheme on a non-uniform mesh are introduced. Subsequently, numerical results for lid-driven cavity flow with* Re* __<__ 65000 are presented. The numerical results of the proposed formulation are also compared with those of earlier studies in the literature. Finally, conclusions are provided.

**GOVERNING EQUATIONS AND NUMERICAL METHODOLOGY**

We consider steady, incompressible and isothermal flow via Navier-Stokes equations in a two-dimensional Cartesian coordinate system (x,y). Stream-wise (x) and cross-stream (y) velocity components are u and v, respectively. Then the equation of continuity and the components of momentum equation in their dimensionless form can be written as follows:

These non-linear PDEs are discretized by using the finite volume method (Patankar, 1980; Versteeg and Malalasekera, 1995) in non-uniform staggered grid arrangement. Continuity and momentum equations can be written in the general form as follows:

where Λ represents the density ρ in the momentum conservation equation and the relaxation time λ in the constitutive equation; Γ is the diffusion coefficient and the source term S_{φ} stands for different flow quantities, depending on the equation, as shown in Table 1. The term φ represents any scalar quantity field, u or v being the x and y-components of the velocity field in this case. Equation (4) contains both diffusion and convective terms. In this study, diffusion terms are approximated by the central difference scheme. On the other hand, the non-uniform version of the QUICK (Yuguo and Baldacchino, 1995; Arampatzis *et al.*, 1994; Freitas *et al.*, 1985) scheme that provides third order spatial accuracy for a uniform grid structure is used for the convective terms. The QUICK scheme is based on the use of the three point upstream weighted quadratic interpolation technique (Hayese *et al.* 1992) for the prediction of the dependent variables at the control volume faces.

To ensure the stability of the higher-order schemes, a well-known and widely used technique, the deferred correction given by Eq. (5) (Khosla and Rubin, 1974), is used to implement the QUICK scheme along with CDS. The brief survey provided in the introduction on CDS demonstrated that the numerical solution with CDS shows no oscillatory solution even with a Peclet number much larger than the critical value of 2. Therefore, CDS can be treated as a lower-order scheme in Eq. (5). Here this approach is adopted for the formulation of the non-uniform version of the QUICK scheme, which is the same formulation proposed by Leonard (1979) for the uniform QUICK.

The superscripts LO, HO and 0 represent lower-order, higher order and values of the previous iteration, respectively.

For the sake of clarity, in the following lines the implementation of QUICK scheme with CDS on a one-dimensional system approximated by a non-uniform mesh is shown in Figure 2.

Here east and west cell-face values of the dependent variable can be plugged into Eq. (5) as:

The other face values of the dependent variables in the cross-stream direction can be obtained in a similar way. For one-dimensional flow, Eq. (4) simplifies to

Integration of Eq. (8) over a control volume, which is depicted in Figure 2, yields

After implementation of CDS for diffusive terms Eq. (9) becomes

Substituting Eqs. (6) and (7) into Eq. (10) and then rearranging leads to the final form of the discretized flow equation as follows:

where the coefficients are expressed through the following relations, provided that the central differences scheme is employed

where b_{S} in the source term is one of the extra deferred correction terms of the QUICK scheme treated explicitly from the values obtained in the previous iteration. In addition, the under-relaxation factor α_{φ} is introduced implicitly in the equations.

In this study, SIMPLE (Patankar, 1980) is employed to solve the coupled system of the continuity and momentum equations. A non-uniform staggered grid arrangement is employed to discretize the flow geometry. The set of linear algebraic equations is then solved by using the Line-by-Line technique based on the TDMA (Thomas algorithm or the tridiagonal matrix algorithm) and the alternative direction implicit (ADI) scheme. To stabilize the calculations, global under-relaxation factors are employed, depending on the values of the Reynolds number. The solution process is reiterated until the maximum relative changes of flow variables (u, v, p) are less then a prescribed tolerance or residual as:

where φ = (u,v,p)^{T.}

**RESULTS AND DISCUSSIONS**

The aforementioned numerical algorithm was applied to incompressible steady 2-D lid-driven cavity flow. In the following parts, mesh refinement analysis is first carried out by employing the series of refined meshes given in Table 2 for the purpose of comparison of the CDS and QUICK schemes at *Re* = 10000. Then the results obtained on highly graded dense mesh with the QUICK scheme are compared with the available results in the literature for various *Re* numbers. After that, a detailed investigation of lid-driven cavity flow is carried out at *Re* numbers up to __<__ 65000.

To ensure stability of the solution, under-relaxation factors of 0.7 and 0.1 are employed for the velocity at *Re* numbers of 1000 and 65000, respectively. Figure 3 demonstrates differences between the minimum values of the stream function, Ψ, computed by the QUICK and CDS schemes for *Re* = 10000 as a function of the minimum grid sizes on a log-log scale. The reference value is obtained by using the QUICK scheme and the finest mesh of M6. Although the CDS scheme provides convergent as well as non-oscillatory solutions for the meshes used, QUICK with CDS enables the same degree of relative error at much coarser mesh structures than the CDS scheme. The slopes of the lines that provides the relative accuracies are 2.51 and 1.86 for QUICK and CDS schemes, respectively.

In Tables 3 and 4 minimum values of the stream function and vorticity obtained at the center of the primary vortex for various *Re *numbers, including 1000, 10000, 15000 and 20000, are compared with those obtained in the earlier studies. Good agreement between the results of different studies is observed, especially at low *Re* numbers (*Re* = 1000), as pointed out by Sahin and Owens (2003) and Erturk *et al.* (2005). At higher *Re* (*Re *= 10000) differences between the results become more pronounced. The order of the scheme used for discretization of the convective transport terms and the mesh sizes used are believed to play an important role in these differences (Erturk *et al.,* 2005). The results included in Table 3 show that, in the case of lower order schemes, to achieve a similar degree of accuracy to that of the higher order one, the number of grids should be increased substantially.

Discretization of the flow using a graded mesh instead of a uniform mesh size throughout the flow domain has been observed to affect the accuracy of the computational results. For example, Hartmann and Peric (1990) used a finite volume multigrid method for the solution of natural convection flow in a square cavity. They presented accurate solutions obtained from both a uniform and non-uniform mesh of 640x640 up to a Rayleigh number of 10^{6}. They reported improved accuracy provided by the use of the non-uniform mesh structure compared to the uniform one having the same number of nodes. Sahin and Owens (2003) introduced a novel implicit cell vertex finite volume method for the solution of both the steady and unsteady Navier-Stokes equations. They applied their method to the lid-driven cavity flow up to *Re* = 10000. Their results indicate that the positive impact of adopting a non-uniform mesh structure or higher resolution near the lid boundary becomes more pronounced at higher *Re*. The advantage of a non-uniform grid structure is also observed in Table 4. Both results compare well with each other. Therefore, a non-uniform grid structure seems to compensate the disadvantages associated with the lower order schemes, especially at high *Re*. Hence, using a non-uniform mesh on the lid and cavity walls not only leads to improved accuracy of the solution, but also to non-oscillatory solutions at high *Re*.

In the following parts of this paper, the effect of increasing *Re *on the flow field is investigated in detail. The finest mesh M6 was used unless otherwise stated. Furthermore, the incremental continuation technique (Schreiber and Keller, 1983) was employed to obtain a convergent solution at *Re* higher than 10000. The incremental value of *Re* was 5000.

Figures 4 and 5 show profiles of the horizontal (u) and vertical (v) velocity components as a function of *Re* along the centerlines x=0.5 and y=0.5, respectively. In addition, the maximum, positive direction, and minimum, negative direction, values of the velocity components at the corresponding centerlines and intensities of the primary vortex at its center with their corresponding locations are tabulated in Table 5. Here subscripts denote values pertaining to the minimum and maximum velocities. One could expect that, as the inertial effects increase, the movement of the lid should be convected deeper and deeper in the flow, resulting in velocity components of higher magnitude and boundary layers of smaller thickness, in addition to the continuous movement of primary vortex center position towards the cavity lid. Up to a *Re* value of 20000 this trend is obtained in the flow, as shown in Figures 4 and 5 and Table 5. At higher *Re* no appreciable change in the magnitudes of the velocity components is observed. Furthermore, the upward shift in the primary vortex center position halts after *Re* of 25000, while the stream function at the center of the primary vortex starts to decrease after *Re* of 10000, which was also observed in an earlier study (Erturk *et al.*, 2005).

Figures 6 and 7 depict computed streamlines and vorticity contours, respectively. In the figures, the stream function contour levels are set with values similar to those of Sahin and Owens (2003) and Barragy and Carey (1997). Additional vortices increasing in size and intensity with respect to *Re* are apparent in the figures. The vorticity contours shown in Figure 7 indicate that it is possible to obtain a smooth solution at a relatively high *Re *of 65000 by using the QUICK with CDS scheme for the approximation of the convective terms.

For the assessment of the accuracy of the present results, the vorticity values at the center of the primary vortex in Table 5 are plotted with respect to *Re* in Figure 8. Also shown in the figure is the theoretically obtained value of the corresponding vorticity, 1.886, at the infinitely large *Re* limit (Burggraf, 1966). The results asymptotically approach the limiting value, demonstrating good agreement between the present computations and the theory.

In Figure 9, enlarged views of the stream function close the corners of the cavity are shown. The abbreviations BL, BR and TL refer to bottom left, bottom right and top left corners of the cavity, respectively. The stream function contours are plotted for three different *Re* numbers, including 20000, 50000 and 65000. The secondary, tertiary and quaternary vortices can be identified in the bottom left and right corners for each *Re* number.

At *Re* = 65000 an additional new quinary vortex (BR4) appears at the bottom right corner. However, numerical results tabulated in Table 7 indicate that this new vortex forms after *Re* = 25000 in the bottom left (BL4) and right corners (BR4). This quinary vortex has not been observed in previous numerical studies. Detailed views of the quinary vortex at the bottom left (BL4) and right corners (BR4) are presented in Figures 10 and 11, respectively. The size of this vortex grows with increasing *Re* number.

Results associated with the secondary vortices are listed in Table 6, which also includes relevant results obtained by Erturk (2009) for two different *Re* numbers of 10000 and 20000. It should be pointed out that, in the latter, a uniform mesh size of 1025x1025 was used, as opposed to the non-uniform mesh of 768x768 used in this study. The results of both studies agree with each other well. Furthermore, secondary vortices up to a *Re* number of 65000 with the increment of 5000 are tabulated in Table 7 to enable further comparisons in the future.

**CONCLUSIONS**

In this work, a numerical solution of 2-D steady Navier-Stokes equations for a lid-driven cavity at Reynolds numbers up to 65000 is presented. A non-uniform version of the QUICK scheme is employed for the approximation of the convective transport terms. To implement the QUICK scheme, the deferred correction technique is used. In this approach a lower-order scheme is treated by the central difference scheme. Detailed numerical experiments demonstrate that this formulation not only provides a smooth solution, but also enables convergence for the range of Reynolds number investigated. Furthermore, by using CDS, the compromise in the accuracy of the solution due to UDS can be limited. Discretization of the flow field by means of a non-uniform mesh structure instead of a uniform one leads to both improved accuracy and non-oscillatory solutions at high *Re*.

The present computational results compare well with those published for the investigated *Re* number range. Moreover, close agreement between the present results and theory is obtained as far as the vorticity value at the primary vortex center is concerned. New vortices that have not been reported before appear at the bottom left and right corner of the cavity after *Re* number 25000.

**ACKNOWLEDGEMENTS**

The authors acknowledge the financial support provided by The Scientific and Technological Research Council of Turkey through the research project of 109M012.

**REFERENCES**

Arampatzis, G., Assimacopoulos, D., Mitsoulis, E., Treatment of numerical diffusion in strong convective flows. International Journal for Numerical Methods in Fluids, 18, pp. 313-331 (1994). [ Links ]

Barragy, E., Carey, G. F., Stream function-vorticity driven cavity solutions using p finite elements. Computers and Fluids, 26, pp. 453-468 (1997). [ Links ]

Botella, O., Peyret, R., Benchmark spectral results on the lid-driven cavity flow. Computers & Fluids, 27, pp. 421-433 (1998) [ Links ]

Burggraf, O. R., Analytical and numerical studies of the structure of steady separated flows. J. Fluid Mech., 24, pp. 113-151 (1966). [ Links ]

Erturk, E., Corke, T. C., Gökçel, C., Numerical solution of 2-D steady incompressible driven cavity flow at high Reynolds numbers. International Journal for Numerical Methods in Fluids, 48, pp. 747-774 (2005) [ Links ]

Erturk, E., Discussion on driven cavity flow. International Journal for Numerical Methods in Fluids, 60, pp. 275-294 (2009). [ Links ]

Erturk, E., Gökçel, C., Fourth order compact formulation of Navier-Stokes equations and driven cavity flow at high Reynolds numbers. International Journal for Numerical Methods in Fluids, 50, pp. 421-436 (2006). [ Links ]

Freitas, C. J., Street, R. L., Findikakis, A. N., Koseff, J. R., Numerical simulation of three-dimensional flow in a cavity. International Journal for Numerical Methods in Fluids, 5, pp. 561-575 (1985). [ Links ]

Ghia, U., Ghia, K. N., Shin, C. T., High-Re solutions for incompressible flow using Navier-Stokes equations and multigrid method. Journal of Computational Physics, 48, pp. 387-411 (1982). [ Links ]

Hartmann, M., Perić, M., Scheuerer, G., Finite volume multigrid prediction of laminar natural convection: Bench-mark solutions. International Journal for Numerical Methods in Fluids, 11, pp. 189-207 (1990). [ Links ]

Hayese, T., Humphrey, J. A. C., Greif, R., A consistently formulated QUICK scheme for fast and stable convergence using finite-volume iterative calculation procedure. Journal of Computational Physics, 98, pp. 108-118 (1992). [ Links ]

Khosla, P. K., Rubin, S. G., A diagonally dominant second order accurate implicit scheme. Computers and Fluids, 2, pp. 207-209 (1974). [ Links ]

Leonard, B. P., A stable accurate convective modeling procedure based on quadratic upstream interpolation. Comput. Methods Appl. Mech. Eng., 19, pp. 59-98 (1979). [ Links ]

Lilek, Z., Muzaferija, S., Perić, M., Efficiency and accuracy aspects of a full-multigrid simple algorithm for three-dimensional flows. Numerical Heat Transfer, Part B, 31, pp. 23-42 (1997). [ Links ]

Lilek, Z., Perić, M., A fourth-order finite volume method with colocated variable arrangement. Computers and Fluids, 24, pp. 239-252 (1995). [ Links ]

Marchi, C. H., Suero, R., Araki, L. K., The lid-driven square cavity flow: Numerical solution with a 1024x1024 grid. J. of the Braz. Soc. of Mech. Sci & Eng., 31, pp. 186-198 (2009). [ Links ]

Nacer, B., David, L., Pascal, B., Gėrard, J., Contribution to the improvement of the QUICK scheme the resolution of the convection-diffusion problems. Heat Mass Transfer, 43, pp. 1075-1085 (2007). [ Links ]

Patankar, S. V., Numerical Heat Transfer and Fluid Flow. McGraw-Hill, New York (1980). [ Links ]

Ramšak, M., Škerget, L., A subdomain boundary method for high-Reynolds laminar flow using stream function-vorticity formulation. International Journal for Numerical Methods in Fluids, 46, pp. 815-847 (2004). [ Links ]

Sahin, M., Owens, R. G., A noverl fully implicit finite volume method applied to the lid-driven cavity problem-Part I: High Reynolds number flow calculations. International Journal for Numerical Methods in Fluids, 42, pp. 57-77 (2003). [ Links ]

Schreiber, R., Keller, H. B., Driven cavity flows by efficient numerical techniques. Journal of Computational Physics, 49, pp. 310-333 (1983). [ Links ]

Shankar, P. N., Deshpande, M. D., Fluid mechanics in the driven cavity. Annual Review of Fluid Mechanics, 32, pp. 93-136 (2000). [ Links ]

Shyy, W., Thakur, S., Wright, J., Second-order upwind and central difference schemes for recirculating flow computation. AIAA Journal, 30, pp. 923-932 (1992). [ Links ]

Versteeg, H. K., Malalasekera, W., An Introduction to Computational Fluid Dynamics: The Finite Volume Method. Prentice Hall (1995). [ Links ]

Wang, X., Yang, Z. F., Huang, G. H., High-order compact difference scheme for convection-diffusion problems on nonuniform grids. Journal of Engineering Mechanics, 131, pp. 1221-1228 (2005). [ Links ]

Yuguo, L, Baldacchino, L., Implementation of some higher-order convection schemes on non-uniform grids. International Journal for Numerical Methods in Fluids, 21, pp. 1201-1220 (1995). [ Links ]

*(Submitted: November 6, 2012 ; Revised: March 7, 2013 ; Accepted: March 7, 2013)*

* To whom correspondence should be addressed