• 沒有找到結果。

GDQ Method for Natural Convection in a Cubic Cavity Using Velocity-vorticity Formulation

N/A
N/A
Protected

Academic year: 2021

Share "GDQ Method for Natural Convection in a Cubic Cavity Using Velocity-vorticity Formulation"

Copied!
25
0
0

加載中.... (立即查看全文)

全文

(1)

GDQ METHOD FOR NATURAL CONVECTION IN A CUBIC

CAVITY USING VELOCITY–VORTICITY FORMULATION

D. C. Lo, D. L. Young, and K. Murugesan

Department of Civil Engineering and Hydrotech Research Institute, National Taiwan University, Taipei, Taiwan

Natural convection in a differentially heated cubic enclosure is studied by solving the velocity–vorticity form of the Navier–Stokes equations by a generalized differential quad-rature (GDQ) method. The governing equations in the form of velocity Poisson equations, vorticity transport equations, and energy equation are solved using a coupled numerical scheme via a single global matrix for velocities, vorticities, and temperature. Vorticity and velocity coupling at the solid boundaries is enforced through a higher-order approxi-mation by the GDQ method, thus assuring accurate satisfaction of the continuity equation. Nusselt numbers computed for Ra¼ 103

, 104, 105, and 106show good agreement with the benchmark results. A mesh independence study indicates that the present numerical procedure requires much coarse mesh compared to other numerical schemes to produce the benchmark solutions of the flow and heat transfer problems.

1. INTRODUCTION

Vortex-dominant flow is an important characteristic of natural convection in differentially heated cavities. In situations where the main concern is heat transfer rather than the pressure field, the velocity–vorticity form of the momentum equa-tions is a better choice compared to the primitive-variable form to achieve the divergence-free velocity field constraint for the incompressible Navier–Stokes equa-tions in three dimensions. In addition, the vortex flow-dominated natural-convection flow field can be simulated directly by the velocity–vorticity formulation, without the need to handle the pressure term. For the case of incompressible fluid flows, if the incompressibility condition imposed by the continuity equation is satisfied by some means, then a divergence-free flow field can be computed by solving the velocity– vorticity equations [1–5]. Mallinson and de Vahl Davis [6] were the first to propose the vorticity–stream function formulation without the pressure term, to study natu-ral convection in a rectangular box. Though this formulation is easy to implement for the case of two-dimensional problems, an extension to three-dimensional pro-blems is not straightforward. Further, the velocity, which is required for computing

Received 24 August 2004; accepted 5 February 2005.

We acknowledge financial support from the National Science Council of Taiwan through Grant NSC-93-E-002-017, and it is gratefully appreciated.

Address correspondence to D. L. Young, Department of Civil Engineering and Hydrotech Research Institute, National Taiwan University, Taipei, Taiwan. E-mail: dlyoung@ntu.edu.tw

363 Copyright # Taylor & Francis Inc.

ISSN: 1040-7790 print=1521-0626 online DOI: 10.1080/10407790591009251

(2)

the Nusselt number, has to be computed as a derived variable. Wong and Baker [5] used the velocity–vorticity equations to study natural convection in a cubic enclos-ure. However, the two main issues encountered with the velocity–vorticity formu-lation are (1) that the number of variables is increased from four to six as compared to the primitive-variable form for three-dimensional problems, and (2) enforcing the vorticity definition at the solid boundaries to satisfy the continuity equation [7]. As a solution to the first problem, Davis and Carpenter [4] solved only three governing equations by considering two velocities and one vorticity as the primitive variables and computed the remaining three field variables as secondary variables. However, they handled the convective part of the governing equations using a predictor–corrector scheme, thus deviating not much from the existing algo-rithms for treating the pressure term.

With regard to the second problem, enforcing the vorticity definition at the wall, a higher-order scheme is essential to compute the boundary vorticity values in terms of the velocity gradients. Wong and Baker [5] used a second-order-accurate Taylor’s series expansion to compute the vorticity values at the bound-aries. Davis and Carpenter [4] used an integral approach for vorticity definition at the boundary, as followed by Guevremont et al. [10]. Such low-order schemes have been commonly employed with the finite-element methods [4, 5, 8, 9] and the finite-difference methods [2, 3] in the solution of velocity–vorticity equations. Since the finite-difference method and the finite-element method use only low-order approximations for spatial discretization of the differential operators, a large number of grid points have to be used to achieve benchmark solutions. Though higher-order elements can be used in the case of the finite-element method, the implementation becomes unwieldy for the case of three-dimensional flow problems. One way to overcome both the problems discussed above is to use a higher-order approximating scheme such as the generalized differential quadrature (GDQ) method.

The GDQ method approximates a partial differential equation with a higher-order polynomial by using all the grid points in the coordinate direction [11]. Hence, even with a coarse mesh, the vorticity boundary values can be evaluated with

NOMENCLATURE

g acceleration due to gravity

L length

Numean mean Nusselt number throughout the

cavity

Nuoverall overall Nusselt number on the

boundary at x¼ 0.5 or x ¼ 0.5 Pr Prandtl number Ra Rayleigh number Re Reynolds number T temperature Ta reference temperature

u; v; w dimensionless velocity in the x; y; z directions

x; y; z dimensionless Cartesian coordinate directions

a diffusion coefficient b thermal expansion coefficient C boundary of the closed domain m dynamic viscosity

t kinematic viscosity q mass density

n; g; f dimensionless vorticity in the x; y; z directions

(3)

higher-order accuracy naturally, without the need for a separate approximation scheme. Though the GDQ method has been used to solve the Navier–Stokes equa-tions in two and three dimensions, works are available only for the stream function– vorticity and primitive-variable forms of the Navier–Stokes equations [12–14]. Further, most of the works on 3-D natural convection, except Mallinson and de Vahl Davis [6] and Wong and Baker [5], used the primitive-variable form of the Navier–Stokes equations. Numerical solutions for natural convection in a cubic cavity have also been reported using numerical schemes based on a pseudo-spectral Chebyshev algorithm [15], the lattice Boltzmann model [16], the finite-difference method [17], the finite-element method [5], etc. Wakashima and Saitoh [18] used a vorticity–vector potential formulation to produce benchmark solutions for natural convection in a cubic cavity using uniform grids of size 1203. All the above numerical procedures employed finer meshes of the order of 513and above in order to achieve benchmark solutions for 103 Ra  106. Wong and Baker [5] used a mesh of size

483to solve the velocity–vorticity equations by the finite-element method using a parallel computational algorithm. Recently, Lo et al. [19] reported results for natural convection in a square cavity using a simple numerical algorithm based on the velocity–vorticity formulation and the GDQ method for Ra¼ 103

–106. They demonstrated that the use of the GDQ method in combination with the velocity–vorticity formulation requires only a much coarser mesh compared to other numerical schemes to achieve the benchmark solutions. In the present work, a numerical scheme based on the GDQ method is proposed to solve the velocity–vorticity form of the Navier–Stokes equations with the aim of offseting the increased computa-tional effort by the increased number of field variables for three-dimensional problems in the velocity–vorticity formulation. The methodology underlying the 2-D model of the velocity–vorticity formulation using the GDQ method developed by Lo et al. [19] is now extended to 3-D work.

All the flow variables and the temperature scalar are coupled together through their coefficient matrices resulting from the approximation of the governing equa-tions by the GDQ method. The time derivatives in the governing equaequa-tions are discretized using a second-order-accurate, implicit, three-time-level scheme. The final simultaneous equations are solved using a BiCG iterative scheme by storing only the nonzero entries of the coefficient matrix. The efficiency of the proposed numerical scheme is demonstrated by computing flow and heat transfer results for natural convection in a cubic cavity for the Rayleigh number range of 103–106, as discussed in the following sections.

2. GENERALIZED DIFFERENTIAL QUADRATURE METHOD

The method of GDQ replaces a given partial space derivative of a function f(x) by a linear weighted sum of the function values at the discrete sample points considered along a coordinate direction, resulting in a set of algebraic equations. Hence the GDQ method can be used to obtain numerical solution of partial differ-ential equations with higher-order accuracy. The details about this method can be obtained elsewhere [12, 20]. For a function of three variables f(x,y,z), the pth-order derivatives, qth-order derivatives, and rth-order derivatives of the function with

(4)

respect to x, y, and z can be obtained as fxðpÞðxi; yj; zkÞ ¼ XL l¼1 AðpÞi;lfðxl; yj; zkÞ p¼ 1; 2; . . . ; L  1 ð1aÞ fxðqÞðxi; yj; zkÞ ¼ XM m¼1 BðqÞj;mfðxl; ym; zkÞ q¼ 1; 2; . . . ; M  1 ð1bÞ fxðrÞðxi; yj; zkÞ ¼ XN n¼1 Ck;nðrÞfðxi; yj; znÞ r¼ 1; 2; . . . ; N  1 for i¼ 1; 2; . . . ; L; j ¼ 1; 2; . . . ; M; k ¼ 1; 2; . . . ; N ð1cÞ

where l, m, n represent the indices for the grid points in the x, y, and z coordinates, respectively; L, M, N are the number of grid points in the x, y, z directions, respect-ively; and AðpÞi;l; BðqÞj;m; Ck;nðrÞ are the weighting coefficients. The first-order weighting coefficients Að1Þi;l; Bð1Þj;m; Cð1Þk;ncan be determined as follows:

Að1Þi; j ¼ L ð1Þðx iÞ ðxi xjÞLð1ÞðxjÞ i; j¼ 1; 2; . . . ; L; but j 6¼ i ð2aÞ Bð1Þi; j ¼ M ð1Þðy iÞ ðyi yjÞMð1ÞðyjÞ i; j¼ 1; 2; . . . ; M; but j 6¼ i ð2bÞ Ci; jð1Þ¼ N ð1Þðz iÞ ðzi zjÞNð1ÞðzjÞ i; j¼ 1; 2; . . . ; N; but j 6¼ i ð2cÞ in which Lð1ÞðxiÞ ¼ YL j¼1; j6¼i ðxi xjÞ Mð1ÞðyiÞ ¼ YM j¼1; j6¼i ðyi yjÞ Nð1ÞðziÞ ¼ YN j¼1; j6¼i ðzi zjÞ ð3Þ

Similarly, the weighting coefficients for the second- and higher-order derivatives can be obtained as AðpÞi; j ¼ p Aðp1Þi;i Að1Þi; j A ðp1Þ i; j xi xj ! for i; j¼ 1; 2; . . . ; L; but j 6¼ i; l ¼ 2; 3; . . . ; L  1 ð4aÞ

(5)

BðqÞi; j ¼ q Bðq1Þi;i Bð1Þi; j B ðq1Þ i; j yi yj ! for i; j¼ 1; 2; . . . ; M; but j 6¼ i; m ¼ 2; 3; . . . ; M  1 ð4bÞ Ci; jðrÞ¼ r Ci;iðr1ÞCi; jð1ÞC ðr1Þ i; j zi zj ! for i; j¼ 1; 2; . . . ; N; but j 6¼ i; n ¼ 2; 3; . . . ; N  1 ð4cÞ When j¼ i, the weighting coefficients are written as

AðpÞi;i ¼  X L j¼1; j6¼i AðpÞi; j i¼ 1; 2; . . . ; L; p ¼ 1; 2; . . . ; L  1 ð5aÞ BðqÞi;i ¼  X M j¼1; j6¼i BðqÞi; j i¼ 1; 2; . . . ; M; q ¼ 1; 2; . . . ; M  1 ð5bÞ Ci;iðrÞ ¼  X N j¼1; j6¼i CðrÞi; j i¼ 1; 2; . . . ; N; r ¼ 1; 2; . . . ; N  1 ð5cÞ

It should be noted from the above equations that the weighting coefficients of the second- and higher-order derivatives can be computed from the first-order deriva-tives themselves.

3. GOVERNING EQUATIONS

The governing equations for natural convection can be described by the incompressible Navier–Stokes equations and the energy equation. Assuming the Boussinesq approximation, the primitive-variable forms of these equations can be written in dimensional form as follows.

Continuity equation: r  ~uu¼ 0 ð6Þ Momentum equations: q~uu qtþ ð~uu   rÞ~uu¼ 1 qrP þ tðr2~uu Þ þ gbðT TaÞ~kk ð7Þ Energy equation: qT qt þ ~uu  ðrT Þ ¼ a r2T ð8Þ on a Cartesian coordinate frame in which x–y represents the horizontal plane and z refers to the vertical direction. In the velocity–vorticity form of the Navier–Stokes

(6)

equations, the vorticity vector is defined as ~ x

x¼ r  ~uu ð9Þ By taking the curl on both sides of Eq. (7) and using the vorticity definition, the momentum conservation equation (7) can be rewritten as

q~xx qt þ ~uu  r~ x x ¼ ~xx r~uuþ tr2~xx þ r  ½gbðT TaÞ~kk ð10Þ

where ~uu¼ ðu; v; wÞ and ~xx¼ ðn;

g;1Þ are the velocity and the vorticity vectors in dimensional form representing the components in the x, y, and z directions, respectively; ~kk is the unit vector in the z direction; T¼ dimensional form of temperature; and Tais the reference temperature. The vorticity transport equation

(10) is a single vector equation with two unknown vectors, ~uu and ~xx. A second vector equation can be obtained by taking curl of the vorticity definition given by Eq. (9). After using the vorticity definition and the continuity equation, the second vector equation is obtained as

r2~uu¼ r  ~xx ð11Þ The above equation is the velocity Poisson equation, which expresses the kinematic relationship between the velocity vector and the vorticity vector. By standard nondi-mensionalization of Eq. (10), the vorticity transport equation in nondimensional form can be written as

q~xx

qt þ ~uu r~xx¼ ~xx r~uuþ Prr

2~

x

xþ Ra PrrðTÞ~kk ð12Þ The nondimensional parameters are defined as Prandtl number Pr¼ t=a and Rayleigh number Ra¼ gb DT L3=at.

Similarly the nondimensional form of the velocity Poisson equation (11) can be rewritten as

r2~uu¼ r  ~xx ð13Þ

Then the energy equation (8) can be written in nondimensional form as qT

qt þ ~uu ðrTÞ ¼ r

2T ð14Þ

Equations (12)–(14) are the final form of the governing equations that characterize the flow and heat transfer during a natural-convection process. These equations have to be solved in a computational domain X which is enclosed by a solid boundary C. For the case of natural convection in a differentially heated cubic cavity, no-slip velocity boundary conditions are assumed on all the boundary walls. The Dirichlet boundary conditions for velocity on all the walls can be imposed as

(7)

The boundary conditions for the vorticity transport equations are computed from its definition given by Eq. (9). The energy equation is solved by assuming Dirichlet temperature boundary conditions equal to 0.5 and 0.5, respectively, on the left and the right walls of the cavity. Other walls of the cavity are assumed to be adia-batic for heat transport. The Dirichlet and Neumann boundary conditions for the energy equation can be written as

T¼ Tb ð16aÞ qT qy ¼ 0 qT qz ¼ 0 ð16bÞ 4. NUMERICAL SOLUTION

In the classical numerical solution algorithms [5, 10], an iterative solution pro-cedure is adopted to resolve the coupling between Eqs. (12)–(14). This will result in huge computational effort in terms of computer memory and computational time for 3-D problems. In order to achieve a significant saving in computational effort, a coupled solution scheme via a single global matrix for all the field variables is adopted in the present work. Since all the field variables are solved using a coupled algorithm, the boundary vorticity values are computed implicitly, without the need to compute the boundary vorticity values externally using a separate scheme such as the Taylor’s series expansion scheme used by Wong and Baker [5]. The coupled algorithm has also enabled the implicit enforcement of the kinematic Poisson equa-tion as well as the coupling of the velocity and the vorticity at the wall. In addiequa-tion to this, the use of the GDQ method enables the approximation of the vorticity defi-nition at the boundaries with higher-order accuracy. Hence the continuity constraint and the conservation of the solenoidality of vorticity field are easily satisfied. The time derivatives of the vorticity transport equations and energy equation are discre-tized using a second-order-accurate, three-time-level implicit scheme. For example, the time derivative of a variable u can be approximated as

qu qt ¼

3utþ1 4utþ ut1

2 Dt

4.1. Approximation of the Governing Equations Using the GDQ Method

Application of the GDQ method to spatial discretization of the governing equations results in a set of algebraic equations. The GDQ form of the governing equations can be expressed as follows.

Velocity–Poisson equations in the three coordinate directions: XL l¼1 Að2Þi;lul; j;kþ XM m¼1 Bð2Þj;mui;m;kþ XN n¼1 Cð2Þk;nui; j;nþ XM m¼1 Bð1Þj;mfi;m;kX N n¼1 Ck;nð1Þgi; j;n¼ 0 ð17aÞ

(8)

XL l¼1 Að2Þi;lvl; j;kþ XM m¼1 Bð2Þj;mvi;m;kþ XN n¼1 Cð2Þk;nvi; j;n XL l¼1 Að1Þi;lfl; j;kþX N n¼1 Cð1Þk;nni; j;n¼ 0 ð17bÞ XL l¼1 Að2Þi;lwl; j;kþ XM m¼1 Bð2Þj;mwi;m;kþ XN n¼1 Ck;nð2Þwi; j;nþ XL l¼1 Að1Þi;lgl; j;k XM m¼1 Bð1Þj;mni;m;k¼ 0 ð17cÞ Vorticity transport equations in the three coordinate directions:

3 2 Dt   ðni; j;kÞ tþ1 þ uqi; j;kX L l¼1 Að1Þi;lnl; j;kþ v q i; j;k XM m¼1 Bð1Þj;mni;m;kþ w q i; j;k XN n¼1 Ck;nð1Þni; j;n !tþ1  nqi; j;kX L l¼1 Að1Þi;lul; j;kþ gqi; j;k XM m¼1 Bð1Þj;mui;m;kþ fqi; j;k XN n¼1 Ck;nð1Þui; j;n !tþ1  Pr X L l¼1 Að2Þi;lnl; j;kþX M m¼1 Bð2Þj;mni;m;kþX N n¼1 Cð2Þk;nni; j;n !tþ1  Ra Pr X M m¼1 Bð1Þj;mTi;m;k !tþ1 ¼ 4 2 Dt   ntijkþ 1 2 Dt   nt1ijk ð18aÞ 3 2 Dt   ðgi; j;kÞ tþ1 þ uqi; j;kX L l¼1 Að1Þi;lgl; j;kþ vqi; j;kX M m¼1 Bð1Þj;mgi;m;kþ wqi; j;kX N n¼1 Ck;nð1Þgi; j;n !tþ1  nqi; j;kX L l¼1 Að1Þi;lvl; j;kþ gqi; j;k XM m¼1 Bð1Þj;mvi;m;kþ fqi; j;k XN n¼1 Ck;nð1Þvi; j;n !tþ1  Pr X L l¼1 Að2Þi;lgl; j;kþX M m¼1 Bð2Þj;mgi;m;kþX N n¼1 Cð2Þk;ngi; j;n !tþ1 þ Ra Pr X L l¼1 Að1Þi;lTl; j;k !tþ1 ¼ 4 2 Dt   gtijkþ 1 2 Dt   gt1ijk ð18bÞ 3 2 Dt   ðfi; j;kÞ tþ1 þ uqi; j;kX L l¼1 Að1Þi;l fl; j;kþ vqi; j;k XM m¼1 Bð1Þj;mfi;m;kþ wqi; j;kX N n¼1 Cð1Þk;nfi; j;n !tþ1  nqi; j;kX L l¼1 Að1Þi;lwl; j;kþ gqi; j;k XM m¼1 Bð1Þj;mwi;m;kþ fqi; j;k XN n¼1 Ck;nð1Þwi; j;n !tþ1  Pr X L l¼1 Að2Þi;lfl; j;kþX M m¼1 Bð2Þj;mfi;m;kþX N n¼1 Ck;nð2Þfi; j;n !tþ1 ¼ 4 2 Dt   ftijkþ 1 2 Dt   ft1ijk ð18cÞ

(9)

Energy equation: 3 2 Dt   ðTi; j;kÞtþ1þ uqi; j;k XL l¼1 Að1Þi;lTl; j;kþ vqi; j;k XM m¼1 Bð1Þj;mTi;m;kþ wqi; j;k XN n¼1 Ck;nð1ÞTi; j;n !tþ1  X L l¼1 Að2Þi;lTl; j;kþ XM m¼1 Bð2Þj;mTi;m;kþ XN n¼1 Ck;nð2ÞTi; j;n !tþ1 ¼ 4 2 Dt   Tijkt þ 1 2 Dt   Tijkt1 ð19Þ Combining Eqs. (17)–(19), all the seven field variables can be represented by means of a single global matrix as

A 0 0 0 D C 0 0 A 0 D 0 B 0 0 0 A C B 0 0 E 0 0 F 0 0 G 0 E 0 0 F 0 H 0 0 E 0 0 F 0 0 0 0 0 0 0 I 2 6 6 6 6 6 6 6 6 4 3 7 7 7 7 7 7 7 7 5 u v w n g f T 8 > > > > > > > < > > > > > > > : 9 > > > > > > > = > > > > > > > ; ¼ fu fv fw fn fg ff fT 8 > > > > > > > > < > > > > > > > > : 9 > > > > > > > > = > > > > > > > > ; ð20Þ

where A, B, C, D, E, F, G, and H are the influence matrices which represent the vari-ous differential operators that appear in the GDQ approximation of Eqs. (17)–(19), u, v, w, n; g; f, and T are the vectors representing the unknown field variables, and fu; fv; fw; fn; fg; ff; fT are the load vectors. The influence matrices and the load

vec-tors are computed as ½A ¼X L l¼1 Að2Þi;l þX M m¼1 Bð2Þj;mþX N n¼1 Ck;nð2Þ ½B ¼X L l¼1 Að1Þi;l ½C ¼X M m¼1 Bð1Þj;m ½D ¼X N n¼1 Cð1Þk;n ½E ¼  nqi; j;kX L l¼1 Að1Þi;l þ gqi; j;kX M m¼1 Bð1Þj;mþ fqi; j;kX N n¼1 Ck;nð1Þ ! ½F  ¼ 3 2 Dtþ u q i; j;k XL l¼1 Að1Þi;l þ vqi; j;kX M m¼1 Bð1Þj;mþ wqi; j;kX N n¼1 Ck;nð1Þ !  Pr X L l¼1 Að2Þi;l þX M m¼1 Bð2Þj;mþX N n¼1 Cð2Þk;n ! ½G ¼ Ra PrX M m¼1 Bð1Þj;m ½H ¼ Ra PrX L l¼1 Að1Þi;l ½I ¼ 3 2 Dtþ u q i; j;k XL l¼1 Að1Þi;l þ vqi; j;kX M m¼1 Bð1Þj;mþ wqi; j;kX N n¼1 Ck;nð1Þ !  X L l¼1 Að2Þi;l þX M m¼1 Bð2Þj;mþX N n¼1 Cð2Þk;n !

(10)

½ fu ¼ 0 ½fv ¼ 0 ½fw ¼ 0 ½ fn ¼ 4 2 Dtn tþ 1 2 Dtn t1 ½f g ¼ 4 2 Dtg tþ 1 2 Dtg t1 ½ ff ¼ 4 2 Dtf tþ 1 2 Dtf t1 ½ f T ¼ 4 2 DtðTÞ t þ 1 2 DtðTÞ t1

where q is the iterative number index and t is the time-level index.

4.2. Determination of Vorticity Boundary Conditions

The vorticity boundary conditions for all the three vorticity transport equa-tions in the principal coordinate direcequa-tions are computed using the vorticity defi-nition given by Eq. (9). The velocity gradients used in this defidefi-nition can be computed with a higher-order accuracy using the first-order weighting coefficients Að1Þi;l; Bð1Þj;m; Ck;nð1Þ. By applying the GDQ approximation to the vorticity definition given by Eq. (9), the three vorticity components on a boundary can be expressed as

ni; j;kX M m¼1 Bð1Þj;mwi;m;kþ XN n¼1 Ck;nð1Þvi; j;n¼ 0 ð21aÞ gi; j;kþ XL l¼1 Að1Þi;lwl; j;k XN n¼1 Ck;nð1Þui; j;n¼ 0 ð21bÞ fi; j;k XL l¼1 Að1Þi;lvl; j;kþ XM m¼1 Bð1Þj;mui;m;k¼ 0 ð21cÞ

The Dirichlet boundary conditions for the temperature are expressed as TL; j;k¼ 0:5 j¼ 1; . . . ; M; k¼ 1; . . . ; N

T1; j;k¼ 0:5 j¼ 1; . . . ; M; k¼ 1; . . . ; N

ð22Þ The adiabatic boundary conditions can be achieved by computing the normal derivatives of the temperature at the adiabatic walls and equating them to zero. Hence the GDQ forms of the adiabatic boundary conditions can be represented by the following expressions [20]:

Ti;1;k¼ 1 Bð1Þ1;1Bð1ÞM;M Bð1Þ1;MBð1ÞM;1 X M1 m¼2 Bð1Þ1;MBð1ÞM;m Bð1ÞM;MBð1Þ1;m   Ti;m;k " # ð23aÞ Ti;M;k¼ 1 Bð1Þ1;1Bð1ÞM;M Bð1Þ1;MBð1ÞM;1 X M1 m¼2 Bð1ÞM;1Bð1Þ1;m Bð1Þ1;1Bð1ÞM;m   Ti;m;k " # ð23bÞ Ti; j;1¼ 1 C1;1ð1ÞCð1ÞN;N Cð1Þ1;NCN;1ð1Þ X N1 n¼2 Cð1Þ1;NCN;nð1Þ  Cð1ÞN;NCð1Þ1;n   Ti; j;n " # ð23cÞ

(11)

Ti; j;N ¼ 1 C1;1ð1ÞCN;Nð1Þ  C1;Nð1ÞCð1ÞN;1 X N1 n¼2 CN;1ð1ÞC1;nð1Þ C1;1ð1ÞCN;nð1Þ   Ti; j;n " # ð23dÞ

Equations (23a)–(23d) also involve the implicit scheme for the Neumann boundary conditions. The simultaneous equations resulting from the global matrix system of Eq. (20) are solved using a BiCG iterative equation solver. Since the coefficient matrix of the global matrix system is sparse, only the nonzero entries are stored in column storage format. In the present study we employed the following convergence criteria to terminate the iterative process used to resolve the coupling between the field variables at a given time step:

ðutþ1qþ1 utþ1q Þ=utþ1q      106 ðvtþ1qþ1 vtþ1q Þ=vtþ1q      106 ðwtþ1qþ1 wtþ1q Þ=wtþ1q      106 ðntþ1qþ1 n tþ1 q Þ=n tþ1 q      106 ðgtþ1 qþ1 gtþ1q Þ=gtþ1q      10 6 ðftþ1qþ1 fqtþ1Þ=fqtþ1  106 ðTqþ1tþ1 Tqtþ1Þ=Tqtþ1      106

In the successive time step, we used the velocity, vorticity, and temperature compo-nents at the previous time step as the initial guess for the next iteration. The compu-tations are carried out until steady-state conditions are reached. The convergence criteria used in the time loop to achieve steady-state conditions are

ðutþ1 utÞ=ut     106 ðvtþ1 vtÞ=vt  10 6 ðwtþ1 wtÞ=wt  106 ðntþ1 ntÞ=nt     106 ðgtþ1 gtÞ=gt  10 6 ðftþ1 ftÞ=ft  106 ðTtþ1 TtÞ=Tt     106

For the GDQ method, the mesh point distribution in the three spatial coordinates is assumed to be the same and is expressed as

xi¼ cos½p=ð2LÞ  cos½ð2i  1Þp=ð2LÞ cos½p=ð2LÞ  cos½ð2L  1Þp=ð2LÞ i¼ 1; 2; . . . ; L yj¼ cos½p=ð2MÞ  cos½ð2j  1Þp=ð2MÞ cos½p=ð2MÞ  cos½ð2M  1Þp=ð2MÞ j¼ 1; 2; . . . ; M zk¼ cos½p=ð2NÞ  cos½ð2k  1Þp=ð2NÞ cos½p=ð2NÞ  cos½ð2N  1Þp=ð2NÞ k¼ 1; 2; . . . ; N ð24Þ

where L, M, N are the number of grid points in the x, y, and z directions, respectively.

5. RESULTS AND DISCUSSIONS

The accuracy of the present numerical procedure is verified by (1) a grid inde-pendence study and (2) validating the predicted results with the benchmark

(12)

solutions. Air is assumed to be the working fluid, with Pr¼ 0:71 for the natural-convection simulation cases.

5.1. Grid Independence Study

For the grid independence study, numerical results are obtained for natural convection in a cubic cavity (see Figure 1) for Ra¼ 103and 104. The computational grid was refined successively from 153to 193and 253. The following mean and over-all Nusselt numbers are considered as the parameters for comparison purposes. 1. The mean Nusselt number throughout the cavity is

NumeanðyÞ ¼ Z 1 0 qTðy; zÞ qx     x¼0:5; or x¼0:5 dz ð25Þ

2. The overall Nusselt number on the boundary at x¼ 0.5 or x ¼ 0.5 is

Nuoverall¼

Z 1 0

NumeanðyÞ dy ð26Þ

The effect of grid refinement on the values of the mean and the overall Nusselt numbers is shown in Table 1 for Ra¼ 103 and 104. As the mesh is refined successively from 153to 193and 253, there is a consistent improvement in the accu-racy of the predicted values of the Nusselt numbers for both values of the Rayleigh number, with no difference being observed between the grids 193and 253. This indi-cates that the 253grid can produce grid-independent numerical results and hence will be used throughout this work. The grid independence study is also demonstrated by

(13)

plotting the temperature profiles along the centerline of the symmetry plane at y¼ 0.5. Figures 2a and 2b show the effect of grid refinement on the temperature pro-files for Ra¼ 103

and 104, respectively. The temperature profiles of all three meshes almost coincide with each other for both Ra¼ 103

and 104. The effect of grid refine-ment on the velocity profile u–z along the vertical symmetric central line is shown in Figures 3a-1 and 3a-2 for Ra¼ 103

and 104, respectively. Similarly, Figures 3b-1 and 3b-2 represent the x–w plots along the horizontal symmetric central line for Ra¼ 103

and 104, respectively. The exact coincidence of the above velocity profiles for all three meshes indicates that the present numerical algorithm based on the velocity– vorticity formulation and the GDQ method predicts the temperature and flow results consistent with grid refinement.

5.2. Validation of Numerical Results

For the case of a natural-convection problem, the effect of Rayleigh number on the convective heat transfer between the end walls is evaluated in terms of the mean and the overall Nusselt number for Ra¼ 103, 104, 105, and 106. In order to validate the present numerical algorithm, the predicted results are compared with the results obtained by three different numerical schemes: (1) the pseudo-spectral Chebyshev algorithm (Tric et al., 2000), (2) the lattice Boltzmann model (Peng et al., 2003), and (3) the finite-difference method (Fusegi et al., 1991). Table 2 shows the compari-son of the present results with the results of the above three numerical schemes for Ra¼ 103

and 104. For both Rayleigh numbers, the present results are in close agree-ment with the results of the above three numerical schemes, which adopted grids finer than 253, used in the present method. Table 3 shows a similar comparison for Ra¼ 105

and 106. Even for higher Rayleigh number values, the present results obtained using 253grid are closer to the results of Tric et al. [15], who used a grid size of 813. A close look at the tabular values indicates that the present predictions are in better agreement with the results of Tric et al. [15] compared to the results of the other two numerical schemes [16, 17], obtained using finer meshes. For all four values of Rayleigh number considered in the present analysis, the average and over-all Nusselt numbers on the isothermal wover-all have been predicted with less than 0.4% error in comparison to the benchmark results of Tric et al. [15]. Though Wong and Baker [5] also solved the velocity–vorticity equations to study natural convection in a cubic cavity, they had to use a fine mesh of size 483in order to compute the vorticity boundary values using a second-order-accurate Taylor’s series scheme, since they employed the finite-element method. Moreover, the numerical schemes used in

Table 1. Grid independence study results for Ra¼ 103and 104

Ra 103 104 Method GDQ GDQ GDQ PSC [15] GDQ GDQ GDQ PSC [15] Mesh size 153 193 253 323 153 193 253 613 Nuoverall 1.069 1.070 1.070 1.070 2.052 2.054 2.054 2.0542 Numean(y¼ 0.5) 1.085 1.087 1.087 1.0873 2.248 2.251 2.251 2.2505

(14)

[5, 15–17] were implemented on high-speed computers such as Cray computers. Using the present numerical algorithm, the benchmark results could be achieved using a coarse mesh on a Pentium IV personal computer. The reason is that the

Figure 2. Temperature profiles in the centerline of the symmetry plane at y¼ 0.5 for (a) Ra ¼ 103,

(15)

Figure 3. Velocity profiles along the centerline in the symmetry plane at y¼ 0.5 for (a) Ra ¼ 103, (b)

Ra¼ 104.

Table 2. Validation results for Ra¼ 103and 104: Comparison of Nusselt number at the isothermal wall,

x¼ 0.5 Ra 103 104 Method GDQ PSC [15] LBM [16] FDM [17] GDQ PSC [15] LBM [16] FDM [17] Mesh size 253 613 813 323 253 613 61 45  45 623 Nuoverall 1.070 1.070 1.075 1.085 2.054 2.0542 2.085 2.100 Numean(y¼ 0.5) 1.087 1.0873 1.097 1.105 2.251 2.2505 2.304 2.302

(16)

Table 3. Validation results for Ra¼ 105and 106: Comparison of Nusselt number at the isothermal wall, x¼ 0.5 Ra 105 106 Method GDQ PSC [15] LBM [16] FDM [17] GDQ PSC [15] LBM [16] FDM [17] Mesh size 253 813 91  65  65 623 253 813 N=A 623 Nuoverall 4.335 4.3370 4.378 4.361 8.666 8.6407 N=A 8.770 Numean(y¼ 0.5) 4.610 4.6127 4.658 4.646 8.91 8.8771 N=A 9.012

Figure 4. Velocity vectors at y¼ 0.5 plane for (a) Ra ¼ 103, (b) Ra

¼ 104, (c) Ra

¼ 105, (d) Ra

(17)

GDQ method makes use of all the grid points in a given coordinate direction while approximating a differential with respect to that coordinate direction. Hence, while computing the vorticity values at the boundary, there is no need to use very fine mesh near the boundary as required in the finite-element scheme [5].

5.3. Results on Flow and Thermal Characteristics

The ability of the present coupled numerical scheme is further demonstrated by plotting the velocity, vorticity, and temperature contours on various symmetric mid-planes along the principal axes of the cubic cavity. Figures 4a–4d show the velo-city vector distributions on the xz plane at y ¼ 0.5 for Ra ¼ 103, 104, 105, and 106. As the Rayleigh number increases, the flow is characterized by a combination of an

Figure 5. Vorticity contours at y¼ 0.5 plane for (a) Ra ¼ 103, (b) Ra

¼ 104, (c) Ra

¼ 105, (d) Ra

(18)
(19)

almost stagnant interior fluid core and distinct boundary layers near the side walls. Since the vorticity is the primitive variable in the present formulation, the vortical flow patterns due to the natural convection are obtained directly. The vorticity con-tour distributions on the xz plane at y ¼ 0.5 are shown in Figures 5a–5d for Ra¼ 103, 104, 105, and 106. The vorticity values near the isothermal walls increase with increase in the value of Rayleigh number, resulting in a strong circulatory motion. Higher values of vorticity on the isothermal walls indicate the intense convection effect on the fluid structure. The above constant vorticity contours clearly demonstrate again the existence of a near-stagnant interior core along with the distinct boundary layers near the end walls.

Figures 6a and 6b show the variations of u velocity along the vertical centerline and w velocity along the horizontal centerline at the y¼ 0.5 plane for Ra ¼ 103, 104, 105, and 106. The peak values of the horizontal and the vertical velocities increase due to the intensified convective activities with increase in Rayleigh number. The

Figure 7. Contour maps of temperature at y¼ 0.5 plane for (a) Ra ¼ 103, (b) Ra

¼ 104, (c) Ra

¼ 105,

(20)

steep rise in the w-velocity gradient at the points on the hot and the cold walls also confirm the increased convective activity at higher Rayleigh number values as observed in Figure 6b. The uz and xw velocity profiles predicted in the present work are in qualitative agreement with the results of Wong and Baker [5].

The temperature contours on the xz plane at y ¼ 0.5 are shown in Figure 7 for Ra¼ 103, 104, 105, and 106. With increase in Rayleigh number, a high degree of convection is observed such that distinct thermal boundary layers start appearing near the isothermal walls. The thickness of the thermal boundary layer decreases as Rayleigh number increases. The 3-D variations of temperature field can be under-stood by plotting the temperature contours on the xy plane at z ¼ 0:5 as shown in Figures 8a–8d for different Rayleigh numbers. Sharp temperature changes are noted at the left and right edges of the xy plane due to the increase in the Rayleigh num-ber and are clearly indicated in the above figures in combination with temperature contours on the vertical plane at y¼ 1. Figures 9a–9d represent the temperature

Figure 8. Contour map of temperature T at z ¼ 0.5 for Pr ¼ 0.71 and Rayleigh numbers (a) Ra ¼ 103, (b)

Ra¼ 104, (c) Ra

¼ 105, (d) Ra

(21)

contours on the yz plane at x ¼ 0 for different Rayleigh numbers. The vertical stratification of the temperature field becomes distinct with increase in Rayleigh number, as indicated in these figures. The accurate predictions of the coincidence of the temperature profiles on the x¼ 0 plane with those on the y ¼ 1 plane demon-strate that the present numerical scheme is capable of predicting the 3-D variation of the temperature field accurately. In order to understand the 3-D convective effect, the variation of the mean Nusselt number along the y direction are plotted in Figures 10a–10d for Ra¼ 103, 104, 105, and 106. As the Rayleigh number increases, the mean Nusselt number increases as the symmetry plane is approached, and the peak of the mean Nusselt number occurs at y¼ 0:5. Due to the presence of intensive convective flow in the y direction, the convective heat transfer is enhanced, resulting in the appearance of two minor peaks at y¼ 0:2 and 0.8.

Figure 9. Contour map of temperature T at x¼ 0 for Pr ¼ 0.71 and Rayleigh numbers (a) Ra ¼ 103,

(b) Ra¼ 104, (c) Ra

¼ 105, (d) Ra

(22)

6. CONCLUSIONS

In the present numerical study, natural convection in a differentially heated cubic cavity has been studied using the Navier–Stokes equations in velocity–vorticity form and the energy equation. The governing equations have been solved by a coupled numerical algorithm based on the GDQ method. A grid independence study conducted for Ra¼ 103and 104indicated that the accuracy of the numerical predic-tions of heat transfer parameters increases consistently with the mesh refinement. Validation test results obtained for Ra¼ 103, 104, 105, and 106show that the bench-mark results could be achieved using only a 253mesh, which is much coarser com-pared to the finer grids used in other benchmark algorithms. As the Rayleigh number increases, the convective flow field is intensified near the end walls in the y direction, which was confirmed by the appearance of two minor peaks of the mean Nusselt number. The thickness of the boundary layer near the isothermal walls

Figure 10. Distribution of the mean Nusselt number along the y direction for (a) Ra¼ 103, (b) Ra

¼ 104,

(c) Ra¼ 105, (d) Ra

(23)

decreases with increase in Rayleigh number. The results obtained demonstrate that the combination of the GDQ method and the velocity–vorticity form of the Navier– Stokes equations is an efficient numerical procedure to study flow and heat transfer in a differentially heated cubic enclosure.

REFERENCES

1. H. Fasel, Investigation of the Stability of Boundary Layers by a Finite-Difference Model of the Navier-Stokes Equations, J. Fluid Mech., vol. 78, pp. 355–383, 1976.

2. M. Napolitano and L. A. Catalano, A Multigrid Solver for the Vorticity–Velocity Navier– Stokes Equations, Int. J. Numer. Meth. Fluids, vol. 13, pp. 49–59, 1993.

3. G. Guj and F. Stella, A Vorticity–Velocity Method for the Numerical Solution of 3D Incompressible Flows, J. Comput. Phys., vol. 106, pp. 286–298, 1993.

4. C. Davis and P. W. Carpenter, A Novel Velocity–Vorticity Formulation of the Navier-Stokes Equations with Application to Boundary Layer Disturbance Evolution, J. Comput. Phys., vol. 172, pp. 119–165, 2001.

5. K. L. Wong and A. J. Baker, A 3D Incompressible Navier–Stokes Velocity–Vorticity Weak Form Finite Element Algorithm, Int. J. Numer. Meth. Fluids, vol. 38, pp. 99–123, 2002.

6. G. D. Mallinson and G. de Vahl Davis, Three-Dimensional Natural Convection in a Box: A Numerical Study, J. Fluid Mech., vol. 83, pp. 1–31, 1977.

7. O. Daube, Resolution of the 2D Navier-Stokes Equations in Velocity–Vorticity Form by Means of an Influence Matrix Technique, J. Comput. Phys., vol. 103, pp. 402–414, 1992. 8. M. S. Ingber and S. N. Kempka, A Galerkin Implementation of the Generalized Helmholtz Decomposition for Vorticity Formulations, J. Comput. Phys., vol. 169, pp. 215–237, 2001.

9. D. C. Lo and D. L. Young, Arbitrary Lagrangian-Eulerian Finite Element Analysis of Free Surface Flow Using a Velocity-Vorticity Formulation, J. Comput. Phys., vol. 195, pp. 175–201, 2004.

10. G. Guevremont, W. G. Habashi, and M. M. Hafez, Finite Element Solution of the Navier-Stokes Equations by a Velocity–Vorticity Method, Int. J. Numer. Meth. Fluids, vol. 10, pp. 461–475, 1990.

11. R. E. Bellman, B. G. Kashef, and J. Casti, Differential Quadrature: A Technique for the Rapid Solution of Nonlinear Partial Differential Equations, J. Comput. Phys., vol. 10, pp. 40–52, 1972.

12. C. Shu and B. E. Richards, Application of Generalized Differential Quadrature to Solve 2-Dimensional Incompressible Navier–Stokes Equations, Int. J. Numer. Meth. Fluids, vol. 15, pp. 791–798, 1992.

13. C. Shu, K. S. Yeo, and Q. Yao, An Efficient Approach to Simulate Natural Convection in Arbitrarily Eccentric Annuli by Vorticity-Stream Function Formulation, Numer. Heat Transfer A, vol. 38, pp. 739–756, 2000.

14. C. Shu, L. Wang, and Y. T. Chew, Numerical Computation of Three-Dimensional Incom-pressible Navier-Stokes Equations in Primitive Variable Form by DQ Method, Int. J. Numer. Meth. Fluids, vol. 43, pp. 345–368, 2003.

15. E. Tric, G. Labrosse, and M. Betrouni, A First Incursion into the 3D Structure of Natural Convection of Air in a Differentially Heated Cubic Cavity, from Accurate Numerical Solutions, Int. J. Heat Mass Transfer, vol. 43, pp. 4043–4056, 2000.

16. Y. Peng, C. Shu, and Y. T. Chew, A 3D Incompressible Thermal Lattice Boltzmann Model and its Application to Simulate Natural Convection in a Cubic Cavity, J. Comput. Phys., vol. 193, pp. 260–274, 2003.

(24)

17. T. Fusegi, J. M. Hyun, K. Kuwahara, and B. Farouk, A Numerical Study of Three-Dimensional Natural Convection in a Differentially Heated Cubical Enclosure, Int. J. Heat Mass Transfer, vol. 34, pp. 1543–1557, 1991.

18. S. Wakashima and T. S. Saitoh, Benchmark Solutions for Natural Convection in a Cubic Cavity Using the High-Order Time-Space Method, Int. J. Heat Mass Transfer, vol. 47, pp. 853–864, 2004.

19. D. C. Lo, D. L. Young, and K. Murugesan, GDQ Method for Natural Convection in a Square Cavity Using a Velocity-Vorticity Formulation, Numer. Heat Transfer B, vol. 47, pp. 321–341, 2005.

20. C. Shu, Differential Quadrature and Its Application in Engineering, Springer-Verlag, London, UK, 2000.

(25)

數據

Figure 1. Geometry and boundary conditions for the buoyancy-driven cavity problem.
Table 1. Grid independence study results for Ra ¼ 10 3 and 10 4
Figure 2. Temperature profiles in the centerline of the symmetry plane at y ¼ 0.5 for (a) Ra ¼ 10 3 , (b) Ra ¼ 10 4 .
Figure 3. Velocity profiles along the centerline in the symmetry plane at y ¼ 0.5 for (a) Ra ¼ 10 3 , (b) Ra ¼ 10 4 .
+7

參考文獻

相關文件

Mie–Gr¨uneisen equa- tion of state (1), we want to use an Eulerian formulation of the equations as in the form described in (2), and to employ a state-of-the-art shock capturing

The underlying idea was to use the power of sampling, in a fashion similar to the way it is used in empirical samples from large universes of data, in order to approximate the

Success in establishing, and then comprehending, Dal Ferro’s formula for the solution of the general cubic equation, and success in discovering a similar equation – the solution

Now, nearly all of the current flows through wire S since it has a much lower resistance than the light bulb. The light bulb does not glow because the current flowing through it

In this study, we compute the band structures for three types of photonic structures. The first one is a modified simple cubic lattice consisting of dielectric spheres on the

wavy wall in the Y-Z plane, amplitude A=2.54mm· · · 43 Fig.3-27 Distributions of instantaneous spanwise velocity contours for flat boundary.. in the Y-Z planes, amplitude A=0mm · ·

Ching , “Investigation of a large top wall temperature on the natural convection plume along a heated vertical wall in a square

Detailed information about the mean velocity, pressure, and turbulent kinetic energy is provided to illustrate how the porous hedge affects the cross ventilation in