A11107

NBSIR 84-2932

Finite Difference Calculations of

Buoyant Convection in an

Enclosure

Part II: Verification of the

Nonlinear Algorithm

Ronald G. Rehm* Howard Baum**

P. Darcy Barnett* and Daniel M. Corley*

* Center for Applied Mathematics

**Center for Fire Research

* U S. DEPARTMENT OF COMMERCE

National Bureau of Standards

National Engineering Laboratory

Center for Applied Mathematics

Gaithersburg, MD 20899

**U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

National Engineering Laboratory

Center for Fire Research

Gaithersburg, MD 20899

September 1984

Final Report

* QC

100

,U56

34-2932

1934

U S. DEPARTMENT OF COMMERCE

National Bureau of Standards

National Engineering Laboratory

Center for Applied Mathematics

Gaithersburg, MD 20899

NBSIR 84-2932

OF £T

•i'Y

FINITE DIFFERENCE CALCULATIONS OF

BUOYANT CONVECTION IN AN

ENCLOSURE

PART II: VERIFICATION OF THE

NONLINEAR ALGORITHM

Ronald G. Rehm*, Howard Baum**

P. Darcy Barnett* and Daniel M. Corley*

* Center for Applied Mathematics

**Center for Fire Research

* U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

National Engineering Laboratory

Center for Applied Mathematics

Gaithersburg, MD 20899

**U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

National Engineering Laboratory

Center for Fire Research

Gaithersburg, MD 20899

September 1984

Final Report

U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

National Engineering Laboratory

Center for Applied Mathematics

Gaithersburg, MD 20899

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

NATIONAL BUREAU OF STANDARDS. Ernest Ambler, Director

Finite Difference Calculations of Buoyant

Convection in an Enclosure

Part II: Verification of the Nonlinear

A1 gorithm

Ronald G. Rehm*, Howard R. Baumt,

P. Darcy Barnett* and Daniel M. Corleyt

National Bureau of Standards

Washington, D. C. 20234

* Center for Applied Mathematics

t Center for Fire Research
/
ABSTRACT

Earlier, a novel mathematical model of buoyant convection in an enclosure was

developed. The nonlinear equations constituting this model have recently been

solved by finite difference methods in two dimensions.

In this paper two solutions, obtained in special cases, to the model equations

are presented. For both cases the solutions to the partial differential equations and to the finite difference equations used to approximate the differential equations are obtained by combinations of analytical and numerical techniques. Agreement between the exact solutions to the difference equations

described in this paper and the independently obtained numerical solutions

was found nearly to the accuracy specified (usually 10“ 6

) for an iterative

procedure used in the computational scheme.

The first solution is for a time-dependent irrotational , incompressible flow in

an enclosure driven by sources and sinks of fluid as specified by the heat

source. This problem arises from the full nonlinear equations with boundary

conditions, in continuous or discrete form, by requiring that the velocity

field be irrotational and the density remain constant.

The second set of solutions arises when several other simplifications are made

to the equations. The density is taken to be constant, the heating is assumed

to be zero, the velocity field is taken to be two dimensional and derivable

from a stream function only, the vorticity is taken to be a constant, and the

flow is independent of time. These solutions are used to determine the accuracy with which the code described in Reference 2 solves the nonlinear finite

difference equations in special cases.
- J

• ?
1. Introduction

Over the past few years, the National Bureau of Standards has been engaged in a

research project to develop, starting from basic conservation laws, a mathematical model of fire development within a room. Large scale convection is an

essential component of such a model because this fluid motion governs the smoke

and hot gas transport within a room and also supplies fresh oxygen to the fuel

to sustain combustion. Therefore, development of a mathematical model of

buoyant convection was begun as a first step toward a more complete room-fire

model, which would include effects of combustion chemistry, radiation and smoke

dynamics. The mathematical model for convection, the partial differential

equations and boundary conditions, are derived in Reference 1.

Because fluids admit a rich variety of phenomena and because the equations of

fluid dynamics are nonlinear, it is difficult to obtain solutions to these

equations except in very special cases. It is the nonlinear nature of the

equations of fluid dynamics that makes them difficult to solve: analytical

tools for the solution of nonlinear problems are very limited. For this reason

numerical techniques have become very widely used; numerical solution of nonlinear equations is the only systematic method of solving them under a wide

variety of conditions. Yet, if these numerical solutions are to be trusted, it

is essential that they be carefully checked to determine their accuracy.

It is the purpose of this paper to present two solutions, obtained in special

cases, to the equations derived in References 1 and 2. Solutions to both the

partial differential equations and to the finite difference equations are

obtained by a combination of analytical and numerical techniques entirely independent of the purely numerical procedure described in Reference 2. Agreement

1
between the exact solution to the difference equations described in this paper

and the numerical solutions described in Reference 2 was found nearly to the

accuracy specified for the iterative procedure used to solve the pressure equation described in Reference 3 (usually 10" 6

). On the other hand, comparison

of the solution to the difference equations and of the corresponding solution

to the partial differential equation determines the magnitude of the discretization error made by replacing the partial differential equations by finite

difference equations. For the number of grid points used in most of the calculations performed to date, the discretization error was a few percent. In

Reference 4 an additional analytical solution is presented, which was also

used to assess stability and accuracy of the numerical scheme.

In Section 2 the partial differential equations are presented. These equations

already are restricted from the more general ones presented in References 1 and

2. The grid is also shown upon which the discretized variables are defined.

The first solution, presented in Section 3, is for a fully three-dimensional

time-dependent irrotational , incompressible flow in an enclosure driven by

sources and sinks of fluid as specified by the heat source. This problem

arises from the full nonlinear equations with boundary conditions, in continuous or discrete form, by requiring that the velocity field be irrotational

and the density remain constant. The resulting nonlinear partial differential

equations and difference equations can then be solved. The solution to the

difference equations is used to determine, in detail, the accuracy with which

the code described in Reference 2 solves the nonlinear finite difference equations in this special case. The solution to the corresponding continuous

problem can be used to assess the discretization errors involved as a function

of mesh size.

2
The second set of solutions is presented and discussed in Section 4. This

problem arises when the density is taken to be constant, the heating is assumed

to be zero, the velocity field is taken to be two dimensional and derivable

from a stream function only, the vorticity is taken to be a constant, and the

flow is independent of time. The discretized problem can be solved exactly,

using a combination of analytical and numerical techniques and this solution has

been used to check the code described in Reference 2.

2. Equations

The general problem considered is that of the fluid motion in a rectangular enclosure produced by a specified heat source. The equations, derived in

Reference 1 and rewritten in more convenient form for numerical integration in

Reference 2, are for an inviscid, non-heat-conducting perfect gas. In the

solutions presented in this paper, density variations are taken to be small,

and buoyancy effects are neglected. Then the effects of density variations do

not alter the flow to a first approximation.

The equations describing such a flow fluid are presented below. They are a

simplified version of Equations (11) of Reference 2.

3p 3p 3p 3p — + U — + V — + w —

3t 3 X ay 3 Z

P

0 D(x,y,t) (la)

3U 2 + + 1 —- + 1/2 vq - u x to = - - vp

3t p 0 db)

3
V* (1c)

1 2 ~ V P

• Po

3D- + 1/2 V

3

1

2 2

q

- (u X to)

v

.u = D(x,y,z,t)

where u = (u,v,w) is the vector velocity and

(Id)

2 2 , q

= u + v

2 + w (2a)

-»• >

to = V

X u (2b)

dp o — Y_i

= y / v Q(x,y,z,t) dV (2c)

D(x,y,t)

1 YP (t)

0

(y- 1) Q(x,y,z,t)

dp 0

"

dt

Q(x,y,z,t) = Q

q X ,p){T) tanh At exp [-$ x (x-xc

)‘

- 3 y (y-yc ) “ Xz ^

(2d)

(2e)

Po is a constant, reference density, Q(x,y,z,t) is a heat source of form specified in (2e) with constants 3 x »3y and x used to determine its spatial extent and

A its temporal scale, Q

q the heat source magnitude, and P 0 (t) is a mean pressure

in the enclosure.

4
The boundary conditions specify that there will be no outflow or inflow at

boundaries: the normal velocity is zero. Hence

u = 0 or v = 0 at vertical boundaries

(3a) w = 0 at horizontal boundaries

and consequently the normal derivative of pressure at boundaries is zero

i-P

. 0 (3b)

3n

The equations can be made dimensionless by introducing a length scale equal to the

height of the enclosure, a velocity scale determined by the magnitude of the divergence (or the heat source) specified in Equation (2d), the ambient density p

Q and

the ambient or initial pressure, Poo. When this is done, the enclosure becomes one

unit high, and in the equations above, p , p* and Q

q can be taken as unity.

In Figure 1, a two-dimensional rectangular enclosure in dimensionless variables is

is shown together with a schematic two-dimensional representation of the spatial

grids used for the finite difference scheme. The grid formed from solid lines represents the basic mesh into which the enclosure is divided: in general there are

I mesh cells in the x-direction, J mesh cells in the y-direction, and K mesh cells

in the z-direction.

Upon this basic mesh, the components of the vector (u,v,w) and the component of

the vector vorticity are defined: the location of the velocity components for a

sample mesh cell are also shown in Figure 1.

5
A second grid, formed by joining the center points of the basic grid cells, is

that upon which scalar quantities such as pressure p are defined. In Figure 1

the pressure at the center of the sample mesh cell is shown.

The following discretely evaluated functions will denote approximations to the

corresponding solutions to Equations (9) and (10):

u*

1 = u(i5x, (j-l/2)sy, (k-l/2)sz, nst)

i jk

v" = v((i-l/2)sx, jSy, (k-l/2)sz, nst)

i jk

w" = w((i-l/2)6x, (j-l/2)6y, k6z, nst)

i jk

( 4 )

p

n s p( ( i -1/2 )6x, ( j-1/2 )6y , (k-l/2)sz, nst)

i jk

D° s D (

( i - 1/2 ) 6 x , (j-l/2)5y, (k-l/2)«z, n6t) , i jk

<j> n

= <j> ( (i-l/2)6x, (j-l/2)5y, (k-1/2 )6 z , nfit)

l jk

where 6x = 1/(I-AR), 6y = 1/ ( J • BR ) and Sz = 1/K are the mesh cell sizes in the

x- , y- and z-directions respectively and where St is the time-step size. Such a

staggered grid is commonly used for multidimensional finite difference integrations [5].

6
In the following two sections, solutions to nonlinear equations (1) - (3) under

additional restrictions are given. Then the discretized approximation to each

restricted set of partial differential equations is also presented and solved.

It should be noted that each of the discretizations chosen to approximate the restricted set of partial differential equations is compatible with the discretizations chosen in Reference 2 for the more general partial differential equations

describing buoyant convection due to heating. Because of this compatibility, the

solutions given in Sections 3 and 4 have been found to be very useful for comparison with the numerical results obtained from the general algorithm described in

Reference 2.

3. Irrotational Flow Driven by a Divergence

The nonlinear problem is one in which the density is taken to be constant and the

velocity is assumed to be derivable from a velocity potential. In this section,

first the solution to the partial differential equations under these assumptions

will be presented, and then the corresponding solution for the difference equations

is given.

A. Continuous Solution

We start the analysis by noting generally that the fluid velocity can be written

as the a gradient of a scalar potential, plus the curl of a vector potential

+ +

u = V 4> +

v

x ^ (5)

For an irrotational flow, ^

= 0, and Equation (Id) becomes:

V% = D(x,y ,z,t)

7
on 0 < x < 1/AR, 0 < y < 1/DR, 0 < z < 1. The boundary conditions are found from

j\ Equation (3) to be — = 0 on all boundaries. When D(x,y,z,t) is separable, the

equation for <j> is separable, and <j> can be determined either by Fourier series or

by Green's functions and complex-variable techniques in the two-dimensional case.

Then the spatial distribution of <p is specified by the spatial distribution of D

and the time dependence of <p is equal to that of D. The velocity field is found

from

The pressure field is found from the dimensionless version of Equation (lb) by

noting that the vorticity is zero. With the density constant and the above definitions for the velocity components in terms of the potential. Equation (lb) becomes in component form

3n

u - (7)

(8a)

(8b)

(8c)

where the density p

Q in dimensionless form is unity and

(9)

Integration of Equations (8) yields

( 10 )

8
The function g(t) is determined from the condition

1 1/AR 1/BR , / dz / dx / dyp(x,y,z,t) = 0 (11)

0

J 0 o

See the discussion following Equation (8) of Reference 2 concerning this condition. Integration of Equation (10) over the total volume 0 < x < 1/AR,

0 < y < 1/BR, 0 < z < 1 and requiring that

/* dz

/J

/AR

dx /^dy<t>(x,y,z,t) = 0 (12)

yields for g(t)

1 1/AR 1/BR

g(t ) = (AR-BR/2) !q

dz l

Q

dx f

Q

dy

The pressure field is then found from the potential field through Equation (10)

where g(t) is determined from Equation (13). Comparison of the solution obtained

above for the continuous case with that obtained below for the discretized equations then provides an estimate of the discretization error for a specified mesh

size.

B. Discrete Solution

When the density is constant and the discrete velocity components are derivable

from a discretized potential, a solution to the difference equations can be obtained by a procedure analogous to that used above to obtain a solution to the

partial differential equations. First, reference should be made to Figure 1

where the discretized potential is defined at the center of a grid cell. Then,

the discretized velocity components and the potential are related, to second

order accuracy, by the relations

9
n 1 / n n u

= — <p - <(> ijk 5x v i+l,j,k ijk

(14a)

n 1 / n n

v = — <t> - ijk Sy Vi,j+l,k ijk

(14b)

n 1 / n n w = -— ( <t> -<()

ijk 6z \ ij,k+l ijk

(14c)

Substitution of these relations into the discrete version of Equation (Id),

yields for l<i<I,l<j<J,l<k<K,

1/n n n\l/n n n

~—-

( d> - 2d> +d> ] + —- d> -2d> + d> Sx2 \ i+l,jk ijk i-l,jk/ 6y2 V i,j+l,k ijk i,j-l,k

1/n n n \ n + —- / <p - 2<p + <p | = D 3_

2

\ i J» k+1 ij k i j»k-l/ ijk

(15)

where D is the discretized version Equations (2d) and (2e),

ijk

The discrete boundary conditions are

n n n n

d> = <t> and d. = d> for OcjcJ ,0<k<K

0,jk

v l,jk i+I

,

jk

v I,jk

(16a)

n n n n

4> ^

= <b and d> = d> for 0<i<I,0<k<K

i,0, k i,l, k i,J+lk

y i,J,k

(16b)

n n n n and <j> . = <j>. for 0 < i < I , 0 < j < J ij,0 i j ,1 i J »K+1 l j ,K

(16c)

10
n n n n n n where d> ,

d> ,

<b , <p and <j>. . ^ , <j> 0,jk I+l,jk i,0,k i,J+l,k ij,0 ij,K+l

are values defined at the center of fictitious cells adjacent to and outside the

region considered. These values are eliminated from the problem when Equations (16)

are substituted into Equation (15). Equation (15) and boundary conditions (16) can

be solved using fast direct methods such as a three-dimensional generalization of the

n package of Swarztrauber and Sweet, Reference 6. When D is a separable function of

ijk

i, j, k and n as is the case for the forms chosen, the spatial portion of <p for all

ijk

n can be determined by only one solution of the linear algebraic system represented by

Equations (15) and (16). The discretized velocity field is then determined from

Equations (14).

The discretized pressure field is derived from the discretized version of Equations

(lb). Since the vorticity is zero and p 0

= 1, these equations can be written

n+1 n

u = U + 6 ijk ijk

1

26 x i + l > jk.

( n \

kq

i jki

2_

1 /n n _ n

^ ' ~x V

i+1 * j

' 1jk

J (Ha)

n+1 n (1 v = v +6

ijk ijk ( 2Sy

( n

Iq

i ,j+l,k.

( n

Vq ijk

1 / n n

<$y \P i,j+l,k

P ijk

(17b)

w

n+1

/ ijk

n

= W +6

ijk

1 26z V,.

2 —

1 / n n

'

6z \P ij,k+l

P ijk

(17c)

where

6 = St

2st

for n = 1

for n > 1 (18a)

11
for n = 1

i jk

i jk

n-1

i i jk

i jk

n

i i jk

n-1

f i jk

w

i jk w

i jk

n-1

i i jk

(18b)

for n > 1

1

1

and where q is defined by Equation (2a) in discretized form,

i jk

~n and U

i jk

1/2

i-1 » jk

,

etc.

The upper bracketed values in Equations (18) are the ones to be chosen during a firstorder, explicit starting time step whereas the lower value is to be chosen for a leapfrog time step. In this form, therefore. Equations (17) and (18) can be used to represent succinctly both first-order and leap-frog time steps. It is important to note

that the time step may be reduced during the computation to avoid violating the

stability criterion.

n 6t _< max

1<A<1 l<j<J

1 <k<K

+

'liiljkl IVljlcI liiljkl + +

- 6x Sy Sz

-1

(19)

When this is done, the time-marching algorithm is restarted with the values of the

dependent variables at the last successful time step used as initial conditions. A

first-order time step is then taken using the new time step size to obtain values of

the dependent variables at two levels, and then the leap-frog is resumed.

12
We substitute definitions (14) into Equations (17) and define a potential

i jk

n

<t> . .. 1 jk

n-1

. .. ljk

( 20 )

as in Equation (18) to take account of the different form of Equation (17)

during a first-order, starting time step. Then Equations (17) become

16x

n+1 n / n i+ljk ' $ i+l,jk \ i+1 , jk/ n + i— + n

6 2 i+1 , jk

r n+1 n / n \

<Mjk " ^ i jk

j vi j k/ t _n

6 2

P

i jk

= 0

(21a)

1«y

n+1 n / n

1 »j + l»k ' $ i ,j+l ,k \

q i,j+l,k/ n — + + p

6 2 i,j+l,k

r n+1 n / n jk " $ ijk \ H ijk) n

6 2 i jk

= 0

16Z

n+1 n / n \

2

^ i j » k+1 ’ $ ij,k+l VQij.k+l) n + A-.- /. + pn

- ,• * 2

MlJ Pi j, k+1 ‘

’ n+1 n <Mjk " $ ijk

n ^ijk/ n +

Pi jk = 0

(21b)

(21c)

13
The difference equations are satisfied by a function independent of i, j and k

( 22 )

Equation (22) is the discretized version of Equation (10). To determine g

11 , we

must apply the discrete analogue of Equation (11)

I J K n «x «y 62 l l l p

= 0

i=l j=l k=l 'Jk

and require

I J K n

I l I = 0. i=l j=l k=l

Then

1 I J K IM-j-jk

I6x J6y Kfiz g

n = gn =

j J \ 6x 6y 6z

v AR-BR i=l j=l k=l ^

or

n AR-BR l J K

g

n = __ J l l 6x 6y 6z

2 i=l j=l k=l

^i+1, jk

~ *i-l,jk

N

V 2sx /

2 2

*1»j +1 .k ^i>j-l,k\ +

/ *i,jk+l ~ ^i j ,k-l\

26y 26 z

(23)

(24)

(25)

(26)

14
where we have used the relation

2 2

(

q.

i jk

n )

2

l+l, jk

~ i-l.Jk

26 x

( 27 ) 2

+

*i j ,k+l ~ ( Hj,k-l

^ 26 z

n n

Therefore, once the discretized potential has been determined from Equation

The solution procedure outlined above has been carried out, namely Equation

(15) was solved numerically using a three-dimensional Poisson solver and then the

velocity field was determined from Eqs. (14) and the pressure field from Equation (22). The results in the two dimensional case were compared with those

computed from the algorithm described in Reference 2 when, in this algorithm, ~n n the perturbation density p and the vorticity u were specified to be zero.

i J i j

The parameter in the pressure solver described in Reference 2, which determines

the accuracy to which the linear algebraic system for the discretized pressure is

solved, was taken to be e = 10" 6

. This value was chosen as a reasonable compromise given the machine accuracy (36 bits) of the NBS UNIVAC 1100 and the

iterative nature of the pressure solver. Agreement up to a few tens of time

steps, between dependent variables (both components of velocity and pressure)

determined by the two different procedures, was found to a few parts in 10“ 6

.

Accuracy of the numerical procedure described in Reference 2 was found to

degrade slowly, being a few parts in 10“5 after a few hundred time steps.

It is interesting to describe the physical content of the problem discussed in

this section. When the density is taken to be constant, the driving function

D(x,y,z,t) can no longer be assumed to result from heating and cannot be inter-

(15), the discretized pressure field is obtained from Equation (22) using

Equation (26) to determine g

n and using Equation (27) for

15
preted as a source of specific volume. Rather, D(x,y,z,t) must be interpreted

as a distribution of sources and sinks of fluid in an enclosure. With this

interpretation the potential flow, with no outflow or inflow at the boundaries

of the enclosure, has some interesting features, which are described briefly

bel ow.

The form of D(x,y,z,t) used for these calculations is defined by Equations (2c)

and (2d) with Q(x,y,z,t) given by Equation (2e). (For the solutions shown below, Q is centered along floor of the enclosure.) Since the continuous potential is given by Equation (6), with the corresponding velocities given by

Equations (7), and, since Q and therefore D is separable in time, the spatial

dependence of the potential and the velocity components remains the same

throughout the history of the problem. The function D(x,y,z,t) represents

sources of fluid of greatest strength at the center of the "heat source" Q and

decreasing in strength as the distance from the center increases. At some

distance from the center, the sources become sinks so that the integral of the

sources (and sinks) over the enclosure is zero. Therefore, the velocity field

is one describing flow from the "heat source" to the remainder of the enclosure

as shown in Figure 2 for all times during the history of the flow.

In contrast, the pressure changes dramatically during the flow. During the

early portion of the flow, all nonlinear manifestations can be neglected, and

/ \ 8 <p in Equation (10), p

= - — . The potential <t> will behave qualitatively like

3

1

the negative of the source function D; i.e., where the source is maximum (at

the center of the "heat source") the potential is minimum, and the pressure is

highest, as shown in Figure 3. Hence the pressure is maximum where D is and

decreases to negative values away from the "heat source" Q.

16
Later in the flow, when the source is asymptoting to its final value, the nonlinear terms dominate and — becomes relatively small in Equation (10). At

3

1

any particular time g(t) is a constant, and Equation (10) becomes approximately

P = (28)

The term in square brackets is the magnitude of the velocity squared. This

equation states that the pressure is minimum where the velocity magnitude is

largest, and this velocity magnitude is largest where the gradient of the "heat

source" Q is maximum. Figure 4 shows contours of constant pressure late in the

flow when the heat source is asymptoting to its final value.

4. Two-Dimensional, Steady-State Flow with Vorticity

A. Continuous Solution

When the density is constant, no heat source is present and we restrict the flow

to be two dimensional, then the continuity equation, Equation (Id), becomes — + - = 0

3x 3y

and the velocity components

ponent, the stream function

3 ^

u

= — , v

= - ay

In the steady state, the xEquation (lb), become

(29)

are derivable from a vector potential of one comip, (see Equation (5))

dip

3X

and y- components of the momentum equation.

(30)

17
( 31 )

1 3 1 ap

- — (u z + V

2

) - Vw = - — —

2 3X p 08

X

1 d

r ?

1 3P

- — (u 2 + V

2

) + Uw = - — —

2 ay p 0 ay

where p 0 is constant and the scalar vorticity is defined by

av a u ..

OJ = — - — = - V

2]P ax ay

(32)

Define

H = p/ Po + 1/2 (u 2 + v

2

) (33)

and substitute for u and v from Equation (30) into Equations (31). Then

aH ai|> — + (o — = 0 ax ax

aH dip — + o) — = 0

ay ay

If we now assume that u is a constant

(34)

0) = + b (35)

then Equations (34) become

(H + b\p ) = 0 ax

— ( H + bip ) = 0

3y

18
or

H + bi|i = C = constant (36)

where C is evaluated, for example, by requiring that the mean value of p,

integrated over the enclosure be zero.

Equation (32) becomes

V

2 *

= - b (37)

subject to the impervious wall boundary conditions, which become

4

> = 0 at x = 0, 1

and y

= 0, 1 (38)

Therefore, procedurally, one can solve the Helmholtz equation (37) subject to

boundary conditions (38) to obtain the stream function 4

>. Then Equations (30)

yield the velocities, and Equations (33) and (36) the pressure.

B. Discrete Solution

As in the continuous case, when the density is constant, no heat source is present and the flow is two dimensional, the continuity equation becomes

u n

. . u

n

,

. v

n

. _ v

n

. . 1J 1-1,

J

+

1J 1.J-1 _ 0

fix <5y

and the discretized velocity components are obtained from the stream function

as follows

(39)

n u

i j

'P n

i ,j-l

5y

n v ij

n1-1 > J

fix (40)

19
The momentum equations. Equation (31), become in discretized form

1 1

P 0

<5x

(41a)

(41b)

where

n

O)

i j Sx sy

(42a)

ij

1

2

(42b)

(42c)

Equations (40) imply that

n 1 / n w . = —r (

6X

2 \ 1+l.J

n n \ 1/n n n \ , % 2*. . + ,

. ) + - 2*.. + . . )( 43 ) ij i-l, J/ sy \ !»J +1 iJ i,J-l/

20
and

n 1 /

n n

u

= — U - ^ u 25y \i,j+l i,J-l

i j

1/n n 26x \i+l,j ’N-IJ

ij

(44)

Equations (41) can be written

1 6x

1 + -2

1 V

1 / n \

2

n P.., ./Pn “ « q. . - P. .Ip

2 \ i+1 , j ) i+l» j

0 2 l ij ij

nl/n n

\ n 1/n n

to —• I \b -ib i+io | ib - ib ij 26 x V i+l,j i-1 , j / i,j-l 26 x \ i+l,j-l i-l,j-l

(45a)

= 0

1

1 + -2

1/n \

2 n 1

2 (

q T,j+l ) + P

i ,j+l

/P ° ' 2

- P. ./P, IJ ' nl/n n\n 1 /n

(0 I - lb

) + 0) (\h

ij 25 y V i , j +1 i,j-l/ i- 1 , j 26 y V i- 1 ,

j

+1

n1-l.j-l

(45b)

= 0

As in the continuous case, we take

n

to = + b

i j

(46)

and note that

1 /

n \ n , - q. . + P. 7p

0 2 \

ij / ij 0

(47)

b / n n n n + r I’P. . + 4'. , . + 4> + ) = C = constant

4 \

ij i-l, J i,j-l i-1 , j -

1

21
The constant in this equation is evaluated by requiring that the mean value of

n

p over all mesh points be zero (to make it compatible with the general dei j

fi nit ion of this quantity in the algorithm described in Reference 2).

To determine a solution to the nonlinear difference equations, therefore, the

discretized Helmholtz equation obtained by combining Equations (43) and (46)

with boundary condi tons that i|».

. = 0 on the boundaries, i = 0, I and j = 0, J,

is solved using a software package such as the one developed by Swarztrauber

and Sweet in Reference 6. From this solution for the stream function, the

velocities are determined from Equations (40) and the pressure from Equation

(47). This solution procedure involves analytical and numerical techniques

totally independent of the numerical algorithm described in Reference 2.

The steady-state solution to the difference equations was used to test the

algorithm of Reference 2 by using the values as initial conditions to the code

(for both initial time levels) with the heat source set to zero and the perturbation density also set to zero. The computer code was then allowed to take

several hundred time steps, and the values after different numbers of time

steps were compared. It was found that there was no instability in the algorithm. The solution replicated itself very well, with a slow degradation in

the number of significant figures with which the solution was able to duplicate

itself as the number of time steps increased. Initially, the solution after a

few time steps agreed with the exact solution to a few parts in the sixth

significant figure, as expected from the tolerance set on the pressure solver.

(49)

22
After four hundred time steps, the agreement was to a few parts in the fourth

significant figure, representing a degradation which we felt was reasonable.

In Figure 5 we present plots of the horizontal velocity as a function of vertical location at three horizontal positions and of the vertical velocity as a

function of horizontal location at three vertical positions. This figure shows

a flow field which is rotating in a clockwise direction with the greatest

velocities at the edges and smaller velocities in the interior. The flow field

represents "stirring in a rectangular tea cup" with no dissipative effects of

viscosity. Horizontal velocity plots at a specified horizontal position do not

agree with vertical velocity plots at the corresponding vertical location due

to details in the graphics package used to display the data.

23
REFERENCES

1. Rehm, R. G. and Baum, H. R., The Equations of Motion for Thermally

Driven, Buoyant Flows, J. Research of the National Bureau of

Standards 83, No. 3, pp. 297-308 (May-June 1978).

2. Baum, H. R., Rehm, R. G. , Barnett, P. D. and Corley, D. M. , Finite

Difference Calculations of Buoyant Convection in an Enclosure,

Part I: The Basic Algorithm, SIAM JSSC. 4_, No. 1, pp. 117-135

(March 1983).

3. Lewis, J. and Rehm, R. G. , The Numerical Solution of a Nonseparable

Elliptic Partial Differential Equation by Preconditioned Conjugate

Gradients, NBS Journal of Research, 85, No. 5, pp. 367-370

(Sept. -Oct. 1980).

4. Baum, H. R. and Rehm, R. G. , Finite Difference Solutions for

Internal Waves in Enclosures, National Bureau of Standards Report in

preparation.

5. Harlow, F. H. and Amsden, A. A., Fluid Dynamics, A LASL Monograph , Los Alamos Scientific Laboratory Report LA 4700, Los Alamos,

New Mexico (June 1971).

6. Swarztrauber, P. and Sweet, R., Efficient FORTRAN Subprograms for

the Solution of Elliptic Partial Differential Equations, NCAR-TN/IA109, July (1975).
Figure Captions

Figure 1: A rectangular enclosure in dimensionless variables and the a schematic representation of the basic mesh into which the enclosure

is divided.

There are I cells in the x-direction, J cells in the y-direction

and K cells in the z-di recti on. Also shown is a single mesh cell

and the location within this mesh cell of the discretely defined

dependent variables, velocity components u, v and w, pressure p and velocity potential <f>.

Figure 2: The velocity field in the two-dimensional case. The horizontal

velocity is plotted at three horizontal locations (in each case as

a function of the vertical coordinate). The vertical velocity is

shown at three vertical locations as a function in each case of the

horizontal coordinate. The flow field is from the sources to the

sinks; the spatial pattern does not change with time for this

potential flow.

Figure 3: Contours of constant pressure early in the calculation of the

two-dimensional potential flow. Contours greater than the mean pressure are shown as solid lines, while those below the mean are dashed. The pressure is maximum where the distribution of sources

is and decreases to negative values away from the "heat source"

(where the sinks are).

Figure 4: Contours of constant pressure late in the calculation of the

two-dimensional potential flow. Contours greater than the mean pressure are shown as solid lines, while those below mean are dashed. Note that the pressure is negative where the sources are and positive away from the sources in contrast to Figure 3. The

change in character results from the Bernoulli effect (see the text

for an explanation).

Figure 5: Horizontal velocity as a function of vertical location at three

horizontal positions and vertical velocity as a function of horizontal location at three vertical positons. Flow field is rotating

in the clockwise direction with greatest velocities at the edges

and smaller velocities in the interior, representing a stirring in

the rectangular region.

CO

I

r

'/(I-*)
a
DENSITY

C0NT0U*

NRP DT=

S.CCXC-O*. 3
DENSITY

C0NT0UR

HfiPflT I=

7.30000*00 4
[ T

TT

:> •

S'
NES-114A 1'rev. 2-80

U.5. DEPT. OF COMM. 1. PUBLICATION OR 2. Performing Organ. Report No. 3. Publication Date BIBLIOGRAPHIC DATA

SHEET (See instructions)

REPORT NO.

NBSIR 84-2932 September 1984

4. TITLE AND SUBTITLE

Finite Difference Calculations of Buoyant Convection in an Enclosure, Part II: Verification of the Nonlinear Algorithm

5, AUTHOR(S)

R. G. Rehm, H. R. Baum, P. D. Barnett, and D. M. Corley

6. PERFORMING ORGANIZATION (If joint or other than NBS, see instructions) 7. Contract/Grant No. NATIONAL BUREAU OF STANDARDS

DEPARTMENT OF COMMERCE 8. Type of Report & Period Covered

WASHINGTON, D.C. 20234

9.

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

10.

SUPPLEMENTARY NOTES

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

11.

ABSTRACT (A 200-word or less factual summary of most significant Information. If document Includes a significant

bibliography or literature survey , mention it here)

Earlier, a novel mathematical model of buoyant convection in an enclosure was developed. The nonlinear equations constituting this model have recently been solved

by finite difference methods in two dimensions.

In this paper two solutions, obtained in special cases, to the model equations are presented. For both cases the solutions to the partial differential equations and to

the finite difference equations used to approximate the differential equations are obtained by combinations of analytical and numerical techniques. Agreement between

the exact solutions to the difference equations described in this paper and the independently obtained numerical solutions was found nearly to the accuracy specified

(usually 10"6) for an iterative procedure used in the computational scheme.

.2. KEY WORDS (Six to twelve entries; alphabetical order; capitalize only proper names; and separate key words by semicolon s) Accuracy; Analytical solutions; Euler equations; Finite difference equations; Fluid

dynamics; Numerical solutions; Partial differential equations; Stability

3. AVAILABILITY 14. NO. OF

PRINTED PAGES

EX) Unlimited

E71 For Official Distribution. Do Not Release to NTIS

[ ~\ Order From Superintendent of Documents, U.S. Government Printing Office, Washington, " 20402. D.C.

35

15. Price

[~X] Order From National Technical Information Service (NTIS), Springfield, VA. 22161 $8.50

USCOMM-DC b04J-P»0