Ref

IM BSIR 79-1754

NAT'L INST. OF STAND & TECHi R.I.C,

A1110S DME12D

A Linearized Finite Difference

Computation of Fluid Heating

in an Enclosure

Ronald G. Rehm and Martin R. Cordes

Mathematical Analysis Division

Center for Applied Mathematics

National Engineering Laboratory

Howard R Baum

Center for Fire Research

National Engineering Laboratory

John Lewis

Boeing Computer Services, Inc. Mail Stop 9C-01

P.O. Box 24346

Seattle, Washington 98124

May 3, 1979

Interim Report

"QC

100

U56

79-1754

U S. DEPARTMENT OF COMMERCE

NATIONAL BUREAU OF STANDARDS

TS at M c " NBSIR 79-1754 QCjqD

US&

-VH ^4

A LINEARIZED FINITE DIFFERENCE

COMPUTATION OF FLUID HEATING

IN AN ENCLOSURE

Ronald G Rehm and Martin R. Cordes

Mathematical Analysis Division

Center for Applied Mathematics

National Engineering Laboratory

Howard R. Baum

Center for Fire Research

National Engineering Laboratory

John Lewis

Boeing Computer Services, Inc. Mail Stop 9C-01

P.0. Box 24346

Seattle, Washington 98124

May 3, 1979

Interim Report

U.S. DEPARTMENT OF COMMERCE, Juanita M. Kreps, Secretary

Jordan J. Baruch, Assistant Secretary for Science and Technology

NATIONAL BUREAU OF STANDARDS, Ernest Ambler, Director
, . * .-a I '
A Linearized Finite-Difference Computation of

Fluid Heating in an Enclosure

Ronald G.

Howard R.

John Lewi

Martin R.

Rehm

* Baum

-H- s

Cordes"*"

+ Center for Applied Mathematics, National Engineering

Laboratory, National Bureau of Standards

* Center for Fire Research, National Engineering

Laboratory, National Bureau of Standards

++ Boeing Computer Services, Inc., Mail Stop 9C-01 , P.0.

Box 24346, Seattle, Washington 98124
Abstract

In an earlier paper, approximate equations of motion were derived

which are applicable to nondissipative , very nonadiabatic, buoyant flows

of a perfect gas. In the present paper, these approximate equations are recast in a form in which they can be integrated by finite difference

techniques. These nonlinear equations are made dimensionless, specialized

to two dimensions and linearized. Finite difference approximations, using

central differences in space and leapfrog in time, are made to the linearized equations. To start the computation, a first-order scheme is used

for the initial time step. The dependent variables are density, the

horizontal and vertical velocity components and pressure, and these

variables are defined at various positions within a grid cell (see Figure

lb) . The computer time required for various size computations and the

accuracies obtained are reported.

At each time step, a large sparse system of linear algebraic equations,

the finite difference approximation to the elliptic equation for the pressure, must be solved. A hybrid method, combining the interative algorithm

called conjugate gradients with a fast direct method for solving Poisson's

equation, is used to solve the algebraic system efficiently and accurately.

(Typically, the solution of the algebraic equations arising on a 31 x 31

mesh is determined to five significant figures after between two and five

iterations of the conjugate gradients algorithm using from two to five

seconds of CPU time on the N.B.S. UNIVAC 1108.)

Finally, results of some computations are presented using computergenerated plots. A discussion is given of the software required to generate

these plots. A very brief discussion is presented of the analytical procedures required to obtain a solution to the differential equations describing

heating of a homogeneous fluid. Computations of the solution to the difference equations are found to confirm the conclusions of the analysis.

Also, these computations are found to yield flow fields which behave

qualitatively the same as flow fields observed in a room fire. Examples of

computations of the flow fields induced by heating in a stratified fluid

are also presented and discussed.

ii
Contents

Pg.

I. Introduction* 1

II. Formulation 4

A. Full Nonlinear Equations 1+

B. Linear Equations 7

III. Finite Difference Equations n

IV. The Solution of the Pressure Equation j_8

V. Computational Procedure for Solving Linear Finite Difference Equations. .. 2k

VI. Computer-Generated Results 29

A. Graphics 29

B. Linear Heating in an Ambient Homogeneous Fluid 31

C. Linear Heating in a Stratified Fluid 45

C.l. Heating in a Stratified Fluid by a Centered Heat Source k6

C.2. Heating in a Stratified Fluid by a Noncentered Heat Source 57

C.3. Heating in a Stratified, Two- Layer Fluid by a Noncentered Heat

Source 69

VIII. Summary 67

References 70

iii
r'*t»
1

I . Introduction

Convection in fluids heated from below has been of interest for over

one hundred years. Much of this interest has arisen from astrophysical

1-3

and geophysical research. More recently, technological applications in

areas such as nuclear reactor design, liquid natural gas transport, and

fire safety have led to the analysis of phenomena occurring injvery different ranges of geometrical and physical parameters than those considered

in earlier studies.

The principal application of interest to the authors is buoyant

convection induced in room fires. The National Bureau of Standards,

acting through the Center for Fire Research, is involved in a wide rante

of investigations into various pjiysicaland chemical phenomena arising in

unwanted fires. The buoyant convection studies reported here are part of

that effort. The heavy mathematical and computational emphasis in this

project, together with its intrinsic interest, have led to its sponsorship as a joint activity with the Center for Applied Mathematics. This

report is meant to record the current status of the project. As such, it

should be considered an interim working document, rather than an account

of a completed piece of work.

The approximate equations governing the large scale motion of a gas

in an enclosure driven by a prescribed volumetric heat source were derived

4

in an earlier paper. The resulting mathematical model, a ’’thermally expandable fluid," is rather different from other idealized fluid equations,

such as the Boussinesq model, usually employed in buoyant convection

studies. The differences arise from the need to consider large density

variations in a closed geometry of limited vertical extent. This
2

combination of parameters also arises in some boiling water nuclear reactor

studies, where a single fluid approach to the two phase flow results in a 22 21 thermally expandable fluid model. *

Due to the complexity of the equations of interest and the limited

experience in dealing with them, considerable care was necessary in developing numerical computation techniques. The thermally expanding fluid

motions share some common features with both incompressible fluid flow

and with the hydrodynamic models used in numerical weather prediction.

A sampling of the techniques used in the analysis of these related

equations may be found in references 5-9. A study of these methods led

the authors to the conviction that the computational procedures should

be developed on a related system of linear equations before nonlinear

calculations were attempted.

Three features of the thermally expanding fluid model which were

felt to be of particular importance are amenable to analysis within the

context of a linear problem. They are the numerical computation of

internal wave oscillations over many cycles, efficient calculation of

the self adjoint elliptic equation governing the pressure distribution,

and the use of a distributed, volumetric heat source of prescribed strength

(as opposed to boundary conditions) to drive the motion in a rectangular

enclosure. The linear system was chosen by formally linearizing the thermally

expanding fluid model about a rest state stably stratified in the vertical

direction. The linearized motion is driven by the prescribed heat source.

The homogeneous solutions to the problem posed above were studied by

the authors in reference 10. These solutions are internal waves which

arise due to the stable stratification. These waves were used to test
3

finite difference schemes for accuracy and stability. Second order accurate

spatial differences were used in conjunction with both first and second

order accurate time differencing to derive candidate finite difference approximations to the homogeneous (internal wave) equations. The finite

difference equations were solved exactly by analytical means and compared

with the known solutions to the continuous problem. The finite difference

scheme used in the present study to solve the inhomogeneous equations was

selected on the basis of the analysis reported in Reference 10.

This report describes the numerical implementation of the finite

difference scheme selected on the basis of the analysis outlined above.

The motions are driven by a prescribed distributed heat source located at

the bottom of a rectangular enclosure. Section II presents the non-linear

continuous equations derived in reference 4 and describes the linearization.

The finite difference equations, incorporating second order central differences in space and leapfrog in time are shown in section III. The

solution of the pressure equation, based on fast direct methods for solving

11 12 13 Poissons equation combined with an iterative conjugate gradient technique ’

is discussed in section IV. Section V outlines the solution procedure for

the coupled set of linear finite difference equations. Several examples

of computational results, together with the computer graphics used to

display them, are described in Section VI. Finally, a summary of the work

presented in this report is given in Section VII.
h

II. FORMULATION OF CONTINUOUS PROBLEM

A. Full Nonlinear Equations

4

In an earlier paper nonlinear equations describing three

dimensional thermally-driven buoyant flows in an enclosure had been

derived and their properties discussed. In this section these

nonlinear equations will be rewritten in a form appropriate for

numerical integration by finite difference techniques. This recast

set of equations is then linearized and the boundary conditions for

the linear equations are presented.

As in Reference 4 we consider an inviscid, non-heat-conducting

perfect gas. The magnitude and the spatial variation of the heat

source (representing the exothermic reaction in a fire) are taken

as known; justification for such a model is given in Reference 4

.

The fluid and the fire source are assumed confined in a closed

rectangular room with the center of the source along the floor. In

contrast to Reference 4, we consider only a completely enclosed room (no

leaks), and we confine attention in the linear problem to the two dimensional evolution of the flow.

Equations (ll) of Reference 4 are

fo ft) = p£T

( i

)
Here is density, the velocity in the i— coordinate direction

(i » 1, 2, 3), p is the pressure excess above the mean pressure

PQ

(t) in the room, T the temperature, the constant pressure

specific heat, R the gas constant, is the gravitational

acceleration and Q(x^, t) the specified volumetric heat source.

The spatially uniform mean pressure p

Q

(t) depends only upon time

and increases because of the heating within the room. It is

determined in a completely enclosed room by the equation

^

= Vjv

<?(*••>« 4V M

where ^is the ratio of specific heats, V is the volume of the room

and the integration is performed over this entire volume. Equation

(2) is a thermodynamic statement that the mean pressure rise as a

function of time is determined by the total heat added to the room.

(Heat can only be added or removed volumetrically and not through

the walls because thermal conduction and radiative transport have

been ignored in this model.)

Equations (l), derived in Reference 4, are an approximate set

of nonlinear equations which are applicable to highly nonadiabatic,

nondissipative buoyant flows of a perfect gas in three dimensions.

These equations were derived under the assumption that heat is added

slowly compared with the time required to equilibrate the pressure

over the spatial extent of the heat source, an assumption appropriate

to flows induced by fires. They are characterized by the fact that

the spatially uniform mean pressure appears in the energy and state

equations while the spatially variable portion of the pressure appears

only in the momentum equation. The equations admit buoyant or internal-wave
6

motions while "filtering out" high-frequency, acoustic waves. They

reduce to the Boussinesq equations when heating is mild, total density

variations are small, and variations in the mean background pressure can

be neglected (as would be the case if the room considerated here were open

or if the mean pressure variation were comparable to the spatial pressure

perturbation.)

Equations (l), fully nonlinear and three-dimensional, will be

recast into a form suitable for numerical computation. To do this,

we take the substantial derivative of the equation of state and use

this with the energy equation to eliminate the temperature. The

resulting equation describes the evolution of the density under heating

where

=

Tpo - iijr]

( 3)

(+)

Equation (3) and the continuity equation identify D(x^

t) as the divergence

2Ui _

Finally, as in Reference 4, the equation for the spatially variable

D (xz , ir) (sj

portion of the pressure is obtained by dividing the momentum equations

by density and taking the divergence of these equations. The resulting

equation is r~ ^k) =~[^r(a

The boundary conditions on these equations are that velocity normal

to any (impermeable) wall vanish.

Hi /ni = o

'dx-i V P

^3 U-cf)

A / Mr u

rn

where are the normal components of a vector describing the boundary

walls. From Eqs. (l) and these conditions, the appropriate boundary

conditions on the pressure equation are obtained
7

*‘Tfc-

“ P3'n: *

i

The complete set of recast nonlinear equations are gathered and rewritten below

%f, =-p

(

. _ -J— S»t> -- !3xrK§^)+

4e ~ -V

1- Q( k^V 11/ = Y^fcf QOr-D <?(*:,<?-

with boundary conditions

uc /n-i = o , ^

B. Linear Equations

In this subsection the equations derived above are linearized

and specialized to flow in two spatial dimensions (x = x^, y = Xg)

in a rectangular enclosure, 0 £ x <_ L, 0 £ y <_ H. The fluid in the

enclosure initially is assumed to be arbitrarily stably stratified

with density ^(y).

For convenience we will rewrite variables in terms of dimensionless

quantities. This procedure can be applied to the three-dimensional,

fully nonlinear equations also, but we will limit ourselves here to

presenting only the equations in two dimensions in linearized form.

Let u be the velocity in the horizontal or x-direction and v the

velocity in the vertical or y-direction. Define the dimensionless

variables, denoted by hats, as follows:

(*)

0°)
8

u.= si^T1 f- 0. v= e Lr

<J> = C (^ -+• e. § Q

p-o fo (t) -+ ^ £P 0 (|) + & f J (II

)

where

and

* = H £ , t~ (

H

/f)'^ t

<?- Q0 $C**?J)

f© ft) = I + £ £, ft)

£ = (H/g.)*. <?.(dr-ij

Here ^,is the ambient density and the pressure at a specified

altitude, here taken to he y=o. Qq is the magnitude of the volumetric

heat source, and e is a parameter relating the magnitude of the heat

source to ambient pressure (or the energy per unit volume in the

ambient gas )

.

For a heat source characteristic of the rate of heat release

within a small room fire, the parameter e is often small. Provided that

the time does not become too large, then, so that the cumulative effects

of the heat addition remain unimportant, e can be used as a small parameter

and Eqs . (10) can be formally linearized about the stably stratified rest

state. The resulting two dimensional equations are

(§) £

(13)

I
9

p

i r\ !& _ 8(*,p) - T

- * 1 ^ I 0(£,$,t)]

-^-

=*T«<f i*« q(*\§£)

The aspect ratio AR of the room is the ratio of the height H to the

length L of the rectangular room.

The boundary conditions corresponding to Eqs. (13) are that no

mass flow from the enclosure:

tl- o

II 0 o^-dl )c — , A

ir = o ^=o cu»^X =- !

These conditions, Eqs. (I*f), can be transformed into corresponding

boundary conditions on the elliptic equation for the pressure, the

fourth of Eqs. (l3), using the two momentum equations in Eqs. (13)

y ~ o & = )hz

_ as, _ _ 6

^ ^

3 ir » = (

(W

Os)

In the computations presented later, two forms for the ambient

stratification have been selected. For most calculations an exponential

stratification have been used:

£><£) = (~%/ys )

This form includes the special case when Y

g is very large so that is

practically constant. For some of the computations a "two layer" ambient

stratification was used: ^ 9o(§l-«f C-3/%)

P.<3> = >/ys„- y‘/y, ] <5 £ i

The heat source was assumed to have a simple, separable form:

<y*o
- 10 -

where f (t

)

= tanK a£, A is a constant determining the rate at which

A the heat source is "turned on," and xc is the location along the

"floor" (y=o) at which the center of the heat source is located.

A This form for the heat source Q was chosen as representative of a

fire for analytical convenience only; another form could equally have

"been chosen.

Several important observations must be made about Eqs. (l3)-(l6).

First, we note that D(x, y, t), from the form given in Eqs. (l3),

integrated over the rectangular region is zero. This fact, together

with the boundary conditions (IS), imply that the pressure equation

the third of Eqs. (|3), when integrated over the rectangular region

is identically satisfied. Such a condition is necessary if the pressure

equation is to have a solution at all since this elliptic equation

has Neumann boundary conditions. When the pressure equation has a

solution, this solution is only determined to within a constant, and the

constant is determined to be consistent with the initial condition, or

constant, chosen in the integration of the mean pressure equation, the

last of Eqs . (13)

.

Second, it should be noted that the pressure equation in this

linearized formulation is in fact separable: multiplication of the

pressure equation by ^Q,(y) shows that the coefficients of the

derivatives of y in this equation are only functions of y (the coefficients

of derivatives of x are unity). However, in the nonlinear formulation

of this problem, the pressure equation is nonseparable . Since the

linearized equations are being used to test solution techniques, the

solution procedure for the pressure equation discussed in Section IV

is for nonseparable elliptic equations.
11

III. FINITE DIFFERENCE EQUATIONS

In this section the finite difference equations used to approximate

the linear partial differential equations, Eqs. (13), and the boundary

conditions, Eqs. ( IH) and (15), will be presented and discussed.

These difference equations are written on a "staggered mesh" in the style

of reference 6 and were chosen to be approximations to the differential

equations which are second-order accurate in both space and time.

In Figure la, the rectangular enclosure in dimensionless variables

is shown together with a schematic 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 axe M mesh cells in the x-direction and

N mesh cells in the y-direction. (The hat designation on dimensionless

variables will now be dropped for convenience. Subsequently all quantities

will be understood to be dimensionless.') Upon this basic mesh, the two

components of the vector velocity (u, v) and single component of the vector

A ^ , 2vu-_ BU vorticity 0> = -jx ^are defined.

The second grid, formed by joining the center points of the basic

grid cells and denoted by dashed lines, is that upon which scalar

quantities such as density (0 and pressure >ys>are defined. In Figure la

the densities in the left-hand column of cells and in the bottom row

of cells are shown to indicate how they are enumerated for the numerical

computation.

In Figure lb a typical mesh cell is shown, illustrating where all

of the dependent variables in the finite difference scheme are defined

relative to the cell.
12

N

^\xXXX\ 'vV V \ N\ \ > V\ \X‘x \XVVv vVs

1

1

, 5#- r

i

i

1

1

1

|CU

!

\

i

I

i

i

1

1

1

1

l

1!__

i

i[c£

1

1

i

i

I

1

I

I

t

i

i

i

i

1

1 —J—

1

1

l

1

1

1

1

1 L__

P«.

1

1l__

ol

!

1

!ro a:

: v ^ «... "- —err

1

1

1nr o-.

4-vv VV

i

l —

L

61

'V V "V V T

*

: ^

\

II

«s

o

\\ £

a) £

c

o

T3 •H

CO OJ •H 3 O rH nj

i—! c o c X QJ CO 03 cO 14

•rl ai 44 CO

|4 144 o rH

CO 44 rH > •H T3 03 TO •H CJ

CO M

CO 0) bO X 3 X CO

i-4 44 3 0) 3 B o X CO

• H o 3 CJ "d

CO •r4 •H 03 3 X /—

\

CO £

QJ s -3 03 u

B X o

•H 3 1! M-J TO O Q) H

a. 'X x.

03 3 3 4-) Q-

*H II -3 4-1 CO

QJ CO S o

•H

!-i 3

*—

3 B CO C

CO M 4-1 O

o 3 O C3 •H

r*4o £

44 •Ho

4J

03 3 cx

4J

Cl) i-4 Z)

• t-4 M a

(4 iH 3 03 £

cO CJ 4J o

i-4 V| •H c a

3 44 03 bO >s 3 o

03 3 B u

cO V | 3 03 d

44 r- .c

CO O o o 4J CO

QJ 3 03 p4 r» CO U

od a

04

bO

•H

P4

V|

X

v I o

o

CO

"J3

CO

-.CO X

3

6

a) x

u

CO

TO3X

CO

cO X)

X

CO

a)B

3X•u 4=o

4=&3

oa

3

t3

•r4

1-1 bO

a) 4=

i
'f

cr13

Fig.

lbA

typical

mesh

cell,

with

center

located

atx=

(i1/2)

6x

andy=

(j1/2)

Sy,

illustrating

where

all

dependent

variables

for

the

finite

difference

scheme

are

defined.
lb

The following discretely evaluated functions will denote approximations

to the corresponding analytical solutions to Eqs . (13):

U.T; = U.(lcU, /"<*)

d

fij^ jo > (j-fcO-hj-

,

DT;^ , p*,j =

60 ^ = 60 (

* f*., £

/7t cT-6y)

where 5 x = AR/M and 6 y = l/N are the mesh cell sizes in the x- and ydirections respectively and where 6t is the time-step size. Such a

staggered grid is commonly used for multidimensional finite difference

integrations^

With this notation, the following set of finite difference equations

was used to approximate Eqs. (13):

U

m+ 1 M ~ | P J

<Jr 2 eft -

<f-ty

/ /K> ( it- - V "Jr

/-A+- 4

L

• • 2

/f^ 1 a-j - _ Scft 1 <dx

r

rt^ l /r--»

cr-

- P<ft

Oh

4-

- f-i J- - ^ ij

5 v

rr

+l - t-r-'=

^

t

xi

~l

U D - =

" \ i

/KV-

(£)', A O'- ^

9

-0^ 7

NJL (n d )

( /<uj
15

k

= mn

M A) H ZL Qa i ~ i 4~i

^ ^ rr (/“I s* <ft)

t*o)

Spatial derivatives in Eqs . (13) are approximated by central differences

and leapfrog is used for time derivatives. These equations have second

order truncation errors in both space and time; i.e., the local discretization errors made in approximating the differential equations (13) by

these difference equations is 0(5x ), 0(6y ), OCSt^) „ In the first three

of Eqs. (19), the time derivatives have been replaced by central differences which give explicit formulae for computing dependent variables

at a new time level, a difference method commonly referred to as the

leapfrog method. With this method dependent variables at two previous

time levels are required to compute the quantities at a new time level;

hence the scheme is said to be a three-time-level scheme. A starting

procedure is required to generate the dependent variables at the first

time level, given these quantities initially (at t=0)

.

To obtain dependent variables at the first time level, an explicit

first-order (in time) procedure is used:

i icu

I £ I fr)

i a

Sj, it po,^ P

L i)

\

li ~ ^ = ~ fr'V' - f~j)

+

-jfr (p/,^. + ?'j)] (a I c)

\

P

a

»i ( t

:+, 'i ~ 5 p° + fi-hj.)

+ _ s&Zk-r

><*-/ ^ i

<rt

(au)

I O n

p. — - <rt kf fa/e.
1

where all supplementary quantities are defined in Eqs. (20).

The computational procedure by which Eqs. ( 19 ) and (2l) are solved

is very important. At first inspection, it may appear that Eqs. (21) are implicit for example. However, when Eg. (21a) is solved first for

density. Eg. (21d) second for pressure and the remaining equations are solved subsequently, Eqs. (2l) are found to be explicit. Also, the time

level associated with the discretely evaluated variable representing presoire may not be the expected one. These considerations are discussed in

Section V.

The two forms for the ambient density stratification, given by Eqs.

( 16 a) and (l6b) in the continuous case, have been used in discretized form

in all of the computations presented in Section VI. For most computations

n exponential stratification was used

C~ (

i~

&.) ,] (=? * a)

We note that > when Y

s is taken to be large (generally 10

6

in our computations), the ambient density is essentially constant. For some computations, a discretized "two-layer" model of the ambient stratification

was used;

fpOM

- = ^

L - Q

- y*)fy/ys ] l-ji j~

s

In all cases, the initial velocities Ufj , UT

ij

and the density deviation

from ambient were taken to be zero. Boundary conditions for these difference equations are obtained from

Eqs. (l4) and ( 15 ) ; they are

i r

- /yv- /V1 - ^7, m — cJ if^

1
IT

and

$sNThe pressures and densities denoted by a tilde are values defined beyond

jthe grid defined in Figure 1. These values can be regarded therefore as

quantities defined in fictitious cells surrounding the rectangular region

of interest. However, in the procedure for solving the linear algebraic

system of equations for the pressure, the boundary conditions are incorporated directly into the system: the pressure and density values in the

fictitious cells cancel and never appear in the final system of equations

for the pressure in the rectangular region.

The procedure used to solve the discretized pressure equation is

discussed in Section IV. Before this discussion, however, we note that

the discrete pressure equation together with its boundary condition

retains the important property of the partial differential equation for

the pressure with its boundary conditions noted at the end' of Section II.

The definition of D<^ presented in Eqs. (20) shows that

This fact, plus that concerning the boundary conditions noted above, allows

one to demonstrate that, for the discretized pressure equation, the left

hand side and the right hand side, when summed over all grid points.,

each vanish identically. This important condition must hold for a solution

to exist to the linear algebraic system for the discretized pressure field;

it is discussed further in the next section (see Eqs. (29) and (30)).

M Nl = O
18

IV. THE SOLUTION OF THE PRESSURE EQUATION

At each time step in the fully nonlinear problem it is necessary

to calculate the solution to the non-separable elliptic equation

on the rectangle R, subject to the Neumann boundary conditions

= ?(*>$)

Here, the density £>, the forcing function f , and the boundary function g

all depend on time, but are known to us. (in the linear computation

reported here, <p(x, y) is replaced by ^po(y), which yields a separable

elliptic equation as noted in Section II. However, the solution

procedure is developed for the fully nonlinear problem. ) In the

discretization, Equations (25) are replaced by finite difference

equations, Eqs. ( l^d) or (2ld) and (23) in Section 3. The finite

difference equations (Hd) or (5tJd) then directly incorporate the

boundary conditions (23). Hence, the calculation of the pressure

requires only the solution of the linear equations represented by

these combined equations.

Physically Equations (25) produce a unique solution only when the

pressure is prescribed at some particular point. Equivalently, there is

a unique solution with mean pressure zero. Given any particular solution

we can generate a new solution by adding an arbitrary constant to the

pressure. The linear equations (l9d), (21d) behave analogously; the

linear system has a family of solutions such that there is a unique

solution with mean pressure zero, and the family of solutions can be

generated by adding constants to any particular solution.

The calculation of the pressure requires us to successfully meet

the constraints of the larger system of equations:

(a. 5 o.)
19

the solution method must take into account non-uniqueness ; the solution

method must he able to solve large linear systems accurately, since

there are NM equations (l9d); and it is very important that the s lution

be obtained quickly since the calculation is made at each time step.

The finite difference equations (l9d) can be conveniently represented

in matrix notation as

Ap = z

where p (p13_, P^, P-^, ... P^, P^, P

22 , ... P2M>

... ... P

ffl

)

is a vector of dimension NM, as it z, the vector of right sides from Eqs

.

(l9d) or (21d) (modified by boundary conditions (23)). The matrix A gives

the coefficients from the density and the finite differences. Here the

subscript T denotes transpose. Although A is of dimension NM by NM, it

has at most 5 non-zero elements in any row.

As noted above, the solution to Eqs. (25) is determined only to within

an additive constant. Similarly, the solution to the system of linear

equations (2o) is also determined only to within an additive constant.

The non-uniqueness of the solution is represented by

A (p + ae) = z T where e = (l, 1, 1, 1 ... l) ; a is any scalar and p any solution of ( 26 ).

This can be represented also as

|

Ae = 0.

| which is the statement that e is an eigenvector of A corresponding to zero

eigenvalue. The solution p with mean pressure zero is the solution of (26)

which satisfies

-T

p e = 0.

The existence of a solution requires the following equation to hold

z e = 0.

( 26 )

(27)

( 28 )

(29)

(30)
20

This condition on the difference equations (l9d) or (21d) , incorporating

boundary conditions (23b), is simply the statement that the right hand

sides sum to zero. As noted at the end of the last section, see particularly

Eq. (24) and the surrounding discussion, this condition is guaranteed by the

formulation of the linearized difference equations. It also must be

guaranteed when the nonlinear difference equations are developed. We assume

equation (30) and proceed to calculate the solution to the linear system

(26) which satisfies the mean pressure requirement (29).

Standard direct factorization methods (Gaussian elimination) for

2 calculating the solution to (2 6) are infeasible, because 0 ((NM) M)

operations at best are necessary. Iterative methods are often used to

solve sparse linear systems: the method of successive over- relaxation

(SOR) is probably best known. However, a successful implementation of SOR

requires at least a good estimate of the over- re taxation parameter. The

convergence rate even with an optimal parameter still yields 0 ((NM)N log N)

22

operations for a solution of (26).

The specific structure of the finite difference operation on the

rectangle R makes possible special fast direct methods for the solution

of separable elliptic equations. The algorithms discussed in Schwartztrauber

and Sweet"*''

1' are the most current extensions of the ideas behind the fast

finite Fourier transform. However, even these general algorithms require

separability (the coefficients of the derivatives with respect to x involve

only x; the coefficients of the derivatives with respect to y are functions

only of y), and the solution procedure which will be presented here is for

the case of interest in nonlinear problem. The pressure equation can be

transformed to separable form only if is known analytically and has two

continuous derivatives. Since ^

is known only on the staggered grid, these

conditions are not met.
21

The solution method we have adopted is a hybrid method which combines

an iteractive algorithm, conjugate gradients (denoted c.g.), with a fast

direct Poisson solver. The conjugate gradients algorithm provides an

iterative technique for solving Ap = z , provided A is a positive definite

symmetric matrix (equivalently, all eigenvalues of A are positive). Each

iteration of conjugate gradients requires one matrix-vector product, A

times x or roughly 5 MN operations. However, the convergence rate depends

critically on the ratio, k = X^/X^, the largest eigenvalue divided by

the smallest. The ratio is large for the linear system (26), and convergence is fast only if k is close to unity. To improve (drastically)

12 1 the rate of convergence, we write, following Concus and Golub,

A = A + A

s n where A and A are matrices chosen so that A ^z can be calculated quickly

s n s with a fast Poisson solver and the entries in A are small compared to A . n s The best choice of A

g is not known theoretically. Our solution is

to take

A = D LD

s where L is the finite difference matrix using a five-point difference

1 2 formula for the Poisson equation V P = f with Newmann boundary conditions.

Here D is a positive diagonal matrix chosen so that diagonal (A ) = diagonal

s

(A). The remainder A^ is a matrix with zero main-di agonal , and non-zero

entries present on only four off-diagonals. Such a choice is justified

on empirical grounds only if the conjugate gradient convergence rate is

I sufficiently fast. For our test problems this has certainly been the case.

There is also some theoretical basis for this choice; see Concus & Golub

^

.

Then we formally replace the equation

Ap = z

by

(D AD) a = Dz

(31)

(34)

(32)
22

where DAD is a new "scaled" matrix and Dz is the corresponding scaled

right hand side. We then solve the scaled problem using conjugate gradients.

The solution p of (26) is then Da.

If we have made a suitable choice for A , the eigenvalue ratio X /X

for the matrix D AD will be very close to the optimum value, unity. The

solution of (32) requires one matrix-vector product Ax, and the solution

of

A

s x = y

at each iteration. (The matrix D does not appear explicitly in the actual

algorithm. ) The details of this general approach are given in Concus , 13 Golub and O'Leary. Although the solution of (33) requires 0 (NM log N)

operations, the entire process will be efficient if we need to solve (33)

only a very few times. Typical calculations for our code require no more

than three to five Poisson solutions to obtain the solution of (26).

The only difficulty with the approach outlined above is the nonuniqueness of the solution. The matrix A which represents the finite

difference equations is only positive-semi definite , since e is an

eigenvector corresponding to the eigenvalue 0. However, all other

eigenvalues are positive. The iterations in the conjugate gradients

algorithm will generate a solution orthogonal to e if and only if the

iterations are begun with a vector orthogonal to e.

The natural solution produced by the method is the zero mean pressure

solution; any other solution can be produced from this particular solution

by adding a constant.

There is an additional complication in the hybrid algorithm since

both A and L, the discrete Poisson operator, are singular; Ae = Le = 0.

The iteration matrix D AD has a single zero eigenvalue with eigenvector

D

-1

e. In we show that for a more general splitting of matrix A,
23

the iteration matrix corresponding to D AD has the desired single zero

eigenvalue with both A and the matrix corresponding to D singular if and

only if the corresponding eigenvectors of A and the more general D are not

T -1 1

orthogonal. This can he assumed in our calculation since e D e = —

i

+ -— + . . . + —— and all the d. . are positive.

d

12

dNM

It is shown in [l4] that a choice of A

g with the same eigenvector of zero

eigenvalue as A is optimal in some respects, hut such a choice seems

difficult to realize. Our simpler choice, as outlined above, proves

entirely satisfactory. However, we must note that the search direction

in the c.g. iteration must now he orthogonal to D ^e, not e. This requires

that the solution to the discrete Poisson equation Lx = y he normalized

to mean pressure zero. This normalization must he done by the user of the

Schwarztrauher-Sweet routines^

1

and probably is also necessary for other

fast Poisson solvers.
24

V. Computational Procedure for Solving Linear Finite Difference Equations

In Section III the finite difference equations describing linearized

heating, Eqs. (I9)-(23), were presented, and in Section IV the procedure

used to solve the pressure equation was described. In this section the solution procedure for the complete set of linearized difference equations will

be presented and some aspects of the computation will be discussed.

Since the finite difference scheme, Eqs. ( 19 ) is a three level scheme,

values of velocities and density at two levels of time are required to

initiate the computation. Equations (21), a scheme first order in time.

p _ given initial data. Equations (.21; are solved in the following order.

First the densities p^. are obtained from Eq. (2la) . Then the pressures

u_^_. , v_ are computed from ( 21b) and ( 2 lc).

One point in this process is noteworthy. The initial pressure data.

pressure p are obtained using p_. Such a procedure was adopted because

our analytical studies of the finite difference scheme first order in time

demonstrated that it was most appropriate for reproducing the test problem

(see Reference 10 for details)

.

In subsequent time steps the order in which the difference equations

ij

p?. are obtained, using the p^., from Eq. (21d). Finally the velocities

1 1

provided to the computational scheme is used only as an initial "guess" for

the iterative pressure-solving

0

procedure: the initial values of the

1 „ , , , _ , t

are solved is changed somewhat. First, the pressure equations, Eq. (l 9 e)

for the mean pressure and Eq. (19 d) for the pressure variation, are solved.

Then the density and the two velocities are updated using Eqs. ( 19 a)-

(

19 c).
25

It should be noted that in both the first time step and subsequent

ones, the solution to the pressure equation assures that the velocity components will satisfy the proper divergence condition. Equation (5) is the

imposed condition on the divergence of the velocity in the continuous case;

this equation also applies in the linear problem, with the dimensionless

D(x,y,t) being specified by the fourth of Eqs. (9). The discrete version

of Eq. (5) is —1 f n+1 n+1 x 1 , n+1 n+1 x _n+l

(u. . - u. .)+ — (v.. -v. . .) = D

fix i,j x-1, j 6y ij 1,3-1 il

(35)

where ]>” is given by Eqs. (2 0 ) * This condition is very important since

it specifies, in the discrete (and linear) case, what effect the heat source

imposes through continuity upon the velocity field. Substituting for u?*\

v?T^ from Eqs. (19b) and (19c), and using Eq. (35) again, with n-1 replacing

n+1, yields the equation for the pressure, Eq. (25d) for all steps after the

first. Hence proper solution of this pressure equation assures that the

velocities at time step n+1 will satify the proper divergence condition,

Eq. (35). An analogous argument for the first time step, using Eq. (2lb),

(21 ) and (35), shows that the pressure equation (21d) assures the proper

divergence condition at the first time step.

When there is no imposed heating, so that = 0, Eqs. (19)

and (21) are homogeneous (as noted earlier). In Reference 10 the homogeneous

equations (19) with the first-order starting procedure implied by the homogeneous form of Eqs. (21) were solved analytically for internal waves in a

rectangular enclosure. Comparisons were made between this analytical solution to the finite difference equations and one generated by computational
26

procedures. These comparisons showed that the accuracy with which the

difference equations were solved computationally was determined by the

accuracy specified for the solution of the pressure equations (19d) (or

(2ld) on the first time step) . When the accuracy specified for the pressure

equation was 10 the analytically and computationally generated solutions

agreed to 1 or 2 parts in 10^.

Further comparisons were made between the analytical solution to the

finite difference equations and the analytical solution to the continuous

problem, the homogeneous form of Eqs. (13). These solutions differ by the

truncation error resulting from replacing differential equations by difference equations. The truncation errors in turn depend upon the spatial grid

size relative to characteristic spatial variations in the solution and upon

the time step size relative to characteristic temporal variations in the

solution. Internal waves in a rectangular enclosure were examined analytically

to test the finite difference schemes and computational procedure. Such

waves are expected once the fire has stratified the gas in the room.

A single internal-wave mode was specified initially (in either

continuous or discretized form) and its time history was calculated by way

of the differential or difference equations. The mode was specified by

two integers, the number of half wave lengths in the x-direction and the

number of half wavelengths in the y-direction filling the rectangular

enclosure and satisfying boundary conditions. Specification of the mode

then determined the eigenf requency , through the equations of motion, with

which this characteristic mode oscillated. The stability of the finite

difference equations was found to require that
27

rff- *

(‘fef /V*~. V ^ / ^

<%/ys

This condition requires that the time step size be a fraction of the BruntVaisala period for the ambient fluid, eigenf requencies being inversely

proportional to this period. Accuracy of the finite difference equations

was found to require that the mode number relative to the mesh-cell

number M in the x-direction be small, the mode number relative to the

mesh-cell number N in the y-direction be small and the time step size

relative to the Brunt-Vaisala period be small.

As an example, a 5 by 2 mode internal wave was calculated for an

ambient stratification parameter (see equation 16a) equal to unity.

For this stratification, the Brunt-Vaisala period is 2 tt in dimensionless

time units. A 31 x 30 mesh was employed, with 6t = 0.25. At a dimensionless

time t=10, the discrete and continuous solutions differed by about 3%.

During this time the discrete solution oscillated about the continuous

result in an irregular fashion. This behavior in the discretization error

is characteristic of the spurious computational mode introduced by the

leapfrog time differencing procedure. The continuous eigenvalue a which

determines the actual period of the mode under consideration is a=. 92807

for this example. The corresponding discrete mode eigenvalue o* for this

problem may be analytically calculated'*’ to be a*-. 923955. Thus, the

computational scheme should be able to follow oscillations sufficiently

accurately to avoid errors due to phase shifts between modes for times

t 0(100). Throughout the time interval 0 £ t < 10, the finite difference

equations were solved computationally to a few parts in 10^. This estimate

is based on direct comparison between the analytically obtained and numerically

computed results.
28

Last, we mention briefly computational times required. All computations reported in the next section (and those performed for results

reported in Reference 2) were carried out on the N.B.S. UNIVAC 1108. Most

computations were performed using M = 31 and N = 31 for the mesh and were

run for up to two hundred time steps. For a calculation performed

starting with a stratified fluid, between three and four iterations on

average were taken by the conjugate gradients portion of the pressure

solver and the computer time required per time step was between 4 and 5 seconds. For example, one computation of 70 time steps took 302 seconds

of CPU time and another took 3l6 seconds. For a calculation in an

essentially unstratified fluid (in which the elliptic equation for the

pressure becomes Poisson's equation), only one iteration was required and

about 2 seconds per time step was taken. In this case, a computation run

for 200 time steps required 402 seconds of CPU time.

Computations on smaller meshes required substantially less time. For

example, a computation performed for a stratified fluid using M = 7 and N = 8

required 66 seconds of CPU time to run 200 time steps. Another calculation,

performed for an unstratified fluid with M = 15 and N = l6 and also run

for 200 time steps required 122 seconds.

The computational times quoted do not include CPU time required to

produce the computer-generated plots discussed in the next section. When

plots were generated, the results were read to file during the computation

and processed later to produce the graphics. Computer time required for ,

the graphics will be discussed in the next section.
29

VI. Computer-Generated Results

Because Egs. (13) with boundary conditions (14) and (15) are

linear, they can be solved analytically (and these expressions evaluated

numerically) or by a combination of analytical and computational techniques. Unfortunately, neither of these procedures has been implemented

as a computer code at this time so that comparison of the results reported

in this section with results obtained by an alternate procedure cannot be

made yet. Presently, we are implementing one approach using a mixture of

analytical and numerical techniques, and we will make such comparisons

later

.

A. Graphics

From the finite difference equations (19), we see density, horizontal

and vertical velocity components and pressure are all computed dependent

variables at each time step. Typically in the computations we have run,

we have used a grid composed of 31 x 31 mesh boxes: hence at each time

we have 4 dependent variables, each being calculated at approximetely 1000

points. Clearly with such a large number of values computed, graphical

display is required to digest the results. Generally we have used two

computer plotting routines to allow us to visualize results; both routines

were developed by the Mathematical Analysis Division specifically for this

application and have been found to be essential for rapid assimilation of

results. Most plots shown in this report are direct reproductions of these

computer-generated plots.

Both plotting routines display the dependent variables on the rectangular

region of computation. One routine displays profiles of the horizontal and

vertical velocity components. The other plotting routine displays contours
30

of constant density (isopycnic contours).

The horizontal velocity is plotted as a function of vertical position

at prescribed horizontal locations, together with plots of vertical

velocity as a function of horizontal position at prescribed vertical

locations. (The locations at which horizontal or vertical velocity plots

are desired can be specified.) These superimposed plots give a clear

picture of the velocity field at any instant of time.

At a horizontal location, this plotting routine takes the values of

the computed horizontal velocity components at the vertical mesh points

and fits this data with cubic splines using a package "splines under

tension" written by A.K. ClineV * iince the vertical mesh points at which

the horizontal velocity components are computed is staggered with respect

to the grid lines covering the computational region, horizontal velocities

are not computed at the top or at the bottom boundaries. The last computed

horizontal velocity is one half mesh point from the boundary, and the

velocity profiles are extrapolated to the boundaries using the fitted cubic

splines (not under tension)

.

Similarly at the vertical location, the plotting routine fits the

computed vertical velocity values as a function of horizontal position with

cubic splines. However, for the vertical velocity, the normal derivative

at the left and right walls can be shown to be zero for the linearized

equations

.

Using this condition and the values of the vertical velocity at the

two points adjacent to each wall, a quadratic equation was used to extrapolate the value of the vertical velocity at each wall. Then the package
31

written by A.K. Cline was used to fit the vertical velocity data and the

end point information by cubic splines (not under tension) . An earlier

version of this package is described in references 15 and 16.

Finally, the velocity plots shown below were prepared with an interactive Fortran program using Tektronix Plot 10 Advanced Graphing II software and Terminal Control System (TCS) software^’ Hie display device is a

Tektronix 4014 display terminal.

The second plotting routine gives constant-density contours. It first

creates an intermediate code with an interactive Fortran using the NCAR

19

graphics software CONREC and DASHCHAR together with the system plot package

and the Tektronix Plot 10 TCS. Then the intermediate form is translated

19 using the NCAR portable graphics translator written by T. Wright and the

Tektronix Plot 10 TCS software. The plots are displayed on a Tektronix

40l4 display terminal. Typical times required to generate a velocity or

a density plot are five to fifteen seconds of CPU time.

B. Linear Heating in an Ambient Homogeneous Fluid

When the ambient fluid is homogeneous, (y) is a constant and Eqs.

(13) simplify considerably. For this case it is valuable to examine the

analytical form of the solution even though the solution is not calculated

for comparison with the numerical results. Several observations about the

nature of the flow field can be made from the analytical solution, and

these observations are useful discussing the numerical results.

When pQ is constant, the equation for the density in Eqs. (13) can be
32

solved explicitly using the separable form for Q (x,y,t) given in Eq. (17).

Acj.] - flU A

where

t>*Uvd '' Sr'

A vrr

and where the density perturbation from is initially taken to be zero.

The density is found to be a function of time multiplying a function of

space, and density can be expected to display exactly the same spatial

behavior at different times with only the amplitude changing from one

time to the next.

The velocity field can be decomposed into two fields, one derivable

from a potential the other from a stream function:

u. = £>Al

"&x

’ P'V'

where <p is the potential and ip stream function. The boundary conditions

on these two functions arise from the condition that the normal velocity

at a boundary be zero (no outflow at a boundary), Eqs . (10): they are

^

i t — ^ = 0 and i|) = 0 on the boundary, x = 0, and y

= 0, 1 (where — means

the derivative normal to the boundary).

The continuity equation, Eq . (5), implies

v a+ = h (x.,#., t) = fft) b* (*,#)

If we denote by <£ (*^ y J the solution to the equation

V

7 - D*( x >y)

with — O at and ^

~Oi i , then
33

We define the vorticity to be

C

v

=r a>r

Then

The equation for the time evolution of the vorticity is obtained by

taking the curl of the momentum equations in Eqs . (13)

1W - -J- \

L

The solution is

3 X

^

7

The initial conditions that u (x, y, 0) = v (x, y, 0) = 0 have been used

to determine that the initial vorticity distribution is zero.

If, as before, we denote by J (x, y) the solution to the equation

with |

= 0 at x = 0 , ^^5 and y

= 0 , 1 , then

The point to note from the analytical forms derived above is that

the temporal behavior of the velocity field obtained from the potential

is different from that obtained from the stream function. In particular,

when f (t) grows linearly,

(-t) oC t
34

and

Hence, initially the potential contribution to the velocity field will

dominate and the velocity field will grow approximately linearly with

time, whereas later the contribution from the stream function, or rather

due' to the generation of vorticity, will dominate and the velocities will

grow roughly as the cube of time.

In Figure 2-6 plots of horizontal and vertical velocity are displayed for a case of heating in an unstratified fluid. The heat source

is of the form given in Eq (17) with

Values of the parameters for the heat source shown are 3 = 24.5, X = 4,

x^ = 0.5 (cf Eq. (17)) and A= 0.1: hence the source is centered along

the bottom of the floor, its strength diminishing more rapidly laterally

than vertically.

Figure 2 shows, at time t = 0.1, horizontal velocity as a function

of vertical position at three locations (x = 0 . 25 , 0.48 and 0.74) and

vertical velocity as function of horizontal position at three locations

(y = 0.25, 0.48, 0.74). (Note that the heat source and the flow field

are symmetric about the centerline, but the horizontal positions used

for plotting are not.

)

Figure 3 shows the same horizontal and vertical

velocity plots at t = 0.2. Comparison of these two figures shows that

the velocities in Fig. 3 have the same spatial form, but are nearly

double those in Fig. 2. At "early time" during the linear calculation,
35

when these two plots were made, the velocities are dominated by the

potential flow as the analysis above indicates. The potential flow

represents an expansion of the gas due to heating by the source. At this

time the velocities grow linearly with time, and Figs. 2 and 3 confirm

this behavior.
36-

0.00

0.02

0.04

0.06

0.08
9.92

0.04

0.06

0.08
38

Figures 4 and 5 show these horizontal and vertical velocity plots

at t = 2.0 and t = 4.0 respectively. Note the change in the velocity

scales from Figs. 2 and 3. At these times, which are "late time" in

the linear flow, the flow field is dominated by the vorticity field for

which the physical significance is discussed below. The velocity plots

in Fig. 5 have the same spatial form but are approximately eight times

larger than those given in Fig. 4, a fact consistent with the analysis

above, which shows that the stream function (produced by the vorticity

field) grows as the cube of time.

are contributions to the velocities from both the potential and the stream

function; hence t = 0.7 is an "intermediate time" for the linear calculation

in which the velocity field is neither dominated by the potential nor by

the vorticity. (Note that Fig. 6 has the same scales for velocity as Figs.

2 and 3 .

)

An observation should be made at this point regarding the qualitative

behavior of the linear flow field exhibited in Figs. 2-6. The overall

nature of the flow field is what one would expect (and, for the later

times, qualitatively what people find in a room fire). During the "early

time" of the linear calculation as noted before (when the velocities are

dominated by the potential function), the flow is strictly outward from

the heated region as the gas is heated and expands. At "later time"

during the linear calculation, the flow is upwardly directed above the

heated region and, by continuity, down along the sides of the enclosure.

Qualitatively, this is the behavior observed in an enclosure fire. In

addition, in the region above the most intensive heating (approximately

y £ 0.17 here), the horizontal flow is mildly outwardly directed from

the region above the heat source (the plume region). Near the bottom of
39

0.9

9.2

9.4

9.6

U
40

9.9

9.3

0.4

9.6

9.3
0.C a

0.02

0.04

0.06

0.03
42

the enclosure, however, there is a relatively strong inflow. In room

fires a strong inflow or "floor jet" is often observed near the floor.

Such a qualitatively good picture from a linear calculation is very

encouraging.

Finally, in Figs. 7 and 8 plots of constant- density contours are

shown at two different times, t = 0.5 and t = 2.0 respectively. Two points

should be noted. First, as observed above from the analytical solution,

the spatial behavior of the density at the two different times is the same,

only the magnitude of the density perturbation has changed between the two

times. Second, we note that the spatial behavior of the density perturbation

is proportional to D (x, y) and hence directly reflects the nature of the

specified heat source Q (x, y)

.
C0NT01K

HW

>.80000-01

T«

.50000-02

C&Tt*

MEW*.

V

.SOCOO-W

1*T(3.S)*

.7«e»-02

l*»a3

5C«.t0

IT

45

C. Linear Heating in a Stratified Fluid

When the ambient fluid is stratified with an exponential variation

as given in Eq. (16a), a formal solution to Eqs . (13) can still be obtained

(because the equations are linear) . However, because of the complexity of

this solution, it is not very helpful for gaining an understanding of the

behavior of the flow field. Therefore, in this section results of some

computations solving the complete difference equations, Eqs. (19) and (21),

will be presented and discussed, but no additional analysis will be presented.

In all of the computations reported in this section, a reasonably

refined mesh, I = J = 31, was chosen. The enclosure was taken to be a

square, /'i = 1.0, and the gas constant = 1.4. Parameters describing the

heat source were 6 = 24.5, \ = 4.0 in Eqs. (19) and A = 0.1 in Eqs. (20).

The parameter e, specifying the relative error to which the pressure equation

is to be solved, was taken to be 10 ^

, and the dimensionless time step was

taken to be 6t = 0.2. All other parameters relevant to the individual

computations reported will be stated when the results are discussed.

An interesting observation can be made at this point regarding the

relative effort required to obtain numerical solutions (in the form of plots

say) as the complexity of the (linear) problem increases. In a previous paper,

Ref. 10, the authors obtained solutions to the equations for internal waves

in enclosures. In the continuous case, these waves are solutions to the

homogeneous form of Eqs. (13), and in the discrete case, they are solutions

to the homogeneous forms of Eqs. (19) and (21). Analytical solutions to

both the continuous and the discrete cases were obtained, and a numerical

solution was also obtained by the time-stepping procedure described in this

report. To generate a code and compute numerical values for the homogeneous
46

problem was very simple using the analytical solutions, but required much

effort for the direct numerical time-stepping procedure. Once the timestepping procedure was implemented as a code, it was a, easy task to modify

it to calculate the inhomogeneous problems with constant or variable ambient

density (described in the last section and this section respectively). To

calculate values from the analytical solution to the inhomogeneous problem

with constant background density in the continuous case requires some

effort, but not too much. However, to perform a similar task when the

background density is variable is a rather big effort. In short, as the

complexity of the problem, even though linear, increases, the effort required to obtain and evaluate analytical results may become much greater

than straightforward numerical integration by difference methods. In the

nonlinear case, of course, analytical solution will not even be possible.

Cil. Heating in a Stratified Fluid by a Centered Heat Source

Internal waves can be sustained in a stratified fluid. As noted

above these waves are solutions to the homogeneous form of Eqs. (13). The

heating source is the inhomogeneous, driving term for these equations, and

this source excites internal waves during heating. One of the most interesting features of the plots shown in this and the subsequent two examples

in this section are the internal waves generated by the heat source. With

each of the examples other observations will be madp.

In this first example the stratification parameter Y = 1.0. Hence

/ / ja

the buoyancy period T = 2tt {-g ~ = 2tt, and to observe anY effects

of internal waves, the calculations must be performed for at least a period

or two. This computation was run for 14 dimensionless time units with

6t = 0.2 so that 70 time steps were taken.
In Fig. 9 plots of horizontal and vertical velocities at dimensionless

time t = 4.0 for this case are presented. These plots should be compared

with those presented in Fig. 5: qualitatively, the plots are the same.

However, some quantitative differences can be observed. First, all velocities

are somewhat smaller in Fig. 9 (the stratified fluid) than in Fig. 5; the

difference is particularly noticable in the vertical velocities. This

effect is directly traceable to the stable background stratification which

tends to suppress buoyant vertical velocities induced by the heating.

Another difference appears in the horizontal velocities where the transition

between outflow from the center of the heated region and inflow occurs lower

in the enclosure due to the background stratification.

The next three figures, Figs. 10, 11 and 12, show the subsequent development with time (t =8.0, 12.0 and 14.0 respectively) of the velocity profiles. Note the change in the velocity scales from Fig. 5. It appears from

these plots that the vertical velocity changes little in magnitude over the

time period represented. However, the horizontal velocities change substantially during this period. The maximum horizontal velocity increases

with time. and the transition between outflow from the center of the heat

source to inflow occurs lower in the enclosure at each time. Other

distortions of the profiles as functions of time are also apparent, and

these changes are a result of the background stratification.

The nature of the internal-wave effects can be seen more clearly in

plots of the density perturbation (from the background exponentially varying

density) at times t = 2.0, 4.0, 8.0 and 12.0 in Figs. 13 - 16. In Fig. 13

constant perturbation-density contours at t = 2.0 are shown: clearly these

are qualitatively very similar to the perturbation density contours for the

unstratified case shown in Figs. 7 and 8. However, at later times the
1*8

density contours become qualitatively quite different, the most dominant

feature at any single time being that these contours are dramatically

flattened. Internal- wave effects become apparent by comparing plots at

different times. Such comparison shows that the point (or points) of

highest density, denoted by the letter H, sloshes back and forth within

the enclosure. Furthermore, examination of the zero-perturbation contour

shows undulations with increasing time.
0.2

0.4

0.6
I

- 50 -

0.0
51

X 3

rr

-ry

T0

e-IP

—

TO

EXIT

FROH

THIS

PASS

TYPE

f—•

T*

.CTDCO-Ol

CiNT0UR

INTEKTflL

IF

.80000-01

PTl3.JI =

.10339

LPTELS

SCHLfO

AT

1PO.CO
'ITr.R

Fp-t

- 1

.

000 0

TV

.30000

OWTII*

WTCTVHL

fF

.rr'OO

PT13-3I=

.inis
55

A —

M

rttrR

F~1

-16.000

TO

t.OCOO

CPfTOl*

INTERVAL

P*

I.OCCO

PTI3.3I:
F(&O^G

{&

56 .

CWTft*

Fm

-30.000

Tl

2.0001

CWTHH

IiTCTVflL W

2.0000

PTI3.3J:

-6.5539
57

C.2. Heating in a Stratified Fluid by a Noncentered Heat Source

In this calculation the center of the heat source was shifted to

x = 0.306 (compared with x = 0.5 in previous calculations). All other

parameters remained the same as in previous calculations, and the computation was run for 70 time steps, or to t = 14.0.

Figures 17, 18 and 19 show velocity profiles at t = 4.0, 8.0 and 12.0

respectively. (Note the change of scale ,of the plots between Fig. 17 and

Fig. 18.) The most distinguishing feature of these plots (and of this

calculation) is the asymmetry introduced by this noncentered heat source.

As in the previous case, the magnitude of the vertical velocities change

little, but the shape of these plots change somewhat as functions of time.

The horizontal velocity profiles change both in magnitude and in shape as

functions of time. The locations of the transitions between outflow from

the center of the source to inflow change as functions of time, occuring

lower in the enclosure as time increases.

Figures 20 - 23 show constant-perturbation density contours at t = 2.0,

4.0, 8.0 and 14.0 respectively. At t = 2.0, these contours are qualitative

similar, although displaced in the enclosure because of the noncentered heat

source, to those shown in Fig. 13. At later times the contours become flattened,

and undulations, as described in the previous case, occur.
58 ?

i.6 ®8. 1

••
59

A

«.•

«.2

t.4

o.r,

•.«
0.2

0.4

0.6
'I

I

I

- 61 -

A

trn

-

r—

-i.2?ot

ts

.kjooo-ci

intow*.

r

.icoao-oi

ptis-ji*

-.inn

mu

sc^ro

vr

m.ca
62

A

i

cpTM

H.etno

tv

.area

c*ttvi*

htowl

wr

.srao

-i.cn
63

r-i

-is.tro

t»

j.ccct

cmtw*

htct/*.

r

t.crrxj

rt

rstai*
C-.7-?

t

—1

-M.OOO

If

2.0000

CPU**

UrTETVT.

r

2.

con

PT13.3H
65

C.3. Heating in a Stratified, Two-Layer Fluid by a Noncentered Heat

Source

Q

In the final example all parameters are as described in the

previous case except that the background density stratification changes

abruptly at y

= 0.4, the fluid being of uniform density below this

location. Figure 24 shows plots of constant-perturbation-density at t = 6.0

in this calculation. This figure is shown only to demonstrate the qualitativelv different behavior which occurs when the background density is

composed of two layers, the upper one of variable density and the lower

one of constant density (homogeneous). In the lower layer the profiles

are qualitatively similar to what would be expected for heating in a

homogeneous fluid. The structure in the upper layer is qualitatively

quite different, showing what appears to be much detailed structure

superimposed upon some mild undulations. The detailed structure should

not be seriously considered because the magnitude of these variations is

small. On the basis of previous studies. Reference 10, the discretization

errors expected for the mesh M = N = 31 is of the order of a few percent.

(it could be somewhat larger because the background density has a discontinuity in its gradient between the upper and lower layers. ) The magnitude of the variations of the detailed structure is smaller than the

expected discretization errors.

67

VII. Summary

4

In an earlier paper the equations of motion for thermally driven

buoyant flows were derived. These approximate equations are applicable

to nondissipative, highly nonadiabatic , buoyant flows of a perfect gas in

which the heat is added slowly by a prescribed source. The equations

were derived by assuming that the time scale associated with increase

of the heat source is long compared with the transit time of an acoustic

signal across the spatial extent of the source. Such an assumption is

appropriate for buoyant flows produced by room fires. The resulting

equations admit internal-wave motions while filtering out high-frequency

acoustic waves. They are a generalization of the Boussinesq equations

valid for nearly isobaric flows with large density and temperature differences.

This combination of phenomena has not been extensively examined in studies

of fluid flows.

In the present paper these nonlinear equations were recast in a form

in which they could be integrated by finite difference techniques. In

__ /j 7 o

'< three dimensions the recast equations consist of three evolution equations,

for density and the three components of velocity, along with a linear,

nonseparable elliptic equation for the spatially variable portion of the

pressure. (The mean pressure as a function of time is obtained by solving

an ordinary differential equation) . The formulation in which these fluid

quantities are used as dependent variables is often refered to as a primitive

variable or primitive equation formulation.

The equations were then made dimensionless, specialized to two

dimensions and linearized. Finite difference approximations, second

order accurate in both space and time, were made to the linearized equations.

The dependent variables are density, the horizontal and vertical velocity

components and pressure, and these variables are defined at various positions
68

within a grid cell (see Figure lb)

The leapfrog procedure was used to determine the time evolution of

the difference equations. To start the scheme, a first-order accurate

scheme is used for the first time step. The computer time required for

various size computations "was reported.

At each time step, a system of linear algebraic equations, the finite

difference approximation to the elliptic equation for the pressure, must

be solved. The linear algebraic system of equations can be written in

matrix form, in which case the coefficient matrix for the vector of unknown

pressure values at grid points is large and sparse. It is important to

utilize storage carefully and solve this algebraic system efficiently and

accurately since the system is embedded in the larger fluid calculation and

it must be solved repetitively. A hybrid method, combining the interactive algorithm called conjugate gradients with a fast direct method

for solving Poisson's equation, was used to solve the algebraic system.

(it was found typically that the solution of the algebraic equations

arising on a 31 x 31 mesh could be determined to five significant figures

after between two and five iterations of the conjugate gradients algorithm

in between two to five seconds of CPU time on the N.B.S. UNIVAC 1108.

)

Finally, results of some computations were presented using computergenerated plots. A discussion was given of the software required to

generate these plots. A very brief discussion was presented of the

analytical procedures required to obtain a solution to the differential

equations describing heating of a homogeneous fluid. Computations of the
69

solution to the difference equations were found to confirm the conclusions

of the analysis. Also these computations are found to yield flow fields

which behave qualitatively the same as flow fields observed in a room

fire. Examples of computations of the flow fields induced by heating in

a stratified fluid were also presented and discussed.
70

References

1. Chandrasekhar, S., Hydrodynamic and Hydromagnetic Stability , Oxford

University Press, 1961.

2. Swarzschild, M. , Structure and Evolution of the Stars , Princeton

University Press, Princeton, N.J., 1958.

3. Turner, J.S., Buoyancy Effects in Fluids , Cambridge University Press,

Cambridge England, 1973.

4. Rehm, R.G. and Baum, II. R., The Equations of Motion for Thermally

Driven, Buoyant Flows , Journal of Research of the National

Bureau of Standards, Vol. 8j3, No. 3, May-June 1978, pp . 297-308.

5. Richtmyer, R.D. and Morton, K.W., Difference Methods for InitialValue Problems , Wiley-Interscience , New York, 1967.

6. Harlow, F.H. and Amsden, A. A., Fluid Dynamics, ALASL Monograph , Los

Alamos Scientific Laboratory Report LA4700, June 1971, Los

Alamos, New Mexico.

7. Haltiner, G.J., Numerical Weather Prediction , John Wiley and Sons,

Inc., New York, 1971.

8. Roache, P.J., Computational Fluid Dynamics , Hormosa Publishers,

Albuquerque, N.M., 1976.

9. Kreiss, H. and Oliger, J., "Methods for the Approximate Solution

of Time Dependent Problems", Global Atmospheric Research Programme

(GARP) Publication Series No. 10, February 1973.

10. Baum, H.R. and Rehm, R.G., "Finite Difference Solutions for Internal

Waves in Enclosures", National Bureau of Standards Report in

preparation.

11. Schwarztrauber , P. and Sweet R. , "Efficient FORTRAN Subprograms for

the Solution of Elliptic Partial Differential Equations", NCARTN/IA-109, July, 1975.

12. Concus, P. and Golub, G., "Use of Fast Direct Methods for the

Efficient Numerical Solution of Nonseperable Elliptic Equation",

SIAM J. Numer. Anal., 10, 6, Dec. 1973.

13. Concus, P., Golub, G. and O'Leary, D.P., "A Generalized Conjugate

Gradient Method for the Numerical Solution of Elliptic Partial

Differential Equations", in Sparse Matrix Computations , J.

Bunch & D. Rose, ed., pp . 309-332, Academic Press, New York, 1976.

14. Lewis, J. "The Numerical Solution of Pressure Equation by Hybrid

Conjugate Gradients",
71

References

15. Cline, A.K. , "Six subprograms for curve fitting using splines under

tension'). Comm. A.C.M. L7, 4 (April 1974) pp . 220-223.

16. Cline, A.K., "Scalar- and Planar-valued curve fitting using splines

under tension". Comm. A.C.M. 17 , 4 (April 1974) pp . 218-220.

17. Tektronix Plot 10 Advanced Graphing II User's Manual, Tektronix, Inc.,

September 1975.

18. Tektronix Plot 10 Terminal Control System User Manual, Tektronix, Inc.,

February 1977.

19. NCAR Graphics Software, edited by Tom Wright, NCAR Technical Note (preliminary edition)', April 1977.

20. Meijer, J.E., "Hydrodynamic Models for the Treatment of Reactor Thermal

Transients", Nuclear Science and Engineering, K) (1961), pp . 269-277.

21. Porsching, T.A., "A Finite Difference Method for Thermally Expandable

Fluid Transients", Nuclear Science and Engineering, 64_, (1977)

pp . 177-186; also Stewart, C.W. and George, T.L., "A Eulerian

Computation Method for Fluid Flows with Large Density Gradients

at all Speeds", Nuclear Science and Engineering, 6-4, (1977)

pp . 237-243.

22. Dorr, F.W., "The Direct Solution of the Discrete Poisson Equation on a Rectangle", SIAM Review, 12_, (April 1970), pp . 248-263.
NlJS-1 14A (REV. B-7BI

U.S. DEPT. OF COMM. BIBLIOGRAPHIC DATA

SHEET

1. PUBLICATION OR REPORT NO.

NBSIR 79-1754

m

4. TITLE AND SUBTITLE

A Linearized Finite-Difference Computation of Fluid Heating

in an Enclosure.

5. Publication Date May 1979

7. AUTHOR(S) Ronald G. Rehm, Howard R. Baum, John Lewis, and

Martin R. Cordes.

8. Performing Organ. Report No.

9. PERFORMING ORGANIZATION NAME AND ADDRESS

NATIONAL BUREAU OF STANDARDS

DEPARTMENT OF COMMERCE

WASHINGTON, DC 20234

11. Contract/Grant No.

12. SPONSORING ORGANIZATION NAME AND COMPLETE ADDRESS (Street. City. state, ZIP)

SAME AS ITEM 9

13. Type of Report & Period Covered

Interim

15. SUPPLEMENTARY NOTES

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

It. ABSTRACT (A 200-word or less factual summary of moat significant information . It document includes a significant bibliography or

literature surveyt mention it here.) In an earlier paper, approximate equations of motion were derived which are applicable to nondissipative, very nonadiabatic, buoyant flows of

a perfect gas. In the present paper, these approximate equations are rec-rt in a form in which they can be integrated by finite difference techniques. The nonlinear

equations are made dimensionless, specialized to two dimensions and linearized.

Finite difference approximations, using central differences in space and leapfrog in

time, are made to the linearized equations. To start the computation, a first-order

scheme is used for the initial time step. The dependent variables are density, the

horizontal and vertical velocity components and pressure, and these variables are defined at various positions within a grid cell. The computer time required for

various size computations and the accuracies obtained are reported.

At each time step, a large sparse system of linear algebraic equations, the finite

difference approximation to the elliptic equation for the pressure, must be solved.

A hybrid method, combining the interactive algorithm called conjugate gradients with

a fast direct method for solving Poisson's equation, is used to solve the algebraic

system efficiently and accurately. Examples of computations of the flow fields

induced by heating are presented and discussed.

17. KEY WORDS (six to twelve entries; alphabetical order; capitalize only the firat letter ol the liret key word unteea a proper name;

separated by eemicoiona) Buoyant convection; finite — difference computations; fire in an enclosure: heat addition: hyperbolic, elliptic equations; initial, boundary - value

problem; linearized equations jpartial differential equations.

18. AVAILABILITY [X] Unlimited

I I For Official Distribution. Do Not Release to NTIS

Q j Order From Sup. of Doc., U.S. Government Printing Office, Washington, DC

20402, SD Stock No. SN003-003-

\X] Order From National Technical Information Service (NTIS), Springfield, VA. 22161

19. SECURITY CLASS

(THIS REPORT)

UNCLASSIFIED

20. SECURITY CLASS

(THIS PAGE)

UNCLASSIFIED

21. NO. OF

PRINTED PAGES

76

22. Price

$$6. 00

USCOMM-DC