The set of coupled partial differential equations and the boundary conditions described in the previous section are converted to algebraic equations by means of finite-volume techniques using rectangular meshes on a staggered arrangement. Diffusive terms at the boundaries of the control volumes are evaluated by means of second-order central differences, while the convective terms are approximated by means of the high-order SMART scheme [13], with a theoretical order of accuracy between 1 and 3.

The domain where the computations are performed and a schematic of the mesh adopted is shown in Fig. 1b. The mesh is represented by the parameter n. According to Fig.1b, 1.5n*n* n control volumes are used (for example, when the problem is solved on 30 * 20 * 20 control volumes, it means n = 20). The numerical solutions have been calculated adopting a global /г-refinement criterion. That is, all the numerical parameters (numerical scheme, numerical boundary conditions, etc.) are fixed, and the mesh is refined to yield a set of numerical solutions. This set of numerical solutions have been post-processed by means of a tool based on the Richardson extrapolation theory and on the concept of the Grid Convergence Index (GCI) [4][5]. When the numerical solutions are free of programming errors, convergence errors and round-off errors, the computational error is only due to the discretization. The tool processes a set of three consecutive solutions in the /і-refinement. The main output is an estimate of the uncertainty of the numerical solutions due to discretization, the. Only solutions with order of accuracy between 1 and 3 GCI less than 1percent and Richardson nodes higher than 60 percent are considered. The mesh is intensified near the two plates using a tanh-like function with a concentration factor of 2 [10], in order to properly solve the boundary layer. This aspect is indicated in the Fig. 1b with solid triangles at the boundaries intensified.

The resulting algebraic equation system was solved using SIMPLE algorithm [2], and the iterative convergence procedure was truncated when relative increments of the computed Nusselt number in the last 50 iterations dropped below 0.0001%.