Numerical modelling of enclosure convection

Acta Astronautica Vol. 22, pp. 36%374, 1990

0094-5765/90 $3.00 + 0.0 Pergamon Press plc

Printed in Great Britain

Numerical Modelling of Enclosure Convection

J. C. Duh Sverdrup Technology, Inc. NASA Lewis Research Center Cleveland, Ohio 44135, U.S.A.

Abstract

v

velocity component in y-direction

A numerical study on the steady and unsteady natural convection in two-dimensional rectangular enclosures has been performed by a time-accurate ADI finite difference scheme. The study covered a range of Rayleigh numbers (Ra) between 103 and 107, aspect ratios (Ar) between 0.2 and 10.0, and tilt angles

o~ 13 8 v 0 co

thermal diffusivity coefficient of thermal expansion residual momentum diffusivity tilt angle under-relaxation parameter

(O) between -90° (heating from bottom) and +90° (heating from top). Various Prandtl numbers (Pr) have been studied, but only the results of water (Pr = 7.0) are reported here due to space limitations. The physics revealed, however, includes the convection phenomena and the Rayleigh-Benard stability, as well as the combined mechanism of these two. The onset of secondary cells is determined by using a velocity map, which is simpler and cleaner, instead of a slreamline plot The critical Ra number for the occurence of these secondary cells is shown to be lower than can be detected by experimental studies. On the Rayleigh-Benard stability part, a second transition from stable single-cell convection to periodic multi-cellniar convection is disclosed. This paper also briefly discusses various issues concerning numerically treating the enclosure convection, including under-relaxadon and the effect of grid size and grid stretching. Nomenclature Ar CX CY H h k L Nu Pr

enclosure aspect ratio, H/L grid stretching constant in x-direction grid stretching constant in y-direction length of the heated/cooled side walls heat transfer coefficient thermal conductivity length of the insulated side walls Nusselt number, h.L/k Prandd number, v/ix

Ra T t u

Rayleigh number, [g.~*(Th-Te)*H3]/(ct.v ) temperatme tin'e velocity component in x-direction

vonicity stream function

Introduction Natural convection in enclosures induced by thermal and/or concentration gradients has attracted extensive research interest in the past because of its significance in many technological applications such as thermal insulation engineering [ 1], geophysics [2], solar collectors [3], materials processing [4], etc. In recent years, a particular emphasis has been placed on natural convection under reduced gravity conditions because of its importance in a number of space-based systems (e.g., thermal storage/control system and cryogenic fluid management system) and applications such as microgravity materials processing. These convective flows also present great challenges to the fundameatal science community because the fluid motion that results is fully coupled to the driving mechanism originated from the density variation. The convective process can be further complicated if the driving potential (temperature and/or concentration gradient) is applied at various directions relative to the gravitational vector. Due to the richness of these phenomena and the vast technological needs, a great deal of experimental work has been conducted in the past 30 years. Carton [5] reviewed the works on high-aspect-ratio enclosures; while most recently, Ostrach [6] presented a comprehensive review which covered various geometries, aspect ratios, boundary conditions, as well as driving mechanisms. A complete survey of the experimental, analytical, and numerical studies can also be found in Gebhart, et al. [7]

This pape~is declareda work of the U.S. Governmentand is not subject to copyrightprotection in the UnitedStates 367

368

J.C. DUH

On the analytical studies of this problem, the potential for further future progress appears to be limited. In the past, most of the analytical works provided only approximate solutions for the limiting values of governing parameters [8,9,10], and the driving potential for convection was limited to the horizontalheating configurations. Bejan and Tien [11] have predicted the Nusselt number in a more general way, however, their results overpredict the Nu values in the intermediate Ra region when Ar is small, as is pointed out by Kamotani, Wang, and Ostrach [12].

( 1 CosO) ~-~q + v*" Ar z

~=-V*

oy

re) ~}

2

2 ~2 where, V* = A t ~ + Ox

~2 1 Ar Oy2

~4)

The nondimensionalization was carried out by using Xl=LX, Zl=HZ, Ul=(Cffl-I)u, Vl=(Ot/L)v, tl=(LH/~t)t, Vl=mg,

On the numerical side, the growth of computational capability has been significant. In the past three decades, cost of performing a given calculation has been reduced by a factor of 10 every 8 years [13], and when improvements in numerical algorithms are taken into account, the reduction in cost is even greater. As a result, extensive numerical solutions have been obtained for wide ranges of the governing parameters, and a good collection of these works can be found in Torrance and Catton [ 14]. As the momentum of computation capability keeps growing, so will the numerical effort on this complicated problem, and an important direction the effort should be towards is that of relating the numerics and the physics more closely, in addition to improving the accuracy, efficiency, and flexibility of the computation.

~I=[0ff(LH)]~, and TI=Tc+T(Th-Tc), where the subscript 1 denotes dimensional variables. The relevant governing parameters are defined as Ar = H/L, Pr = v/a, Ra = [gl3(ThTc)H3]/(otv), Ra L = Ra/Ar 3, Nu = (hL)/k, and NUll = NuAr. A no-slip condition was assumed for velocities on all boundaries. Two isothermal end wails, i.e., T(x=0)= 0, T(x=l)=l, connected by two adiabatic side walls, i.e., VT--0 at y=0 and y=l, were adopted for the temperature boundary conditions. In an effort to minimize the CPU time required to obtain solutions, calculation always started with the results of a previous calculation having similar values of Ra, Pr, ~, and At.

As pointed out by Ostrach [15], some important issues have remained unsolved, or at least confusing; for example, the differentiation between boundary layer and core flows, the dimensional scaling, and the proper description of flow fields such as the formation of secondary cells. This paper is intended to help clarify some of these issues by discussing briefly the effect of grid stretching, as well as the relevance of relaxation parameters for calculation. The use of velocity maps is proposed to resolve the issue of the onset of secondary cells as opposed to the most popular way of plotting the streamline contours. Heat transfer rate and the flow and temperature fields are examined as functions of Ra, Ar, and 0. By tilting the enclosure from -90° to +90° , the two basic modes of transport - convection and stability and their combined mechanism are also disclosed, lq~rnerical Analvsis Thermal convection inside a tilted two-dimensional enclosure, as shown in Figure 1, is considered in this paper. The incompressible, full Navier-Stokes equations that govern this problem can be expressed non-dimensionally in a stream function-vorticity form as follows:

T + ( ~ . V) T -- V *z T "O"T

-J7

+ ( ~ . V) ~ = Pr {Ra [(At SinO) ~

(1) +

L

Figure 1 Coordinate system for a tilted 2-D enclosure The governing equations were solved by using a secondorder time-accurate alternating-direction-implicit (ADI) scheme. Spatial derivatives were approximated with a second-order accurate finite-differencing form on a nonuniform grid system. Time increments (At) were determined such that diagonal dominance was attained. To ensure the full implicitness of the scheme, iteration on vorticity boundary values was performed within each time step until the difference between successive iterations dropped below 10-7 . Steady-state solutions were considered reached when INun+l - Nunl < 10"6, and e~ < 10-3, where e~ is defined as:

e~=

M N

• where the superscript n represents the number of the time step. When oscillatory flow patterns occurred, Nu versus time and phase diagrams were plotted.

369

40th IAF Congress Under-relaxation is necessary for the inner-time step

iterations. With this set of governing equations, the underrelaxation parameter ¢0depends largely on the product of Ra and Pr. With Pr fixed, lower Ra calls for a ¢orange between 0.5 and 0.6. For higher Ra cases, the required o) drops to around 0.2 on a 61x61 grid system. A possible explanation for this is that large Ra numbers amplify any error associated with the calculated temperature field (this can be seen from equation (2)) and pass the enlarged error along to the velocity field, which, in turn, is iterated back to calculate the temperature again. Under-relaxation is required to reduce this error enlargement and to keep the iterations stable. Adequate grid resolution is important in obtaining accurate solutions. Too coarse a grid mesh can yield inaccurate or totallyerroneous answers. Alternatively,the calculation can simply failto converge. A good grid system must be fine enough to resolve the solutions in high-gradient area. Failure to meet the necessary resolution is the main cause for numerical problems and errors. Since high gradients exist near the walls in thisproblem, a stretched grid system which provides finer grids near the wall and coarser grids near the center is employed as follows: A X i + l = C X i ~[x i fori=l to(M-l)/2,

(6)

Axi+t =~X i Ax i fori= (M-I)/2+I toM-l;

(7)

31x31 u n i f o r m - - - ~ 41x41 uniform,......f...'w"~N~ 5 Ix51 u n f f o n n ~ 61x61 u n f f o r m ~ - ~ 71x71unit~

3 lx31 strtched 41x41 stretched

[

~

(b) Enlarged portion Figure 2 Effect of grid resolution on vertical velocity component

Grid Uniform

CPU Nu 3 I x 31 274 5.4791 41 × 41 673 5.2005 51 x 5l 1325 5.0517 61 x 61

2345

4.9643

71 × 71 4453 4.9087 Stretched 31× 31 730 5.0660 4£×4]. 2093 4.8801 Table 1 CPU time required and the calculated Nu (Ra=105, Ar=l, 0= 0o)

where CX is the stretching constant in x-direction; similar expressions can also be written for Ay.

azsa=.=aa~z~am~ A comparison between a uniform grid system with various densities and the stretched grid system is summarized in Figure 2 and Table1. Figure 2 plots the distributionofv along y=0.5. The summary shows that the accuracy of the stretched grid solution is comparwise to that of the solutions obtained from a much denser, uniform grid. This plot also shows that the "grid-independent" stateis reached asymptotically. Table I further demonstrates the computational efficiency by applying the grid stretchingnear the end walls. The initialstatesfor

these calculations all start with a linear temperature distribution between end walls and zero velocities everywhere. 90

The results presented in this paper were obtained using a 61x61 mesh. The stretching parameters CX and CY used vary from 1.0 to 1.1, and depend on the magnitude of Ra and Ar. Higher Ra and/or smaller Ar demand more stretching. All the calculations were performed on a Cray-YMP computer. A wide range of cases have been studied by systematically varying the parameters Ra, Ar, and 0. The heat transferrate is determined by

[ ~T(x----0,y) j Nu = ]0 ~ oy

(8)

80

and Nu H = Nu Ar. The calculated Nussclt numbers of the horizontally-heating configurations (0=0 °) have bccn compared favorably with the experimental data from Kamotani, Wang, and Ostrach [12], and has been presented in [16].

7O ~o

o= ~' 5o 40 3o

In this paper, three issues are addressed: the driving and resistingmechanisms in the inclined heating-from-above configuration, the second transitionof the transport mode for the Rayleigh-Bernard system, and the onset of secondary cells for the horizontally-heating cases.

2o l0

0 x-ooordlnst, e

-10 0.0

0.1

0.2

0.3

0.4

(a) Raffil05, At=l, 0= 0o, y=0.5

0.5

H e a t i n 2 - f r o m - a b o v e Configuration

370

J.C. DUH

When 8 increases from 0 ° to +90 °, i.e., heating from above, the velocity field and the temperature field interact in a way that the stronger the driving force, the more stabilized thermal stratification it will result in the core region. As is shown in Figure 3 for Ra=103, the driving potential is relatively small, the isotherms are only slightly perturbed from the conduction profile which offers little resistance to the flow induced, and only single cell is convected in the enclosure. As Ra number increases, stronger driving potential induces stronger convection currents near the end walls. This velocity field repositions the temperature distribution to a more stabilized thermal stratification near the core region, as can be seen from Figure 4(b). Such a temperature field, in return, prevents the flows from penetrating through the core area, and forms two recirculating corner cells, as is shown in Figure 4(a).

cases, only single cell is convected within the enclosure. As 0 moves to - 9 0 °, however, the enclosure becomes a classical Rayleigh-Benard system, and the stability criteria become the important issue. The first transition of the transport mode associated with such a system is the onset of convection. Chandrasekhar [17] extended the fundamental paper by Rayleigh [18] and solved for the critical Ra numbers for convection to occur in an infinite fluid layer. Catton [19] included rigid side walls into the mathematical model and solved the linearized purturbation equations for critical Rayleigh numbers at various aspect ratios with a unit Pr number. Ozoe, et al. [20] has also numerically studied the heat transfer characteristics of this problem at various Ra numbers up to 8x103. The results agreed with Catton's prediction, but did not reveal the transition of the transport mode at higher Ra numbers. In the current study, various Ra numbers from 103 to 107 have been investigated for the Rayleigh-Bemard system (0=-90°), and an oscillatory flow was observed when Ra reached certain magnitude. Figure 5 depicts the oscillation of Nu as a function of time with Ra=105 and Ar=-0.28. The associated variation of streamlines and isotherms within one period of oscillation are plotted in Figures 6 and 7~ Results of a similar nature have been reported recently by Goldhirsch, et al. [21], but their analysis opted for perfect conducting side walls instead of adiabatic ones.

(a) Streamlines

(b) Isotherms

Figure 3 Streamlines and isotherms at Ra=103, Ar=0.75, and 0--40°

When the Ra number is lowered to 104, only uni-cellular convection is present. This indicates the existence of a second critical Rayleigh number for transition from stable motion to multi-cellular periodic oscillation. Moreover, this instabiliwinduced oscillation can be observed for 0 > - 9 0 ° up to a critical 0 c which is a function of Ra, Pr, and Ar. When 0c
O=Tt

* 0=4r~/3 •=57t/3 0=2M3 . O=2~

(a) Streamlines

(b) Isotherms o=~/3

Figure 4 Streamlines and isotherms at Ra=107, Ar=0.75, and 0--40° Ravlei~h-Bernard Svstem When the enclosure is tilted towards the other direction, a heating-from-below configuration is present where in most

Figure 5 Nu versus time (Ra=105, At=0.28, 0 =-90°)

40th IAF Congress

371

P

114/

m II

I

• IHUlll

>Jj

%

f//4"~ } | I

~NII]

I~lfll i V/I///-/~

I Iltll

II |fl ~ ' ~

I~l w

(a) ~---~/3

(b) ¢=2H3

(c) ,~

(a) ~----~/3

(b) 0--2~/3

(c) # ~

r-: ~

"' f

S///I,O

J

II

I [ G,'~//~_

fJz

.

i ~ . . ~~

Ira\\ %/~

k,,dd (d) ~=4~/3

(e) ¢=5n/3

(f) ¢=2~

Figure 6 Time variation of streamlines within one period of oscillation (Ra=105, Ar=0.28, 0=--90 o)

(d) ¢=4H3

(e) ¢=5~13

(f) ¢=2~

Figure 7 Time variation of isothea-ms within one period of oscillation (Ra=105, At---0.28, 0=--90 o)

In Figure 8, Nu is plotted as a function of 0 with Ra=105 and At---0.28. The vertical bars represent the range of

H e l t i n l from I b ~

_ :

Hemlinll lroln below

J

/ ~mblmOcl e~vmlal~

na I~bma~,

oscillation. In this case, the 0 c is around -75 o, and both

im~Uly

convection and instability influence the range of 0 between 0 c and --913° with time period of the oscillation increasing as 0 approaches 0 c. The numerical calculation thus requires timeaccurate schemes in order to capture the unsteady flows. If the Ra number is further lowered to 103, the flow subsides and the system returns to conduction. This corff'mns the existence of critical Ra number for onset of convection, i.e., the first transition of fluid transport in this Rayleigh-Bemard system.

m ~hm

e

Figure 8 Effect of inclination on heat transfer

(Ra=105,Ar=0.28) Wirtz and Tseng [22] have also studied the effect of inclination on enclosure convection at high Ra numbers. Steady-state governing equations were used in their formulation and the solution failed to converge when 0 dropped below 0c. Wirtz and Tseng attributed the non-convergency to the possible formation of oblique rolls in the third dimension without examining the consequences from thermal instability. Based on

this study, the non-convergence of their calculation appears to be caused by the time-oscillation that the steady-state formulation is Simply unable to capture. The critical Rayleigh number for the second transition as a function of Pr and At, as well as the critical Oc, which

372

J.C. Dun

delineates the oscillation from stable motion, as a function of Ra, Pr and Ar will be presented elsewhere. Onset of Secondary Cells In the case of heating-from-the-side, a small Ra number induces slow convection, only one cell is formed within the enclosure, and the isotherms are only slightly perturbed from the conduction temperature profile. As the Ra number increases, parallel flow forms at the core region, boundary layer-like flow forms near the end walls, and the temperature field is stratified at the core region and has high gradients near the ends.

•~--~"~\~i (a) 10 contours

When the Ra number increases to a certain critical number, secondary cells form within the main cell. Flow visualization has been conducted by Vest and Arpaci [23] and Elder [24] for high-aspect-ratio enclosures, and by Ostrach, Loka and Kumar [25], and Simpkins and Dudderar [26] for low-aspect-ratio ones. In the high-Ar cases, the critical Ra L for the onset of secondary ceils was determined to be 3x105 + 30% for Ar>20. For low-Ar enclosures, Ostrach, et al. [25], observed secondary cells when (Ra At)2 = 1013 and Ar--0.2 while Simpkins and Dudderar proposed the criterion as Ra At3/4 = 6.4x105 + 10%.

(b) 30 contours

(c) 50 contours

Figure 9 Streamline contours for Ra=105, At=l, 0=0 °

Eider [24] states, however, that determining the onset of secondary cells is extremely difficult because the secondary motion near the critical Ra number is very weak and hard to detect. To this end, it is believed that "numerical experiments" can help by capturing these "weak" cells and determining the critical Ra number in a cleaner way. A number of numerical studies have been conducted in the past, and streamline contours have been plotted to show the calculated secondary ceils [27, 28] although no major effort has been reported in terms of defining the critical Ra number for the onset of secondary ceils as a function of Ar and Pr.

(iv)

C~)

01i)

6)

(a) Ra=10 3

The use of a streamline plot to determine the critical Ra, however, has a drawback; it is graphic-resolution dependent. This point can be better illustrated, as shown in Figure 9, by plotting 10, 30, and 50 streamlines equally spaced between •rnax and ~gminfor a Ra=105 and Ar=l case. With 10 contour lines, it appears that only one nicely-shaped cell exists. However, as the number of contour lines increases, a dumbell shaped cell eventually turns into two imbedded secondary ceils. The resolution issue presents some difficulties near the critical point because the secondary motion is so weak that it requires a large number of contour lines to resolve the secondary ceils.

In light of this difficulty, a simple, clean method can be obtained by using velocity maps to determine the onset point of secondary ceils. The map is constructed by plotting zero u and zero v lines; these lines divide the domain into four possible velocity arrangements, i.e., (a) u>o, v>0, (b) u>o, v<0, (c)

D

o~:)

;¢) ~:; i

I

(b) Ra=105

i::) t

(c) Ra=107

Figure 10 Velocity map for Ar=l; where region (i) shows u>0, v>0, region (ii) shows u>0, v<0, region (iii) shows u<0, v>0, and region (iv) shows u<0, v<0

40th IAF Congress

373

u<0, v<0, and (d) u<0, v>0. In Figure I0, the velocitymap is constructodfor Ra=103, 105,and 107,rospoctively.Clearly, only one cellexistsin the enclosureat Ra=103 while secondary cells are formed at Ra=105. When Ra reaches 107, the secondary cells move close~ to the end walls and the flow field gets more complicated.

that the velocity map isa more accurate way to define the onset point. The critical Rayleigh number thus determined is shown to be lower than has been detected by experimental studies. The relationship between the critical Ra and the aspect ratio also is shown to be different from what have been reported previously.

The difficulty encountered in experimentally detea'minm"g the onset point of the secondary cells can also be illustrated by plotting the v-component of velocity across the midheight of the enclosure as shown in Figure 11. At Ra--105, the IvmaxI/lvminl is 30, and this ratio increases significantly as Ra gets closer to Rao which makes it exttmnely difficult to resolve

~clcnowledt,ments

Based on thisstudy,the criticalRayleigh number for the onset of secondary cellsis 1.2x104 with At=l, while for Ar--~.2case,secondary cellscan occur as earlyas Ra=103. The currentfindingshows a dramaticallydifferenttrendfrom what have been suggested by Ostrach et al.[25],and Simpkins and Dudderar [26],where they both predictedhigherRa t for lower As.

[

I"

(a) Ra=103

(b) Ra=105

Figure 11 Distribution of vertical velocity-component along y=0.5.

Conclusion A time-accurate numerical study has been conducted and reported on natural convection in two-dimensionalrectangular enclosures. The effect of grid stretching has been discussed, as well as the grid size and under-relaxation. A wide range of Rayleigh numbers, aspect ratios, and tilt angles have been investigated for a Prandtl numer at 7. Various transport modes have been discussed, including convection, instability and the combined mechanism of these two. On the Rayleigh-Bemard instability, a second transition from stable uni-celiularconvection to oscillatory convection has been found. A simple and straightforward method of using the velocity map is employed to determine the onset of secondary cells in the heating-from-side cases. It is demonstrated in this paper

The author is grateful to his colleague and good friend, Mr. Mark Weislogcl, for his help in editing the graphs and table. The author also thanks Dr. Dave Jacqmin and Mr. Bill Devol for proofreading this paper. The encouragement and guidance receivedfrom Mr. Jack Salznmn, Chief of the Microgravity Science and Technology Branch, during the course of thisstudyis deeply appreciated.

References 1.

Sernas, V., Lee, E. I., "Heat Transfer in Air Enclosures of Aspect Ratio Less than One," J. Heat Transfer, Vol. 103, pp. 617-622, 1981. 2. Turner, J. S., "Double-Diffusive Phenomena, "Ann. Rev. Fluid Mech., Vol. 6, pp.37-56, 1974. 3. Buchberg, H., Catton, I., and Edwards, D. K., "Natural Convection in Enclosed Space - A Review of Application to Solar Energy Collection," J. Heat Transfer, Vol. 98, pp. 182-188, 1976. 4. Ostrach, S., "Fluid Mechanics in Crystal Growth - the 1982 Freeman Scholar l.e,cture," J. Fluids Engng., Vol. 105, pp. 5-20, 1983. 5. Catton, I., "Natural Convection in Enclosures," Heat Transfer 1978, Vol. 6, pp. 13-31, National Research Council of Canada, 1979. 6. Ostrach, S., "Natural Convection in Enclosures," J. Heat Transfer, Vol. 110, pp. 1175-1190, 1988. 7. Gebhart, B., Jaluria, Y., Mahajan, R. L., and Sammakia, B., Buoyancy-Induced Flows And Transport, Hemisphere, 1988. 8. Cormack, D. E., Leal, L. G., Imberger, J., "Natural Convection in a Shallow Cavity with Differentialiy Heated End Walls. Part 1. Asymptotic Theory," J. Fluid Mech., Vol. 65, pp. 209-229, 1974. 9. Tichy, J., Gadgil, A., "High Rayleigh Number Laminar Convection in Low Aspect Ratio Enclosures with AdiabaticHorizontalWalls and DifferentiallyHeated VerticalWalls," J.Heat Transfer,Vol. I04, pp. 103-I 10, 1982. 10. Kassemi, S. A., "High Rayleigh Number Convection in RectangularEnclosureswith DifferentiallyHeated Vertical Walls and Aspect Ratios Between Zero and Unity,"N A S A TM- I00277, 1988. 11. Bejan, A., and "lien,C. L., "Laminar Natural Convection Heat Transfer in a Horizontal Cavity with Different End Temperatures," J. Heat Transfer, Vol. 100, pp. 641-647,

374

J.C. DUH

1978. 12. Kamotoni, Y., Wang, L. W., Ostrach, S., "Experiments on Natural Convection Heat Transfer in Low Aspect Ratio Enclosures," AIAA Journal, Vol. 21, No. 2, pp. 93-97, 1983. 13. Anderson, D. A., Tannehill, J. C., and Pletcher, R. H , Computational Fluid Mechanics and Heat Transfer, Hemisphere, 1984. 14. Torrance, K. E., and Carton, I., Eds., Natural Convection in Enclosures, ASME HTD Symposium, Vol. 8, 1980. 15. Ostrach, S., "Natural Convection Heat Transfer in Cavities and Cells," Proceeding 7th International Heat Transfer Conference, VoL 1, Hemsphere, 1982. 16. Duh, J. C., Wang, L. W., Rashidnia, N., and Tan, C.H., "Natural Convection Inside Inclined Low-AspectRatio Enclosures," 7th International PCH Conference, Cambridge, Massachusetts, June, 1989. 17. Chandrasekhar, S., Hydrodynamic and Hydromagnetic Stability, Dover, 1981. 18. Rayleigh, L., "On Convective Currents in a Horizontal Layer of Liquid when the Higher Temperature is on the Under Side," Phil. Mag., Vol. 32, pp. 529-546, 1916. 19. Catton, I., "The effect of Insulating Vertical Walls on the Onset of Motion in a Fluid Heated from Below," Int. J. Heat Mass Transfer, Vol. 15, pp. 665-672, 1972. 20. Ozoe, H., Yamamoto, K., Churchill, S. W., and Sayama, H., "Three-Dimensional, Numerical Analysis of Laminar Natural Convection in a Confined Fluid Heated from Below," J. Heat Transfer, Vol. 98, pp. 202-207, 1976. 21. Goldhirsch, I., Pelz, R. B., and Orszag, S. A., "Numerical Simulation of Thermal Convection in a TwoDimensional Finite Box," J. Fluid Mech., Vol. 199, pp. 1-28, 1989. 22. Wirtz, R. A., and Tseng, W.-F., "Natural Convection Across Tilted, Rectangular Enclosures of Small Aspect Ratios," in Natural Convection in Enclosures, eds. Torrance and Catton, ASME HTD - Vol. 8, pp. 47-54, 1980. 23. Vest, C. M., and Arpaci, V. S., "Stability of Natural Convection in a Vertical Slot," J. Fluid Mech., Vol. 36, pp. 1-15, 1969. 24. Elder, J. W., "Laminar Free Convection in a Vertical Slot," J. Fluid Mech., Vol. 23, pp. 77-98, 1965. 25. Ostrach, S., Loka, R. R., and Kumar, A., "Natural Convection in Low Aspect Ratio Rectangular Enclosures," in Natural Convection in Enclosures, eds. Torrance and Catton, ASME HTD - Vol. 8, 1980. 26. Simpkins, P. G., Dudderar, T. D., "Convection in Rectangular Cavities with Differentially Heated End Wails," J. Fluid Mech., Vol. 110, pp.433-456, 1981. 27. De Vahl Davis, G., and Mallinson, G. D., "A Note on Natural Convection in a Vertical Slot," J. Fluid Mech., Vol. 72, pp. 87-93, 1975. 28. Chikhaoui, A., Maslanik, M. K., and Sani, R. L., "Steady Three-Dimensional Convection in a Vertical Rectangular Enclosure," The fifth International Conference

on Numerical Methods in Thermal Problems, Montreal, Canada, 1987.