AlllDS Reference

NBS

PUBLICATIONS NAri INST OF STANDARDS & TECH R.I.C. NBSIR 86-332S All 102490992

Jones, Walter W/Buoyancy driven fiow as QC100 .U56 NO.86-3329 1986 V19 C.1 NBS-P Buoyancy Driven Fiow As the

Forcing Function of Smoke

Transport Models

Walter W. Jones

Xavier Bodart

U S. DEPARTMENT OF COMMERCE

National Bureau of Standards

National Engineering Laboratory

Center for Fire Research

Gaithersburg, MD 20899

May 1986

U S. DEPARTMENT OF COMMERCE

*‘'*CAU O*

QC

100

.056

»NAL BUREAU OF STANDARDS

86-3329

1986

NBSIR 86-3329

' • ir BUOYANCY DRIVEN FLOW AS THE

FORCING FUNCTION OF SMOKE

TRANSPORT MODELS

RBS

RESEARCH

INFORMATION CENTER ^ -

Walter W. Jones

U S. DEPARTMENT OF COMMERCE

National Bureau of Standards

National Engineering Laboratory

Center for Fire Research

Gaithersburg, MD 20899

USA

and

Xavier Bodart

Centre Scientifique et Technique du Batiment

Champs Sur Marne

France

May 1986

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

NATIONAL BUREAU OF STANDARDS, Ernest Ambler, Director
AT^?n

T

otet at RittM

^3«?JA WOJ1 W?V>na YOMAYOU8

ijioMaSs woifr>titu3 omioro^

aj^OOM THO<?ai^Af]T

.I

m:

Affuji tt hJSbyV

t , i

-f *'^

3r>«-3M^<r»l li!) rt’^*6n»{t^ 'y »l&r>fi'4Si’i< ywfffcrjdKJ itiinjt»iitf( rUri»rti^''tin »!/»tiO'‘;{ f:*'>9V -HM ,||l« 4'‘W’»0'rl'' a

'w

pj

: f. •,

*»li8 *on.«ijrt$^

' «cia9i^

'Si' >'?.'V/:|’.'

,30fiaMM<)a %D it u

^n«»«ria \c9imA :»Mnil ,KCJIAaM^rt ^

IA^O^AM m
TABLE OF CONTENTS

Page

LIST OF FIGURES iv

LIST OF TABLES v

Abstract 1

1. INTRODUCTION 1

2. PROBLEM TO BE SOLVED 2

3. SOURCE TERM 5

4. A SAMPLE PROBLEM 10

5. CONCLUSIONS 14

6. REFERENCES 14

-ill-
LIST OF FIGUEIES

Page

Figure 1. Typical flow patterns through vents which are subjected

to fires 13

Figure 2. Pictoral of the integration scheme used in eqn. (7) 16

Figure 3. Schematic of the notation used in the source term

calculation 17

Figure 4. A subset of the flows of class I; see Table I 18

Figure 5. An example of the effect which change in the interface

height will cause for one of the diagrams in Fig. 4 19

Figure 6. Output from the pressure function evaluation at 0, 20,

40 and 60 seconds. The surface plot shows the actual

values and the contours are: yellow: df/dt = 0; green

and red: df/dt = ± 6. Pressure variation is left to

right on the diagram and temperature variation is into

the picture (upper left figure) 20

Figure 7 . Data capture “on-the-f ly" from the function evaluation

from 49 to 51 seconds (upper right figure) 20

Figure 8. A composite of the output of the function evaluation

from 0 to 60 seconds (bottom of page) 20

-Iv-
LIST OF TABLES

Page

Table 1. Classes of flow diagrams 9

Table 2. Condition used to separate and order the flow diagrams ... 11

-V

BUOYANCY DRIVEN FLOW AS THE FORCING FUNCTION OF SMOKE TRANSPORT MODELS

Walter W. Jones

Xavier Bodart

Abstract

Flow at vents is the major driving force in smoke transport models. The

precision with which we can calculate these flows determines to a great extent

how accurately we can model buoyant flow and the inherent speed of the models.

This report describes some of the problems encountered in calculating these

flows, and gives a general algorithm for their calculation.

Key words: vent flow, smoke transport model, fire modeling, differential

equation.

1 . INTRODUCTION

Over the past few years modeling fire growth and smoke transport has

become an important aspect of fire research. As computers have become cheaper

and faster, the ability to handle the equations associated with such phenomena

has improved to the point where numerical experiments can replace some

physical experiments. This is useful in avoiding the expense and effort of

actual fire tests. To this end there are two aspects of modeling which are

important. The first is that the physical algorithm must be correct and the

second is that the numerical routines which are utilized must be fast and able

to handle the wide range of values which occur in natural phenomena. Typical

time scales range from microseconds for chemical kinetics to tens of seconds

for heat conduction.

- 1 -
The predictive equations which are solved arise from the conservation of

mass, momentum and energy. When integrated over a finite sized control

volume, they are first order, nonlinear, ordinary differential equations, of

the form

^

= f(x), (1)

where x represents a vector of unknown dependent variables. Much of the

modeling has assumed a simplification which is to solve the pseudo steadystate as the transient term vanishes

f -£(x) >0.

In either case a solution is found by varying x in the phase space about some

initial value x and looking for a minimum in the value of the function f(x)

o or actually - f(x)}. The thesis of this paper is the difference in the

physical meaning of these two types of equations and to demonstrate that

solving the original ODE's is superior to solving their algebraic counterpart.

The underlying point is to demonstrate the rationale that the apparently more

complex form of the equations is actually easier to solve, and in addition to

show that the algebraic form may not always yield the correct solution.

2. PROBLEM TO BE SOLVED

The difficulty which arises in the solution of eqn. (1) is the bumpiness

in the multidimensional phase space allowed by these equations (one for each

zone or control volume). The term "phase space" is used in this context to

-2 -
mean the manifold of pressure, temperature, etc., in which the solution of the

relevant set of equations lie, and the possible physical variations which can

occur. Even for a single equation, e.g. the pressure equation, the topology

is not simple. There is an assumed simplification of considering only relaxed

dx

states — = 0 which can actually make achieving the final state more diffidt

cult. No direction is provided in finding a minimum and it gives no sense of

whether a specific bump in the manifold is an absolute minimum. By retaining

the transient term some relief is obtained for this situation. The most

important improvement in finding a solution is that search direction in phase

space is given by the transient term. Thus the basic search algorithm is

considerably simplified. The second bit of help comes from the observation

that the derivative terms do not vanish at a pseudo minimum but only at

absolute minimum. Further, the point for which the transient term vanishes

may not be the correct solution. The correct solution actually minimizes the

difference between the left and right side of eqn. (1).

Several examples should serve to illustrate the problem, a possible

solution and the means by which this translates into a technique for a general

solution of the equations.

The actual equations which we solve are somewhat more complex than

eqn. (1); however, it is useful to use a simplified example to test the

techniques. Then one can test the heuristic without getting bogged down in

technical details of numerical programming. The set of predictive equations

which we use for each compartment in a typical zone model [1], are

dt (3-l)V

s

( 2 )

-3 -
dT

dt

1®

/ +—S—s

“ ( 6- 1 )

dt 3 \PV^j(^ A

dV

u

dt (h)

V

c m T + E

- —^ s p u u u ^

where

s = c m T + c m.T. + E + E. puu p Z Z u Z and

3 = c /R = y/(y-1)

P

(3)

(A)

(5)

where T is the temperature and E^ and E^^ are the energy release rates into the

upper and lower zones, respectively. The primary source term is the s which

depends mainly on fluid transport, that is enthalpy flux driven by density and

pressure gradients. It also varies the most rapidly. Actually, the microscopic chemical kinetics have a shorter time constant but the spatially

averaged zone properties (T, p) do not vary on such a time scale and thus it

is reasonable to ignore these very short times. In addition to improving the

solution technique itself, we have developed a better method of finding the

dominant portion of the source term, m, the enthalpy flux from mass flow. Due

to its importance, we will discuss it in the next section before giving

examples of the solution of this type of problem.

The conservation equations which we solve are for energy and mass.

Momentum is not defined within a control volume but only at the boundaries

where zones are connected. So we do not solve the momentum equation directly

-4 -
but rather use the integral form known as Bernoulli’s equation which yields

mass flux as a function of the pressure and density differences.

3 . SOURCE TERM

Equation (2) is the one on which we shall focus as it is this one that is

most often "simplified" in the as3nnptotic sense. In order to explain some of

the difficulty encountered in solving this equation, we digress for a moment

to discuss the primary source term, namely fluid transport through vents.

This fluid flow phenomenon connects the control volumes and is dominant

because this term fluctuates most rapidly of all the source terms in response

to changes in the environment. One of the inq)rovements which we have incorporated into our current models is a means of calculating these flow fields with

the correct number of neutral planes (up to three) and without discontinuties

in the function. This latter feature implies no discontinuities in the first

order derivatives for the ODE's.

Typical types of flows which can occur in fires are illustrated in

fig. 1. The notation for the flow is

SS = upper layer to upper layer

SA = upper layer to lower layer

AS = lower layer to upper layer

AA = lower layer to lower layer

- 5 -
This notation was originated earlier [2] for the single neutral plane case but

is useful for a physical description of the general problem. One generally

uses the Bernoulli equation to calculate the flow velocities between two

compartments which are connected by an opening. Indeed this is the solution

for the momentum equation which allows us to exclude it specifically when

solving the conservation equations in general.

The general form is

•m. = lO

C • S • /2p (kg/sec).

where m = mass flow rate

C = orifice coefficient

S opening area (m^)

P = gas density on side "i"

P .1 = pressure on side "i"

P

o = pressure on side "o".

The implication of using this equation is that the pressure at a stagnation

point is used. That is, the flow velocity vanishes where the pressure Is

measured. As the pressure always appears as a difference, in principle it

does not matter whether absolute pressures or pressure defects are used. More

accuracy is obtained by using pressure defect, however. This avoids, to some

extent, the problems of the small difference of large numbers.

The pressure is always calculated with respect to the base of a compartment. With this in mind we can express the pressure on the other side of a

-6 -
partition as a function of the variable (y), the height above the base:

implicit in our use of equation (6) is that the opening is rectangular, so

that the area integral of the flow term will allow us to remove the width from

the integral.

That is

Z

2 flow = / /, . , . pV dzdb width / pVdz. (7) •'width ’height Z'

Thus the total flow becomes

m_ = C W Z S, /2 • p. , |p.(z) - P (z)| . i->-o ^

k l»k I 1 o '

The pressure term will be reversed if the flow is o -> i. Thus we have the

integral over the area shown in fig. 2. Terms are as before except that

p. , is the average inlet mass density within area "k". The simplest way to

1 , K

define the limits of integration is with neutral planes, that is the height at

which flow reversal occurs, P. (z) = P (z).

1 o

Each side of an opening is assumed to consist of two homogeneous gas

layers (zones) of uniform density and temperature. There is an apparent

inconsistency in that the equation of state dictates P = pRT, and we assume P

varies but p and T remain constant, at least within a zone. This pressure

fluctuation is so small compared to the magnitude of the base pressure that

ignoring it for all calculations except the flow field is reasonable. For

simplicity of notation we will use a slightly different means of Identifying

the zones. The correspondence is

-7 -
upper "i" zone = 1

lower "i" zone = 2

upper "o" zone = 3

lower "o” zone = 4

The former terra using "SA" is physically more understandable and the results

can be put into these terms, but the derivations are more compact with the

numerical indicies. Figure 3 shows a schematic of the notation. The corresponding densities are , P^ p^. and are the height of

a sill, soffit, hot/cold interface in the "i" compartment and hot/cold interface in the "o" compartment. With the base (reference) pressures P^(0) and

Pq( 0), the internal pressure on each side is given by

P^ (z) = P^(0) - min (y, Z^) • p^g - max (z-Z^, 0) • p^g

P (z) = P (0) - min (y, Z )

• p, g

- max (z-Z , 0) • p g

o o o 4 o J

The function F (z) = P. (z) - P (z) can be considered a family of functions of

i e one variable, z. In all cases the pressure appears only in this form. In

principle, this family of curves can lead to an inordinate number of possible

flow fields. By Imposing the restrictions found in fire scenarios, however,

we end up with only a few possiblities . For stratification to occur, the

following restrictions can be imposed;
We can also require

without loss of

reverse the "i"

These cases and

generality, since for the reverse, it is only necessary to

and "o" compartments. We are left with five different cases,

their restrictions are shown in Table I.

TABLE I

Class Restrictions Max // of neutral planes Figure 4

1 ^2 Q. z

.

1 < z — 0

1 a

II ^2 < a. z

.

1

< z — 0

2 b

III ^3 <

^2 < P, , Z . > z

4 1 0

3 c IV ^2 < 0K Cl z,

1 < z

0

2 d

V ^2 < P3. z

.

1 < Z

0

1 e

Classification of types of flow which occur in a vent based on the relative

densities and interface heights. If there exist soffits or sills, then the

number of neutral planes within a vent can be less than the number indicated.

If there were no soffits or sills to consider, then the calculation would be

fairly straightforward. However, the possibility of soffit/sill combinations

requires many numerical tests in the calculation. Class I is the basis of the

analysis of classes 11-V. It can have 44 different flow combinations, depending on the relative position of H^, B^, and Z^. It may contain at most a

single neutral plane (flow reversal). Twenty four of these combinations are

without a neutral plane and twenty with a neutral plane. Figure 5 shows the

effect of the "o" compartment interface height on the flow field from one of

the diagrams shown in Figure 4 .

- 9 -
Table II shows the criteria hsed for solution in classes II-V. The

interval H^] can subsequently be partitioned to contain at most a single

neutral plane so that the logic used for class I can be utilized. The only

caveat is to be sure that the equalities and inequalities are treated

correctly.

4. A SAMPLE PROBLEM

We limit our discussion to eqn. (2) with the appropriate source terms.

It is most often this equation to which the asymptotic approximation is

applied, namely

dt - 0 . ( 8 )

This approximation is appropriate for a steady state situation, for example,

where a fire is fully developed and the flow fields have been established.

For a transient problem such as a developing fire or the case when a window Is

broken, such is not appropriate. Most models of fire imply this approximation

by using large initial time steps, perhaps of one or two seconds. However,

when faced with a true transient phenomenon, all sub-divide the time step,

generally to a value below the pressure relaxation time.

Furthermore, a fundamental problem arises in that the solution for

eqn. (2) may not be identical to that of eqn. (8). The requirement Is that

the difference between the left- and right-hand sides (LHS and RHS) of the

conservation equations should be minimal. This is truly equivalent to the

statement of eqn. (8) only in a steady state regime. So this approximation

- 10-
TABLE II

F(z) = P^(z) - Pg(z)

Condition Max. no. of neutral planes

in [Bj, H^]

II. j< B^ 0 , 1

Zi > Hf

Bf < Z. < F(Zi) • max (F(Bf), F(Hf))20

F(Z^) • max (F(Bj), F(H^) < 0

F(Zj^) • min (F(Bf), F(H^)) 2 0

F(Z^) • min (F(Bf), F(Hf) ^0

III. 2 Hf

Zo <. Bf, Z^ 2 Hf

2 Bf, Zi < Bf

0 , 1

0

1

2

0 , 1

0

0

Zo < Bf* Bf < Zf < Hf

F(Zf) • max (F(Bf), F(Hf)) 2.0 0

F(Z^) . max (F(Bf), F(Hf)) > 0

F(Zf) • min (F(Bf), F(Hf)) 20 1

F(Z^) • min (F(Bf), F(Hf)) 2 2

B

f < Zo < Hf

Zi > Hf

F(Zq) • min (F(Bf), F(H.))2 0 0

F(Z^) • min (F(Bf), F(Hf)) < 0

F(Zq) max (F(Bf). F(H.)) 1

F(Z )

• max (F(B.), FtHf)) 2 0 2

Z^ < H^ [Bj < Z < Z, < Hj^]

min (F,(Bf1, F(Z.)) • max (F(Hf), F(Z )) >00

min Fj(B^), F(Z^7) * max (F(Hf), F(zJ) >_ 0 1,2 or 3

IV. Z . 2 Hf 0,1

Zf 2 Bf 0,1

Bf < Z. < Hf

FCZf) • max (F(Bf), F(Hf))2 0 0

F(Zf) • max (F(Bf), F(Hf)) < 0

F(Zf) * min (F(Bf), F(Hf)) > 0 1 F(Z.) • min (F(Bf), F(Hf)) 2 0 2

V. F is monotonic over [B^, Hf] 0,1

- 11 -
should only be made if there is a real gain, such as reduced computing

I requirements or equations which are simpler to manipulate. An example should

Illustrate some of the difficulties.

The context of a problem which we will examine is a four compartment

calculation using the FAST [1] model to find a solution. The point of

Interest is to examine the solution of eqn. (2) in the second compartment.

There are three ways to examine the equation; first, a pseudo analytic technique is to start with the converged solution of temperature and pressure in

all compartments and then form a perturbation expansion of the RHS of

eqn. (2), in terms of pressure and temperature for compartment //2 only. The

result is shown in fig. (6) as a surface plot of the value of the PHS of eqn.

(2) for a variation in pressure (left to right) and temperature (into the

picture). The total calculation lasts 60 seconds and fig. (6a, 6b, 6c, and

6d) are done at 0, 20, 40, and 60 seconds, respectively. The variation in

pressure was temperature T^ 100 K. If this equation were

solved by Itself, then the solution would come fairly easily and eqn. (8)

might be appropriate; second, we can capture the results from the ODE solver

which is used in FAST. This was done at 50 ± 1 seconds and the results are

shown in fig. (7). The axes are the same as in fig. (6). The latter figure

shows the somewhat more complex interaction of all sixteen equations, since we

capture the intermediate results with variations in all parameters. From the

figure it can be seen that it may not be possible to get from a pseudo minimum

on the left hand side to the correct answer on the right hand side using

eqn. (8). Equation (2) helps by providing a direction and distance for subsequent solutions, namely the dP/dt term, the term which has been dropped In thiasymptotic approximation. Figure (8) shows a composite of the result of

- 12 -
applying the second technique to the entire time history of the pressure

equation.

Finally we can look at a simple piece of this problem analytically. We

can extend earlier w/ork on pressure relaxation. Attention is directed to fig.

(1) and eqn. (28) in the paper by Rehm and Baum [3]. Two similar examples can

be given, the first being the forced flow in the initial stages of a fire.

For this case, expansion in the room of fire origin forces air into an

adjacent compartment. For this case, eqn. (2) becomes (P' ~ ~ P

3

)

P'

1_ ^

dt

1 (3 - 1)V C C, 6 T p P d

3/2

with P^ = P^(t=o). If the asymptotic approximation were correct then the

right hand side would be constant for an e-foldlng time for P'. Integrating

this equation, we find a characteristic time for the pressure equation is

V

where V is the volume of the compartment of interest, A is the vent opening

and 6T is the temperature difference between the two compartments. For

typical initial conditions we obtain a time constant about 1 second, but it

can vary over a wide range. This is certainly not less than the typical time

step used for the algebraic solver in fire models. The other extreme would be

a near steady state regime of equal inflow and outflow. In this case, the

equation reduces to

dt

= ,0^ A |.l/2 (6T)^'^^

(RT)^

,3/2

- 13-
where R is the gas constant. For typical values of the parameters, this

equation has a time constant of about 0.10 seconds, once again being within

the range of times used in solving the other conservation equations.

5. CONCLUSIONS

We have tried to indicate where the problems arise in solving the

conservation equation used in predicting fire growth and smoke spread. The

type of predictive equations used and the primary driving term have been

discussed in detail. Finally, a sample calculation has been presented to show

where the problems actually exist for equation solvers and the reason that an

ODE solver, in general, has an easier job in threading its way through the

phase space thicket than does an algebraic solver. Experience has shown [4]

that this method yields a reduction in computing time of a factor of two to an

order of magnitude over the time required for a model wherein the asymptotic

approximation has been made.

6 . REFERENCES

[1] Jones, W.W., A Multicompartment Model for the Spread of Fire, Smoke and

Toxic Gases, Fire Safety Journal^, 55 (1985).

[2] Tanaka, T., A Model on Fire Spread in Small Scale Buildings, BRI (Japan)

Research Paper 79 (1978).

[3] Rehm, R.G. and Baum, H.R., The Equations of Motion for Thermally Driven

Buoyant Flows, Journal of Research 8^ 297 (1978).

[4] Jones, W.W., An Improved Numerical Method for Fire Models, Fall Technical

Meeting of the Combustion Institute, Atlantic City (1982).

- 14 -
15

Figure
c

CN

16

Figure
'////A

I

CM

n

17

Figure
0

Q.

<

Q.

<

Q.

<

a.<

D

O

n

0

<r

18

Figure
t f

1

1

CM 1

1

1

1

1

1

1

1 o

1

1

•

N

kVOs

1

1

1

1 •N19

Figure

20

NBS-H4A (REV. 2-80

U.S. DEPT. OF COMM. BIBLIOGRAPHIC DATA

SHEET (See instructions)

1. PUBLICATION OR

REPORT NO.

NBSIR-86/3329

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

4. TITLE AND SUBTITLE

Buoyancy Driven Flow as the Forcing Function of Smoke Transport Models

5. AUTHOR(S) y^iter W. Jones, Xavier Bodart

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

national bureau of standards

DEPARTMENT OF COMMERCE

Gaithersburg, Maryland 20899

7. Contract/Grant No.

8 . Type of Report & Period Covered

IR Final

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

10.

SUPPLEMENTARY NOTES

(321 Document describes a computer program; SF-185, FlPS 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)

Flow at vents is the major driving force in smoke transport models. The precision

with which we can calculate these flows determines to a great extent how accurately

we can model buoyant flow and the inherent speed of the models. ihis report

describes some of the problems encountered in calculating these flows, and gives

a general algorithm for their calculation.

12.

KEY WORDS (Six to twelve entries; alphabetical order; capitalize only proper names; and separate key words by semicolon s) fire models; numerical methods; smoke movement; zone models

13. AVAILABILITY

^31 Unl imited

|_J For Official Distribution. Do Not Release to NTIS

(33 Order From Superintendent of Documents, U.S. Government Printing Office, Washington, D.C.

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

14. NO. OF

PRINTED PAGES

26

15. Price

$9.95

USCOMM-DC 0043-P80
AIMoi 41^ 4R4> ATM atTMlIMliSM

»-Wi'n--ntwt

x0JTtTHKJii f3MA >* Jt»T#

*V V - •; iT^ gotajon *(t» » woJt oaygnt toBoxauS

'I

1 .1

. , F" 'i'^wp

%1 - <lWWg5^/JW\'Ti

'TAlTfflX xfstoot .W imitlM

%.R ’^‘^^

lp» I ^„;

.

, '1/ iiyi, 1,1 _ ^»i»oyts }M^< \5) MOttAllw*#|*OOw»HK‘' i«» 4

' MIB! ^

WAf»«« JAHilTAi-

^ ll!li®I lOAJMHIOa ^^TwailTAA-Hf

b«r#lrf»»H »«iud«*t»43lPO

MTOW ^Afiimrnj^kjrit

’• *-* ,<r*»y>»Mt i>^<w^oa 4>H-'^» ciTMn(|ei< iii;<?»_-t<«h<»j < tXto—t wp»«wvj^_r t

!i- a ’^^. ”‘ ^ ’* Mun^wH1 1» ‘ -^- ' a

”‘ ^ ’* '«*’•*'"* Mwniy^l

IT^IWBU -^aim at .3-iaX gntvni to(<« wli li c^>»y i2 «,«

^

****^‘“^!^®* ^

l»7..xu,B,a.

.iK454Hi..U5

waido,,

70’i«ti 5o

^

U»or,9|j-'l

r»l Hl' l l^l H I ?l^l V*# I'** pa^e '»A)r«««ik44tr «f taiMASB' ' :iT

L^ .. _®

^ {SPSillsiVO* AldOM ;»&CM|hiJ0: lc.3ii(^HiHnt 79:©iji*-« *il : A • *«•'' >*H <*0 .•«M«(U»**0 [ i

*^i^

:

^^^^A*** t/\.90iy^ rMpa<nr%«yAO 4>0 ,»xi'—<—<^{| >* in' ~ ^aintnani fl '-^ "T

Hr*;'

'.' - Av .'Wat^jtrXt^ .letTKV aMnOf^

I ^

<aJi.. .. 1. #