The lid-driven square cavity flow: numerical solution with a 1024 x 1024 grid

TECHNICAL PAPERS

The lid-driven square cavity flow: numerical solution with a 1024 x 1024 grid

Carlos Henrique Marchi^I; Roberta Suero^II; Luciano Kiyoshi Araki^III

^IEmeritus Member, ABCM, [e-mail], Federal University of Paraná - UFPR, Department of Mechanical Engineering, Caixa postal 19040, 81531-980 Curitiba, PR, Brazil

^[e-mail], Federal University of Paraná - UFPR, Post-Graduate Program in Numerical Methods in Engineering

^[e-mail], Federal University of Paraná - UFPR, Department of Mechanical Engineering

ABSTRACT

The problem of flow inside a square cavity whose lid has constant velocity is solved. This problem is modeled by the Navier-Stokes equations. The numerical model is based on the finite volume method with numerical approximations of second-order accuracy and multiple Richardson extrapolations (MRE). The iterative process was repeated until the machine round-off error achievement. This work presents results for 42 variables of interest, and their discretization errors estimates, on the 1024 x 1024 nodes grid and the following Reynolds numbers: 0.01, 10, 100, 400 and 1000. These results are compared with sixteen sources in literature. The numerical solutions of this work are the most accurate obtained for this problem to date. The use of multiple Richardson extrapolations reduces the discretization error significantly.

Keywords: discretization error, error estimate, CFD, Richardson extrapolation, finite volume method

Introduction

This work addresses the classical problem (Kawaguti, 1961; Burggraf, 1966; Rubin and Khosla, 1977; Benjamin and Denny, 1979; Ghia, Ghia and Shin, 1982) of laminar flow inside a square cavity of which lid moves at constant velocity: Fig. 1; where u and v are the components of the velocity vector in x and y directions, ρ and µ are fluid density and viscosity. This problem is widely employed to evaluate numerical methods and to validate codes for solving the Navier-Stokes equations (Botella and Peyret, 1998). In the works cited in Table 1, the problem was solved for 11 x 11 up to 2048 x 2048 node grids, and for Reynolds numbers (Re) from zero to 21,000.

As can be seen in Table 1, several numerical methods have been used, including finite difference method (FDM), finite volume method (FVM), finite elements method (FEM), lattice Boltzmann (LB), and the spectral method (Spectral). In addition, a variety of mathematical formulations has been used, including stream function and vorticity (ψ -ω); stream function and velocity (ψ -V); lattice Boltzmann equation (LBE), and the Navier-Stokes equations (u-v-p). The problem considered here is also known as "singular driven cavity" (Botella and Peyret, 1998), because there are two discontinuities in the boundary condition of u, at lid corners: 0 in one side and 1 in another. In contrast, there is a problem called "regularized driven cavity" (Botella and Peyret, 1998), which does not present discontinuities.

The main objective of this work is to present the most accurate numerical solutions found to date for the problem of "singular driven cavity" with Re = 0.01, 10, 100, 400 and 1000. To achieve this aim, we use the Navier-Stokes equations; the finite volume method; co-located arrangement of variables; segregated solution for the three conservation equations; numerical approximations of second-order accuracy; 1024 x 1024 control volumes uniform grid; the iterative process repeated until the machine round-off error achievement; double precision in calculations; and multiple Richardson extrapolations (Richardson, 1910). Solutions are presented for 42 variables of interest, which involve velocity profiles, mass flow rate, minimum value of the stream function, minimum and maximum velocities (and their coordinates), and wall forces on the fluid.

Other objectives of this work are: (1) propose an error estimator for use with numerical solutions obtained through multiple Richardson extrapolations; (2) verify (Roache, 1998) if the proposed estimator provides reliable error estimates for a problem of which analytical solution is known (Shih, Tan and Hwang, 1989); (3) apply the proposed estimator to each of the 42 variables of interest and five values of Re, presenting the estimated discretization error for each numerical solution; (4) confirm the order of accuracy (p_L) of the numerical solutions; and (5) compare the results with sixteen sources in literature. This work does not have as aim to present an optimized numerical model neither for CPU time nor for computational memory consumption. It is emphasized that the main objective is to provide the most accurate results to date for literature.

Although there is extensive literature on the problem considered here, this work is justified by the following reasons:

No work appears to have been developed to date to estimate the numerical error involved in the solution of each variable of interest (U_h in Table 2). This is important, however, in order to know the reliability of numerical solutions, allowing more careful comparisons to be made. Some authors have presented the solution variation for some variables for two or three grids; they, however, did not estimate numerical errors with any discretization error estimator.

Only Bruneau and Saad (2006) and Wright and Gaskell (1995) present solutions on grids as fine as those of the present work, but only for Re = 1000, few variables and without using multiple Richardson extrapolations. In Table 2, the column RE indicates whether the authors have used Richardson extrapolation or not; if the answer is positive, it is presented how many times it was used for each variable. In this work, it is going to be shown that RE reduces significantly the discretization error.

Most authors have stopped their iterative process (Ui) based on the residue criterion (R) or on the variation of any variable (Ferziger and Peric, 1999), with the tolerance value between 1.0e-3 and 1.0e-12. In the present work, the iterative process is taken until the machine round-off error achievement.

In Table 2, "digits" represents the number of significant figures for each solution, which has the highest value in the presented work.

In Table 2, p_L is the theoretical accuracy of discretization error of the approximations employed by each author. In the present work, the practical value obtained for this order is shown, for each variable of interest, confirming or not the theoretical value. In the next sections, the following subjects are discussed: the mathematical and numerical models; the theory and equations used to calculate effective and apparent orders of error to perform multiple Richardson extrapolations and the discretization error estimator; the results for the problem with known analytical solution; the classical problem results; and conclusions of this work.

Nomenclature

E = true discretization error F = viscous drag force of the cavity's lid or wall on the fluid in x direction (N) h = size of the control volumes (m) m = number of Richardson extrapolation M = mass flow rate (kg/s) nm = number of grids p = pressure (Pa) p_E = effective order of the discretization error p_L = asymptotic order of the discretization error p_U = apparent order of the discretization error r = grid refinement ratio Re = Reynolds number S = source term in Eq. (3) u = component of the velocity vector in the x direction (m/s) U = estimated discretization error v = component of the velocity vector in the y direction (m/s) x = Cartesian coordinate in the horizontal direction (m) y = Cartesian coordinate in the vertical direction (m) z = depth of the cavity (m) Greek symbols = numerical solution of the variable of interest Φ = exact analytical solution of the variable of interest µ = viscosity (Pa.s) ρ = density (kg/m3) ψ = stream function (m2/s) Subscripts 1 fine grid 2 coarse grid 3 supercoarse grid

Mathematical Model

The mathematical model of the problem consists of the Conservation of Mass and the Conservation of Momentum laws (the Navier-Stokes equations). Simplifications considered for the problem are: steady state; two-dimensional laminar flow in x and y directions; incompressible fluid; µ as constant; and absence of other effects. Thus, the resulting mathematical model is:

where p is the pressure and S is a source term, which is null at the classical problem (without known analytical solution) and given by Shih, Tan and Hwang (1989) in the case of a manufactured solution problem. The domain is a square of unitary side with the origin of the system of coordinates, as shown in Fig. 1.

The variables of interest of the problem involve the primitive variables themselves (u and v) and integrations of u and v, which are:

Profile of u in x = ½ at 15 selected points of y.

Profile of v in y = ½ at 15 selected points of x.

The minimum value (u_min) of profile of u in x = ½ and its respective y coordinate.

The minimum (v_min) and the maximum (v_max) values of profile of v in y = ½ and their respective x coordinates.

The minimum value of the stream function (ψ _min) and its coordinates x and y.

The mass flow rate (M) that flows through y = ½ line between x = 0 and ½ , i.e.,

where z is the cavity depth, which is considered unitary.

The viscous drag force (F) in direction x (Hou et al., 1995) is the force exerted by the fluid boundary surface, calculated by

where F_n is F in y = 1 (cavity movable lid) and F_s is F in y = 0 (the lower cavity wall).

Numerical Model

Briefly, the numerical model adopted to solve the mathematical model described by Eqs. (1) to (3) has the following characteristics: (1) finite volume method (Ferziger and Peric, 1999); (2) central difference scheme (CDS) (Ferziger and Peric, 1999) for diffusive and pressure terms; (3) CDS scheme with deferred correction (Ferziger and Peric, 1999; Khosla and Rubin, 1974) on upstream difference scheme (UDS) for advective terms; (4) Eqs. (1) to (3) are solved sequentially using the MSI (Modified Strongly Implicit) method (Schneider and Zedan, 1981); (5) SIMPLEC (Semi IMPlicit Linked Equations Consistent) method (Van Doormaal and Raithby, 1984) to treat the pressure-velocity coupling; (6) uniform grids; (7) the boundary conditions for u and v, Fig. 1, are applied employing ghost volumes (Ferziger and Peric, 1999); (8) Eqs. (1) to (3) are written for unsteady state, aiming the use of time as a relaxation parameter in the iterative solution process of the discretized mathematical model; and (9) a co-located arrangement of variables is used (Marchi and Maliska, 1994). This numerical model does not require boundary conditions for pressure (Ferziger and Peric, 1999).

The numerical solution of the variables of interest is obtained as follows:

• The numerical solution of the profile of u in x = ½ is obtained by the mathematical average of u stored at the east face of the two adjacent volumes to each desired y coordinate. This u at each control volume east face is that one of the co-located arrangement of variables of Marchi and Maliska (1994). This is necessary because the number of volumes used in each coordinate direction is even, so that no control volume center coincides with the line x = ½.

• The numerical solution of the profile of v in y = ½ is obtained analogously to the profile of u, through the mathematical average of v stored at the north face of the two adjacent volumes to each desired x coordinate.

• u_min is the minimum value of the solution of u stored at east faces among all control volumes of the grid with x = ½. And its y coordinate is the one in the east face center of the volume corresponding to u_min.

• v_min and v_max are the minimum and the maximum values of the solution of v stored at the north face among all control volumes of the grid with y = ½. And their x coordinates are the ones in the north face center of the volumes corresponding to v_min and v_max.

• For each vertical line, numerical solution of the stream function field (ψ) is obtained through integration of the product of u, stored at each control volume east face, by the height of each control volume (Δy), starting from the lower wall, in y = 0. The vertical lines coincide with x coordinates of the faces of each control volume. The numerical integration used here is based on the rectangular rule (Kreyszig, 1999). The minimum value of the stream function (ψ_min) is obtained directly from the ψ field.

• The numerical solution of the mass flow rate (M) defined by Eq. (4) is obtained by numerical integration employing the rectangular rule, through

where i represents the number of the control volume in x direction; i = 1 is the real control volume at the left-hand wall of the cavity; N_x is the total number of real control volumes in x direction; Δ x is the width of each control volume; and v_n is v at the north face of each control volume.

• The numerical solution of the force (

  Fn), which is defined by Eq. (5), is obtained with one upstream point (UDS) (Ferziger and Peric, 1999) and the numerical integration used here is the rectangular rule, which results in

where u_T,i is the velocity of the cavity lid at the x coordinate of each control volume i center; and is the nodal velocity u at each real i control volume center, whose volume north face coincides with the cavity lid.

• The numerical solution of

  Fs, defined by Eq. (5), is obtained analogously to

  Fn but with one point downstream (DDS) (Ferziger and Peric, 1999), resulting in

where u_i,1 is the nodal velocity u at each real i control volume center, whose volume south face coincides with the cavity lower wall, whose u = 0.

The Stokes 1.5 computational code was implemented in Fortran 95, using Intel Fortran 9.1 compiler and double precision. The iterative process is repeated until the machine round-off error achievement. This is verified by monitoring the l₁-norm (Kreyszig, 1999), along the iterations, of the residue (R) sum of the three solved systems of Eqs. (1) to (3). The residue sum value of the three systems, in each outer iteration, is nondimensionalized by its value at the end of the first outer iteration (R1).

Discretization Error

The numerical error (E) can be defined as the difference between the exact analytical solution (Φ) of a variable of interest and its numerical solution (), i.e.,

The sources that cause numerical error can be split into four types (Marchi and Silva, 2002): truncation, iteration, round-off, and programming errors. When other sources of error are inexistent or much smaller in relation to truncation ones, the numerical error can also be called as discretization error.

Considering the numerical model described in the previous section, the predicted asymptotic order (p_L) of the discretization error is equal to two (Ferziger and Peric, 1999; Schneider, 2007) for all variables of interest, except for x and y coordinates of ψ _min, u_min, v_min and v_max, of which p_L are unknown. In literature (Roache, 1998, 1994), p_L is also called formal order or accuracy order.

In theory (Marchi, 2001), it is expected that p_E (effective order) and p_U (apparent order) → p_L for h → 0. In other words, it is expected that the practical orders (p_E and p_U), which are calculated with the numerical solutions of each variable of interest, tend toward the asymptotic order (p_L), foreseen a priori, when the size of the control volumes (h) tends toward zero.

The effective order (p_E) of the true error is defined by (Marchi, 2001)

where E(₁) and E(₂) are true discretization errors of the numerical solutions ₁ and ₂ obtained, respectively, with fine (h₁) and coarse (h₂) grids; h = size of the control volumes (in this work, h = Δx = Δy); and r = h₂/h₁ (grid refinement ratio).

According to Eq. (10), the effective order (p_E) is function of the true discretization error of a variable of interest. Thus, for problems which analytical solution is known, it can be used to verify a posteriori if, as h → 0, one obtains p_L. When E is unknown, (p_E) cannot be calculated. In this case, one can use the concept of observed or apparent order (p_U) defined by (De Vahl Davis, 1983; Marchi and Silva, 2002)

where

₁, ₂ and ₃ = numerical solutions obtained, respectively, with fine (h₁), coarse (h₂) and supercoarse (h₃) grids, and r = h₃/h₂ = h₂/h₁.

Some studies (Benjamin and Denny, 1979; Schreiber and Keller, 1983; Erturk, Corke and Gökçöl, 2005) achieved excellent results when employing multiple Richardson extrapolations (MRE) to reduce the discretization error of ψ _min. However, these authors used this process with at most four grids, resulting in up to three extrapolations for the finest grid they used. In the present work, this procedure was utilized with up to ten grids, resulting in up to nine extrapolations for the finest grid used (1024 x 1024), and applied to almost all variables of interest. This was done by means of

where

_1,m∞ is the numerical solution of the variable of interest () with m extrapolations on the fine grid (h₁); _1,(m-1)∞ and _2,(m-1)∞ are numerical solutions with (m-1) extrapolations on the fine (h₁) and coarse (h₂) grids; r = h₂/h₁ (grid refinement ratio); m = number of Richardson extrapolations, with m = 0 being the numerical solution obtained in grid h without any extrapolation; nm = number of different grids with numerical solutions of without any extrapolation; p_V(m) = true orders (Marchi and Silva, 2002) of the discretization error, with p_V(1) = p_L. For the numerical model used in this work, p_V = 2, 4, 6 ... for all variables of interest, except for x and y coordinates of ψ_min, u_min, v_min and v_max, of which values of p_V are unknown.

In this work, Eq. (12) is applied to all variables of interest to reduce the discretization error, except to x and y coordinates of ψ_min, u_min, v_min and v_max, of which results are the ones obtained with the finest grid, of 1024 x 1024 nodes, without any extrapolation. In theory, the accuracy order of the results of u(0.5;05), v(0.5;0.5), M, Fs, Fn, v_min, v_max, u_min and ψ_min is 20, since they are obtained with nine extrapolations through Eq. (12). And, in the case of the profiles of u and v, it is 14, because they are obtained with six extrapolations.

In practical situations, a numerical solution is obtained because the analytical solution is unknown. Hence, the true value of the numerical error is also unknown. Therefore, the numerical error must be estimated. The estimated discretization error (U) of _1,(nm-1)∞, i.e., of the numerical solution with the highest possible number of extrapolations in the finest grid, will be considered equal to

which is the module of the difference with the highest number of extrapolations that can be calculated between the two finest grids. In the case of the x and y coordinates of ψ _min, u_min, v_min and v_max, one adopts

where

_1024x1024 and _512x512 are the numerical solutions obtained without extrapolation on 1024 x 1024 and 512 x 512 nodes grids.

Literature (Roache, 1998) offers several discretization error estimators. The use of Eqs. (13) and (14) is justified based on an analysis of the problem presented in the next section. This problem is similar to the classical problem of square cavity with movable lid, but its analytical solution is known. Thus, it was possible to evaluate the performance of Eqs. (13) and (14), which resulted in reliable error estimates, i.e., U/|E| > 1 for all variables of interest of the present work.

Problem with Known Analytical Solution

There is a variant of the classical problem of which analytical solution is known and is given by Shih, Tan and Hwang (1989). In this case, the source term (S) of Eq. (3) is different from zero, and is presented in Shih, Tan and Hwang (1989). The analytical solution of u and v is (Shih, Tan and Hwang, 1989)

The lid velocity varies with x, according to Eq. (15) for y = 1. The other boundary conditions are shown in Fig. 1. The Reynolds number (Re) is one for Eq. (17).

All numerical solutions in this work were obtained with ten different grids: 2 x 2, 4 x 4, 8 x 8 and so on up to 1024 x 1024 real control volumes. All simulations of this study were made in a Xeon Quad Core X5355 Intel processor, 2.66 GHz, using one core. The maximum RAM memory employed was 242 MB, for the simulations with 1024 x 1024 nodes grid. To obtain the numerical solutions, the initial estimate used for u, v and p was the analytical one for the problem given, respectively, by Eqs. (15) and (16) and by Shih, Tan and Hwang (1989).

The column named as Shih, in Table 3, presents the main parameters of the iterative process involved in the problem solution for the 1024 x 1024 nodes grid. At this table: Δt is the time interval used to further the iterative process; Itmax is the total number of outer iterations made; It(Eπ) is the approximate number of outer iterations made when the round-off error level is achieved; R1 is the l₁-norm of the residue sum of the three solved systems after the first iteration; Rf/R1 is the l₁-norm of the residue sum of the three solved systems at Itmax, made dimensionless based on the first iteration, showing how much the residue was reduced along the iterative process; Alg is the number of significant figures of the solution which does not have round-off errors after Itmax iterations; and tCPU is the CPU time needed to make Itmax iterations. The typical behavior of the residue drop along the iterations and the machine round-off error achievement are shown in Fig. 2 for the 256 x 256 nodes grid. In that figure, Shih denotes the problem from this section, and Ghia the problem treated in the next section.

Ideally, for the model used in this work, the value of the stream function (ψ) in y = 1 should be null for each of the 1024 control volumes at the cavity lid. Its absolute maximum value, which represents the mass conservation error, resulted in 1.4 x 10^-14; which is very near the null value and is at the level of double precision used in this work.

The velocity profiles in the two directions in the cavity center are shown in Fig. 3. The congruence between the analytical solution of Shih, Tan and Hwang (1989) and the numerical solution of the present work, with the 1024 x 1024 grid, can be considered excellent.

For each variable of interest, Tab. 4 presents: (a) the effective order (p_E) calculated by Eq. (10), based on the true discretization error (E) of the numerical solutions ( ) without extrapolation, on the 1024 x 1024 and 512 x 512 nodes grids, i.e., with r = 2; (b) the numerical solution () on the 1024 x 1024 nodes grid, with extrapolation calculated by Eq. (12), except for x and y variables for which none extrapolation was employed; (c) the value of E calculated by Eq. (9); and (d) the ratio of the estimated discretization error (U) to the module of E. In this table and in the following ones, the notation 1.0e-3 and Ind. represent, respectively, 1.0 x 10^-3 and indefinite.

Table 4

In Tab. 4, it can be noted that for the velocity profiles of u and v, p_E varies from 1.995 to 2.000, confirming the value of p_L = 2, which was predicted a priori. For coordinate type variables, p_E is indefinite or assumes values of null and close to unity or to two. For other variables, p_E varies from 1.639 to 2.584, i.e., around p_L.

For the velocity profiles u and v, M and Fs, the |Eh/E| ratio varies from 1.6 x 10³ to 3.8 x 10⁶. This ratio represents the extension to which the discretization error of the solution without extrapolation (Eh), obtained on the 1024 x 1024 grid, is reduced with the use of multiple extrapolations through Eq. (12). This reduction was not so effective for the variables v_min, v_max, u_min, Fn and ψ _min, which reductions were, respectively, of 1.9, 2.0, 2.6, 4.6 and 75.

The magnitude of E and U may vary considerably along the velocity profiles. The ratio between the maximum and the minimum values of E is 652, and for U it is 478.

Summarizing, in this section, for a two-dimensional flow problem with known analytical solution, we showed: (1) the importance of using multiple extrapolations to reduce the true discretization error (E); and (2) the discretization error estimated (U) with Eq. (13) or (14) reliability, i.e., U/|E| > 1. In the next section, the same procedure is applied to the classical cavity flow problem of which analytical solution is unknown.

Classical Problem with Unknown Analytical Solution

In the classical problem (Kawaguti, 1961; Burggraf, 1966; Rubin and Khosla, 1977; Benjamin and Denny, 1979; Ghia, Ghia and Shin, 1982) of laminar flow inside a square cavity, the lid velocity (U_T) is constant and has unitary value. The other boundary conditions are shown in Fig. 1. At lid corners, u = 0 on one side and u = 1 on the other. The source term (S) of Eq. (3) is null. The Reynolds number (Re) is defined by

where L = 1 m, dimension of the side of the square cavity; ρ = 1 kg/m3, density; and µ is the viscosity in Pa.s, obtained from Eq. (17) for a given Re. Numerical solutions were obtained for Re = 0.01, 10, 100, 400 and 1000. The initial estimate used was u = v = p = 0.

Table 3 shows the main parameters of the iterative process involved in the problem solution achievement for the 1024 x 1024 nodes grid. It is noted that the typical behavior of the residue drop along iterations and the achievement of the machine round-off error are shown in Fig. 2 for the 256 x 256 nodes grid. For these five values of Re and the 1024 x 1024 grid, CPU time varied from 2 days and 9 hours to 9 days and 1 hour. The l₁-norm of the residue sum of the three solved systems varied from 1.6 x 10^-15 to 9.5 x 10^-13. The obtained solutions present from 10 to 13 significant figures without machine round-off error.

The stream function value (ψ) in y = 1, which should be null for each of the 1024 control volumes at the cavity lid, resulted in the following maximum values: 5.9 x 10^-16, 1.7 x 10^-15, 5.4 x 10^-16, 1.0 x 10^-15 and 2.3 x 10^-15, respectively, for Re = 0.01, 10, 100, 400 and 1000. These values are very close to the null one, being at the level of double precision employed in this work.

The velocity profiles in the two directions at cavity center are shown in Fig. 4. The congruence between the numerical solutions of Ghia, Ghia and Shin (1982), Botella and Peyret (1998) and Bruneau and Saad (2006), and the numerical solution of this work using the 1024 x 1024 grid can be considered very good.

Table 5 presents the apparent order (p_U) for Re = 0.01, 10, 100, 400 and 1000. This order was calculated using Eq. (11), based on the numerical solutions () of each variable of interest, without extrapolation, on the 1024 x 1024, 512 x 512 and 256 x 256 nodes grids, i.e., with r = 2. For the u and v velocity profiles, p_U varies from 1.760 to 2.345, with most results very close to the value of p_L = 2, which was predicted a priori. For coordinate type variables, p_U is indefinite or assumes a unitary value. For other variables, p_U can differ substantially from p_L, varying from 1.102 to 2.191.

An unexpected behavior of p_U, shown in Table 5, occurred for the variable Fn. For all five values of Re, p_U varies between -3.50 x 10^-4 and 3.06 x 10^-2. Based on p_U versus h curves, it was found that p_U → 0 for h → 0. This means that the value of Fn diverges with grid refinement. In Nie, Robbins and Chen (2006), for the same problem studied in this work, it was shown that Fn cannot be obtained solely through continuum mechanics, i.e., with the Navier-Stokes equations. It is also necessary to consider the movement on microscopic scale. This problem is due to the discontinuity existing in the boundary condition of u, at lid corners: 0 on one side and 1 on the other. In Bruneau and Saad (2006), two other variables were found to diverge when h → 0, due to the discontinuity existing in the boundary condition of u. It should be noted that for Fs: (a) in all five values of Re, its value converges to h → 0; (b) it is calculated with the same types of numerical approximations as Fn; (c) its p_U varies from 1.322 to 2.029; and (d) its calculation does not involve discontinuities.

Table 6 presents, for Re = 0.01, 10, 100, 400 and 1000, the numerical solution () of each variable of interest. In the case of coordinate type variables, the solution is presented without extrapolations, obtained on the 1024 x 1024 nodes grid. For the remaining variables, the solution is the one obtained on the 1024 x 1024 nodes grid with extrapolations through Eq. (12). No results are shown for Fn due to its divergence. The number of significant figures presented for each variable is defined by its respective U value from Table 7.

Table 7

Tables 8 to 13 list the results of this work and those of several other authors for the variables of interest, where Ref. indicates the works cited in Table 1. Among all results of the sixteen works reported in literature and cited here, those of Botella and Peyret (1998) are probably the most accurate. However, considering the estimated error (U) reported by Botella and Peyret (1998) and the tolerance they adopted in the iterative process, the results of the present work are probably more accurate than those of Botella and Peyret (1998). Keeping in mind that, in the present work, U is presented in Tab. 7 and its value probably overestimates the true error; moreover, the iterative process was repeated until the achievement of machine round-off error. Among all variables of interest compared in Tables 8 to 13, the results of Botella and Peyret (1998) are probably more accurate than those of the present work only for the following variables: v_min for Re = 100 and 1000; and v_max for Re = 1000.

It is worth noting the congruence among all results of the present work, which are compared to those of Botella and Peyret (1998). The results of Botella and Peyret (1998) lie within the interval comprised between ± U of this work results. For example, the result of Botella and Peyret (1998) for v_max in Re = 1000 is 0.3769447, which is between 0.3769398 and 0.3769544, given in the present work. An exception is ψ _min, for which Botella and Peyret (1998) report the result of -0.1189366, which is not comprised the interval of -0.118936739 and -0.118936677 of this work, presenting a very slight difference of 7.7 x 10^-8.

Conclusion

In this work, numerical solutions were obtained for laminar flow inside a square cavity of which lid moves at variable velocity and analytical solution is known (Shih, Tan and Hwang, 1989). Results were presented for 42 variables of interest () in the 1024 x 1024 nodes grid. It was found that:

1. For all variables of interest, the discretization error estimated (U) with Eqs. (13) and (14), proposed here, is reliable. In other words, U/|E| > 1, where E is the true discretization error.

2. The use of multiple Richardson extrapolations (MRE) with Eq. (12) reduced E between 1.6 x 103 and 3.8 x 10⁶ times for velocity profiles u and v, M and Fs. This reduction was not so effective for variables v_min, v_max, u_min, Fn and ψ _min, which reductions were of 1.9, 2.0, 2.6, 4.6 and 75 times, respectively. For coordinate type variables, this procedure does not apply.

3. For 34 variables, the effective order value (p_E) is very close (1.96 to 2.08) to the theoretical asymptotic order (p_L) = 2 predicted a priori. For coordinate type variables, p_E seems to tend towards unity. For other variables, p_E varies from 1.64 to 2.58, i.e., around p_L.

The main focus of this work was to solve the problem of laminar flow inside a square cavity of which lid moves at a constant velocity and analytical solution is unknown (Kawaguti, 1961; Burggraf, 1966; Rubin and Khosla, 1977; Benjamin and Denny, 1979; Ghia, Ghia and Shin, 1982). Results were presented for 42 variables of interest (), and their estimated discretization errors (U) on a grid of 1024 x 1024 nodes and Reynolds numbers (Re) = 0.01, 10, 100, 400 and 1000. It was found that:

• Among all results of the sixteen works reported in literature and cited here, those of Botella and Peyret (1998) are probably the most accurate. However, considering the estimated error (U) reported by Botella and Peyret (1998) and the tolerance they adopted in the iterative process, the results of the present work are probably more accurate than those of Botella and Peyret (1998). Among all variables of interest compared in Tables 8 to 13, the results of Botella and Peyret (1998) are probably more accurate than those of the present work only for the following variables: v_min for Re = 100 and 1000; and v_max for Re = 1000. There is a notable consistency among all results of the present work, comparing them with those of Botella and Peyret (1998): the results of Botella and Peyret (1998) fall inside the interval comprised between ± U of the results of the present work.

• For velocity profiles u and v, the apparent order (p_U) varies from 1.76 to 2.34, with most of the results very close to the theoretical value of p_L = 2 which was predicted a priori. For coordinate type variables, p_U seems to tend towards unity. In only five cases out of more than 200, the value of p_U varied from 1.10 to 1.63, remaining distant from p_L. An exception is the variable Fn: its value does not converge with the grid refinement, causing p_U to tend towards zero. This is apparently due to the discontinuity in the boundary condition (B.C.) of u at lid corners. For Fs, which does not present discontinuities in the B.C., the solution converges with the grid refinement, with p_U varying from 1.33 to 2.03 for the five values of Re.

Acknowledgments

The second author gratefully acknowledges the support provided by the Laboratory of Numerical Experimentation (LENA), at Federal University of Paraná (UFPR), and the guidance of Professors Carlos Henrique Marchi and Fábio Alencar Schneider, as well as CAPES (Coordenação de Aperfeiçoamento de Pessoal de Nível Superior) - Brazil for granting a scholarship. And the first author thanks the support of CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico) - Brazil for granting him a scholarship. The authors acknowledge The "UNIESPAÇO Program" of The Brazilian Space Agency (AEB) and thank the two referees for their suggestions and comments.

Paper accepted December, 2008.

Technical Editor: Francisco R. Cunha

• Barragy, E. and Carey, G.F., 1997, "Stream function-vorticity driven cavity solution using p finite elements", Computers & Fluids, Vol. 26, pp. 453-468.
• Benjamin, A.S. and Denny, V.E., 1979, "On the convergence of numerical solutions for 2-D dlows in a cavity at large Re", Journal of Computational Physics, Vol. 33, pp. 340-358.
• Botella, O. and Peyret, R., 1998, "Benchmark spectral results on the lid-driven cavity flow", Computers & Fluids, Vol. 27, pp. 421-433.
• Bruneau, C.H. and Saad, M., 2006, "The 2D lid-driven cavity problem revisited", Computers & Fluids, Vol. 35, pp. 326-348.
• Burggraf, O.R., 1966, "Analytical and numerical studies of the structure of steady separated flows", Journal of Fluid Mechanics, Vol. 24, pp. 113-151.
• De Vahl Davis, G., 1983, "Natural convection of air in a square cavity: a bench mark numerical solution", International Journal for Numerical Methods in Fluids, Vol. 3, pp. 249-264.
• Erturk, E., Corke, T.C. and Gökçöl, C, 2005, "Numerical solutions of 2D steady incompressible driven cavity flow at high Reynolds numbers", International Journal for Numerical Methods in Fluids, Vol. 48, pp. 747-774.
• Ferziger, J.H. and Peric, M., 1999, "Computational Methods for Fluid Dynamics", Springer-Verlag, Berlin, Germany, 2 ed., 389 p.
• 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, Vol. 48, pp. 387-411.
• Goyon, O., 1996, "High-Reynolds number solutions of Navier-Stokes equations using incremental unknowns", Computer Methods in Applied Mechanics and Engineering, Vol. 130, pp. 319-335.
• Gupta, M.M. and Kalita, J.C., 2005, "A new paradigm for solving Navier-Stokes equations: streamfunction-velocity formulation", Journal of Computational Physics, Vol. 207, pp. 52-68.
• Hayase, T., Humphrey, J.A. and Greif, R., 1992, "A consistently formulated QUICK scheme for fast and stable convergence using finite-volume iterative calculation procedures", Journal of Computational Physics, Vol. 98, pp. 108-118.
• Hou, S., Zou, Q., Chen, S., Doolen, G. and Cogley, A.C., 1995, "Simulation of cavity flow by lattice Boltzmann method", Journal of Computational Physics, Vol.118, pp. 329-347.
• 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, Vol. 16, pp. 2307-2315.
• Khosla, P.K. and Rubin, S.G., 1974, "A diagonally dominant second-order accurate implicit scheme", Computers & Fluids, Vol. 2, pp. 207-209.
• Kreyszig, E., 1999, "Advanced Engineering Mathematics", Wiley, New York, United States of America, 8 ed., 1156 p.
• Marchi, C.H., 2001, "Verification of Unidimensional Numerical Solutions in Fluid Dynamics" (in Portuguese), Ph.D. thesis, Federal University of Santa Catarina, Florianópolis, Brazil, 289 p.
• Marchi, C.H. and Maliska, C.R., 1994, "A nonorthogonal finite-volume method for the solution of all speed flows using co-located variables", Numerical Heat transfer, Part B, Vol.26, pp. 293-311.
• Marchi, C. H. and Silva, A.F.C., 2002, "Unidimensional numerical solution error estimation for convergent apparent order", Numerical Heat Transfer, Part B, Vol. 42, pp. 167-188.
• Nie, X., Robbins, M.O. and Chen, S., 2006, "Resolving singular forces in cavity flow: multiscale modeling from atomic to millimeter scales", Physical Review Letters, Vol.96 (134501), pp. 1-4.
• Nishida, H. and Satofuka, N., 1992, "Higher-order solutions of square driven cavity flow using a variable-order multi-grid method", International Journal for Numerical Methods in Engineering, Vol. 34, pp. 637-653.
• Richardson, L. F., 1910, "The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam", Phylosophical Proceedings of the Royal Society of London Serial A, Vol. 210, pp. 307-357.
• Roache, P.J., 1994, "Perspective: a method for uniform reporting of grid refinement studies", ASME Journal of Fluids Engineering, Vol. 116, pp. 405-413.
• Roache, P.J., 1998, "Verification and Validation in Computational Science and Engineering", Hermosa, Albuquerque, United States of America, 446 p.
• Rubin, S.G. and Khosla, P.K., 1977, "Polynomial interpolation methods for viscous flow calculations", Journal of Computational Physics, Vol. 24, pp. 217-244.
• Schneider, F.A., 2007, "Verification of Numerical Solutions in Diffusive and Advective Problems with Nonuniform Grids" (in Portuguese), Ph.D. thesis, Federal University of Paraná, Curitiba, Brazil, 223 p.
• Schneider, G.E. and Zedan, M., 1981, "A modified strongly implicit procedure for the numerical solution of field problems", Numerical Heat Transfer, Vol. 4, pp. 1-19.
• Schreiber, R. and Keller, H.B., 1983, "Driven cavity flows by efficient numerical techniques", Journal of Computational Physics, Vol. 49, pp. 310-333.
• 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, Vol. 9, pp. 193-212.
• Van Doormaal, J.P. and Raithby, G.D., 1984, "Enhancements of the SIMPLE method for predicting incompressible fluid flow", Numerical Heat Transfer, Vol. 7, pp. 147-163.
• Vanka, S.P., 1986, "Block-implicit multigrid solution of Navier-Stokes equations in primitive variables", Journal of Computational Physics, Vol. 65, pp. 138-158.
• Wright, N.G. and Gaskell, P.H., 1995, "An efficient multigrid approach to solving highly recirculating flows", Computers & Fluids, Vol. 24, pp. 63-79.
• Zhang, J., 2003, "Numerical simulation of 2D square driven cavity using fourth-order compact finite difference schemes", Computers and Mathematics with Applications, Vol. 45, pp. 43-52.

Publication Dates

• Publication in this collection
  04 Dec 2009
• Date of issue
  Sept 2009