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

PROCESS SYSTEMS ENGINEERING

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

K. YapiciI, ^*; Y. Uludag^II

^ICumhuriyet University, Department of Chemical Engineering, Phone: + 90.346.219 1010/1338, Fax: + 90.346.219 1179, 58140 Sivas, Turkey. E-mail: [e-mail]

^IIMiddle 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⁶. 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.

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

• 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).
• Barragy, E., Carey, G. F., Stream function-vorticity driven cavity solutions using p finite elements. Computers and Fluids, 26, pp. 453-468 (1997).
• Botella, O., Peyret, R., Benchmark spectral results on the lid-driven cavity flow. Computers & Fluids, 27, pp. 421-433 (1998)
• Burggraf, O. R., Analytical and numerical studies of the structure of steady separated flows. J. Fluid Mech., 24, pp. 113-151 (1966).
• 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)
• Erturk, E., Discussion on driven cavity flow. International Journal for Numerical Methods in Fluids, 60, pp. 275-294 (2009).
• 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).
• 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).
• 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).
• 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).
• 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).
• Khosla, P. K., Rubin, S. G., A diagonally dominant second order accurate implicit scheme. Computers and Fluids, 2, pp. 207-209 (1974).
• Leonard, B. P., A stable accurate convective modeling procedure based on quadratic upstream interpolation. Comput. Methods Appl. Mech. Eng., 19, pp. 59-98 (1979).
• 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).
• Lilek, Z., Perić, M., A fourth-order finite volume method with colocated variable arrangement. Computers and Fluids, 24, pp. 239-252 (1995).
• 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).
• 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).
• Patankar, S. V., Numerical Heat Transfer and Fluid Flow. McGraw-Hill, New York (1980).
• 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).
• 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).
• Schreiber, R., Keller, H. B., Driven cavity flows by efficient numerical techniques. Journal of Computational Physics, 49, pp. 310-333 (1983).
• Shankar, P. N., Deshpande, M. D., Fluid mechanics in the driven cavity. Annual Review of Fluid Mechanics, 32, pp. 93-136 (2000).
• Shyy, W., Thakur, S., Wright, J., Second-order upwind and central difference schemes for recirculating flow computation. AIAA Journal, 30, pp. 923-932 (1992).
• Versteeg, H. K., Malalasekera, W., An Introduction to Computational Fluid Dynamics: The Finite Volume Method. Prentice Hall (1995).
• 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).
• 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).

*

To whom correspondence should be addressed

Publication Dates

• Publication in this collection
  10 Jan 2014
• Date of issue
  Dec 2013

History

• Received
  06 Nov 2012
• Accepted
  07 Mar 2013
• Reviewed
  07 Mar 2013