Computational Mechanics (1989) 4, 4 I7M27

Computational Mechanics © Springer-Verlag 1989

High-accuracy bench mark solutions to natural convection in a square cavity T. Saitoh Department of Mechanical Engineering II, Tohoku University, Sendai 980, Japan

K. Hirose Yonezawa Technical High School 9, Yonezawa 992, Japan

Abstract. A fourth-order high-accuracy finite difference method is presented for the bouyancy-driven flow in a square cavity with differentially heated vertical walls. The two bench mark solutions against which other solutions can be compared were obtained. The present solution is seemed to be accurate up to fifth decimal. The proposed scheme is stable and convergent for high Rayleigh number, and will be applicable to general problems involving flow and heat transfer, especially in three dimensions. List of symbols a

g h k L Pr Ra t

Tc, Th T AT bl~ V x,y

zXx, A y

thermal diffusivity gravitational acceleration mesh length in the x-direction time step side length of cavity Prandtl number Rayleigh number = g fl A TL3/v a time surface temperatures (see Fig. t a) temperature temperature difference = Th - TC velocity in the x and y directions, respectively coordinates mesh lengths in the x- and y-directions, respectively

Greek symbols fl v q) co

volumetric expansion coefficient kinematic viscosity stream function vorticity

1 Introduction

Natural convection heat transfer in an enclosed cavity has been regarded as a standard comparison problem for checking the computational methodologies pertinent to general Navier-Stokes equations, The papers by Wilkes and Churchill (1966), Cormack et al. (1974), Saitoh (1977, 1984), Jones (1979), Mallinson and de Vahl Davis (1973), de Vahl Davis (1983) and Markatos and Pericleous (1984) are typical examples of such work. Following the suggestion by Jones, the contributions were reported and compared with a bench mark solution prepared by de Vahl Davis and Jonens (1983). Later, their solutions were refined in greater detail (de Vahl Davis 1983). The "double glazing" problem is one of the typical applications of the above problems. Other include ventilation of green house, convection in the thermal energy storage reservoirs, and convection in various types of solar collectors. Recent literature is reviewed by Ostrach (1982). Effort has

418

Computational Mechanics 4 (1989)

been lately devoted to the high Rayleigh number problem which has a complicated nature compared with that at low Rayleigh numbers. At very high Rayleigh numbers, the entire natural convection flow is divided into two domains, i.e. the boundary layer and the core flow. Besides this, the boundary layer gets very thin as the Rayleigh number increases. Therefore, in order to attain sufficient resolution of the numerical solution, the amount of grids should be greatly augmented. However, this may result in requiring tremendous computing time. The appearance of the supercomputers, such as NEAC SX-2, which has a computing performance of more than 1300 MFLOPS, enables us to solve multi-dimensional problems of various types. Nonetheless, the above example poses a time-consuming problem even if the SX-2 were employed as far as the ordinary computational procedure with low order truncation error is adopted. Another difficulty is the degree of accuracy of the solution obtained. The coarse grid-spacing leads to less accurate results. To resolve the above problems, Markatos and Pericleous (1984) have recently presented an interesting method for laminar and turbulent natural convection in an enclosed cavity. They claim that their method is cost-effective and strongly convergent up to high Rayleigh numbers. They used the laminar formulation in the Rayleigh numbers less than or equal to 106, and the turbulent formulation using a two-equation model for Ra > 106. The solution procedure used is the same proposed by Spalding (1981), in which finite-domain equation are derived by integration of the differential equations over an elementary control volume or cell surrounding a grid node. To enhance numerical stability, use was made of upwind differencing. To detect the "false diffusion" griddependence was examined by changing the grid size. They showed that the computer time to get one set of solution may be one order of magnitude less than that of the existing methods. Streamlines, isotherms, horizontal and vertical velocity contours, vector fields, secondary vortices, and the Nu - Ra correlations were clarified in detail for a Ra number range from 103 to 1016. It seems to the present authors that their method is the only one that makes it possible to solve the natural convection problems at high Ra numbers up to 1016. However, as is well known, the truncation error of the upwind differencing is of an order of O (h), h being the mesh length. Therefore, the degree of accuracy of their solution, especially at high Ra numbers, should be examined with great care. In this paper, the fourth order high accuracy finite difference method is proposed for solving the bouyancy-driven flow in a square cavity with differentially heated vertical walls. Two bench mark solutions for Ra numbers of 104 and 106 against which other solutions can be compared are presented by making the best use of one of the characteristic of the truncation error of the fourth order. Sample solutions were also obtained for high Ra numbers up to Ra = 108. Though the present results are restricted to square cavities of limited aspect ratio, Prandtl number, and Rayleigh numbers, the procedure is general and can be equally applied to various problems including turbulent flow. In particular, this method will be cost-effective for three-dimensional problems, which are getting more and more important in recent days.

2 The comparison problem The problem being considered is the two-dimensional natural convection flow and heat transfer in an upright square cavity of side L, which is the same as that treated everywhere, e.g. de Vahl Davis (1983). The horizontal two walls are insulated and the vertical sides are at constant temperatures Th and To. The Prandtl number is fixed to be 0.71. The problem is then to calculate the flow and thermal field for Rayleigh numbers of 103, 104, 105, and 106. When the Nusselt number Nu along any line parallel to the y-axis is defined as Nux= S u r - ~ 7 x

dy,

o

the following results should be supplied: average Nusselt number Nu = ~Nuxdxlx_Oo~x=l

419

T. Saitoh and K. Hirose: High-accuracy bench mark solutions to natural convection in a square cavity adiabatic /

....

.7-h

a

]'cx

x, ~

b

adiabatic

Fig.1. a A schematic description of the problem considered; b an example of scaled-grid mesh

x

maximum and minimum Nusselt numbers on the hot wall (x = 0), and their location

{

N1Amax = U T-- ~x

maxatx=0

NUmin= {biT-- ~}minatx=O average Nusselt number on the vertical mid-plane of the cavity

NUll2 = { Nux}x= 1/2 average Nusselt number on the vertical boundary of the cavity at x = 0

Nuo = { Nux}x=o maximum vertical velocity on the vertical mid-plane and its location Umax

maximum vertical velocity on the horizontal mid-plane and its location /;max

stream function at the mid-point of the cavity

I~l..d A schematic is shown in Fig. 1 a with a coordinate system.

3

Mathematical

formulation

The stream function and vorticity formulation was adopted. The same dimensionless variables as those employed in the past studies were used, i.e. 1

{x,y}= v{x',y'}, L,

L {u,v}=-{u',v'}, a

L2 ~o =

--

a

11)' ~o',

W -

T ' - Tc T-

a'

.

(1)

Th-T c

The governing equations are then 5oo 5t5T

+ u

&9

5(o 5T + v ~ = Pr v2~o + Pr Ra U 5x 5T

5T

8~- + U ~ x + V 0y

co = ve,p

V2T

(2) (3) (4)

420

Computational Mechanics4 (1989)

where, Laplacian V2 means 02 82 V2_= + -8x 2

Oy2"

In the above equations, the false unsteady terms were added. The solution is then obtained as an asymptotic steady state solution. The boundary conditions are listed below. x = 0;

T = 1,

u = v = 0,

x=l;

T=0,

u=v=0,

y = 0, 1;

6T6y- 0,

&P - 0,

co = ©2q3

ax

Ux l =0

&¢ - 0, ;-ax

u = v = 0,

co = 62~p Sx21x=

~P - &P 6y - 0,

(5)

~21~ co = -~y2 - [y=0,1

4 Solution methodology In order to cope with the high Rayleigh number problem under moderate cost performance and allowable accuracy, two major tools are prepared; i.e. high degree of accuracy finite difference scheme and scaled-grid spacing with new transformation functions. These are indispensable factors when dealing with multi-dimensional, especially 3-dimensional, problems since the computer running time (CPU time) is roughly proportional to h ~+2, n being the number of dimensions. Therefore, to decrease the CPU time, it is most efficient to use a coarse mesh length. However, this will result in considerable degradation of resolution of the numerical results. 4.1 Multi-point finite difference scheme

The problem described in the previous section was solved by using the multi-point explicit finite difference method ( M E F D M ) which was first presented by Saitoh (1974, 1977, 1984). Four multipoint schemes, i.e. multi-point explicit, multi-point implicit, multi-point Adams-Bashforth, and multi-point DuFort-Frankel, have been proposed. The multi-point approximations of the first and second derivatives, for example, are described as follows.

~XXi,j ~ bli--2'J-- 8Ui--l'J @ 8Ui+l'J-- Hi+2'J -~- 4h4 o5bl 12h

5!

6x 5"'"

(6)

06U6 ~X21i,j~2bl --~'~ Ui--a'J+ 16Ui--l'J 12h 230ui'j + 16ui+l'J- ui+2'J -}- --6!]i4 --'"~x --

--

(7)

8

The above central difference was used for the advections terms of Navier-Stokes equations. The truncation error T;,j of M E F D M for the one-dimensional heat conduction equation is given by the next equation.

ri,j -- Hi'j+lk-- bli'j k ~2U

h4

-- bli-2'J "~- 16Ui-l'J --12h a30ui'j + 16Ui+l'J- Ui+2,j ~6U

k 2 ~3U

h6

.(~bl~ ~X2]i,

(8)

~Su

= 2 ~6t + 90 ~6x + 6- ~6t + 1008 6x 8 + ""

(9)

The multi-point finite difference representation for the vorticity equation (2) is written as, ca!Z+1) = (,~(.k)+ At [ -

•i,j

1

-- < , j ~ t ~'y + Pr {A 2

+ A2

} + P r R a A 1T~)].

(10)

T. Saitoh and K. Hirose: High-accuracy bench mark solutions to natural convection in a square cavity

421

Here, AI(_D(xk) = O i - 2,J -- 8(.0i_1, j q- 8 O i + 1 , j -

0i+2, j

12Ax A109~ k) 7- (Di, J - 2 -- 8(.Oi, j _ 1 °c 8(Di, j + 1 -- (Di,j+ 2

12Ay A2co~) _- - Col_z, j + 16coi_l,j - 30coi, ) + 16o)i+1,j - coi+2, j 12(Ax) 2 A2 (D(yk).~-- -- 03i,j_ 2 q- 1 6 c o i , j _ 1 -- 30coi, j +

16coi, j+1

- - (Di,j+ 2

12(Ay)2 Another fourth order method for two-dimensional problems is presented recently by Gupta et al. (1983). The stability range of M E F D M is k 3 _ < _ (11) h 2 = 8' The stability analysis for the M E F D M is shown in Appendix. It is important that all boundary conditions are expressed by the fully fourth order equations. For example, the fully fourth order expressions for the typical boundary conditions are shown below. (i) Adiabatic condition: bli,j_ 1 ~ b/i,j+l;

bli,j--2 ~ bli,j+ 2

(ii) Constant temperature conditon (ui = C): ui=

C;

ui_ 1 = 2C-

ui+ l

(iii) Vorticity at the wall in Poisson's equation: coi=V2,#

when

1{1

'Pi=~xx = 0 . &P i

coi,j - (Ax) 2 ~ i + l , j -

s

3~+2,7+ ~ i + 3 , j -

,

gt&+4,j

}

O.)i_l, j = 6~0i, j -- 15C0i+1,j + 20 COi+2,j -- 15C0i+3, j + 6 ~ 0 i + 4 , j -

03i+5, j

4.2 Scaled-grid spacing

In order to improve the resolution in the boundary layer settled near the walls, finer mesh is needed. For instance, Markatos and Pericleous (1984) have used at least 15 grid points in the boundary layer. However, if the usual uniform mesh were employed, the total number of grid points amounts to a formidable quantity, e.g. in their paper, at most a 160 x 100 grid arrangement was used. At present, it seems that mesh division more than 100 in the one-direction is realistic only for two dimensions. For example, 100 x 100 x 100 grids are almost prohibitive in 3-dimensional problems. The idea of the scaled-grid spacing (SGS) was first introduced by Newell and Schmidt (1970) and improved later by Cormack et al. (1974). Since it seems that the previous functions for transformation including Newell-Schmidt polynomials are not suitable for the present purpose, a new transformation was introduced. One of our transformations is given by the next equation: X-

1 2 ( n - 1)

((c - 1) (2x - 1)n + (n - c) (2x - 1) + n - 1)

(12)

422

Computational Mechanics 4 (1989)

Here, X means the transformed coordinate, and c and n designate gradient at x -- 0 and the degree of the polynomial, respectively. The SGS is particularly effective for high Rayleigh numbers in which the boundary layer gets very thin, thereby requiring high resolution. If R a is 108, the thickness of the boundary layer is estimated to be 0.05 (i.e. 5% of side length). The SGS with the number of division M -- 40, c = 5, and n = 7 may generate at least 7 grid points within the boundary layer, which gurantees moderate accuracy. In parenthesis, this is equivalent to a 200 x 200 grids arrangement for uniform mesh. A sample of the scaled-grid spacing using the above equation is shown in Fig. 1 b.

5 The bench mark solution

De Vahl Davis (1983) has already presented four sets of bench mark solutions for Rayleigh numbers ranging from 103 to 106. He at first obtained the original solutions with 20 x 20, 40 x 40, 60 x 60, and 80 x 80 uniform grids. Then, bench mark solutions X t were extrapolated from the two original solutions by using the following equation: X 1 -- (hl/h2)nX2

Xt =

1 - (hi~h2) n

(13)

Here, X1 and X2 are the computed original values at h = hi and h = h2, respectively. The order of discretization error n is assumed to be constant and equal to 2. From his results for R a = 1 0 4, the discretization error ~ is represented by ~- 40h 2 •

(14)

Here, e has been defined as e -- z7- u

(15)

with fi and u being the numerical and exact values, respectively. From the above equation, in order to get solutions which are accurate up to four significant figures, at least 630 x 630 grids are required. This implies that the computer time will be some 3500 times as large as that for 80 x 80 grids, which is almost impossible judging from the cost performance of the present computers. For this reason, a high-order method is most suitable for this purpose. The present original fourth order results for R a = 1 0 4 with uniform 30 x 30 and 50 x 50 grids are presented in Table 1. The problem solved and the definitions of Nuo, Nul/2, Numin, Numax, ... are the same as those used by de Vahl Davis. Table 1 also includes extrapolated values by using Eq. (13) with n = 4. These values are to be called the new "bench mark solution". A full run at R a = 1 0 4 using a 50 x 50 grids would require 7200 sweeps and take some 4500 seconds on N E A C ACOS 1000 at Tohoku University. The time per grid point per sweep per variable was 5 x 10- 5 sCPU. The computing speed of ACOS 1000 is 7-14 MFLOPS. As seen from a comparison between the original and bench mark solutions, it may be confirmed that the degree of accuracy of original solutions (e.g. ~Pmid)is approximately 0.002%. The results obtained originally are accurate to four significant figures. Umax, Vmax, Numax, ... are extrapolated by a similar procedure as that described by de Vahl Davis and others. Although we believe that our bench mark solution is reliable, we recommend the reader to compare his results with ~Pmidand Numi n. The next recommendable values are Umax, Vmax, and Numax. Since an integration procedure is involved to obtain N u values, the results obtained may include some errors. Therefore, it is not guaranteed that these values are accurate to the last decimal. Table 1 also contains our original solutions for R a = l06 which were obtained with the scaledgrid spacing (40 x 30). The bench mark solution was not computed in this case because the results were obtained by virtue of the scaled-grid spacing. Figures 2 and 3 show the streamlines, isotherms, and vorticity contours for the original solutions at R a - - 1 0 4 and 106, respectively. General trends are the same as shown in the de Vahl Davis' (1983) paper.

423

T. Saitoh and K. Hirose: High-accuracy bench mark solutions to natural convection in a square cavity Table 1. The original and bench mark solutions; the results by de Vahl Davis are also included

Ra

=

Ra

10 4

=

10 6

Bench Mark Solution 30 x 30

50 x 50

40 x 30

de Vahl Davis' bins

50 x 50

16.379 64.886 0.8505 220.47 0.03783 8.7956 8.7989 8.8487 17.140 0.0473 1.015 1

16.32 64.63 0.850 219.36 0.0379 8.800 8.799 8.817 17.925 0.0378 0.989 1

16.245 64.389 0.8512 216.76 0.03943 8.7126 8.7218 8.4843 17.451 0.0392 1.011 1

Present de Vahl Davis tPmid Uma* y Vmax X Nu NUll2 Nuo Numax y Numin y

5.0722 16.1819 0.8232 19.5906 0.1196 2.2356 2.2377 2.2227 3.4791 0.1497 0.5893 1

J

1

L

I

5.0730 16.1838 0.8232 19.6165 0.1191 2.2415 2.2423 2.2385 3.5158 0.1465 0.5864 1

~

I

,

I

5.0731 16.1841 0.8232 19.6204 0.1190 2.2424 2.2430 2.2409 3.5213 0.1460 0.5860 1

,

-,

0

0.20

0.40

0.60

0.80

a

5.071 16.178 0.823 19.617 0.119 2.243 2.243 2.238 3.528 0.143 0.586 1

1.00 0

.

,

0.20

•

0.40

0.60

0.80

1.00 0

b

0.20

0.40

0.60

0.80

1.00

0.20

0.40

0.60

0.80

1.00

e

2

1.00

0.80

0.60

0.40

0.20

0 0 a

0.20

O.LO

0.60

0.80

1,00 0

0.20

0./-,0

0.60

,~ 1.00 0

0.80

b

e

3

Figs. 2 and 3. a Streamlines, b isotherms, and e vorticity contours for 2 Ra

=

1 0 4,

and 3 Ra

=

10 6

Figure 4 designates the velocity profile for Ra = 106. For a comparison, the scaled mesh results (30 x 20) is also shown. The agreement between the two is excellent. The minimum mesh length for SGS of 40 x 30 is 0.00867 and this is equivalent to uniform 117 x 117. Figure 5 shows the vertical and horizontal velocity profiles at Ra = 106.

424

Computational Mechanics 4 (1989)

20

Pr =0.71 Ra =10 4

.l \

l

• Equal mesh (50x50) o Scaled grid ( 30 x 20)

:

tO Pr = 0.71 Ra = 104

~'10

VaN Davis' result

/

n=2

]~mid I

o

0

0.1

0.2

0.3

0.4

0.5

X

4

U

-1

0

1

xlO2

o e3

Preseot result 0.1

g

xlO 2 ,

U

3 5.

=

c

2

Pr = 0.71 Ra = 10 e

$ G_

1

lo

0.01 ,

'~n = 4

/

-1 -2

0.001 0.01 5

6

I

I

i

I

I

I

h

•

I

,

i ,

1.0

0.1 "

Figs. 4 - 6 . 4 Comparison of scaled-grid and equal grid results for vertical volocity profile. 5 Vertical and horizontal velocities for Ra = 106. 6 Rates of convergence of mid-point streamline value; also shown is the result by de Vahl Davis (1983)

Figure 6 shows a comparison of the rate of convergence n with that of de Vahl Davis. The n = 4 dependence of our method in the finer mesh region is affirmed. The ordinates shows the discretization error in percentage. Figures 7 and 8 indicate the results for high Rayleigh numbers, i.e. Ra = 1 0 7 and 108. These are only of reference value because some researchers have indicated that the flow will become turbulent when Ra exceeds 106. However, in the experiment performed by Torrance et al. (1969), the flow was laminar up to Ra = 109, though different geometry and boundary conditions were employed in the experiment. Judging from the analogy with natural convection around a vertical plate in an infinite space, it seems to the present authors that the flow is laminar at least up to Ra = 108. It is not our purpose to examine at which Ra number the flow transition may occur but to show that our method is stable and convergent for high Rayleigh numbers. It is noted here that stability of a numeric scheme is independent of physical stability of the convection flow as far as the laminar equations are concerned. The high-order method which combines the scaled-grid spacing technique described in this paper is especially efficient for applications to three-dimensional problems in which the reduction of computer running time is of primary importance. The SGS is found to be very efficacious against reduction of grid points. However, as an adverse side effect, the mesh arrangement in the central region gets coarse, thereby introducing large discretization errors when use is made of usual second order methods. In this sence, a combination of fourth-order method with SGS brings a optimum result.

T. Saitoh and K. Hirose: High-accuracy bench mark solutions to natural convection in a square cavity 1.00

1.00.

0.80

0.80.

0.60

0.60-

0.40

0.40-

425

0.20

1.00t/~ ,

I

1.00 I

]

I

,

I

,

I

,

I

,

Hill--

0.80

0.60-

0.60

0.40-

0.40

o2o-

o2o

00

0 0.20

0.40

0.60

7

0.80

1.00

0

,

8

I

0.20

,

I

0.40

,

1

0.60

,

0.80

1.00

Figs. 7 and 8. Streamlines and isotherms for 7 Ra = 107, and 8 Ra = 108

6 Concluding remarks

A new bench mark solution was presented for the two-dimensional bouyancy driven flow of air with Prandtl number 0.71 in an upright square cavity. The following conclusions may be drawn from the present study. (i) The fourth-order finite difference method supported by the scaled-grid spacing technique is suitable means to obtain the bench mark solution against which other solutions can be compared. (ii) The present method utilizes the central difference representation for advection terms and is stable up to high Rayleigh numbers (at least Ra = 108). (iii) The proposed method will be particularly efficient for 3-dimensional problems in which the number of grid points becomes a formidable amount. In closing, it seems very interesting to set a 3-dimensional comparison problem just like the one described in this paper. Appendix: Stability of solution

Consider the stability of a finite difference solution of the one-dimensional energy equation with a convection term: 8u 8u ~ -~- /) ~X

82u OX2 -- 0

(A1)

426

Computational Mechanics 4 (1989)

As is well known, upwind differences have been used when one solves the above equation by the explicit scheme. The discretization error then reduces to the order of O (h). The multi-point explicit finite difference expression is Ui,j+ I =

n=i--2

as,,, Unj

(A2)

Here, r

as, i_2 = - ]~(1 + vh), 2Y

as, i-1 = ~ - (2 + vh), 5

(A3)

as, i = 1 -- ~ r, 2F

as, i+ 1 = - f (2 -- v h ) , p

as'i+2=-- -- 1 2 ( 1 - v h ) '

In a matrix representation, Eq. (A1) becomes [K] ( b/j} = £b/j+l}.

(A4)

If the eigenvalue of [K] is represented by )~i,the following equation can be obtained by using Brauer's theorem. I ) ~ i - as, sl <--[as, i-21 + las, i-~ l + 0 + [as, i+ll + [as, i+21 = Ps

(A5)

For the present FDM, it becomes ~i--1+

5

< r (17+91vhl)

2 r =6

(A6)

from which the stability range satisfying 12il < 1 is as follows 12 r _< . - 32+3[vh[

(A7)

It is seen that the M E F D M is stable even for problems involving convection terms. In the special case of v = 0 (i. e., heat conduction alone), stability range considerations yield r =< 3/8.

Acknowledgements The authors extend their sincere thanks to the Computer Center of Tohoku University for use of NEAC 2200 ACOS 1000 and Time Sharing System. Thanks are also due to Professor H. Ozoe of Kyushu University for valuable discussions on the flow transition near Ra = 108.

References Cormack, D.E.; Leal, L.G.; Seinfeld, J.H. (1974): Natural convection in a shallow cavity with differentially heated end wails. 2. Numerical solutions. J. Fluid Mech. 65, 231-246 Gupta, M.M.; Manohar, R.P.; Stephenson, J.W. (1983): A fourth-order, cost effective and stable finite difference scheme for the convection-diffusion equation. Proc. 2nd Nat. Symp., pp. 201-209. Washington/DC: Hemisphere

T. Saitoh and K. Hirose: High-accuracy bench mark solutions to natural convection in a square cavity

427

Jones, I.P. (1979): A comparison problem for numerical methods in fluid dynamics: the double glazing problem. In: Lewis, R.W.; Morgan, K. (eds.): Numerical methods in thermal problems, pp. 338-348. Swansea/UK: Pineridge Press Mallinson, G. D.; de Vahl Davis, G. (1973): The method for the false transient for the solution of coupled partial differential equations. J. Comput. Phys. 12, 435~445 Markatos, N. C.; Pericleous, K.A. (1984): Laminar and turbulent natural convection in an square cavity. Int. J. Heat Mass Transfer 27, 755-772 Newell, M.E.; Schmidt, F.W. (1970): Heat transfer by laminar natural convection within rectangular enclosures. J. Heat Transfer 92, 159-168 Ostrach, S. (1982): Natural convection Heat Transfer in cavities and cells. Proc. Int. Heat Transfer Conf., pp. 365-379. Washington/DC: Hemisphere Saitoh, T. (1977): A numerical method for two-dimensional Navier-Stokes equation by multi-point finite differences. Int. J. Numer. Methods Eng. 11, 1439-1454 Saitoh, T. (1984): Numerical analysis of multi-dimensional heat transfer problems via multi-point finite differences. Trans. JSME 50, 2704-2708 Spalding, D. B. (1981): A general purpose computer program for multi-dimensional one- and two-phase flow. Math. Comput. Simul. 23, Amsterdam: North Holland 267576 Saitoh, T. (1974): A numerical method for transient heat conduction equation by 7-point finite differences. Trans. JSME 40, 2617-2631 Saitoh, T. (1984): Numerical methods in heat transfer. Sci. Mach. 36, 1051-1056 Torrance, K.E.; Orloff, L.; Rockett, J.A. (1969): Experiments on natural convection in enclosures with localized heating from below, J. Fluid Mech. 36, 21-31 Vahl Davis, G. de; Jones, I.P. (1983): Natural convection in a square cavity: a comparison exercise. Int. J. Numer. Methods Fluids 3,227-248 Vahl Davis, G. de (1983): Natural convection of air in a square cavity: A bench mark numerical solution. Int. J. Numer. Methods Fluids 3,249-264 Wilkes, J.O.; Churchill, S.W. (1966): The finite-difference computation of natural convection in a rectangular enclosure. AIChE J. 12, 161-171

Communicated by G. Yagawa, August 18, 1987