NAT'L INST. OF STAND & TECH R.I.C. Reference

AlllDl 77bflTD

AlllOS t.lfl701

NBSIR 81-2385

Finite Difference Calculations of

Buoyant Convection in an

Enclosure

Part I: The Basic Algorithm

Center for Applied Mathematics

Center for Fire Research

U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

Washington, DC 20234

100

.U56

81-2385

1981

U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

National Engineering Laboratory

Center for Applied Mathematics

Washington, DC 20234

HATIONAL BVWEJlD

09 nAMDASM

LIB1IART

NBSIR 81-2385 JWI 2 5 1982

no^

FINITE DIFFERENCE CALCULATIONS OF

BUOYANT CONVECTION IN AN

ENCLOSURE

PARTI: THE BASIC ALGORITHM

‘^:i

Howard Baum**, Ronald G. Rehm*,

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

Washington, DC 20234

**U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

National Engineering Laboratory

Center for Fire Research

Washington, DC 20234

December 1 981

Final Report

U.S. DEPARTMENT OF COMMERCE

National Bureau of Standards

National Engineering Laboratory

Center for Applied Mathematics

Washington, DC 20234

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

NATIONAL BUREAU OF STANDARDS, Ernest Ambler, Director
1 ' ••fjf •» : f^ ^'r
Finite Difference Calculations of Buoyant

Convection in an Enclosure

Part I; The Basic Algorithm

Howard R. Baumt , Ronald G. Rehm , 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

A novel mathematical model of buoyant convection in an enclosure, developed

earlier, is solved by finite difference techniques in the two-dimensional

case. This model has been developed as a principal analytical tool for the

prediction of the movement of smoke and hot gases in fires. Effects of large

density variations caused by substantial heating are retained while acoustic

(high-frequency) waves, which are unimportant to buoyant convection, are

analytically filtered out. No viscous or thermal conduction effects are included in the model. These two characteristics (filtering and no dissipative

effects) distinguish the model from all others describing buoyant convection.

The mathematical model consists of a mixed hyperbolic and elliptic set of nonlinear partial differential equations: the problem is a mixed initial, boundary value one. An explicit time-marching algorithm, second-order accurate in

both space and time, is used to solve the equations. The computational procedure uses a software package for solving a nonseparable elliptic equation

developed especially for this problem. The finite difference solutions have

been carefully compared with analytical solutions obtained in special cases to

determine the stability and accuracy of the numerical solutions. The computer

model has been used to compute the buoyant convection produced in an enclosure

by a spatially distributed heat source simulating a fire. The computed results

show qualitative agreement with experimentally observed buoyant convection in

enclosure fires.

i
TABLE OF CONTENTS

Page

Abstract i

I. Introduction I

II. Formulation 4

A. Continuous Problem 4

B. Finite Difference Equations 12

1. Basis for the Selection of the Difference Equations 12

2. The Equations 14

III. Example Calculations 25

A. Square Enclosure with Centered Heat Source 26

B. Rectangular Enclosure with Off-Center Heat Source 30

Appendix 31

References 34

Figure Captions 36

ii
I. Introduction

This paper presents the first results for a finite-difference integration of

an approximate set of equations describing buoyant convection in an enclosure.

The work represents a continuation of the research reported in Reference 1,

where the set of approximate equations was derived. The primary application

of interest to the authors is the movement of smoke and hot gases caused by a

fire in a room.

The research presented here is distinguished from previous numerical computations of buoyant convection in three respects. First, in this model the fluid

is taken to be an inviscid, non-heat-conducting perfect gas, and the spatial

and temporal magnitude and variation of the heat source, which simulates a

fire and drives the flow, are taken as known. These approximations are justified because under conditions characteristic of even a small room fire, the

Grashof numbers (representing the ratio of the inertial to viscous forces for

natural convection) are large enough for molecular transport phenomena to be

important only in wall boundary layers and in the highly convoluted flame

sheets which constitute the region of intense heat addition. The study of the

detailed flame structure of real fires is an extraordinarily complicated subject in its own right, and is bypassed here by specifying the heat source.

Wall boundary layers represent a local refinement to be considered separately

at a later date. Batchelor [2] gives a brief but relevant discussion of the

applicability of the inviscid equations in the context of atmospheric motions.

It should be noted that such simplifications do not preclude a description of

turbulence; but no turbulence model is explicitly included in this study. We

note, however, that any turbulence model appended to the present equations

must be of the "sub grid" variety, since no spatial or temporal averaging is

1
Implied in the equations derived in Reference 1.

The omission of any turbulence model is based on the observation, quantified by

McCaffrey [20], that most of the energy containing fluctuations in buoyant

plumes induced by laboratory scale diffusion flames are of low frequency and

large spatial extent. Such fluctuations, with frequencies typically in the

range 2 - 5 Hz and length scales comparable to the local plume width are

directly resolved by the computational procedure. A knowledge of the behavior

of higher frequencies and smaller length scales is of course crucial to an

understanding of combustion related phenomena. However, as noted above, such

questions have been by-passed in the present study.

Simple models of smoke and hot gas transport which neglect molecular transport

phenomena have been reasonably successful in predicting global properties of

flow fields [3] . The present work is intended as a first step towards more

detailed studies along these lines.

Second, the approximate set of equations integrated in this paper are characterized by the fact that large density variations due to temperature changes

are admitted, but compressibility effects are suppressed. Such a fluid has

been termed thermally expandable in other contexts [4]

.

In the fire setting,

allowance for density variations due to temperature increases during combustion is essential. It is common for temperature in a flame to exceed 1000°C,

implying that the density decreases locally to less than one-quarter its ambient value in the nearly constant-pressure process. In Reference [1], a set of

equations of motion was derived formally which permit description of large

density variations in a flow while ignoring acoustic oscillations arising because of the elastic properties of the fluid. Such model equations include

2
the important features of buoyant flows without requiring excessive computer

time necessary to determine high-frequency sound waves when numerically integrated. In this sense the equations "filter out" the sound waves while describing the lower frequency, organized motions due to buoyant effects such as

internal waves.

Finally, in previous efforts to compute flow fields, produced by buoyancy or

by any other mechanism, there are few reported detailed checks on the quality

of the numerical solutions to the finite difference equations. By contrast,

in References 8 and 22 detailed comparisons are made with analytical/numerical

solutions obtained to the general difference equations in simple, special

cases. Through these comparisons confidence in both the algorithm and its

implementation as a computer code was gained. Such tests showed that the

algorithm is stable and for the example cases performed, the error made in

solving the difference equations at any time step, even in the nonlinear case,

was over two orders of magnitude less than the discretization errors made in

approximating the partial differential equations by finite difference equations. For buoyant flows of the type considered in this paper, it is especially important to have confidence in both the stability and the accurancy

of the algorithm so that the real physical instability represented by the

fluid motion can be distinguished from any computational instabilities.

In Section II the equations derived in Reference 1 are recast into the form in

which they are solved, and the finite difference approximations to the equations are presented. In Section III solutions determined by this model, for

the buoyant flow produced by a heat source in two rectangular enclosures, are

presented and discussed.

3
II. Formulation

A. Continuous Problem

In Reference 1 the authors derived a set of nonlinear equations describing

nondissipative , buoyancy driven flows of a perfect gas. The flows were assumed

to be generated by a localized heat source in which the heat is added slowly so that the time scale associated with the heat-source grov;th and resultant

fluid motion is long compared with the transit time of an acoustic signal

across the spatial extent of the source. Flows Induced by a room fire generally satisfy this assumption. Properties of the equations were discussed in

Reference 1. In this section these equations are rewritten in a form appropriate for numerical integration by finite difference techniques, and the

boundary conditions for the equations are presented.

As in Reference 1 we consider an inviscld, 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 1. 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 1, we consider only a completely enclosed

room (no leaks), and when difference equations are Introduced, we confine

attention to the two dimensional evolution of the flow.

The continuity, momentum, energy, and state equations are given respectively

in Reference 1 as:

A
8p 3 — + (puj[) = 0

9t 3xi

3ujl 9u^ 3(p-pQ(t))

+ U4 ) + pk^g = 0 3t'’3xj/ 3xj^ (1)

/ 3T 3T \

( — + ) Q(xi,t)

^ 3t ^ 3xj/ dt

Po(t) = pRT

th

Here p is density, Uj^ the velocity in the i coordinate direction

(i = 1, 2, 3), p is the pressure excess above the mean pressure Po(t) in the

room, T the temperature, Cp the constant-pressure specific heat, R the gas

constant, k^g is the gravitational acceleration (of magnitude g) and Q(xj^,t)

the specified volumetric heat source* The spatially uniform mean pressure

PqCp) depends only upon time and increases because of the heating within the

the room. It is determined in a completely enclosed room by the equation

^Po y-1

=

/ 0(xi,t)dV (2)

dt V V

where y 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.) It will also turn out to be

a mathematical consistency condition required if a solution for the pressure

field is to exist.

5
Equations (1) and (2), are the approximate set of nonlinear equations solved

by finite difference techniques in two spatial dimensions in this paper. The

equations admit buoyant or internal-wave motions while "filtering out" highfrequency, 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 considered here were open or if the mean pressure variation were comparable to

the spatial pressure perturbation) . To recast the equations into a form more

suitable for numerical computation, 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

3p 3p 3^i — + U£ = - p = - pD(x-;,t) (3)

3t 3xj[^ 3x^

where

D(xj ,t)

1 YPo(t) [

(Y-l)Q(xj ,t) - (4)

Equation (3) and the continuity equation identify D(x^,t) as the divergence

3uj[ = D(x4 ,t) (5)

3xj^ ^

Finally, as in Reference 1, the equation for the spatially variable portion of

the pressure is obtained by dividing the momentum equations by density and

taking the divergence of these equations. The resulting equation is

3 3x£

3 / 1 3p

Sxj’ 'p 3 x£

3ui , 3D(x-j ,t) + —__ ( 6 ) 3x4 3t

6
Equation (6) is the generalization for a "thermally expandable" fluid, which

we consider here, of the well-known incompressibility condition in Boussinesq

fluids. When the density is constant, then D(x,t) = 0, and Equation (6) reduces to Poisson's equation. The boundary conditions on these equations are

that velocity normal to any (impermeable) wall vanish.

u^nji^ = 0 (7)

where n;j^ are the normal components of a vector describing the boundary walls.

From Equations (1) and these conditions, the appropriate boundary conditions

on the pressure equation are obtained

3p

n^ = pg n^kj^ (8)

8ni

An important observation must be made concerning Equation (6) and Neumann

boundary conditions (8). When Equation (6) is integrated over the total volume of the room, both the left hand side of the equation and the first terra on

the right are divergence forms and can be converted into integrals over the

boundaries of the room. Application of boundary conditions (7) and (8) show

3D

that each of the terms is zero; therefore, the integral over the volume of —

3t

must also be zero. The requirement that the integral of Equation (4) for D

over the volume be zero produces Equation (2), the condition for the spatially

uniform background pressure. Therefore, the elliptic equation. Equation (6),

for the pressure, with the Neumann boundary conditions, is seen to produce a

condition which must be satisfied for a solution to the equation to exist.

This condition, Equation (2), determines the time evolution of the spatially

uniform background pressure and demonstrates that the total pressure can be

consistently separated into a spatially uniform background pressure and an

7
inhomogeneous time dependent over-pressure. In the next section describing

the difference equations, exactly analogous considerations are found to apply

to the linear algebraic equations approximating Equation (6) and boundary

conditions (8) .

For selecting a difference scheme, the second of Equations (1), the momentum

equations, are rewritten in vector invariant form, noting that uj[ are coraponents of the velocity vector field u

->•

1

-»• — + 1/2 V (q2) - uxto = -- V p+kg (9)

9t p

where q^ = u*u and w = V x u is the vorticity. The curl of Equation (9)

yields the vorticity-transport equation

^ /lx — - Vx(uxo))=-V (-) xVp

3t ' P

^

( 10 )

Since u is a vector field, it is necessary to calculate correctly both its

divergence and its curl. Equation (5) specifies the divergence of u, and the

equation for the pressure. Equation (6), assures that the divergence is

properly determined at each time. Equation (10) is the equation for the evolution of the curl of the velocity, the vorticity. The difference scheme

selected must assure that Equation (10) is satisfied in difference form.

The complete set of recast nonlinear equations are gathered and rewritten

8
below

9p 9p — + U4 = - pD(x- t)

9t 9xi

where

dpo

dt

Hi

V

/ Q(xi^t) dV

V

D(xi^t)

1 YPo(t)

dPo

(y-l)Q(xi^t)

is the permutation tensor and

9%

o>i = Ei-jii ate the components of the vorticity vector.

9xj

Boundary conditions are

Uini = 0

ni —- = p g ki ni

9xj^

(lla)

(llb)

(llc)

(lid)

(lie)

(12a)

(12b)

9
For numerical integration, Equations (11) and boundary conditions (12) are altered in two additional ways. First, we calculate the density and pressure

differences from their initial values, which may be functions of the vertical

coordinate. This is done to eliminate accuracy problems, since the thermally

induced density and pressure differences can be very small during the early

portion of the heating process. Second, the equations are made dimensionless.

The nondimenslonalization is done both for convenience and to ensure proper

scaling. All dependent quantities are made to be of order unity in magnitude

for purposes of the computation. The remaining dimensionless parameters

characterize the strength and location of the heat source as well as the room

geometry.

In Figure 1 a schematic diagram of a fire evolving a room and a set of coordinate axes are shown. It is assumed that initially the enclosure is filled

with quiescent, stratified fluid of density Po(y)» where we denote xj = x,

X2 = y and X3 = z. We define a density difference from ambient and a pressure

difference as follows:

p (x,y,z,t) = p(x,y,z,t) - Po(y) (13a)

y

p (x,y,z,t) = p(x,y,z,t) - Po(t) + g / Po^y*) ‘^y' (13b)

o

These differences p and p need not be small compared with Po(y) ^ad Po(t)

respectively. Then Equations (lla)-(llc) become

3p 3p dpQ — + Uj[ + V = - (Po(y) + p) D(x:j t) (14a)

3t 8 xi dy ’

10
= " (l/(Po+p)) (lAb)

3ui 3 +1/2 (u.u.) - at axi

9p

3X;{^ - k^gp

(”V 8x£

~ 9p \

(l/(Po+p)) ) ax^/

a ~ ~ 3D 3^ . — (gp/(Po+P)) + — + 1/2 r (q^)

8y . at 3x^2

3

8xj[

where v = U2 and

and the boundary condition (12b) becomes

(14c)

3p

Hi = p gki n^ (15)

8x£

Finally, we form dimensionless equations using the density Pqq = Pq(0), the

height of the room H as the length scale and the free fall time (H/g)^/^ as

the time scale. Then, denoting dimensionless quantities with a hat

A. A »N* P ~ P/PoO » P “ p/poo^n » pQ ~ Po'^PoO

(16)

A A u^ = u^/ZgH , Xji^ = Xj^/H , t = t// (H/g)

Equations (13) remain exactly the same in dimensionless form with g set equal

to one. Subsequently, in this paper all quantities will be understood to be

dimensionless, and the hat notation will be dropped. For the dimensionless

coordinates, we note that 0 ^ x l/AR , 0 ^ y ^ 1 and 0 ^ z _< 1/BR where

AR = H/L and BR = H/W, Also, in the remainder of this paper, the problem will

be specialized to two spatial dimensions so that all quantities will be assumed independent of z.

11
This non-dimenslonallzation, while simplifying the form of the equations, does

not relate the scale of the induced motion to that of the source. An alternative scheme which does have this feature is derived in the Appendix. The

resulting equations are somewhat less convenient for numerical computation,

except in the Boussinesq limit when density differences are small. Hence,

the variables as defined by Equation (16) will be used throughout the remainder of the paper.

B. Finite Difference Equations

1. Basis for the Selection of the Difference Equations

In this section the finite difference equations and the boundary relations for

the solution-algorithm are presented. This algorithm does not represent a new

numerical scheme; rather it represents an adaptation of several well-known and

proven ideas to the equations of interest here. Since these equations are

different than those generally examined in computational fluid dynamics, there

are interesting features which vrill be noted as the difference equations are

presented

.

It is important first, however, to state the requirements which we used in

selecting the scheme. These requirements do not necessarily specify a unique

scheme, but restrict greatly the selection. As noted in the previous subsection, the momentum conservation equations provide relations for determining

the velocity field, which is a vector field, as a function of time. Forming

the difference equations from the vector invariant form of the momentum equations is one criterion. This choice is made to assure that, when the discrete

analogues of the divergence and curl are applied to the difference equations.

12
suitable discretized equations for the pressure and the vorticity transport

are obtained. Thus we can be certain that the divergence and the curl are

calculated correctly in discrete form, the divergence producing the equation

for the pressure and the curl producing an evolution equation for vorticity.

A second criterion imposed is that the difference equations provide secondorder-accurate approximations to the partial differential equations. This

condition is imposed because second order accuracy is necessary to obtain

reasonable spatial and temporal resolution for the hundreds of time steps required to calculate the complete evolution of the flow field in a room fire.

A well-understood buoyant flow field is that produced by internal waves in an

ambient stratified environment. Internal waves have been analyzed and calculated for a long time [6,7]; they are determined as solutions to the linearized partial differential equations of buoyant flow. In addition, internal

waves can be expected to arise naturally in a room fire setting when the

driving fire has heated and stratified the room air; after the fire has

extinguished itself because of vitiated air, for example, pure internal wave

modes will exist. An additional criterion which was imposed on the selection

of the difference scheme was that it accurately reproduce internal-wave modes

in an enclosure: if the difference scheme does not reproduce this rather

simple linear flow field, it is unlikely to reproduce more complicated flows.

In Reference 8 the authors examined internal waves and some second order

linear difference equations that reproduce these waves. This analysis

determined the difference scheme for all but the non-linear convective terms

in the momentum equations.

13
In an important paper for numerical weather forecasting [9], Arakawa discusses

design of computational schemes for long-term numerical integration of the

equations of fluid motion in the case of two-dimensional incompressible flow.

The major thrust of his paper is a derivation of second-order accurate spatial

difference schemes which eliminate the nonlinear computational instability

first noted by Phillips [10]. Arakawa emphasizes that for two-dimensional,

incompressible flow, the discrete approximations to quadratic forms of dependent variables, such as the velocity squared or the vorticity squared (or

both), must be conserved when the continuous variables are, and he uses this

constraint to select three difference approximations which are acceptable.

Another criterion imposed in the selection of our difference scheme is that

it reduce to one of Arakawa* s acceptable, or stable schemes with respect to

the nonlinear computational instability when the flow is incompressible. The

scheme chosen is denoted J3 by Arakawa [11] and conserves energy in the incompressible case. It can be obtained by differencing the vector invariant

form of the momentum equations.

Finally, we v;anted the difference scheme to be easily generalized to threedimensional flow configurations. For a model of buoyant flow driven by a fire

to be successful, it must be able to calculate three-dimensional flows. The

scheme selected and presented below satisfies all of the criteria stated

above

.

2. The Equations

In Figure 2a, the two-dimensional rectangular enclosure in dimensionless

14
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 are I mesh cells in the x-direction and J mesh cells in the y-direction.

Upon this basic mesh, the two components of the vector velocity (u, v) and

8 V 3u

single surviving component of the vector vorticity w = — - — are defined.

3x 8y

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 p and pressure p are defined. In Figure 2a the densities in the lefthand column of cells and in the bottom row of cells are shown to indicate how

they are enumerated for the numerical computation.

In Figure 2b a typical mesh cell is shown. Illustrating where all of the

dependent variables in the finite difference scheme are defined relative to

the cell.

The following discretely evaluated functions will denote approximations to the

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

n

u = u(i6x, ( j-l/2)6y ,n5t)

ij

V = ((l-l/2)6x, j6y, n6t)

Ij

p = p((i-l/2)6x, (j-l/2)6y, n6t)

ij

13
p = p((i-l/2)6x, (j-l/2)6y, n6t)

ij

D = D((i-l/2)6x, (j-l/2)6y, nfit) , ij

n

w = w(i5x, i6y, nSt)

ij

where 6x = 1/(I»AR) and 6y = 1/J 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

[ 12 ]

.

With this notation, the following set of finite difference equations was used

to approximate Equations (11) and boundary conditions (12):

For Equation (14a), l<i<I, l<j<J and n > 2,

-n+1

) ij

1

1+(1/2)D 6t

ij

!

~n-l / \ p (l-(l/2)D 6t)

ij V ij /

~26t

n n n

F + F + (1/2)D Po(j)

pxij pyij ij

where

(18)

n

p = p - Po(j) = the density difference from the initial density,

ij ij

Po(j) = exp[-( j-l/2)6y/Yg] = the prescribed ambient initial stratification,

(19)

Yg = the stratification length scale.

16
n n The flux terras F and F for l<i<I, l<j<J are given by

pxij P7ij

n

F

PXij

n

~n ~n \ / n

“

i,j+l ij i.j-1

pyij 6y

~n ~n \ / n n Po(3+l) - + P “P \ /

^ i,j+l ij l>j-l

26y

(20b)

Equation (18) employs a modification of the second-order accurate central difference (leap frog) temporal discretization. The modification eliminates an

instability that would arise if the leap frog scheme had been applied. It

affects the undifferentiated term pD(x,y,t) in equation (14) that is well

known to lead to a computational instability for ordinary differential equations when leap frog differencing is used [13],

For Equations ( 14b)

n+ 1 n1

u = u - 26t

ij ij

P - P

i+l,j ij

(21a)

.n 2p„( j ) + P + P i+l,j ij

17
and for KKI-l, Kj<J ,

n+1 n-1

' n 2 V = V - 26t < F + —

ij ij

I

•H 6y

(21b)

for Ki<I, Kj<J-l .

n n The fluxes F and F are defined as follows;

xij yij

for KKI-I, Kj<J

n 1 xlj 26x

and for Kj<J-l, Ki<I

n \2 2 1 / n n n n V 03 + V 03 2 \ 03 . , ij w . . , i

ij i,J-l ^

(22a)

n 1 yij 26yn

i.i+1

n n n

U 03 + u

03

n

03 ij ij ‘‘^i-l.j

(22b)

and where

n n n

n

ij

V. , ,

. - V. . U. . , 1 - u. . i+l,j ij i,j+l ij

5x 6y

n In n n ln n V =-(v +v ) ,u =-(u +u)

03 2 ij i+l,i Wj. 2 i»j+l ij

ij ij

(22c)

n 2 (q ) =

ij

n n

u. . + u. ,

. ij i-l,J

n n

V. . + V. ij i.j-1

18
Note that boundary conditions (12a) on the normal velocities imply that

UQ^j = =0 for l<jO and = 0 for Ki^I» These boundary condin n tions are applied formally in the expressions for the fluxes F , F , pxij pyij

n n

F and F in mesh cells adjacent to boundaries,

xij yij

The finite difference analog to Eq, (14C) is, for Ki_<I and KjO

2

p

i+l,j

~n

• P

ij

^n

P

ij

~n

P

i-l,j

6x^ ^n ^n ~n ^n

2Po(j) + P + P 2po(j) + P + P ^ i+1 , j ij ij i-l»j-'

2 +

6y2

~n ^n

P “ P

i,j+l ij ~n ^n

Po(j+l) + PoO) + Pi j+i + Pij

P - P

ij i,j-l

~n

Po^j) + Po(j-l) + Pij + Pi,j -1 -

n+1 n-1 n n n n D-n F -F F F

ij ^°ij x,i-l,j x,ij yi,j-l “ yij

+ +

26t 6x 5y

1

<Sy

^n ^n

P + P i.j+1 ij

~n ^n

Po(j+l) + Po^j) + P + P i,j+l ij

^n ^n

P + P ij i»j-l

~n ^n

PoO) + PoO-1) + P + P ij i,j-l'^

(24)

19
This difference equation for the pressure arises formally by applying the

finite difference analogue of the divergence operator to Equation (21) and

noting that the finite difference divergence of the velocity field satisfies

the equation

n n n n ^Ij

" ^i-l,j

_ ^ij

" ^i,j-l

-J-

' 6x 6y

(This equation is the difference approximation to Equation (5).)

The boundary conditions (15) in discrete form become

= D (25)

ij

~n ^n

P = P

0,j l,j

for KjO

~n ^n

P = P

1+1, j

(26a)

~n

P = 5y

i.O

^n ^n

P + P

. i,l i,0

/2

for Ki<I

~n / ~n ~n \

P - p = - 5y P + P 1 /2

1,J+1 1,J V i,J+l i,J /

(26b)

~n ~n

Although p and p appear in the boundary conditions, Equation 26(b),

i,0 i,J+l

they also appear in Equation (24) with j = 1,J in the same combination. As a

~n ~n

result, p and p never need to be specified to obtain a solution to

1,0 i,J+l

20
Equation (24) and the boundary conditons (26). Equations (24) together with

boundary conditions (26) constitute a singular linear algebraic system of

equations. When Equations (24), with boundary conditions (26) incorporated,

are summed, the left hand side sums to zero, demonstrating that all of the

equations are not linearly independent. Also, the last three terms on the

right hand side sum to zero, producing the requirement that the double sum n+1 n-1

Dij - Dij

over must vanish.

26t

The vanishing of the sura of the left hand sides of Eqs. (24) and (26) is the

discrete analog of the fact that the left hand side of the integral of Eq. (6)

over the room volume vanishes when the boundary conditions u • n = 0 are

applied. The vanishing of the sum of the right hand sides of Eqs. (24) and (26)

is the corresponding discrete analog for the requirement on the integral of the

right hand side of Equation (6). This requirement, that '['l must vanish at

each time level, is the discrete analog of Eq. (lid).

n Examination of Equation (28) for shows that it has been chosen so that its

double sum over all mesh points vanishes, and that the condition which must be

satisfied to allow this choice produces Equation (29) for the mean pressure.

Therefore, the singular linear algebraic system is seen to be consistent and

thus to allow a solution. The solution is made unique by specifying that the ~n double sum over all mesh points of pj^^j is zero. This is tantamount to~n specifying that Pq is literally the mean pressure in the room, with p^j the

perturbation about the mean. Details of the algorithm used to solve Equations

~n

(24) and (26) for pj^j are presented in Reference 17.

21
The heat source has been chosen to have the form

n n

Q = 0 f

ij ij

Q = {-) X exp [-3(xi“Xc)^ ~ Xy^]

ij ^

(27a)

(27b)

xi = (i-l/2)6x , Yj = (j-l/2)6y (27c)

n n

f = Qo tanh At

o

t = 0 nt n-1

I n

* =0

Hence, the discrete divergence of the velocity field becomes

^1 " n

D = [(Y-l)Qii“K]f

ij YPo”

(27d)

(27e)

(28a)

where

K = Y-1

IJ

I

I i=l

J -

I Qij

j=i

(28b)

and the mean background pressure is found from the difference equation

n+ 1 n-1 n n

p = p + Kf 25t

o o (29)

0 1 o with p = p =1 since f = 0 . o o

22
Since difference equations (18) and (21) are three-level, second-order schemes

(leap-frog) in time, a starting procedure is needed. The following firstorder, explicit scheme was used to start the computation and to restart when a

change in the time step has been made:

n+1 n 1

L = u - 6

{ ^

ij ij 1 xij

i+l,j

-n 2Po(j ) + P +

i+l,j

(30a)

n+1 n n 2 = V - 6t F + —

ij ij yij 6y

^n ^n /^n ^n \ 6y

P

iJ+1

-P

ij

+ P +P — \ i,j+l ij/ 2

~n ^n

Po(j+0 + Po^j) P P i,j+l ij

(30b)

~n+l ~n

P = P ij ij

~n

Po(j) + P ij

nD

ij

n + F

Pxij

n + F

Pyij

(31a)

n+1 n n

p = p + Kf 6t

o o (31b)

n+1 ~n+l ~n

When starting. Equations (31) are used to obtain p and p . Then p

° ij ij

n+1 n-1

Dij -

is obtained from Equation (24) with replaced by

25t

n+1 n

i*ij “ i>ij

5t

23
~n n+1

With this solution for p , Equations (30) are used to obtain u and

n+1

V . The starting procedure is completed (and two levels of all dependent

ij

variables have been obtained) by solving Equation (24) with n replaced by

n+1 throughout. Subsequent time steps are taken in a straightforward fashion

with the density and velocity components being advanced through Equations (18)

and (21), and the pressure being updated through Equation (24).

The linear stability of the algorithm is the only other consideration for

discussion. A linear stability analysis of Equation (18) for the density

shows that the time step 6t must satisfy the following condition for

stability

where

n

2 6t £ max +

KKI (V ij/ 6x

Kj<J

n n \

U = (1/2) (xx + u

) ij V ij i-1 »j /

n n

V = (1/2) (v + V

)

ij ^ ij i.j -1^

6y

(32)

When the stability condition, Equation (32), is not satisfied by a time step, n the time step 6t is halved. Then the time-marching algorithm is restarted

using the last time-level values as Initial conditions. A first-order time

step is taken and then leap frog is resumed.

24
When the finite difference analogue of the curl (see Equations (22)) is

applied to the difference equations for the velocity components, Equations

(21), a discretized form of the equation for the vorticity transport. Equation

(10), is formed. A linear stability analysis of the difference equation

yields exactly the same form for the stability criterion as that found above

for the density equation. Reference to Figure 2b shows that the density and

vorticity are evaluated at different points in the mesh, however, and therefore, the divergence D and the velocity components U and V are to be evaluated

at different points than those used in Equation (32), To account for the difference in the stability criterion implied by the different mesh location

points, in all calculations performed using the algorithm described above, the

time step was chosen to be less than or equal to 0.8 the maximum value found for

the right hand side of Equation (32).

III. Example Calculations

The algorithm, described in Section II and tested as discussed in Reference 22,

has been used to compute solutions to the buoyant-flow equations. In this

section results of two computations are presented and discussed. One calculation is for the flow generated by a heat source centered along the floor in a

square enclosure, and the other is for the flow generated by a heat source with

maximum along the floor, but off-center in a rectangular enclosure, twice as

long as it is high. Other calculations, exploring parametrically the features

of this model and examining in detail the numerical results by analyzing the

flow data computed by the model, will be reported in a companion paper.

25
A. Square Enclosure with Centered Heat Source

In Figure 3 contours of constant temperature (isotherms) are shown at dimensionless time 2.0 for a volumetric heat source centered along the floor in a

square room. The rate of heat added per unit volume is largest along the floor

at the center of the room and decreases in a Gaussian fashion with horizontal

distance from the center and exponentially with height above the floor: the

dependence of the heat source upon position in the room is given by Equation

(27) with = 0.5. The heat source is "turned on" slowly according to Equations (27a) and (27d) asymptoting to full strength around t = 10.0. At this

early time the problem is still linear; the flow velocities are sufficiently

small that convection is unimportant, and the temperature increase in the fluid

is directly proportional to the volumetric rate of heat added. Therefore, the

isotherms are also contours along which the volumetric heat-addition rate is

constant. (These contours can be seen to be parabolas by examination of Equation (27b), which describes the spatial dependence of the heat source selected

for these computations.) These computations were performed on a spatial mesh

of I = J = 31; the tick marks along the boundary of the enclosure show the mesh

cell spacing.

In Figure 4 isotherms at dimensionless time t = 10.0 are shov7n. By this time

the flow-field is nonlinear, and the temperature profiles are severely distorted due to buoyancy effects. The temperature has increased and the density has

decreased where heating has occurred. The heated fluid has become lighter than

its surrounding and begins to rise due to buoyancy. By continuity, surrounding fluid begins to be drawn into the region of the heat source, and the

26
isotherms, therefore, appear to be pinched off at the bottom (near the center

of the heat source). Two vortices of equal and opposite strength located on

the two halves of the heat source have developed to the size where their combined flow field is significant.

At time 11.5, Figure 5, a buoyant thermal has developed, giving the appearance of a mushroom cloud. The two vortices mentioned above have begun to rise

with the fluid, being convected out of the region of primary heating. The

buoyant thermal Intensifies in strength as shown in the next two plots at times

12.5 and 13.5, Figures 6 and 7, until the thermal hits the ceiling, as shown

in Figure 8, time 14.5, and begins to spread. Inside the plume a distinctly

periodic structure has begun to develop, as can be seen vividly in Figure 7;

here, progressing up the plume along its centerline, one finds a local low

first, then a periodic sequence of local highs and lows up to about the center

of the head of the thermal.

The heated gases are seen to spread along the ceiling in Figure 9 (time 15.5)

and fill the room from the top down, as shox-m in Figures 10, 11, and 12. This

physical behavior is exactly what is observed in room-fire tests and in other

experimental observations of heating in enclosures. The S 3nnmetry about the

centerline of the room displayed in these computations is some measure of the

accuracy with which they were performed: the heat source is placed S 3mmetrically, but the computations were performed as if no symmetry existed.

To assess the resolution of the computed results shown in Figures 3-12, the

flow field was computed again using a larger number of mesh points I = 63, vT =

64. Selected plots from this larger computation are shown in Figures 13-18.

27
These plots demonstrate that the large-scale features determined by the 31

X 31 computation are correct and agree with those determined from the larger

computation to within about ten percent. The results from the larger computation are characterized by smoother isotherms and more detailed structures because of the greater resolution. We note that the time required for the

buoyant thermal to reach the celling is about ten percent longer in the

63 X 64 results than in the 31 x 31 results, again apparently because of the

greater spatial (and temporal) resolution of the larger computation.

Detail on a length scale of the order of one or two mesh cells must be disregarded because the computations cannot resolve such detail. On the other

hand, features with a larger length scale can be interpreted. The spatially

periodic behavior in the plume noted above is a feature which requires some

discussion. The starting thermal and the plume induced by a heat source in an

enclosure are a result of physical instability of the flow field. In addition, in the introduction, we discussed the fact that this fluid model was one

in which viscosity has been ignored, and, therefore, it could be considered to

result from the Navier-Stokes equations in the limit of very large Grashof

number (roughly Reynolds number squared).

Over the last several years, starting with the pioneering work of Brown and

Roshko,^® there has been a reexamination of the meaning of turbulence in shear

flows. There had been a growing realization that turbulence is not satisfactorily described in terms of velocity correlations and their corresponding

spectra only. Rather there are distinct coherent vortex structures in shear

flovrs, and Brown and Roshko vividly demonstrated the existence of these coher28
ent structures in turbulent shear flows using shadowgraphs to visualize the

flow field. In particular, among many other Interesting features, Brown and

Roshko found that large scale coherent structures of the same type existed in a

shear layer Independent of the value of the Reynolds number provided only that

the Reynolds number is large enough to have turbulence. Subsequent studies in

other flows, see Roshko^^ for some of the references, have shown that organized

structures exist in these flows also. In addition, recent theoretical studies

have shown that many of the features of the large scale coherent structures

observed in high Reynolds number flows can be described by vortex structures

satisfying Euler's equations.

The spatially periodic structures calculated in the starting thermal and the

plume have been found to be vortices of alternating sign produced in the heating region and convected out by buoyancy. These vortices increase in strength

with height above the floor and occur in ant i-S3rmme trie pairs with respect to

the centerline of the room. We interpret these structures as analogous to the

large scale coherent structures observed in turbulent shear flows. In addition, because these vortices are convected with the buoyant flov;, the spatial

periodicity is translated into a temporal frequency: at any point within the

plume, each of the dependent variables oscillates with a frequency related to

the rate at which vortices are generated and convected away. Experiments, both

at and elsewhere, have demonstrated qualitatively the same feature;

namely, buoyant "puffs", or regular upwellings followed by short quiescent

periods, produced at a frequency determined by the experimental arrangement.

The frequency of these puffs is also found to agree with the frequency

predicted by these calculations.

29
B. Rectangular Enclosure with Off-Center Heat Source

In Figure 19 contours of constant temperature are shown at dimensionless time

2.0. As in the previous example, these contours mirror the contours of constant volumetric heat addition rate because the velocity field is very small

and problem is linear yet at this early time. The heat source is centered

along the floor, one quarter of the length of the room from the left wall.

At time 10.0, Figure 20, the Isotherms have the same form as in the previous

example. Figure A, at time 10.0. A buoyant bubble pinched at the bottom by

the inflow is beginning to rise from the heat source. There appears still to

be sjonmetry about the centerline of the heat source even though the source is

not symmetrically placed within the room. At time 13.0 a buoyant thermal has

developed, which is asymmetric with the heated fluid expanding more toward the

center of the room than toward the wall. Figures 21, 22, 23, and 24 show the

thermal rising, growing, hitting the ceiling and spreading. The heated gas

flows across the celling toward the right in a gravity current while the heated

gas at the left moves down the left wall in a fashion similar to that shown in

the first example calculation. As before, there are large scale vortex

structures in the plume, but because of the placement of the heat source within

the room, there is no longer any symmetry.

30
APPENDIX

Alternate Non-dimensional Variables

Consider the dimensional system of equations (1) and (4) written in the form:

3p 3p 1 + Ui + p

I j _

dt dX£

1 , 0

- - / Qdv

V

0

= 0

3ui 3ujL 3p

p :— + :— + — + [p - pAy)] g^i = 0

3t 3x. 3x.

3 1

(Al)

3Ui y_]^ 1 1 = _[Q-- / QdV]

3x, y p V

Here p and uj[ are the density and velocity components as defined in the text.

The quantities pQ(y) and p are respectively the initial density stratification

in the vertical (y) direction and the difference between the local pressure and

the hydrostatic pressure at the height in question. This pressure difference,

which affects the fluid motion, is the quantity p defined in Equation (13b) of

the text. The heat source 0 is prescribed in the form

'*^0f(t)

1 )2 2 _

Q = e c X ^ y (A2)

^x^y^z /it

Note the slight difference in notation from Equation (27b) of the text.

We now seek to introduce non-dimensional variables, denoted with an asterisk

(*), that are close to those defined in the text but which reflect the strength

of the heat source . To this end, we define the following quantities:

Xi = ui = Uui*(xb*, t*)

t (H/U)t*.

> „H

0

Pap *(t*)

0

31
P = Pq(o) p*(xi(.*, t*) (A3)

p = Pq(o) {1 + e p*(X|^*, t)>

pQ^y) = U + 3 pQ*(y*)}

Here, is the undisturbed ambient pressure, Pq(°) the ambient density at the

floor, and H the height of the enclosure. The reference velocity U and the

dimensionless density ratio 3 are as yet undefined. These two quantities are

now determined by requiring that in the Boussinesq limit, when the density

ratio 3^0, all non-georaetric parameters disappear from the problem. This

leads to the following equations for U and 3^

This yields a velocity scale which differs from that employed in the text by a

U^/g = B

(AA)

3U/H = Qq/H^J^z ^a*

/^, and a dimensionless density ratio given by

2/3

(A5)

Finally, the equations of motion (Al) become:

+ (1 + 3 P*) D* = 0

32
(A6)

Note that when 3 is 0(1); i.e., when there are significant density variations,

the non-dimensionalization is for all practical purposes the same as that used

in the text. The most significant feature of this derivation is Equation (A5)

which determines 3, and hence the conditions under which a non-Boussinesq

model is necessary.

33
REFERENCES

1. Rehra, 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. Batchelor, G. K. , The Conditions for Dynamical Similarity of Motions of a Frictionless Perfect-Gas Atmosphere, Quart. J. Roy. Meteor. Soc. pp.

224-235 (1953).

3. Ouintiere, J., Growth of Fire in Building Compartments, Fire Standards and

Safety, Ed. A, F. Robertson, ASTM STP 614, pp. 131-167 (American Society

for Testing and Materials, New York, NY, 1977).

4. Porsching, T. A., A Finite Difference Method for Thermally Expandable

Fluid Transients, Nuclear Sci. and Eng. 64, pp. 177-186 (Sept. 1977).

5. Roache, P, J., Computational Fluid Dynamics , Hermosa Publishers, P. 0. Box

8172, Albuquerque, NM 87108.

6. Lamb, H., Hydrodynamics, Sixth Edition, Dover Publications, New York, pp.

378-380 (1932); also,

Prandtl, L., Essentials of Fluid Dynamics, Blackie and Son Ltd., London,

pp. 374-377 (1952).

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

Chapter 2 (1973); also,

Whitham, G. B., Linear and Nonlinear Waves, John Wiley and Sons, New York,

pp. 421-423 (1974).

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

Waves in Enclosures, National Bureau of Standards Report in preparation.

9. Arakawa, A., Computational Design for Long-Term Numerical Integration of

the Equations of Fluid Motion: Two-Dimensional Incompressible Flow. Part

I, J. Comp. Phys. J_, pp. 119-143 (1966).

10. Phillips, N. A., "An Example of Non-Linear Computational Instability," The

Atmosphere and Sea in Motion , pp. 501-504, Rockefeller Inst. Press in

association with Oxford University Press (1959).

11. Arakawa, A. and Lamb, V. R., Computational Design of the Basic Dynamical

Processes of the UCLA General Circulation Model, Methods in Computational

Physics, Vol. 17 General Circulation Models of the Atmosphere, J. Chang

(ed.). Academic Press, New York, pp. 173-265 (1977).

12. Harlow, F. H. and Amsden, A. A., Fluid Dynamics, A LASL Monograph , Los

Alamos Scientific Laboratory Report LA4700, Los Alamos, New Mexico (June

1971) .

34
13. 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.

14. Schwarz trauber, P, and Sweet, R. , Efficient FORTRAN Subprograms for the

Solution of Elliptic Partial Differential Equations, NCAR-TN/lA-109 , July

(1975).

15. Concus, P, and Golub, G., Use of Fast Direct Methods for the Efficient

Numerical Solution of Nonseparable Elliptic Equation, SIAM J, Numer.

Anal., 10, 6 (Dec. 1973).

16. Corcus, 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).

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

Elliptic Partial Differential Equation by Preconditioned Conjugate

Gradients, NBS Journal of Research, No. 5, pp. 367-390 (SeptemberOctober 1980)

18. Brown, G. L. and Roshko, A., On Density Effects and Large Structure in

Turbulent Mixing Layers, Journal of Fluid Mechanics, Vol. 64, pp. 775-816

(1974).

19. Roshko, A., Structure of Turbulent Shear Flows: A New Look, AIAA Journal,

Vol. 14, No, 10, pp. 1349-1357 (October 1976).

20. McCaffrey, B., Purely Buoyant Diffusion Flames: Some Experimental Results, National Bureau of Standards Report NBSIR 79-1910 (October, 1979).

21. Zukowski, E., Kubota, T, and Cetegen, B., Entainment in Fire Plumes, Fire

Safety Journal^, pp. 107-121 (1980/81).

22. Rehm, R. G. , Baum, H. R., Corley, D. M, and Barnett, P. D., Finite Difference Calculations of Buoyant Convection in an Enclosure, Part II: Verflcation of the Basic Algorithm, National Bureau of Standards Report in preparation.

35
Figure Captions

Figure 1

Figure 2a

Figure 2b

Figure 3

Figure 4

Figure 5

Figure 6

Figure 7

Figure 8

Figure 9

Figure 10

Figure 11

Figure 12

Figure 13

Schematic diagram of an enclosure with a two-dimensional heat

source, plume and region of smoke and hot gases.

1 Rectangular enclosure in dimensionless variables 0 < x < — , AR

0 < y < 1. The mesh upon which the difference scheme is based

as shown schematically for (I = J = 4) as a grid of solid lines.

The mesh of dashed lines joins the center points of the basic

mesh cells and is the grid upon which the pressure computation

is performed.

A typical mesh cell, with center located at x = (1 - 1/2) 6x and

y = (j - 1/2) 6y, Illustrating where all dependent variables for

the finite difference scheme are defined.

Contours of constant temperature at dimensionless time t = 2.0 in

a square enclosure using a 31 x 31 mesh. At this early time con- vection is unimportant, and isotherms reflect contours of constant

volumetric heat addition.

Contours of constant temperature at dimensionless time t = 10.0 in

a square enclosure using a 31 x 31 mesh.

Contours of constant temperature at dimensionless time t = 11.5 in

a square enclosure using a 31 x 31 mesh.

Contours of constant temperature at dimensionless time t = 12.5 in

a square enclosure using a 31 x 31 mesh.

Contours of constant temperature at dimensionless time t = 13.5 in

a square enclosure using a 31 x 31 mesh.

Contours of constant temperature at dimensionless time t = 14.5 in

a square enclosure using a 31 x 31 mesh.

Contours of constant temperature at dimensionless time t = 15.5 in

a square enclosure using a 31 x 31 mesh.

Contours of constant temperature at dimensionless time t = 16.5 in

a square enclosure using a 31 x 31 mesh.

Contours of constant temperature at dimensionless time t = 17.5 in

a square enclosure using a 31 x 31 mesh.

Contours of constant temperature at dimensionless time t = 18.5 in

a square enclosure using a 31 x 31 mesh.

Contours of constant temperature at dimensionless time t = 2.0 in

a square enclosure using a 63 x 64 mesh.

36
Figure 14 Contours of constant temperature at dimensionless time t = 10.0 in

a square enclosure using a 63 x i64 mesh.

Figure 15 Contours of constant temperature at dimensionless time t = 11.225 in

a square enclosure using a 63 x i64 mesh.

Figure 16 Contours of constant temperature at dimensionless time t = 13.225 in

a square enclosure using a 63 x i64 mesh.

Figure 17 Contours of constant temperature at dimensionless time t = 15.225 in

a square enclosure using a 63 x 64 mesh.

Figure 18 Contours of constant temperature at dimensionless time t = 16.187 in

a square enclosure using a 63 x 64 mesh.

Figure 19 Contours of constant temperature at dimensionless time t = 2.0 in a rectangular enclosure using a 62 X 31 mesh.

Figure 20 Contours of constant temperature at dimensionless time t = 10.0 in a rectangular enclosure using a 62 X 31 mesh.

Figure 21 Contours of constant temperature at dimensionless time t = 13.0 in a rectangular enclosure using a 62 X 31 mesh.

Figure 22 Contours of constant temperature at dimensionless time t = 15.45 in

a rectangular enclosure using a 62 X 31 mesh.

Figure 23 Contours of constant temperature at dimensionless time t = 17.95 in

a rectangular enclosure using a 62 X 31 mesh.

Figure 24 Contours of constant temperature at dimensionless time t = 20.45 in

a rectangular enclosure using a 62 x 31 mesh.

37
Smoke

and

X
X

TEMPERfiTURE

C0NT0URS

RT

T

=

2.000
TEMPERATURE

C0NT0URS

ATT=

10.000 4
TEMPERATURE

C0NT0URS

AT

T

=

11.500
TEMPERRTURE

C0NT0URS

AT

T

=

12.500
TEMPERATURE

C0NT0URS

AT

T

=

13.500
TEMPERRTURE

C0NT0URS

ATT=

14.500 4 \f!
TEMPERATURE

C0NT0URS

RT

T

-

15.500

I

j

i
TEMPERRTURE

C0NT0URS

RT

T

=

18.500
TEMPERniURE

C0NT0URS

RT

T

=

2.000

j i 1 1 1

1 rri 1 1 H 1

1 T Ti n rn 1 1 1 1 1 1 1 1 1 1 1 1 1 1 m i 1 1 111 1 111 1 11111 11 1 1.

'»

1 1 1 1 1 1 1 1 1 1 1 1 1 1

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 M 1 1 1 1

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

1 11 1 1 r

1/1

L/ii

r'i.r'i

rvin

iMV

»j

kill

f».i
TEMPERRTURE

C0NT0URS

AT

T

=

10.000
TEMPERfiTURE

C0NT0URS

ftl

J

-

11.225
TEMPERATURE

C0NT0URS

AT

T

=

13.225
TEMPERRTURE

C0NT0URS

AT

T

=

15.225

i

I

'I

t

t.

,t
TEMPERRTURE

C0NT0URS

AT

T

=

16.187
TEMPERATURE

C0NT0URS

AT

T

=

2.000
TEMPERfiTURE

C0NT0URS

AT

T

=

10.000
TEMPERATURE

C0NT0URS

RT

T

=

13.000

_ OT

/ V

'» 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
TEMPERATURE

C0NT0URS

AT

T

=

15.450

.1 I I I I I I rrn i t i ii i itti t m t i i m.

Hrr)

-lJ cn

S' r—

<2

I I I I M I I I 1
TEMPERATURE

C0NT0URS

AT

T

=

17.950
TEMPERRTURE

C0NT0URS

HT

T

=

20.450
NBS-114A (REV. 2-8C) U.S. DEPT. OF COMM. BIBLIOGRAPHIC DATA

SHEET (See instructions)

1. PUBLICATION OR

REPORT NO. NBSIR 81-2385 (NBS)

2. Performing Organ. Report No. 3. PubI ication Date

December, 1981

I. TITLE AND SUBTITLE

Finite Difference Calculation of Buoyant Convection in an Enclosure

Part I: The Basic Algorithm

5. AUTHOR(S)

Howard R. Baum, Ronald G. Rehm, P. Darcy Barnett and Daniel M. Corley

6. PERFORMING ORGANIZATION (If joint or other than NBS. see instructions)

NATIONAL BUREAU OF STANDARDS

DEPARTMENT OF COMMERCE

WASHINGTON, D.C. 20234

7. Contract/Grant No.

8. Type of Report & Period Covered

Final

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

Center for Applied Mathematics

National Engineering Laboratory

National Bureau of Standards

10. SUPPLEMENTARY NOTES

I I 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)

A novel mathematical model of buoyant convection in an enclosure, developed earlier,

is solved by finite difference techniques in the two-dimensional case. This model

has been developed as a principal analytical tool for the prediction of the movement

of smoke and not gases in fires. Effects of large density variations caused by

substantial heating are retained while acoustic (high-frequency) waves, which are unimportant to buoyant convection, are analytically filtered out. No viscous or

thermal conduction effects are included in the model. These two characteristics

(filtering and no dissipative effects) distinguish the model from all others describing buoyant convection. The mathematical model consists of a mixed hyperbolic

and elliptic set of non-linear partial differential equations: the problem is a mixed initial, boundary value one. An explicit time-marching algorithm, second-order

accurate in both space and time, is used to solve the equations. The computational

procedure uses a software package for solving a nonseparable elliptic equation

developed especially for this problem. The finite difference solutions have been

carefully compared with analytical solutions obtained in special cases to

determine the stability and accuracy of the numerical solutions. The computer model

has been used to compute the buoyant convection produced in an enclosure by a spatially

distributed heat source simulating a fire. The computed results show qualitative

agreement with experimentally observed buoyant convection in enclosure fires.

12. KEY WORDS (Six to twelve entries; alphabetical order; capitalize only proper names; and separate key words by semicolons)

buoyant convection; computations-finite difference; Euler equations; finite difference

equations; fire-enclosure; fluid flow; heat source-volumetric; partial differential

equations

13. AVAILABILITY

IXX Unlimited

I I For Official Distribution. Do Not Release to NTIS

r I Order From Superintendent of Documents, U.S. Government Printing Office, Washington, D.C. 20402.

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

14. NO. OF

PRINTED PAGES

67

15. Price

$ 8.00

USCOMM-DC 6043-P80
I

1J