A1110E I'mDM

NBSIR 84-2830

Numerical Solutions for Steady

Natural Convection in A Square

Cavity

U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

Center for Chemical Engineering

Chemical Process Metrology Division

Washington, DC 20234

March 1 984

I

! [ \ i o. [ \ J

«(,•'’fAU

U.S DEPARTMENT OF COMMERCE

NATIONAL BUREAU OF STANDARDS

j*QC

I 100

' .1156 84-2330

1934

(RATIONAL BtTTtrAtf

' OF ST.M’DF.ftwS

1/ Lire

NBSIR 84-2830

NUMERICAL SOLUTIONS FOR STEADY

NATURAL CONVECTION IN A SQUARE

CAVITY

E. F. Moore, R. W. Davis

U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

Center for Chemical Engineering

Chemical Process Metrology Division

Washington, DC 20234

March 1984

U.S. DEPARTMENT OF COMMERCE, Malcolm Baldrige, Secretary

NATIONAL BUREAU OF STANDARDS. Ernest Ambler. Director
Ii
NUMERICAL SOLUTIONS FOR STEADY NATURAL CONVECTION

IN A SQUARE CAVITY

E. F. Moore and R. W. Davis

Chemical Process Metrology Division

Center for Chemical Engineering

National Engineering Laboratory

National Bureau of Standards

Washington, D.C. 20234

Abstract

Numerical solutions have been obtained for steady natural convection in a square cavity. The numerical method used was developed for

unsteady, incompressible, viscous fluid flow. The similarity

parameters were chosen to match those of an international comparison

exercise. Results are presented and compared with those obtained by

other researchers using different methods.

Key Words: Cavity flow; fluid dynamics; natural convection; numerical

methods

.

Introducti on

The idea of carrying out a comparison exercise for numerical

schemes based on the problem of steady natural convection in a square

cavity with vertical sides which are maintained at different

temperatures was proposed in 1979 [1]. The problem is of interest due

to the range of detailed flow patterns found in this simple geometry.

The applicable equations here are the Navier-Stokes and continuity

equations of fluid flow. The Navier-Stokes equations include a body

force term, which requires the solution of an additional equation for

temperature. These equations are all simplified by the Boussinesq

approximation. There are no singularities in the boundary conditions.

.
as occur in the driven cavity problem. The results of this study could

be used as a model problem solution for checking computer codes

developed for practical elliptic fluid flow computations.

Two recent papers summarize the contributions received from

participating researchers [2] and a bench mark solution [3] with which

their results were compared. This report describes a solution obtained

at NBS using a numerical scheme developed for unsteady, incompressible,

viscous fluid flow problems.

In the following section, the comparison problem and similarity

parameters will be presented. The solution procedure will then be

described. Results and comparisons with the bench mark solution and

those of other contributors are also included.

The Comparison Problem

The basic problem considered here is the steady two-dimensional

flow of a Boussinesq fluid in an upright square cavity. The details

are shown in figure 1. Each side of the cavity has length D, and both

velocity components are zero on the walls. The vertical sides are at

temperatures T-| and J

2

> with the normal temperature gradient being

zero on the top and bottom walls. The fluid has thermal diffusivity <

and kinematic viscosity v, with Prandtl number Pr =^. The Rayleigh

3 < number for this flow is defined as Ra = where 3 is the coeffi-

<v

cient of expansion of the fluid, g is gravitational acceleration, and

At = Ti T£- In this study, all lengths are nondimensional ized

K

with respect to D; all velocities with respect to jj ; time with respect

D 2

to — ; and p, the ratio of pressure to reference density, with respect

*2 T

' - T

2

to ^75-. The nondimensional temperature is T = f f— , where

D

2 * 1 "

‘

2
T' is the dimensional value. Therefore, 0 <_ x <_ 1 , 0 <_ z _< 1 , T = 1 at

x = 0 and T = 0 at x = 1

.

The velocity and temperature fields are to be computed with

Prandtl number 0.71 for Rayleigh numbers of 10 3

, 10 4, 10^, and

10^. The comparison also includes Nusselt numbers, or nondimensional

values of heat flux in the horizontal direction.

Numerical Modeling

The two-dimensional Navier-Stokes , continuity, and temperature

equations for a viscous Boussinesq fluid are

3u , 9 (u

2

) , 9 (uw) 9 p 9t 9x

+ =

' 37

9 w +

3 (uw)

+

3 (w

2

) = - 3R

3 1 3 x 3 z 9 z

9U jw . n 3X

+ 37' °

9T

t 3 ( uT) 3 (wT) _ 9 2T

3 1 3 x 3z

9 x

2

a 2 a 2

„ ,o u a Ux Pr ( —£ + —=)

3 2

T

d Z

3 X

: 5 2

v

„ t oX

3 z

33z‘

where u and w are velocity components in the x- and z- directions,

respecti vely.

The basic numerical scheme employed here uses quadratic upwind

differencing for convection and an explicit Leith-type of temporal

differencing [4]. A fast direct method is used to solve the Poisson

equation for pressure at each time step. The solution is accomplished

on a staggered mesh in which pressures and temperatures are defined at

cell centers and normal velocities at cell faces.

3
The boundary conditions on the four cavity walls are that normal

and tangential velocities are zero. The temperatures on the vertical

walls are specified. Any values for grid points just outside the walls

are obtained by quadratic extrapolation using the wall value and two

points just inside. The extrapolated values are used in computing

3T

derivatives at the wall. At the top and bottom walls,^- = 0, so the

points just outside those walls have the same temperatures as those

just inside the walls.

The only mesh employed in this study was a 50 x 50 uniform grid

with Ax = Az = 0.02. The initial conditions for the computations were

either u = 0, w = 0, T = 1-x for Ra = 10 3

, or the results of a

previous computation at a lower Rayleigh number. The time step, At,

was set at a maximum of 7.5 • 10~ 3 so that the diffusion

coefficient in the temperature equation was less than 0.5 and the

Courant number less than one. Computation times on the MBS UNIVAC

1100/82 required to obtain steady-state solutions ranged between 5 and

7 hours. At steady-state, the averages of the absolute values of the

time-derivatives of the velocities [-—

^

and 3— , where n

indicates time level] throughout the cavity were 0(10“ 3

) for the

lower three Rayleigh numbers and 0(10""*) for Ra = 10°. The

maximum change in u or w during the final time step was 0(10“^) for

Ra = 10° and 0(10~ 3

) for the lowest Rayleigh numbers. The

average of the absolute values of the time-derivatives of temperature

was 0(10" 3

) or less. The summation of the absolute values of the

mass flux residuals was 0(10

-3

) for Ra = 10 3 and increased by an

order of magnitude for each tenfold increase in Ra.

4
Results and Discussion

Computations have been carried out for the four Rayleigh numbers

in the comparison exercise. The results are shown in table I. The

quantities listed are those requested for comparison purposes:

u max » z - the maximum u-component of velocity along the line

x = 0.5 and its location;

wmax»

x ” the max imum w-component of velocity along the line

z = 0.5 and its location;

Nu 0- the average Mussel t number on the wall of the cavity at

x = 0;

Mumax , z " the maximum value of the local Mussel t number on

the wall at x = 0 and its location;

Mumi

- n

, z - the minimum value of the local Mussel t number on

the wall at x = 0 and its location;

Mu - the average Mussel t number throughout the cavity;

Mui /2 " the average Mussel t number along the line x = 0.5.

The values for umax , z, wmax , and x are grid point

3T values (no interpolation). The Musselt number, Mu(x,z) = uT was

calculated at each of the mesh points. The quantities Munax

and

Mumi

- n and their locations are values of Nu(0, z) at the z grid

points. The average Nusselt number along a line x = constant, such as

fl Mu or Mu, /? , was obtained by integration: Mu = J Mu(x,z) dz.

o M c x - o The integrals were evaluated by the rectangular rule. Average values

obtained at each value of x were integrated by the trapezoidal rule to

find Mu =

j Mu x dx.

Contour plots of temperature T, velocity components u and w, and

vorticity are presented in figures 2-5 for each of the four Rayleigh

numbers. These illustrate the diagonal symmetry of the problem:

5
T( x , z) = 1 - T( 1 -x, 1-z), u(x,z) = -u(l-x, 1-z), w(x,z) = -w(l-x, 1-z).

Also, the development of the thermal boundary layer as Ra increases is

clearly seen. The same qualitative features appear in the bench mark

solution (figures 4-7 in reference 3).

Figure 6 presents velocity vector plots for the four cases, where

the line segments indicate the local flow direction. They are

analogous to the streamline plots in figure 3 of reference 3.

Secondary recirculating rolls appear for Ra j> 10^ as observed

prevdoti-s-ly [1,5].

The quantitative results presented here in table I can be compared

with those of the bench mark solution given in table I of reference 2.

The percentage differences (with respect to the bench mark solution)

are shown here in table II. The field locations where the parameters

of table II occur have been omitted from the comparisons as was done in

reference 2. Also, stream function was not computed in the present

work. The percentage differences are seen to be small, within about 1%

for the lower three Rayleigh numbers.

It is interesting to compare these results with those in tables

V 1 1 - X of reference 2 which contain the contributed results. For

purposes of comparison, the absolute values of the percentage errors in

Mu, Nu , Mu . , u , and w were averaged for each solution at * max’ min’ max ’ max 3

each Rayleigh number. Mu is either Mu or Mu-]/2» if available,

and has been compared with Mu-j/2 1,1 the bench ^ark solution. The

errors at each Rayleigh number were then averaged to give an overall

average error for each contributed solution.

The average errors in the present work ranged from 0.1% at Ra =

10 3 to 2.1% at Ra = 10^. The error averaged over the four Rayleigh

numbers is 0.3%. Mot all contributors computed all quantities or

submitted solutions for all Rayleigh numbers. Of the 26 solutions from

a
reference 2 that had all quantities, four had a smaller average error

than the present work, placing these results in the top 20% of the

solutions.

Conclusion

Solutions have been presented for laminar natural convection in a

square cavity. The computer code was adapted from one used for

unsteady, incompressible, viscous fluid flow problems. The results

obtained here compare well with the bench mark solution given in

reference 3. The present results also compare well with the

contributed solutions in reference 2. Although the errors in the

present solution increased with increasing Rayleigh number, it should

be noted that all four solutions were obtained using the same grid.

Mesh refinement would have been preferable as Ra increased but was not

possible due to limitations on computer resources.

References

1. Jones, I. P., A comparison problem for numerical methods in fluid

dynamics: the "double-glazing" problem, Numerical Methods in

Thermal Problems , ed. by R. W. Lewis and K. Morgan, Pineridge

Press, Swansea, p. 338, 1979.

2. De Vahl Davis, G. and Jones, I. P., Natural convection in a square

cavity: a comparison exercise, Int. J. Num. Meth. in Fluids , Vol

.

3,

p. 227, 1983.

3. De Vahl Davis, G., Natural convection of air in a square cavity: a

bench mark numerical solution, Int. J. Num. Meth. in Fluids , Vol.

3, p. 249, 1983.

7
4. Davis, R. W. and Moore, E. F., A numerical study of vortex shedding

from rectangles

5. Mallinson, G. D

convection in a

p. 1, 1977.

J. Fluid Mech. , Vol . 116, p. 475, 1982.

and de Vahl Davis, G., Three-dimensional

box: a numerical study, J . Fluid Mech.

natural

Vol. 83,

8
Table I. Results

Ra

ioi IQ

4

10

5

ioi

u max

3.646 16.126 34.62 64.86

z 0.81 0.83 0.85 0.85

wmax3.689 19.507 67.98 211.86

X 0.17 0.11 0.07 0.03

Nu

0

1.118 2.244 4.526 8.902

%ax 1.507 3.538 7.812 18.809

Z 0.09 0.15 0.07 0.03

Nu . min 0.691 0.586 0.732 0.997

z 0.99 0.99 0.99 0.99

ni1.118 2.244 4.523 8.893

Mu

1/2

1.118 2.244 4.524 8.898

Table II. Percentage Differences

Ra

10

3

IQ

4 5 10°

Nu

0

0.1 0.3 0.4 1.0

NlT 0.0 0.0 0.1 1.1

Nu

1/2

0.0 0.0 0.1 i.i

Numax

0.1 0.3 1.2 4.9

Nu . mi n -0.1 0.0 0.4 0.8

umax

-0.1 -0.3 -0.3 0.4

wmax

-0.2 -0.6 -0.9 -3.4
< D >

Fig. 1
RR= 1

000

TEMPERRTURE

Fig•

2a
Rfl=10.000

TEMPERATURE

Fig*

2b
Rfl=

100.000

TEMPERATURE

Fig*

Sc
RA1.000.

000

TEMPERATURE

Fig•

Ed
Rfl=1000

U

VELOCITY o

Fig*

3a
RA=1

0.000

•U

VELOCITY o

Fig*

3b
RA=10Q.000

U

VELOCITY o

Fig*

3c
RR»

1.000.

000

U

VELOCITY o

Fig*

3d
Rfl=1000

h

VELOCITY

Fig.

4a
RA=10.000 \A

VELOCITY

CD to -T r\l

a

C c 2 O 6 a i

|j I

Fig.

4b
Rfl= 1

00-

000

W

VELOCITY Ii !

Fig.

4c
Rfl=

1.000.

000

M

VELOCITY | |

Fig.

Ad
Rfl=1000

VORTICITY

o

o o Z o o o

Fig*

5a
RR=10.000

VORTICITY

o

o d Z d d d

Fig*

5b
RR=100.

000

VORTICITY o

Fig*

5c
Rfi=

1.000.

000

VORTICITY o

Fig*

5d
Rfl-1000

o

Fig*

6a
Rn10,000

o

Fig*

6b
Rfl-100,000

o

Fig.

6c
RH-i,

000,000

Fig*

6d
NBS*1 14A REV. 2*6C ) U.S. DEPT. OF COMM.

BIBLIOGRAPHIC DATA

SHEET (See instructions)

1. PUBLICATION OR

REPORT NO.

84-2830

2. Performing Organ. Report No. 3. Publication Date

March 1984

4. TITLE AND SUBTITLE

Numerical Solutions for Steady Natural Convection in a Square Cavity

5. AUTHOR(S)

E. F. Moore and R. W. Davis

6. PERFORMING ORGANIZATION (If joint or other than N8S. see instructions) 7. Con tract/ Grant No. national bureau of standards

DEPARTMENT OF COMMERCE

WASHINGTON, D.C. 20234

8.

Type of Report & Period Covered

Final

9.

SPONSORING ORGANIZATION NAME AND COMPLETE ADDRESS (Street. City. State. ZIP)

10.

SUPPLEMENTARY NOTES

Document describes a computer program; SF-185, FlPS Software Summary, is attached.

11.

ABSTRACT (A 200-word or less factual summary of most si gnificant information. If document induces a significant

bibliography or literature survey, mention it here)

Numerical solutions have been obtained for steady natural convection in a square

cavity. The numerical method used was developed for unsteady, incompressible,

viscous fluid flow. The similarity parameters were chosen to match those of an international comparison exercise. Results are presented and compared with those

obtained by other researchers using different methods.

12.

KEY WORDS (Six to twelve entries; alphabet: cal order; caoitali ze only prooer names; ana seoarate <ey words zy semicolons

cavity flow; fluid dynamics; natural convection; numerical methods

13. AVAILABILITY

| X i Unlimited

j

! For Official Distribution. Do Not Release to NTIS

Order From Superintendent of Documents, U.S. Government Pointing Office, Washington, D.C.

20402.

~X~ Order From National Technical Information Service (NTIS), Soringfield, VA. 22161

14. NO. OF

PRINTED PAGES

OO

15. Price

58 . 50

jscovm-ac so43-s>so