A NEW BENCHMARK QUALITY SOLUTION FOR THE BUOYANCY-DRIVEN CAVITY BY DISCRETE SINGULAR CONVOLUTION

Description

±Numerical Heat Transfer, Part B, 40: 199 228, 2001Copyright# 2001 Taylor & Francis1040-7790/01 $12.00+ .00A NEW BENCHMARK QUALITY SOLUTIONFOR THE BUOYANCY-DRIVEN CAVITYBY DISCRETESINGULARCONVOLUTIOND. C. Wan, B. S. V. Patnaik, and G. W. WeiDepartment of Computational Science, National University of Singapore,Republic of SingaporeThis article introduces a high-accuracy discrete singular convolution (DSC) for the nu-merical simulation of coupled convective heat transfer problems. The problem of a buoy-ancy-driven cavity is solved by two completely independent numerical procedures. One is aquasi-wavelet-based DSCapproach,whichusestheregularized Shannon skernel, while theother is a standard form of the Galerkin Žnite-element method. The integration of theNavier–Stokes and energy equations is performed by employing velocity correction-based3 8schemes. The entire laminar natural convection range of 10μ Raμ10 is numericallysimulated by both schemes. The reliability and robustness of the present DSC approach isextensively tested and validated by means of grid sensitivity and convergence studies. As aresult, a set of new benchmark quality data is presented. The study emphasizesquantitative,rather than qualitative comparisons.I. INTRODUCTIONConvection by natural means is crucial to ¯ows, in both nature and technol-ogy. There are a variety of real-world applications of natural convection, such asthermal insulation, cooling of electronic equipment, solar ...

Informations

• Publié par: Faut
• Nombre de lectures: 46
• Langue: English

Extrait

Numerical Heat Transfer, Part B, 40: 199± 228, 2001 Copyright # 2001 Taylor & Francis 1040-7790 /01 $12.00 + .00

A NEW BENCHM ARK QUALITY SOLUTION FOR THE BUOYANCY-DRIVEN CAVITY BY DISCRETE SINGULAR CONVOLUTION

D. C. Wan, B. S. V. Patnaik, and G. W. Wei Departmen t of Computational Science, National University of Singapore, Republic of Singapor e

This article introduces a high-accuracy discrete singular convolution (DSC) for the nu- merical simulation of coupled convective heat transfer problems. The problem of a buoy- ancy-driven cavity is solved by two completely independent numerical procedures . One is a quasi-wavelet-based DSC approach , which uses the regularized Shannon’s kernel, while the other is a standard form of the Galerkin  nite-element method. The integration of the Navier–Stokes and energy equations is performed by employing velocity correction-base d schemes. The entire laminar natural convection range of10 3 µRaµ10 8 is numerically simulated by both schemes. The reliability and robustness of the present DSC approac h is extensively tested and validated by means of grid sensitivity and convergence studies. As a result, a set of new benchmark quality data is presented. The study emphasizes quantitative, rather than qualitative comparisons .

I. INTRODUCTION Convection by natural means is crucial to ¯ ows, in both nature and technol- ogy. There are a variety of real-world applications of natural convection, such as thermal insulation, cooling of electronic equipment, solar energy devices, nuclear reactors, heat-recovery systems, room ventilation, crystal growth in liquids, etc. The ¯ uid ¯ ow and heat transfer behavior of such systems can be predicted by the mass, momentum, and energy conservation equations with appropriate boundary condi- tions. The fast-emerging branch of computationa l ¯ uid dynamics (CFD) facilitates the numerical simulation of ¯ uid ¯ ow and heat transfer features. A comprehensive analysis of the ¯ uid ¯ ow and heat transfer patterns in fundamentally simple geo- metries, such as the buoyancy-drive n square cavity, is a necessary precursor to the evolution of better designs for more complex industrial applications. Jones[1]has proposed the problem of a buoyancy-drive cavity as a suitable n vehicle for testing and validation of computer codes for thermal problems. The simplicity of the geometry and the clarity in boundary conditions render this problem

Received 8 December 2000; accepted 1 March 2001. This work was supported by the National University of Singapore and the National Science and Technology Board of the Republic of Singapore. Address correspondence to Dr. G. W. Wei, Department of Computational Science, National University of Singapore, Block S 17, Level 7, 3 Science Drive, Singapore 117543. E-mail: cscweigw@ nus.edu.sg

199

200

D. C. WAN ET AL.

NOM ENCLATURE Bfrequency boundydimensional vertical coordinate C Courant numberathermal di usivity of the ¯ uid Eerrora 1 ;a 2 ;a 3 g Kutta used in the Runge± constantsacceleration due to gravity Gr Grashof number (ˆgbDyL r 3 ef =n 2 ) scheme[ (34) , ( , 35)Eqs. (33)] L ref  width eitherreference length dimension (b 1 ;b 2 ;b 3 or height) constants used in the Runge± Kutta L 2  schemeerror norm[ 35) ( ,Eqs. (33) , (34)] Nu Nusselt numberDgrid spacing, incremental value Nuaverage Nusselt numberd D …x†;d D;s …x† pnondimensional pressure (ˆpL r 2 ef =ra 2 ) convolution kernels pdimensional pressureEa small parameter Pe c cell Peclet number (ˆuDx=n)Zvariable of interest Pr Prandtl number (ˆn=a)ynondimensional temperature rratio ofsoverD[ˆ …y¡y C †=…y H ¡y C †] Ra Rayleigh number (ˆgbDyL 3 ref =na)nkinematic viscosity R s …x†delta regularizerrdensity tdimensional timeswidth of the Gaussian envelope tnondimensional time (ˆta=L r 2 ef )fpotential function unondimensional horizontal velocityorelaxation parameter (ˆuL ref =a)ydimensional temperature udimensional horizontal velocityHgradient Unondimensional velocity vector (u;v) vnondimensional vertical velocity Subscript s (ˆvL ref =a) Ccold wall vdimensional vertical velocity Hhot wall Wbandwidth of support on one side i;jindices in the horizontal and vertical of the grid point  i x  respect vely dimensional horizontal coordinate  direction, xnondimensional horizontal coordinate (ˆx=L ref ) Superscript s x k discrete sampling points around the ¤ value eldintermediate ® point x l;niteration labels ynondimensional vertical coordinate qorder of the derivative (ˆy=L ref )

more appealing for new computationa l algorithms. In the early 1980s, this problem was solved by a number of di erent groups and their results were extensively sum- marized in the standard reference by de Vahl Davis and Jones[2]. De Vahl Davis[3] used forward di erence for the temporal discretization and second-order central di erence for the spatial discretization, to solve the stream function-vorticit y form of the equations. The resulting algebraic equations were solved by an alternating- direction implicit (ADI) algorithm. Although the solution obtained was among the best of those days, the data presented were limited to a Rayleigh number (Ra) of 10 6 . It is pertinent to quote Professor de Vahl Davis from his seminal article[3]: It is hoped that it will lead to further contributions to the search for e cient and high accuracy methods for problems of this type. Even af ter nearly two decades, this statement remains sacrosanct, despite the phe- nomenal progress in numerical methods. A variety of computationa l algorithms has

BENCHMARK SOLUTION FOR THE BUOYANCY-DRIVEN CAVITY201

been tested on this problem and the vast amount of literature is a testimony to this [1± 10]. For example, Hortmann et al.[4] dhave employed a ® nite-volume-base multigrid technique for the simulation of a buoyancy-drive n cavity. Ramaswamy et al.[8]have investigated the performance of two explicit and one semi-implicit projection-base d schemes. They have concluded that the semi-implicit scheme always outperforms the two explicit schemes considered by them. In their numerical in- vestigations, a variety of benchmark problems were solved, including the buoyancy- driven cavity. Shu and Xue[9]employed a global method of generalized di erential quadratur e (GDQ) for solving the stream function-vorticit y form of the Navier- Stokes (N-S) equations. Massarotti et al.[5]used a semi-implicit form of the char- acteristic-base d split scheme (CBS) with equal-order interpolatio n functions for all the variables. Manzari et al.[7] di usion-based algorithm for cialdeveloped an arti® 3-D compressible turbulent ¯ ow problems. This was later extended to 2-D laminar heat transfer problems. The basic idea of Manzari[6, 7]involves modi® of cation continuity equation by employing the concept of arti® cial compressibility. Recently, Mayne et al.[10]employed anh-adaptive ® nite-element method to ensure a very accurate solution for the thermal cavity problem. Despite so much e ort on this problem, there still exist some variations and discrepancies in the available literature (see the detailed discussion in Section III). Inherently, numerical results are ap- proximations . Their accuracy and reliability depends vitally on the underlying computationa l method and the numerical scheme. Hence, further advances in computationa l methodology are crucial to the thermal cavity problem as well as other heat transfer problems. A variety of computationa l methods are available in the literature. In a broad sense they can be classi® ed into global and local methods. Global methods approx- imate a di erentiation at a point by all grid points in the computationa l domain and can be highly accurate. Spectral methods, pseudo-spectra l methods, fast Fourier transforms, and di erential quadratur e come under this category. For example, spectral methods converge exponentially with mesh re® nement for approximatin g an analytical function and thus have the potential to be used in high-precisio n and de- manding large-scale computations . Many global approaches have been successfully applied to the study of ¯ uid ¯ ow and heat transfer in simple geometries, such as cavities, channels, ducts, di users, etc. However, global methods have limited cap- ability in handling irregular geometries and more complex boundary conditions. In- deed, local methods, such as ® nite di erences, ® nite-elements, ® nite strips, and ® nite volumes, are the most popular approache s for solving engineering problems. Local methods utilize information from the nearest neighboring grid points to approximat e the di erentiation at a point and thus are much more ¯ exible. However, local methods converge slowly with respect to mesh re® nement and are not cost-e ective for achieving high precision. Hence, there is a strong demand for a scheme which can exploit the advantage s of both methods. To this end, a high-accuracy quasi-wavelet - based approach, the discrete singular convolution (DSC), was proposed[13]. This method is a promising approach for the numerical realization of singular convolutions [14, 15]. Mathematical foundation for this algorithm stems from the theory of dis- tributions (or generalized functions). Sequence of approximation s to the singular kernels of Hilbert type, Abel type, and delta type can be constructed. Numerical solution to di erential equations are formulated via singular kernels of delta type.

202

D. C. WAN ET AL.

The DSC approach exhibits global methods’ accuracy for integration and local methods’ ¯ exibility in handling irregular geometries and complex boundary conditions, when suitable DSC kernels are chosen. Many DSC kernels, such as (regularized) Shannon’ s delta kernel, (regularized) Dirichlet kernel, (regularized) Lagrange kernel, and (regularized) de la VallÂee Poussin kernel, are constructe  d for a number of applications, such as numerical solution of the Fokker± Planck equation[14, 15]and the Schr Èodinger equation[16]. The DSC algorithm was also utilized for waveguide model analysis, electromagnetic wave propagatio n[17], and structural (plate and beam) analysis[18± 20]with excellent results. Most recently, the DSC algorithm was used to resolve a few numerically challenging problems. The integration of the (nonlinear) sine-Gordon equation[21]with the initial values close to a homoclinic orbit singularity is one such problem for which conventional local methods have encountered great di culty and numerically induced chaos was reported[22]. Another complex problem that has been resolved by using the DSC algorithm is the integration of the (nonlinear) Cahn± Hilliard equation in a circular domain[23] singularity arti® cial r, which is challenging because of the fourth-orde at the origin and complex phase-space geometry. DSC solution ofmachine preci- sionto the Navier± Stokes equations with periodic boundary conditions for the Taylor problem was obtained with 33 grid points in each dimension[15, 18]. Recently, a DSC ® nite-subdomai n method was proposed for the solution of in- compressible viscous ¯ ows under complex geometries[24]. The objectives of the present study are the following: (1) To introduce the highly accurate quasi-wavelet-base d discrete singular convolution[15, 18]for the numerical simulation of coupled convective heat transfer problems; (2) to present benchmark-qualit y data for the entire laminar natural-convectio n range of 10 3 µRaµ10 8 ; and (3) to present a focused and elaborate study on the problem of the buoyancy-drive n cavity, which also looks into some of the discrepancies observed in the literature. The high level of accuracy that could be achieved for the Taylor problem[15, 18]enhances the level of con® and reliability for the dence simulation of the driven-cavity problem. The DSC code is extensively tested and further validated with a second-order-accurat e ® nite-element-base d Galerkin method, against the available numerical simulations of[3± 10]. This article is organized into four sections. A description of the problem under investigation and the methods of solution by both DSC and FEM are presented in Section II. An elaborate and focused study on the results obtained for the buoyancy-drive n cavity is presented in Section III. Conclusions are presented in Section IV.

II. THEORETICAL BACKGROUND AND M ETHODS OF SOLUTION This section de® nes the buoyancy-drive n cavity problem, including the gov- erning partial di erential equations (PDEs) and the boundary conditions. Two completely independent methods of solution, viz., a quasi-wavelet-base d approach and a ® nite-element method, are employed to solve the problem under investigation. The theoretical framework behind the discrete singular convolution is elaborated. The method of solution for the integration of the Navier-Stokes and energy equa- tions by the DSC and FEM is explained.

BENCHMARK SOLUTION FOR THE BUOYANCY-DRIVEN CAVITY203

A. Problem Descr ipt ion A di erentially heated, closed square cavity is depicted in Figure 1, whose left and right vertical walls are maintained aty H andy C , respectively. The horizontal walls are adiabatic (insulated, and there is no transfer of heat through these walls). Fluid is assumed to be viscous, incompressible, Newtonian, and Boussinesq- approximated . The Newtonian assumption guarantees a linear relationship between the shear stress and the velocity gradient. The Boussinesq approximatio n means that the density di erences are con® ned to the buoyanc y term, without violating the assumption of incompressibility. It should be pointed out that there is an additional coupling term in the momentum equations, which indeed dictates the ¯ uid motion within the cavity. 1. Governin g equat ions.The governing PDEs are the coupled mass, mo- mentum, and energy conservation equations, applicable for two dimensions. The equations are given as

Continuity:

q q u x ‡ q qv y ˆ0

Figure 1.Flow domain of interest.

…1†

204

D. C. WAN ET AL. quququqp P ³ q q 2 u x 2 ‡ q q 2 u y 2 ´ xmomentum: qt ‡u qx ‡v qy ˆ ¡ qx ‡r ymomentum: q q v t ‡u q q v ‡ qvqp ³ q 2 vq 2 v ´ ‡Ra Pry x v qy ˆ ¡ qy ‡Pr qx 2 ‡ qy 2

qy qy qy q 2 y q 2 y Energy: qt ‡u qx ‡v qy ˆ qx 2 ‡ qy 2

In the equations above, the following nondimensionalizatio n was employed:

x y xˆ L ef yˆ L ref uˆ uL a ref vˆ vL a ref tˆ L t r 2 a ef r y¡y C pˆ p r L a r 2 e 2 f y y H ¡y C Raˆ gDy na bL r 3 ef Prˆ n a ˆ

…2† …3†

…4†

…5† …6†

In convection, density di erences generate an additional force (popularly known as the buoyanc y force), which competes with the inertial and viscous forces. The ratio of buoyancy to viscous forces is given by the paramete r Grashof number (Gr), which controls natural convection. On the other hand, the ratio of momentum to thermal di usivity, known as Prandtl number (Pr), governs the temperature ® eld and its relationship with the ¯ uid ¯ ow characteristics. Rayleigh number (Ra) , which is the paramete r of interest, is the product of these two dimensionless groups. 2. Boundary condit ions.No-slip velocity boundary condition (uˆvˆ0:0) is applied on all four walls of the square cavity. For temperature , Dirichlet boundary conditions ofy H ˆ1:0 andy C ˆ0:0 are enforced on the left and right vertical walls, respectively. As there is no transfer of heat through the horizontal walls, a Neumann boundary condition (qy=qyˆ0: is applied.0) , 3. The lam inar-flow region.According to Incropera and Dewitt[11], the la- minar natural convection at a local Rayleigh number larger than 10 9 may be pro- moted to turbulent transition in the vertical boundary layer. Undoubtedly, the transition from laminar to turbulence causes an increase in convective heat transfer on the surface of the wall. In their direct numerical simulation (DNS) of air flow in a square cavity, Paolucci and Chenoweth[12]detected the existence of a critical Ray- leigh number between 10 8 and 2£10 8 , where the flow undergoes a Hopf bifurcation into a periodic unsteady flow. Therefore, we confine the present simulations to the lower limit of this Rayleigh number, which is 10 8 . Beyond this value, a Reynolds averaged form of the Navier± Stokes and energy equations would have to be employed together with a workable turbulence model.

B. Discret e Singular Convolu t ion 1. Approxim at ion of singular convo lut ion.For the sake of clarity and integrity in presentation, the section begins with a brief description of the discrete

BENCHMARK SOLUTION FOR THE BUOYANCY-DRIVEN CAVITY205

singular convolution (DSC) algorithm. Elaborate details are available in previous studies[14, 15]. The DSC algorithm concerns with the computer realization of the mathematical distributions. Distributions are not well defined in the usual sense and may not have any value. Particular examples are kernels of the Hilbert type and Abel type. These singular kernels are of crucial importance to a number of fields, such as Radon transform, analytical function theory, linear response theory, etc. For data (surface) interpolatio n and solving partial differential equations, singular ker- nels of delta type are useful. In the DSC algorithm, the functionf…x†and its deriva- tives with respect to the coordinat e at a grid pointxare approximate d by a linear sum of discrete values{f…x k †}in the narrow bandwidth[x¡x W ;x‡x W ]. This can be expressed as

W f … q † …x† º X d … D q ; † s …x¡x k †f…x k † k ˆ¡ W

…7†

where superscriptq…qˆ0;1;2; . . . †denotes theqth-order derivative with respect tox. The{x k }refers to a set of discrete sampling points centered around the pointx. Heresis a regularization parameter,Dis the grid spacing, and 2W‡1 is the total computationa l bandwidth, which is usually much smaller than the computationa l domain. In Eq. (7) ,d D;s …x†is a convolution kernel that approximate the delta dis- s tribution. For band-limited functions, the delta distribution can be replaced by a low-pass ® lter, hence, many wavelet scaling functions can be used as DSC con- volution kernels. One interesting example is Shannon’ s wavelet scaling function,

d D …x† sin…px=D D † ˆ px=

…8†

In fact, Shannon’ s wavelet scaling function forms a sampling basis for the Paley- Wiener reproducing kernel Hilbert space. Shannon’ s wavelet scaling is useful in solving eigenvalue problems with smoothly con® ned potentials[14]and for purely periodic boundary conditions[21]. However, for general computationa l problems, it is important to regularize Shannon’ s delta kernel,

d D;s …x† ˆd D …x†R s …x†

…9†

whereR s …x†is a delta regularizer[13, 14], which can dramatically increase the regularity of Shannon’ s wavelet scaling function. An often-used delta regularizer is the Gaussian

R s …x† ˆexp ³ ¡ 2 x s 2 2 ´ s>0

…10†

wheresthe width of the Gaussian envelop and can be varied in asso-determines ciation with the grid spacing, i.e.,sˆrDbe pointed out that the use. It should of appropriate regularizers can extend the domain of applicability of delta kernels to temporal distributions and even to exponentiall y growing functions[13].

206

D. C. WAN ET AL.

The expression in Eq. (7) provides extremely high computationa l e ciency both on and o a grid. In fact, it can provide exact results when the sampling points are extended to a set of in® nite points for certain band-limitedL 2 functions. Qian and Wei[25]have given a mathematical estimation for the choice ofW,r, andD. For example, if theL 2  anerror for approximatin gL 2 functionf…x†is set to 10 ¡Z , the following relations would give the best accuracy:  r…p¡BD†> p 4:61Zand W r > p 4:61Z  …11† whererˆs=DandBis the frequency bound for the function of interest,f…x†. The ® rst inequality states that for a given grid sizeD, a largeris required for approx- imating the high-frequency component of anL 2 function. The second inequality indicates that if one chooses the ratiorˆ3, then the half-bandwidthWº30 is good enough to achieve the highest accuracy in double-precision computations (Zˆ15) . By appropriately choosingW,r, andD, the resulting approximatio matrix for in- n terpolation and solving the di erential equations has a banded structure. This en- sures that the DSC algorithm attains optimal accuracy and e ciency. It should be pointed out that although regularized Shannon’ s delta kernel is used to illustrate the DSC approximatio n of the delta distribution here, there are a variety of other DSC kernels, such as (regularized) Dirichlet kernels, (regularized) Lagrange kernels, and (regularized) de la VallÂee kernels, which perform equally well [14, 17, 21]grid used in Eq. (7) is uniform because only a single grid spacing is. The prescribed. In computations , Eq. (7) is very e cient since just one kernel is required for the whole computationa l domain[a;b]for a givenDandr. Thus, the kernel ac- quires the property of translationa l invariance. In order to maintain this property near a computationa l boundary, the functionsf…x k †to be located outside the com-have putational domain[a;b] ned. Therefore,, where their values are usually unde® it is necessary to create ® ctitious domains outside the computationa l boundaries. For the DSC algorithm, function values in these ® ctitious domains are generated ac- cording to the boundary conditions and the physical behavior of the solution at the boundaries. For example, when Dirichlet boundary conditions are employed,f…x k †in the ® ctitious domain can be taken to be eitherf…a†orf…b†; for a periodic boundary condition,f…x k †may be enforced by a periodic extension inside the computationa l domain[a;b]to outside. Neumann boundary condition forf…x k †may be obtained by f ’ …a†orf ’ …b†. When the regularized Shannon’ s delta kernel is used, t ed ex res ford … D 0 ; † s …x†,d … D 1 ; † s …x†, andd … D 2 ; † s …x†can be given analyticall  he detail  y as  p sions d … D 0 ; † s …x† ˆ ( sin…p x =D†exp… D ¡ x 2 = 2 s 2 † …x6ˆ0 … † 12† p x = 1…xˆ0†, 8 cos…p x =D†exp…¡ x 2 = 2 s 2 † ¡ sin…p x =D† p e x x 2 p = …¡ x 2 = 2 s 2 † x D d … D 1 ; † s …x† ˆ < > ¡ sin…p x =D† p e s x 2 p = … D ¡ x 2 = 2 s 2 † > : 0

…x6ˆ0† …13† …xˆ0†

BENCHMARK SOLUTION FOR THE BUOYANCY-DRIVEN CAVITY207

and

8 ¡…p=D†sin…p x = x D†exp…¡ x 2 = 2 s 2 † ¡2 cos x 2 …p x =D†exp…¡ x 2 = 2 s 2 † x =D†exp…¡ x 2 = 2 s 2 † > ¡2 cos…p x =D†e s x 2 p…¡ x 2 = 2 s 2 † ‡2 sin… p p x 3 =D d … D 2 ; † s …x† ˆ < ‡ sin…p x = p D† s e 2 x x p = …¡ x 2 = 2 s 2 † ‡ x sin…p x = p D† s e 4 x = p D …¡ x 2 = 2 s 2 † …x6ˆ0 † …14† D

> : ¡ 3 ‡p 3 2 s s 2 2 =D 2 …xˆ0† Once the parameterr cientsis chosen, the coed … D 0 ; † s …x†,d … D 1 ; † s …x†, andd … D 2 ; † s …x† depend only on the grid spacingD. Therefore, when the grid spacing is prescribed, the coe cients need to be computed only once and can be used during the whole computation. The good performance of the present DSC algorithm is due to the unique use of the DSC algorithm both for data interpolation and spatial dis- cretization of the governing equations. Since the computationa l bandwidth is user- de® ned, the approximatio n accuracy is controllable in the present algorithm[15].

2. Solut ion m et hodology.For the convenience of presenting the method of solution, we define the following:

D… † ˆ quqv U q ‡ qy x L…U† ˆF…U† ¡Hp T F…U† ˆ[f;g] T Uˆ[u;v] T Hpˆ µ q q p x ; q q p y ¶ M…y† ˆ[m] f ³ q 2 uq 2 y 2 u ´ ¡ ³ u q q x u ‡v q q u ´ ˆPr qx 2 ‡ qy ³ q 2 zvq 2 v ´ ³ q q x v ‡v q q v y ´ ‡Ra Pry gˆPr qx 2 ‡ qy 2 ¡u mˆ ³ q q 2 x y 2 ‡ q q 2 y y 2 ´ ¡ ³ u q q y x ‡v q q y y ´ Therefore, the system of Eqs. (1)± (4) can be simpli® ed as follows:

…15† …16† …17† …18† …19† …20†

D…U† ˆ0…21† q qt U ˆL…U† ˆF…U† ¡Hp…22† qy ˆM…y† …23† qt A variety of fractional step approaches can be formulated by appropriately combining convective, viscous, and pressure terms of the momentum Eqs. (21)± (23).

208

D. C. WAN ET AL.

When primitive variables are used in the governing PDEs, there is no direct link between the continuity and momentum equations. To bridge this gap, certain rearrangemen t of the momentum equations should be carried out, and thus the popular Poisson equation for pressure is formulated. In the simulation of Navier± Stokes equations, ® rst an intermediate velocity is computed by omitting the pressure gradient, and later, it is corrected by including the same. This seminal idea was originally introduced by Chorin[26]in a ® nite-di erence context. Since then, many special schemes have been designed and developed in this direction. For example, there is a family of implicit and semi-implicit pressure-correctio n methods such as the SIMPLE, consistent SIMPLE (SIMPLEC), and SIMPLER. Moreover, the arti® cial compressibility method, the marker and cell (MAC) method, the fractional-step projection method[26]and its many variants are also commonly used in the literature. In the present investigation, we adopt a fractional time-step and potential- function method (FTSPFM), which is a variant of the MAC method for solving the governing Eqs. (21)± (23). In this approach, an intermediate velocity ® eld and a potential function are introduced and computed to update the velocity, pressure, and temperature in the domain. 3. Spat ial discret izat ion.In the present investigation a staggered grid system is employed. The momentum and energy equations in the horizontal direction is written at the point…i‡ 1 2 ;j†, the momentum equation in the vertical direction is written at the point…i;j‡ 1 2 †, and the pressure and temperature are given at point …i;j†. The continuity equation is approximate d at the point…i;j†. All spatial deriva- tives in Eqs. (21)± (23) are discretized by using the DSC approach. A uniform grid in bothxandydirections is employed. The discretized forms of Eqs. (15)± (20) can be expressed as follows:

…24† …25†

W W D h …U† ˆ X d … D 1 ; † s …kDx†u i ‡ k ; j ‡ X d … D 1 ; † s …kDy†v i ; j ‡ k k ˆ¡ W k ˆ¡ W "  k W W # T H h pˆ X d … D 1 ; † s …kDx†p i ‡ k ; j ; X d … D 1 ; † s …kDy†p i ; j ‡ k ˆ¡ W k ˆ¡ W W W f h ˆPr " X d … D 2 ; † s …kDx†u i ‡ 1 2 ‡ k ; j ‡ X d … D 2 ; † s …kDy†u i ‡ 1 2 ; j ‡ k # k ˆ¡ W k ˆ¡ W W W ¡ " u i ‡ 2 1 ; j X d … D 1 ; † s …kDx†u i ‡ 1 2 ‡ k ; j ‡v i ‡ 1 2 ; j X d … D 1 ; † s …kDy†u i ‡ 1 2 ; j ‡ k # …26† k ˆ¡ W k ˆ¡ W "  j W W ‡ 1 2 ‡ k # g h ˆPr X d … D 2 ; † s …kDx†v i ‡ k ; j ‡ 1 2 ‡ X d … D 2 ; † s …kDy†v i ; k ˆ¡ W k ˆ¡ W W W " k ˆ¡ W k ˆ¡ W ‡ k # ¡u i ; j ‡ 1 2 X d … D 1 ; † s …kDx†v i ‡ k ; j ‡ 1 2 ‡v i ; j ‡ 2 1 X d … D 1 ; † s …kDy†v i ; j ‡ 1 2 ‡Ra Pry i ; j …27†

BENCHMARK SOLUTION FOR THE BUOYANCY-DRIVEN CAVITY m h ˆ " k ˆ X W ¡ W d … D 2 ; † s …kDx†y i ‡ k ; j ‡ k ˆ W X ¡ W d … D 2 ; † s …kDy†y i ; j ‡ k # ¡ " u i ; j k ˆ W X ¡ W d … D 1 ; † s …kDx†y i ‡ k ; j ‡v i ; j k ˆ W X ¡ W d … D 1 ; † s …kDy†y i ; j ‡ k #

209 …28†

L h …U† ˆF h …U† ¡H h p F h …U† ˆ[f h ;g h ] T M h …y† ˆ[m h ]…29† whered … D 1 † s andd … D 2 ; † s  of the regularized Shannon’ sare coe cients delta kernel, given ; in Eqs. (13) and (14), respectively. Here, the detailed labels (i,j) onUare omitted, andDxandDydenote the grid sizes in thexandydirection, respectively. By substituting Eqs. (24)± (28) into Eqs. (21) ± (23), the following semidiscretized approximation s are obtained:

D h …U† ˆ0 d d U t ˆL h …U† ˆF h …U† ¡H h p d d y t ˆM h …y†

…30† …31† …32†

4. Tem poral discret izat ion.A Runge± Kutta scheme is used for temporal dis- cretization. The scheme is of third-orde r accuracy in time and was used by many other authors[34, 35] are formu- (32). In this scheme, ordinary differential Eqs. (31)± lated, which can be solved from the following: U … 1 † ˆa 1 U n ‡b 1 {Dt[F h …U n † ¡H h p … 1 † ]}…33† U … 2 † ˆa 2 U n ‡b 2 {U … 1 † ‡Dt[F h …U … 1 † † ¡H h p … 2 † ]}…34† U n ‡ 1 ˆa 3 U n ‡b 3 {U … 2 † ‡Dt[F h …U … 2 † † ¡H h p n ‡ 1 ]}…35† y n ‡ 1 ˆy n ‡{Dt[M h …y n †]}…36†

where…a 1 ;a 2 ;a 3 † ˆ …1; 43 ; 1 3 †and…b 1 ;b 2 ;b 3 † ˆ …1; 1 2 †. TheU … 1 † ,p … 1 † andU … 2 † ,p … 2 † , 4 ; 3 are their corresponding ® rst and second step values for velocity and pressure, respectively. 5. Treat m ent for t he pressure.Updating the pressure field requires special care. Therefore, a brief description of the treatment of pressure is given. At each step of the Runge± Kutta scheme, the FTSPFM[24]is adopted to solve Eqs. (30) and (31). To illustrate the present scheme, we consider the first step of the Runge± Kutta scheme. Assume that, at timet n , the velocityU n and pressurep n are known, while U … 1 † andp … 1 †  Kutta scheme are unknown. Let us intro-at the first step of the Runge± duce a first-step intermediate velocity fieldU ¤… 1 † :