A Computational Investigation of the Lid-Driven Cavity Flow

Nawal Odah Al-Atawi, Daoud Suleiman Mashat

Department of Mathematics, Numerical Analysis, King Abdulaziz University, Jeddah, Saudi Arabia.

DOI: 10.4236/ajcm.2022.122018                

Abstract

The lid-driven cavity is an important fluid mechanical system that serves as a benchmark for testing numerical methods and for studying fundamental aspects of incompressible flows in confined volumes. These flows are driven by the tangential motion of a bounding wall. The lid-driven cavity serves as a benchmark for testing numerical methods and for studying fundamental aspects of incompressible flows in confined volumes. This article presents a complete study of lid-driven cavity flows, with the primary focus being placed on the development of the flow when the Reynolds number was increased. In order to fully comprehend the physics of flow, it is necessary to take into consideration not only pure two-dimensional flows but also flows that are periodic in one space direction and the whole three-dimensional flow.

Keywords

Steady, Two Dimensional, Navier Stokes, Similarity Transformation, Finite Difference Scheme

1. Introduction

One of the simplest constrained geometries in which fluid motion can be examined is a rectangular or cubic container. Viscosity with constant density can be driven by tangential in plane motion of its enclosing wall, which is the most basic mechanical force [1]. A lid driven cavity (cuboid) in which one of the solid walls moves in an orthogonal manner to the other. The structure of the fluid flow in a hollow driven cavity by a sliding lid is a classical topic that has been the focus of extensive investigation for more than three decades [1] [2]. This problem has been studied by a broad variety of people. There are three explanations that might account for such a widespread interest [1] [2]. To begin, this issue has direct relation to many of the production techniques that are used in the coating sector, which is important to note from a purely practical standpoint. Second, since it is a prototype issue, many of the flow phenomena that are present in flows with closed streamlines may be investigated inside it. These flow features include a shear flow, boundary layers, eddies, and a core vortex [1] [2] [3]. Finally, since it has a straightforward geometry but displays a complicated flow behavior, it is an excellent environment for evaluating the accuracy and effectiveness of novel numerical approaches to the solution of the Navier-Stokes equations. This is because the geometry is straightforward [1] [2] [3]. Burggraf investigated the lid-driven cavity issue (see Figure 1) to examine whether the Prandtl-Batchelor theorem might be used. If a flow has closed lines, it will have a core of constant vorticity at infinite Reynolds numbers, according to the theory. It was Burggraf’s goal to prove this theorem quantitatively using finite differences. Until the Reynolds number of 400, he was able to achieve converging solutions using his numerical technique [4] [5] [6]. The vortices in the corners and a central vortex characterize the lid-driven cavity flow in this Reynolds number range (BR and BL vortex). The downstream eddy is another name for the BR vortex (DSE). Pan and Acrivos4 carried out flow experiments in square and rectangular cavities in order to derive the asymptotic behaviour of the flows. They experimented in the region of 20 to 4000 Reynolds numbers. The following is a summary of their findings about the flow in a cubical box. As the Reynolds number grows, core vortex moves to the geometric center of cavity and consumes a large portion of its volume. When Reynolds number 500 is reached, the downstream eddy’s size begins to decrease rapidly [6]. Their research was predicated on the premise that the observed flow patterns were two-dimensional [4] [5] [6].

Because of the ease with which it may be assembled, the lid-driven cavity has been the subject of a significant amount of research. Both as a numerical benchmark problem and as a test bed for the investigation of specific physical phenomena, it has been used. There are more than 1800 results returned when you search the term “lid-driven” on the Web of Science [7] [8]. A review on lid-driven cavity flows seems to be warranted for the aforementioned reasons, in addition to the rapid development of this area of research. Considering that approximately 20 years have elapsed since Shankar and Deshpande published an introduction of the topic, this seems to be the case. Following the first numerical investigations conducted by Kawaguti and Burggraf, the pursuit of efficiency and accuracy got underway with the work of Ghia [8] [9] [10], and Schreiber, who computed the steady two-dimensional flow for Reynolds numbers up to 10⁴ in a square cavity that was surrounded by three rigid walls and a lid that moved at a constant velocity. This was the beginning of the quest for efficiency and accuracy. Koseff and Street conducted a number of tests on the flow in three-dimensional cavities with varying lengths in the third dimension. A significant number of these studies were described in Koseff and Street. Dedicated three-dimensional test cases have been designed and analyzed numerically by several research groups, with the results gathered in Deville et al. These test cases were stimulated by the experimental results as well as the problems that still needed to be answered. After this joint effort, which did not yield very conclusive results for the targeted Reynolds number of Re = 3200, a new level of accuracy has been reached for two-dimensional flows by Botella and Peyret, who employed spectral methods combined with a dedicated treatment of the singular corners/edges where the moving wall meets with a stationary wall. Botella and Peyret’s work was able to achieve this new level of accuracy for two-dimensional flows by employing a dedicated treatment of the singular their approach produces a highly precise numerical solution for the two-dimensional issue up to the value of 10³, and it does so all the way up [8] [11] [12] [13]. Benchmarks for three-dimensional flows became interesting as computer power continued to advance and the process of routinely computing three-dimensional flows grew more common. Albensoeder and Kuhlmann provided highly accurate three-dimensional flow fields for Re = 10³ by applying the method of Botella and Peyret to three dimensions. This was done for a variety of cavity lengths in the spanwise direction as well as for rigid and periodic boundary conditions at the end walls. In addition to acting as a numerical benchmark, the lid-driven cavity problem is responsible for the emergence of a large number of fundamental fluid mechanics phenomena. The discontinuous boundary conditions along the edges at where moving and stationary walls meet are an important aspect for both an analytical and numerical approach of the issue. This issue is a particular instance of Taylor’s scraping problem, for which he has suggested solutions that are analogous to those previously proposed. At the apex of such an edge, the vorticity and the pressure move in opposite directions. This occurs because the boundary conditions for the velocity that is perpendicular to the edge are discontinuous [8] [13]. Closed-form solutions have been derived in terms of a series expansion of the steady flow at short distances away from the discontinuous corner for flows in two dimensions [8] [13].

Singular velocity and pressure fields are captured even in a truncated series expansion, which can be used to speed up the convergence of the entire numerical task of solving the Navier-Stokes equations [14] [15]. Regularization/smoothing of the discontinuity is an artifice to remove the singularity. Because of the pressure divergence that occurs when approaching the unique corners of the mathematical problem, there is an indication that experimental realizations of the lid-driven cavity will encounter difficulties [14] [15] [16]. Cavitation will occur at the corner when the wall slides away from the edge and the pressure decreases. Pumping and leakage can occur as a result of a gap between a fixed wall and a moving one. The viscous flow around the acute corners formed by two stationary walls is another oddity. A specific asymptotic behavior is observed in the form of an infinite succession of self-similar vortices whose size and intensity decrease geometrically as the edge is approached. Global vortex formations are also affected by these local effects, as well as changes in boundary motion and the cavity’s height to width ratio [14] [15] [16]. This is especially true in two-dimensional Stokes flow. In addition to numerical methods, the bi-orthogonal series method or Greens functions can be used to solve these types of issues. More difficult is the theoretical asymptotic approach of the three-dimensional problem of the local flow near a corner where three rigid walls meet [14] [15] [16].

The development of the two-dimensional flow in response to an increase in wall velocity is yet another essential component of the lid-driven cavity flow [10] [17]. The steady, two-dimensional flow should, in the absence of instabilities and according to the Prandtl-Batchelor theorem, evolve for large Reynolds numbers to a vortex with an inviscid core of uniform vorticity that is surrounded by viscous boundary layers that relate the vortex core to the boundary conditions [10] [17]. Experimentation and numerical analysis both backed up this pattern of behavior. However, the two-dimensional steady flow is not stable at high Reynolds numbers, and when the Reynolds number increases beyond a critical value, smaller-scale vortices are shed into the cavity from the downstream end of the moving wall [10] [17]. This occurs when the Reynolds number exceeds a threshold value. Also numerical investigation of a mixed convection flow in a lid driven cavity is done by using commercial finite volume package fluent to visualize the nature of the flow and estimate the heat transfer inside the cavity for different aspect ratio [18].

2. Governing System

The biological part of the modelling framework is based on the mass balance (continuity) equation for the state variables describing the transport (momentum conservation) and reaction among the species which can be propagated as

$\frac{\partial {\beta }_i}{\partial t}+\Delta \cdot \left({\beta }_i\right)=\Delta \cdot \left(D_e\Delta {\beta }_i\right)+R\left({\beta }_i\right)+S\left({\beta }_i\right)\mathrm{,}$ (1)

where, ${\beta }_i\left({\tau }_{\alpha }\mathrm{,}t\right)$ is a conservative quantities like cell density or concentration, v is the velocity of the flow field and ${\tau }_{\alpha }\in \Omega \subset R^3$ stands for a position vector in a Cartesian coordinate system [10] [19]. The diffusive coefficient is $D_e$ which describes the mixing in the system and it is influenced by the velocity profile and advection of the quantity. The reaction kinetics and temporal changes of the reacting species and described by the reaction term $R\left({\beta }_i\right)$ and the source terms $S\left({\beta }_i\right)$ respectively [10] [19]. The objective of this research article to identify the two numerical techniques with existence literature such that shortcoming and merits of each scheme recommended for incompressible flows. Two dimensional lid driven cavity flow which was originated by Burggraf. He was also unsure of the singularity at the two corners due to moving lid in contact with stationary walls, later on which were rectified that specification of the velocity of either unity or zero at the two corners, alters the numerical results significantly [10] [19]. Dimensionless approach typically generalized the problem for fluid mechanics point of view which stated that dimensionless solution depends on the set of dimensionless parameters, i.e. Reynolds number [10] [19]. Nondimensionalizing the equations of motion is our goal in this part. Since it will allow us to accurately compare the order and magnitudes of the various terms which involved in the aforementioned system [10] [19]. So, the incompressible continuity equation serves as a starting point,

$\nabla \cdot V=0.$ (2)

The Navier Stokes equations for incompressible Newtonian fluid with constant characteristics are;

$\rho \left[\frac{\partial V}{\partial t}+\left(V\cdot \nabla \right)V\right]=-\nabla P+\mu {\nabla }^{2}V+g,$ (3)

where, $g$ represents surface force (gravity), $g=0$ while $\mu$ is the viscosity of the fluid [20] [21] [22]. Also $\rho$ is the density while P is the pressure. In order to nondimensionalize the equations of motion, we now incorporate certain characteristic scaling parameters (reference) [20] [21] [22]. Base on scaling parameters, we then define various nondimensional variables and a nondimsional operator [20] [21] [22]. The scaling parameters for L (characteristic length), V (characteristic speed), f (characteristic frequency), $P_0-P_{\infty }$ (reference pressure difference) and g (gravitational acceleration) in term of primary dimensions are L, Lt⁻¹, t⁻¹, mLt⁻² and Lt⁻² respectively [20] [21] [22].

$t^{*}=ft,x^{*}=\frac{x}{L},V^{*}=\frac{V}{V},P^{*}=\frac{P-P_{\infty }}{P_0-P_{\infty }},{\nabla }^{*}=L\nabla .$ (4)

The Equation (3) can be rewritten as;

$\rho Vf\frac{\partial V^{*}}{\partial t^{*}}+\frac{\rho V^2}{L}\left(V^{*}\cdot {\nabla }^{*}\right)V^{*}=-\frac{P_0-P_{\infty }}{L}{\nabla }^{*}{P}^{*}+\frac{\mu V}{L^2}{\nabla }^{*2}{V}^{*}.$ (5)

To make all the terms dimensionless, we can multiply Equation (5) with $\frac{L}{\rho V^2}$ which lead to the following expression;

$\left[\frac{fL}{V}\right]\frac{\partial V^{\mathrm{<mo>\; *\; </mo>}}}{\partial t^{\mathrm{<mo>\; *\; </mo>}}}+\left(V^{*}\cdot {\nabla }^{*}\right)V^{\mathrm{<mo>\; *\; </mo>}}=-\left[\frac{P_0-P_{\infty }}{\rho V^2}\right]{\nabla }^{\mathrm{<mo>\; *\; </mo>}}{P}^{\mathrm{<mo>\; *\; </mo>}}+\left[\frac{\mu }{\rho VL}\right]{\nabla }^{\mathrm{<mo>\; *\; </mo>2}}{V}^{\mathrm{<mo>\; *\; </mo>}}\mathrm{.}$ (6)

According to the Buckingham pi theorem, each of the components in square brackets in the preceding equation is a nondimensional grouping of parameters. Each term in square bracket in Equation (6) represents nondimensional numbers which can be elaborated as;

$\left[\frac{fL}{V}\right]=St\mathrm{.}=\text{Strouhal}\text{Number}$

$\left[\frac{P_0-P_{\infty }}{\rho V^2}\right]=Eu\mathrm{.}=\text{Euler}\text{Number}$

$\left[\frac{\mu }{\rho VL}\right]=R{e}_L\mathrm{.}=\text{Reynolds}\text{Number}\mathrm{.}$

So, nondimensional Navier Stokes equation is,

$St\frac{\partial V^{*}}{\partial t^{*}}+\left(V^{*}\cdot {\nabla }^{*}\right)V^{*}=-Eu{\nabla }^{*}{P}^{*}+\left[\frac{1}{R{e}_L}\right]{\nabla }^{*2}{V}^{*}.$ (7)

Using special assumption to Equation (7), such as

$St=Eu=1$ (8)

2.1. Assumptions

Some assumptions on the governing Equation (7) of flow can be written as;

● Steady, two dimension such that $\frac{\partial }{\partial z}=0$ and $w=0$.

● Density is constant which implies $\frac{\partial \rho }{\partial t}=0$.

● No mass or surface forces such that $f_x=0$ and $f_y=0$.

2.2. Research System

The governing Equation (7) of flow can be written as in components form so, two dimensional incompressible continuity equation (CE) is

$\frac{\partial \tilde{u}}{\partial \tilde{x}}+\frac{\partial \tilde{v}}{\partial \tilde{y}}=0.$ (9)

Momentum equation along x direction (MOM_x) is,

$\tilde{u}\frac{\partial \tilde{u}}{\partial \tilde{x}}+\tilde{v}\frac{\partial \tilde{u}}{\partial \tilde{y}}=-\frac{\partial \tilde{p}}{\partial \tilde{x}}+\frac{1}{R{e}_L}\left(\frac{{\partial }^2\tilde{u}}{\partial {\tilde{x}}^2}+\frac{{\partial }^2\tilde{u}}{\partial {\tilde{y}}^2}\right).$ (10)

Momentum equation along y direction (MOM_y) is,

$\tilde{u}\frac{\partial \tilde{v}}{\partial \tilde{x}}+\tilde{v}\frac{\partial \tilde{v}}{\partial \tilde{y}}=-\frac{\partial \tilde{p}}{\partial \tilde{y}}+\frac{1}{R{e}_L}\left(\frac{{\partial }^2\tilde{v}}{\partial {\tilde{x}}^2}+\frac{{\partial }^2\tilde{v}}{\partial {\tilde{y}}^2}\right).$ (11)

3. Methodology

We have to discretize these Equations (9) $\mapsto$ (11) to be able to handle in computers. We use the SIMPLE algorithm to solve the research problem [23] [24] [25] [26]. A flow chart of the algorithm is shown in Figure 1. To solve an incompressible flow problem which is complicated because of the lack of an independent equation of pressure. In the SIMPLE algorithm, the momentum equations are solved using a guessed pressure field [25] [26] [27] [28]. Then, the correction part of the pressure is derived by using the equation from continuity equation. After this, the new velocity and pressure fields can be derived. The fields can approach to the exact values as the iteration runs until its convergence. This algorithm is used in many commercial CFD codes. We also use the algorithm because the algorithm is suitable way to solve the fluid governing equations when it comes to the incompressible flow [28] [29]. The discretized equations with non-dimensional form of variables are shown below.

There are many computational techniques that are most often employed in computational fluid dynamics. These include finite difference, finite volume, finite

Figure 1. The SIMPLE algorithm [29].

element and their variations [30]. The current study focuses mostly on the finite difference approach for computing solutions to the aforementioned systems (9) $\mapsto$ (11) in term of variables (u,v & p) [8] [30] [31]. The discretized equations with non-dimensional form of variables are shown below (1).

$CE\mathrm{<mo>\; :\; </mo>}\frac{u_{i\mathrm{,}j}-u_{i-\mathrm{1,}j}}{\Delta x}+\frac{v_{i\mathrm{,}j}-v_{i-\mathrm{1,}j}}{\Delta y}=\mathrm{0,}$ (12)

$MO{M}_x\mathrm{<mo>\; :\; </mo>}{u}_{i\mathrm{,}j}^{n+1}=u_{i\mathrm{,}j}^n+\Delta t{A}_{i\mathrm{,}j}^n-\Delta t\frac{p_{i+\mathrm{1,}j}^n-p_{i\mathrm{,}j}^n}{\Delta x}\mathrm{,}$ (13)

$MO{M}_y\mathrm{<mo>\; :\; </mo>}{v}_{i\mathrm{,}j}^{n+1}=v_{i\mathrm{,}j}^n+\Delta t{B}_{i\mathrm{,}j}^n-\Delta t\frac{p_{i\mathrm{,}j+1}^n-p_{i\mathrm{,}j}^n}{\Delta y},$ (14)

where

$\begin{array}{c}{A}_{i\mathrm{,}j}^n=-\left(u_{i\mathrm{,}j}^n\frac{u_{i+\mathrm{1,}j}^n-u_{i-\mathrm{1,}j}^n}{2\Delta x}+\bar{v}\frac{u_{i\mathrm{,}j-1}^n-u_{i\mathrm{,}j+1}^n}{2\Delta y}\right)\\ +\frac{1}{R{e}_L}\left(\frac{u_{i+\mathrm{1,}j}^n-2u_{i\mathrm{,}j}^n+u_{i-\mathrm{1,}j}^n}{{\left(\Delta x\right)}^2}+\frac{u_{i\mathrm{,}j+1}^n-2u_{i\mathrm{,}j}^n+u_{i\mathrm{,}j-1}^n}{{\left(\Delta y\right)}^2}\right)\end{array}$

$\begin{array}{c}{B}_{i\mathrm{,}j}^n=-\left(\bar{u}\frac{v_{i+\mathrm{1,}j}^n-v_{i-\mathrm{1,}j}^n}{2\Delta x}+v_{i\mathrm{,}j}^n\frac{v_{i\mathrm{,}j-1}^n-v_{i\mathrm{,}j+1}^n}{2\Delta y}\right)\\ +\frac{1}{R{e}_L}\left(\frac{v_{i+\mathrm{1,}j}^n-2v_{i\mathrm{,}j}^n+v_{i-\mathrm{1,}j}^n}{{\left(\Delta x\right)}^2}+\frac{v_{i\mathrm{,}j+1}^n-2v_{i\mathrm{,}j}^n+v_{i\mathrm{,}j-1}^n}{{\left(\Delta y\right)}^2}\right)\end{array}$

$\bar{u}=\frac{1}{4}\left(u_{i\mathrm{,}j}^n+u_{i-\mathrm{1,}j}^n+u_{i\mathrm{,}j+1}^n+u_{i-\mathrm{1,}j+1}^n\right)$

$\bar{v}=\frac{1}{4}\left(v_{i\mathrm{,}j}^n+v_{i+\mathrm{1,}j}^n+v_{i\mathrm{,}j-1}^n+v_{i+\mathrm{1,}j-1}^n\right).$

The corrected values are defined as

$\begin{array}{l}u=u^{*}+u^{\prime },\\ v=v^{*}+v^{\prime },\end{array}$ (15)

and

$p=p^{*}+p^{\prime }.$ (16)

Now correction terms are define as;

$\begin{array}{l}{u^{\prime }}_{i\mathrm{,}j}^{n+1}=-\Delta t\frac{{p^{\prime }}_{i+\mathrm{1,}j}^n-{p^{\prime }}_{i\mathrm{,}j}^n}{\Delta x}\mathrm{,}\\ {v^{\prime }}_{i\mathrm{,}j}^{n+1}=-\Delta t\frac{{p^{\prime }}_{i\mathrm{,}j+1}^n-{p^{\prime }}_{i\mathrm{,}j}^n}{\Delta y}\end{array}$ (17)

${p^{\prime }}_{i,j}=\frac{{\left(\Delta x\right)}^2\cdot {\left(\Delta y\right)}^2}{{\left(\Delta x\right)}^2+{\left(\Delta y\right)}^2}\cdot \left[\frac{{p^{\prime }}_{i+1,j}^n-{p^{\prime }}_{i-1,j}^n}{{\left(\Delta x\right)}^2}+\frac{{p^{\prime }}_{i,j+1}^n-{p^{\prime }}_{i,j-1}^n}{{\left(\Delta y\right)}^2}-D_{i,j}\right]$ (18)

where

$D_{i\mathrm{,}j}=\frac{1}{\Delta t}\left(\frac{u_{i\mathrm{,}j}^{\mathrm{<mo>\; *\; </mo>}}-u_{i-\mathrm{1,}j}^{\mathrm{<mo>\; *\; </mo>}}}{\Delta x}+\frac{v_{i\mathrm{,}j}^{\mathrm{<mo>\; *\; </mo>}}-v_{i\mathrm{,}j-1}^{\mathrm{<mo>\; *\; </mo>}}}{\Delta y}\right)$

4. Results and Discussion

In this article, a finite difference numerical approach has been discussed in which two dimensional coupled system is studied. Validation of the results has been done by using computer program (MATLAB) by making comparison with benchmark. The numerical scheme (Equations (12) to (18)) using finite difference methodolgy is applied on two dimensional steady incompressible system. Figure 2 describes the results at Reynolds number Re = 40 in which the primary vortex can be seen from u-component of the velocity while two corner vortices in v-component of velocity. Figure 3 describes the results at Reynolds number Re = 400 in which the primary vortex can be seen from u-component of the velocity while two corner vortices in v-component of velocity while Figure 4 describes the results at Reynolds number Re = 1000 in which the primary vortex can be seen from u-component of the velocity while two corner vortices in v-component of velocity. Figure 5 represents a pressure profile for incompressible two dimensional flow. Second part of this article is to focus on the accuracy of the numerical scheme which linked with convergence criteria. Since, convergence criteria define the accuracy of the result of the steady state problem [32] [33] [34]. The more converged, the more accurate results we could get. The “Residual” means the imbalance of the equations. Residuals would decrease as the calculations are proceeding and equations are converging [32] [33] [34]. When

Figure 2. Velocity components such as u and v with vortex generation at Re = 40.

Figure 3. Velocity components such as u and v with vortex generation at Re = 400.

Figure 4. Velocity components such as u and v with vortex generation at Re = 1000.

Figure 5. Pressure profile at Re = 50.

the residuals reach a specified tolerance, then the calculation recognized as being converged and when all the residuals reach the tolerance, the calculation will stop. Small tolerance gives accurate results but too small tolerance requires huge number of iteration step. We tried to compute with some tolerances and compare with the results. From this, we found that the velocity distributions are almost same with each case which is 10⁻³, 10⁻⁴ and 10⁻⁵ as the tolerance [32] [33] [34]. Since we can get almost same results, so 10⁻³ as the convergence tolerance because the larger value of the tolerance gives shorter calculation time as mentioned above which always reduce the calculation time [32] [33] [34]. Now the velocity of the moving wall is 1 m/s which indicate that 10⁻³ is at least qualitative convergence. However, even though the literature also says that stricter convergence consideration is required for transport variables, thus 10⁻³ is enough because it gives, almost same results with the 10⁻⁵ in shorter time [32] [33] [34].

5. Conclusion

Fluid dynamics can be studied in great depth using the lid-driven cavity problem. It is widely used for both qualitative and quantitative benchmarking, which has resulted in a wealth of data. However, there are still many unanswered questions that will require further examination. In constrained systems, the stability and transition scenario remains a big concern. There are still a lot of unanswered questions about the flow’s dependence on the spanwise confinement, which is given by the span aspect ratio. Analysis of the global stability of three-dimensional flows in constrained geometries shows promising techniques. Additionally, the Lagrangian topology and flow kinematics have recently attracted a lot of attention. The mixing qualities are determined by these flow features, as well as by the flow’s chaotic and regular sections. For three-dimensional flows, numerical computations are expensive because of the high level of accuracy required. Streamline topology in steady three-dimensional flows has only recently been possible to address.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Schneider, T.M., Gibson, J.F., Lagha, M., Lillo, F.D. and Eckhardt, B. (2008) Laminar Turbulent Boundary in Plane Couette Flow. Physical Review E, 78, Article ID: 037301. https://doi.org/10.1103/PhysRevE.78.037301
[2]	Kawaguti, M. (1961) Numerical Solution of the Navier-Stokes Equations for the Flow in a Two-Dimensional Cavity. Journal of the Physical Society of Japan, 16, 2307-2315. https://doi.org/10.1143/JPSJ.16.2307
[3]	Burggraf, O.R. (1966) Analytical and Numerical Studies of the Structure of Steady Separated Flows. Journal of Fluid Mechanics, 24, 113-151. https://doi.org/10.1017/S0022112066000545
[4]	Ghia, U., Ghia, K.N. and Shin, C.T. (1982) High-Re Solutions for Incompressible Flow Using the Navier-Stokes Equations and a Multigrid Method. Journal of Computational Physics, 48, 387-411. https://doi.org/10.1016/0021-9991(82)90058-4
[5]	Schreiber, R. and Keller, H.B. (1983) Driven Cavity Flows by Efficient Numerical Techniques. Journal of Computational Physics, 49, 310-333. https://doi.org/10.1016/0021-9991(83)90129-8
[6]	Koseff, J.R. and Street, R.L. (1984) The Lid-Driven Cavity Flow: A Synthesis of Qualitative and Quantitative Observations. Journal of Fluids Engineering, 106, 390-398. https://doi.org/10.1115/1.3243136
[7]	Deville, M., Le, T.-H. and Morchoisne, Y. (1992) Numerical Simulation of 3-D Incompressible Unsteady Viscous Laminar Flows. Volume 36 of Notes on Numerical Fluid Mechanics. Vieweg, Braunschweig. https://doi.org/10.1007/978-3-663-00221-5
[8]	Botella, O. and Peyret, R. (1998) Benchmark Spectral Results on the Lid-Driven Cavity Flow. Computers & Fluids, 27, 421-433. https://doi.org/10.1016/S0045-7930(98)00002-4
[9]	Benjamin, A.S. and Denny, V.E. (1979) On the Convergence of Numerical Solutions for 2-D Flows in a Cavity at Large Reynolds Numbers. Journal of Computational Physics, 33, 340-358. https://doi.org/10.1016/0021-9991(79)90160-8
[10]	Erturk, E., Corke, T.C. and Gokcol, C. (2005) Numerical Solutions of 2-d Steady Incompressible Driven Cavity Flow at High Reynolds Numbers. International Journal for Numerical Methods in Fluids, 48, 747-774. https://doi.org/10.1002/fld.953
[11]	Albensoeder, S. and Kuhlmann, H.C. (2005) Accurate Three-Dimensional Lid-Driven Cavity Flow. Journal of Computational Physics, 206, 536-558. https://doi.org/10.1016/j.jcp.2004.12.024
[12]	Taylor, G.I. (1962) On Scraping Viscous Fluid from a Plane Surface. In: Batchelor, G.K., Ed., The Scientific Papers of Sir Geoffrey Ingram Taylor, Volume IV: Mechanics of Fluids—Miscellaneous Papers, Cambridge University Press, 1971, Cambridge, 410-413.
[13]	Gupta, M.M., Manohar, R.P. and Noble, B. (1981) Nature of Viscous Flows near Sharp Corners. Computers & Fluids, 9, 379-388. https://doi.org/10.1016/0045-7930(81)90009-8
[14]	Riedler, J. and Schneider, W. (1983) Viscous Flow in Corner Regions with a Moving Wall and Leakage of Fluid. Acta Mechanica, 48, 95-102. https://doi.org/10.1007/BF01178500
[15]	Batchelor, G.K. (1956) On Steady Laminar Flow with Closed Streamlines at Large Reynolds Numbers. Journal of Fluid Mechanics, 1, 177-190. https://doi.org/10.1017/S0022112056000123
[16]	Shih, T.M., Tan, C.H. and Hwang, B.C. (1989) Effects of Grid Staggering on Numerical Schemes. International Journal for Numerical Methods in Fluids, 9, 193-212. https://doi.org/10.1002/fld.1650090206
[17]	Bozeman and Dalton (1973) Numerical Study of Viscous Flow in Cavity. Journal of Computational Physics, 12, 348-363. https://doi.org/10.1016/0021-9991(73)90157-5
[18]	Omari, R. (2016) Numerical Investigation of a Mixed Convection Flow in a Lid-Driven Cavity. American Journal of Computational Mathematics, 6, 251-258. https://doi.org/10.4236/ajcm.2016.63026
[19]	Goodrich, J.W., Gustafson, K. and Halasi, K. (1990) Hopf Bifurcation in the Driven Cavity. Journal of Computational Physics, 90, 219-261. https://doi.org/10.1016/0021-9991(90)90204-E
[20]	Yates, G.T. (1986) How Microorganisms Move through Water. American Scientist, 74, 358-365.
[21]	Heinsohn, R.J. and Cimbala, J.M. (2003) Indoor Air Quality Engineering. Marcel-Dekker, New York. https://doi.org/10.1201/9780203911693
[22]	Kundu, P.K. (1990) Fluid Mechanics. Academic Press, San Diego.
[23]	Van Dyke, M. (1982) An Album of Fluid Motion. The Parabolic Press, Stanford. https://doi.org/10.1115/1.3241909
[24]	White, F.M. (1991) Viscous Fluid Flow. 2nd Edition, McGraw-Hill, New York.
[25]	Saqib, M., Hasnain, S. and Mashat, D. (2017) Highly Efficient Computational Methods for Two Dimensional Coupled Non-Linear Unsteady Convection-Diffusion Problems. IEEE Access, 5, 7139-7148. https://doi.org/10.1109/ACCESS.2017.2699320
[26]	Hasnain, S., Bashir, S., Linker, P. and Saqib, M. (2019) Efficiency of Numerical Schemes for Two Dimensional Gray Scott Model. AIP Advances, 9, 105-121. https://doi.org/10.1063/1.5095517
[27]	Hasnain, S., Saqib, M., Afzal, M.F. and Hussain, I. (2019) Numerical Study to Coupled Three Dimensional Reaction Diffusion System. IEEE Access, 7, 46695-46705.
[28]	Harlow, F.H. and Welch, J.E. (1965) Numerical Calculation of Time Dependent Viscous Incompressible Flow of Fluid with Free Surface. The Physics of Fluids, 3, 2182-2189. https://doi.org/10.1063/1.1761178
[29]	Chorin, A.J. (1968) Numerical Solution of Navier-Stokes Equations. Mathematics of Computation, 22, 745-762. https://doi.org/10.1090/S0025-5718-1968-0242392-2
[30]	Gupta, M.M. (1991) High Accuracy Solutions of Incompressible Navier-Stokes Equations. Journal of Computational Physics, 93, 343-359. https://doi.org/10.1016/0021-9991(91)90188-Q
[31]	Rubin, S.G. and Khosla, P.K. (1981) Navier-Stokes Calculations with a Coupled Strongly Implicit Method. Computers and Fluids, 9, 163-180. https://doi.org/10.1016/0045-7930(81)90023-2
[32]	Goyon, O. (1996) High-Reynolds Number Solutions of Navier-Stokes Equations Using Incremental Unknowns. Computer Methods in Applied Mechanics and Engineering, 130, 319-335. https://doi.org/10.1016/0045-7825(95)00923-X
[33]	Wan, D., Zhou, Y. and Wei, G. (2002) Numerical Solutions of Incompressible Flows by Discrete Singular Convolution. International Journal for Numerical Methods in Fluids, 38, 789-810. https://doi.org/10.1002/fld.253
[34]	Marchi, C. and Silva, H. (2002) Unidimensional Numerical Solution Error Estimation for Convergent Apparent Order. Numerical Heat Transfer, Part B, 42, 167-188. https://doi.org/10.1080/10407790190053888