NUREG/CR-5233

A Computer Code for

Fire Protection and Risk

Analysis of Nuclear Plants

Prepared by A.K. Singhal, S.D. Habchi, A.J. Przekwas

CFD Research Corporation

Prepared for

U.S. Nuclear Regulatory

Commission
NOTICE

This report was prepared as an account of work sponsored by an agency of the United States

Government. Neither the United States Government nor any agency thereof, or any of their

employees, makes any warranty, expressed or implied, or assumes any legal liability of responsibility for any third party's use, or the results of such use, of any information, apparatus.

product or process disclosed in this report, or represents that its use by such third party would

not infringe privately owned rights.

NOTICE

Availability of Reference Materials Cited in NRC Publications

Most documents cited in NRC publications will be available from one of the following sources:

1. The NRC Public Document Room, 1717 H Street, N.W.

Washington, DC 20555

2. The Superintendent of Documents, U.S. Government Printing Office, Post Office Box 37082,

Washington, DC 20013-7082

3. The National Technical Information Service, Springfield, VA 22161

Although the listing that follows represents the majority of documents cited in NRC publications,

it is not intended to be exhaustive.

Referenced documents available for inspection and copying for a fee from the NRC Public Document Room include NRC correr:oridence and internal NRC memoranda; NRC 3ffice of Inspection

and Enforcement bulletins, circulars, information notices, inspection and investigation notices;

Licensee Event Reports; vendor reports and correspondence; Commission papers; and applicant and

licensee documents and correspondence.

The following documents in the NUREG series are available for purchase from the GPO Sales

Program: formal NRC staff and contractor reports, NRC-sponsored conference proceedings, and

NRC booklets and brochures. Also available are Regulatory Guides, NRC regulations in the Code of

Federal Regulations, and Nuclear Regulatory Commission Issuances.

Documents available from the National Technical Information Service include NUREG series

reports and technical reports prepared by other federal agencies and reports prepared by the Atomic

Energy Commission, forerunner agency to the Nuclear Regulatory Commission.

Documents available from public and special technical libraries include all open literature items,

such as books, journal and periodical articles, and transactions. Federal Register notices, federal and

state legislation, and congressional reports can usually be obtained from these libraries.

Documents such as theses, dissertations, foreign reports and translations, and non-NRC conference

proceedings are available for purchase from the organization sponsoring the publication cited.

Single copies of NRC draft reports are available free, to the extent of supply, upon written

request to the Division of Information Support Services, Distribution Section, U.S. Nuclear

Regulatory Commission, Washington, DC 20555.

Copies of industry codes and standards used in a substantive manner in the NRC regulatory process

are maintained at the NRC Library, 7920 Norfolk Avenue, Bethesda, Maryland, and are available

there for reference use by the public. Codes and standards are usually copyrighted and may be

purchased from the originating organization or, if they are American National Standards, from the

American National Standards Institute, 1430 Broadway, New York, NY 10018.
NUREG/CR-5233

RP

A Computer Code for

Fire Protection and Risk

Analysis of Nuclear Plants

Manuscript Completed: June 1988

Date Published: September 1988

Prepared by

A.K. Singhal, S.D. Habchi, A.J. Przekwas

CFD Research Corporation

3313 Bob Wallace Avenue, Suite 205

Huntsville, AL 35805

Prepared for

Division of Engineering

Office of Nuclear Regulatory Research

U.S. Nuclear Regulatory Commission

Washington, DC 20555

NRC FIN D2041

ABSTRACT

A fire modeling computer code has been developed with sufficient

flexibility for accurate representations of geometry, ventilation

and other conditions as may be present in cable rooms, control

rooms, and other enclosures in nuclear power plants. The

computer code is capable of three-dimensional, transient,

turbulent flow and heat transfer calculations with chemical

reaction and radiation. The code has a modular struc-ture,

specifically designed, for fire problems.

The code employs the latest relevant finite-volume solution

techniques. The code has been applied to a series of benchmark

problems and a recent fire test problem. This has confirmed the

feasibility of the fire code. Considerable further work is

needed to enhance the physical models (to improve the realism of

predicted solutions) and to validate and document the final code.

Specific recommendations are made for Phases II and III of the

Project.

iii

TABLE OF CONTENTS

pacre

ABSTRACT

1. INTRODUCTION AND SUMMARY

1.1 Background

1.2 Objective of Phase I Study

1.3 Results of Phase I Study

1.4 Recommendations for Phase II Work

1.5 Recommendations for Phase III Work

1.6 Outline of the Remainder of the Report

2. MATHEMATICAL BACKGROUND

2.1 Introduction

2.2 Governing Equations

2.3 Description of Physical Models

Turbulence Model

Heat Transfer Model

Combustion Model

Thermal Radiation Model

2.4 Description of Numerical Procedure

Finite-Volume Differencing Scheme

Solution Procedure

3. CODE STRUCTURE

1

1

1

2

4

5

5

4.

5.

CODE VERIFICATION STUDY

4.1 Overview

4.2 Two-Dimensional Driven Cavity Flow

4.3 Two-Dimensional Laminar Flow Over a Back

Facing Step

4.4 Developing Turbulent Pipe Flow

4.5 Two-Dimensional Axisymmetric Flow Over a Back

Facing Step

4.6 Two-Dimensional Natural Convection

4.7 Three-Dimensional Natural Convection

4.8 Shear Layer Mixing (Propane-Air Non-Reacting)

CODE DEMONSTRATION FOR FIRE PROBLEMS

5.1 Problem Selection

5.2 Test Case Selected (FMRC/SNL/NRC)

Experimental Setup

Grid and Boundary Conditions

Fire Source and Combustion Setup

5.3 Results of 3-D Computations

Steady-State Results

Transient Results

5.4 Investigatory 2-D Computations

Grid Size Sensitivity Study

Fire in Room with Obstacles

Effects of a Wider Exhaust Outlet

wardward5

5

6

8

8

10

10

12

12

13

14

15

17

17

19

19

19

25

25

25

31

36

36

38

38

38

43

44

44

44

53

54

54

58

V
TABLE OF CONTENTS (CONT.)

Page

Effects of Simulating Ventilation Air with Finite

Momentum 62

Effects of Increased Ventilation Flow Rate 62

6. CONCLUSIONS AND RECOMMENDATIONS 62

6.1 Achievements of Phase I Study 62

6.2 Recommendations for Phase II Work 62

6.3 Recommendations for Phase III Work 66

7. REFERENCES 67

BIBLIOGRAPHICAL DATA SHEET (NRC Form 335) 70

vi
LIST OF ILLUSTRATIONS

Figure ag

1.1 Main Modules of the Fire Code 2

2.1 Partitioning of the Turbulence Energy Spectral 9

Density in the Multiple Scale Turbulence Model

2.2 Staggered Grid Notation in the x-y Plane 13

3.1 Structure of the Pre-Processor 16

3.2 Structure of the Main Code 18

4.1 Geometry and Computational Grid for Driven 20

Cavity Problem

4.2 Streamline Contours from, a) Present Fire Code; 21

b) Ghia, et. al. [23]

4.3 Magnified Inlet Region of Backward-Facing Step; 22

a) Geometry, and b) Computational Grid

4.4 Comparison of Predicted Axial Velocity Profiles 23

with Experimental Data of Armaly, et. al. [24]

4.5 Geometry of the Developing Pipe Flow Problem 24

4.6 Predicted and Measured Axial Velocity Profile 24

for Fully Developed Turbulent Pipe Flow

4.7 Geometry of the Inlet Section of the 2-D Axisym- 26

metric Backward-Facing Step Problem

4.8 Predicted and Measured Axial Velocity Profiles 26

4.9 Grid and Boundary Conditions for 2-D Natural 27

Convection Problem

4.10 2-D Natural Convection Velocity Vectors for Re=106 27

4.11 Streamlines in a Square Enclosure for a) Ra=10 3 ; 28

b) Ra=10; c) Ra=107

4.12 Temperature Contours in a Square Enclosure for 29

a) Ra=103; b) Ra=106; c) Ra=lO7

4.13 Grid and Geometry for 3-D Natural Convection 30

Problem

4.14 Predicted Vertical (xy) Velocity Distribution for 32

the 3-D Natural Convection Problem

vii
LIST OF ILLUSTRATIONS (CONT.)

Figure Page

4.15 Predicted Horizontal Velocity Distribution for 33

3-D Natural Convection Problem

4.16 Predicted Vertical (yz) Velocity Distribution for 33

3-D Natural Convection Problem

4.17 Calculated Isotherms in Three Mid Planes 34

4.18 Geometry and Computational Grid for Shear Layer 35

Mixing Problem

4.19 Predicted and Measured Mixture Fraction at the 37

Centerline

5.1 Enclosure Ventilation System of FMRC/SNL/NRC 39

(1986/87) Tests

5.2 Location of Fire Source and Ventilation Inlet and 40

Outlet Ports

5.3 Measured and Simulated Heat Release Rate Profile 41

5.4 Computational Grid Chosen for the Demonstration 42

Calculation

5.5 Calculated Velocity Distribution in the Test Room 45

for Initial (Steady-State Pre-Fire) Condition

5.6 Temperature Contours at Near-Fire Vertical y-z 46

Plane at t = 1, 2, 3, and 5 Minutes

5.7 Temperature Contours at Near-Fire Vertical x-y 47

Plane at t = 1, 2, 3, 4, 5, and 6 Minutes

5.8 Velocity Contours at the First Vertical y-z Plane 49

Near Ventilation Exit

5.9 Velocity Vectors at Mid-Vertical Plane Near 50

Ventilation Inlets (yz - 6)

5.10 Velocity Vectors at the Last Vertical y-z Plane 51

Near the Wall

5.11 Comparison of Predicted Vertical Temperature 52

Profiles with Experimental Data

5.12 Computational Grids Used for Grid Sensitivity 55

Study; (2-D Investigatory Problem)

viii
LIST OF ILLUSTRATIONS (CONT.)

Figure Page

5.13 Velocity Vectors for Fine Grid 56

5.14 Temperature Contours for Fine Grid 56

5.15 Vertical Temperature Profiles Near the Fire Source 57

5.16 Fine Grid with Obstacles for Mock-Up 2-D Test Case 58

5.17 Predicted Velocity Vectors for Fine Grid with 59

Obstacles

5.18 Temperature Contours for Flow with Obstacles 60

5.19 Predicted Velocity Distribution for a Wider 61

Ventilation Outlet

5.20 Velocity Distribution for Ventilations with No 63

Diffusers

5.21 Calculated Velocity Vectors for Higher Ventilation 64

Inlet Flow Rates

5.22 Temperature Contours for Higher Ventilation Inlet 65

Flow Rates

ix

LIST OF TABLES

Table Page 1 Exchange Coefficient (r) and Source Term (St) 7 for Different 0 Variables

xi

1. INTRODUCTION AND SUMMARY

1.1 Background

This is the final report of the SBIR Phase I (six-month) study, entitled: "A Computer Code for Fire Protection and Risk Analysis

of Nuclear Power Plants". The study was sponsored by the U.S.

Nuclear Regulatory Commission (NRC Contract No. NRC-04-87-374)

and conducted by CFD Research Corporation (CFDRC Project No.

4025).

The overall objective of the proposed research program (Phases I,

II, and III) is to provide NRC, utilities, and affiliated

organizations, with a computer program capable of predicting

thermal environments, of postulated fires, inside control rooms,

cable rooms, and other enclosures of nuclear power plants. The

nature of the analysis is deterministic (rather than probabilis- tic) and the selected approach is "field modeling" (rather than

"zone modeling"). In general, the field modeling approach is

known as the Computational Fluid Dynamics (CFD) approach.

Typical results of such analyses include three-dimensional,

time-dependent distributions of velocity, temperature, pressure,

and concentrations of selected species. Such flow details are

important for:

1. better understanding of fire environments and their

sensitivity to various geometric and ventilation

parameters of complex enclosures;

2. improved basis of regulatory requirements (such as

those in Appendix R to 10CFR50 and Section 9.5.1 of the

Standard Review Plan [1]); and

3. improved design and operation of complex enclosures and

their ventilation and fire protection systems.

In the past, several field modeling studies (2-8) have been

performed, either with general-purpose proprietary CFD codes or

with research codes at Universities. The purpose of the proposed

project is to develop and validate a CFD code, specifically designed for fire modeling problems, by using the latest relevant

numerical techniques and physical models. The code will represent the current state-of-the-art of field nodeling of fires, and

will be a practical analysis tool.

1.2 Objective of Phase I Study

The general objective of SBIR Phase I studies is to establish the

feasibility of the proposed approach. For the present study,

specific objectives (as stated in the Phase I proposal) are:

1. Design a code architecture to maximize modularity, portability, and ease-of-use.

2. Construct the basic code by meticulous coding of

selected techniques and physical models.

1
3. Perform check-out and verification calculations.

4. Prepare the final report.

1.3 Results of Phase I Study

The Phase I study has met all of the above-stated objectives. A

three-dimensional, transient fire analysis computer program has

been developed and demonstrated up to the planned stage. It

consists of three modules as shown in Figure 1.1.

Pre-Processor Main

Module >" Code

~Post-Processor

Graphic Module

Figure 1.1 Main Modules of the Fire Code

The pre-processor module provides the means to set up a problem

in either "interactive" or "batch" mode. It has provisions for

geometric prescriptions of complex enclosures, multiple internal

solid objects, partitions, and multiple inlets and outlets of the

ventilation system. The main code consists of about 20 major

subroutines and several subprograms, so as to maximize the

modularity of the code. In Phase I work, maximum emphasis has

been placed on the selection and implementation of numerical

schemes, as these form the foundation of the code. The selected

numerical schemes include:

* A pressure-based, finite-volume solution method.

* A fully implicit formulation, for efficient simulation

of transient flows.

* First and second order accurate differencing schemes

(the latter is implemented, but not yet tested).

* Iterative, sequential solutions of dependent variables

(Each variable is solved for over the whole calculation

domain (rather than over a plane, line, or point)

before the solution of the next variable starts].

Use of CFDRC's latest whole-field equation solver with

extended, symmetric links (the use of an efficient and

robust linear equation solver is important for the

overall computational efficiency of solutions).

• Use of modified SIMPLEC algorithm.

* Use of dynamic storage and vectorizable coding

practices.

In addition to the above numerical schemes, the code has

provisions for the following physical models.

2
k-c turbulence model with modifications for buoyancy

effects.

* Multi-scale turbulence models.

* Instantaneous (diffusion-controlled) reaction model

• One-step (global kinetics) reaction model, with

modifications for turbulence interaction, such as eddy

break-up model.

Multiple-step (two or four step) reaction model, with

provisions for production and reaction of intermediate

species such as CO, C2H,,, etc. As compared to one-step

models, a multiple-step model will be more realistic

for fires in small enclosures or in low ventilation

conditions.

Six-flux radiation model; this does not require gases

to be transparent to radiation, and therefore provides

a good basis for radiation in fire/smoke environment.

Transient heat conduction model, based on an integral profile method to calculate heat losses through thick

walls, with proper account of time-dependent temperature profiles and thermal inertia of walls.

The basic k-e turbulence model and the one-step reaction model

have been completed and used in the present demonstration

calculations. Most of the other physical models have been

incorporated into the code, but not yet debugged or tested.

Likewise, the post-processor module of the code is in its initial

stage of development, which sufficed to produce basic graphics capability of plotting velocity vectors and contours of temperature, etc., as presented in Sections 4 and 5 of this report.

To verify the developed code, several bench-mark problems have

been simulated. These problems include:

1. Driven Cavity Problem

2. 2-D Plane Flow Over a Backward Step

3. Developing Turbulent Pipe Flow

4. 2-D Axisymmetric Flow Over a Backward Step

5. 2-D Natural Convection in a Square Enclosure

6. 3-D Natural Convection in a Closed Enclosure

7. Shear Mixing Layer (Propane-Air Non-Reacting)

Each of the above problems includes at least one critical

physical phenomenon of typical fire problems. Results of these

simulations have been found to be satisfactory (i.e. in agreement

with benchmark data).

For the demonstration of the code's capability of simulating

nuclear power plant fire problems, the code was also applied to a

recent (1986/87) fire test problem (selected from the fire tests

conducted at the Factory Mutual Research Corporation, under the

guidance of Sandia National Laboratory on behalf of the U.S.

Nuclear Regulatory Commission). Reference 9 provides sample data

for three tests using propylene gas as fuel. The test enclosure

3
is 60'x40'x20' and has a ventilation system with one exhaust and

six inlet ports. Test 14 has been used as the base case in the

present study. The ventilation rate for this case is one room

change per hour. At present, for demonstration purposes, walls

have been assumed to be adiabatic, and a coarse grid has been

used. The three-dimensional calculations have been made for 0-6

minute durations. Fire has been simulated by prescribing mass

source of fuel (rather than directly prescribing heat release

rate) and using a one-step chemical reaction model. The initial

results show qualitatively plausible flow distribution. However,

the predicted temperature profile, along a vertical line near the

fire source, did not show very good agreement with the measured

data. To determine contributing factors for this discrepancy,

several two-dimensional investigatory calculations have been

performed. It is established that finer computational mesh and

improved boundary conditions (e.g. convective and radiative heat

transfer to walls) are needed to improve the agreement of

predictions with experimental data, to a satisfactory level.

Improvements in other physical models are also important,

however, for the selected test case, these will contribute to a

relatively smaller extent. Several additional two-dimensional

calculations have also been made to examine the sensitivity of

solutions to: (a) internal solid obstacle; (b) wider exhaust

outlet; (c) alternative practice of simulating ventilation air

inlet; and (d) increased ventilation flow rate. Results of these

parametric calculations showed physically correct trends, and

confirmed the basic soundness of the developed fire code.

1.4 Recommendations for Phase II Work

The work completed in Phase I has clearly established the

feasibility of an advanced state-of-the-art fire simulation code.

To convert this code into a reliable and usable code (i.e. fully

validated and well-documented code), the following work items are

recommended for SBIR Phase II study.

1. Completion and testing of the physical models (described in Section 1.3).

2. Systematic parametric studies to establish the

sensitivity of results to computational mesh, and

various physical models.

3. Validation of the code against several (at least six)

test cases of the 1986-87 NRC/SNL/FMRC experiments. If

necessary, additional data from earlier tests (published by SNL, NBS, EPRI, FRS (UK), etc.) will also be

used.

4. Finalization of the code, with emphasis on the ease-ofuse, modularity, and portability.

5. Documentation of the code (including all assumptions,

theoretical models, code structure and user's guide).

A detailed work plan with time and cost estimates will be

submitted to NRC in the SBIR Phase II proposal.

4
1.5 Recommendations for Phase III Work

After the satisfactory completion of Phase II work, to effectively utilize and disseminate the developed code, the following

three activities are recommended to be carried out with commercial and/or non-SBIR government funds.

1. Incorporation of selected data bank, i.e. physical

properties and models for typical fire sources,

pertinent to nuclear power plants.

2. Enhancement of graphical display of results (improved

color graphics and animation of results) so that the

time evolution of the fire environments can be

presented (and recorded) in the form of movies and/or

video cassettes. Such films (real-time results) will

enhance the general use of results, and will also be

very useful in assessing various safety aspects such as

available evacuation time.

3. Development of technical training and user support

material, and a provision for code maintenance and

necessary enhancements.

For items 1 and 2 above, necessary support will be sought from

commercial organizations such as: Factory M4utual, Ove Arup &

Partners, and IBM. In the past, these organizations have shown

positive inclinations, however, due to the uncertainty and delay

in NRCI's request for 1988 SBIR proposals,, no formal agreements

have been established.

For item 3 above, CFDRC plans to provide the necessary support

from its internal funds.

1.6 Outline of the Remainder of the Report

Following this Introduction and Summary section, Section 2

discusses the mathematical background of the numerical and

physical models and solutions procedure. Section 3 describes

salient features of code structure and overall logic. Section 4

describes verification test cases and results. Section 5

describes the results of a fire demonstration problem (selected

from 1986/87 fire tests of FM4RC/SNL/NRC) and several investigatory test cases to support the analysis of results. Conclusions

drawn from the present study and recommendations for future work

are presented in Section 6.

2. MATHEMATICAL BACKGROUND

2.1 Introduction

The fire code predicts fluid flow, heat transfer, and chemical

reaction processes in a three-dimensional cartesian or cylindrical domain by solving a system of governing transient, partial

differential transport equations. objectives and capabilities of

5
the code have been described in Section 1. This chapter

describes the mathematical formulation of the selected equations

and the assumptions in physical models and relevant numerical

techniques.

2.2 Governing Equations

The general form of the conservation equations is as follows:

ao+ div (pvh) = div (rff grad 0) + So (1) at

Where 0 stands for general conserved property like the velocity

components (u, v, w), the enthalpy (h), the turbulence kinetic

energy (k) and its dissipation rate (c), the mass fraction of the

fuel (m..), and various species concentration Ci. All variables

represent time-averaged, rather than instantaneous, values of a

turbulent flow. The variables pi , and So stand for density,

velocity vector, and source term for 0 per unit volume respectively. Finally, r,*f is the effective exchange coefficient for

the transport of property 0. The values for different rff and S

for different O's are listed in Table 1.

In the present study, variable-density forumulation is used, i.e.

the Boussinesq approximation in not used. Local density is used

in all convection and source terms. The buoyancy term g(pp,,f) in the v-momentum equation implies that the reduced (or

relative) pressure practice is used. For buoyancy-dominated

flows, this practice is known to be superior to the practice of

using absolute pressure (or pg as the buoyancy term). The local

density, p, is calculated from the equation of state as a

function of local temperature and local species composition.

The following three equations apply to Table 1 and define the

effective and turbulent viscosities and the turbulence generation

rates Gk and Gb:

Atf =t + is (2)

t = CA p k2 /e and (3)

2 2 2

Gk " t 22 )2 + faV2+ . r i] r__xi +ay rr... ir

2 2 2 a

+ + 2(4)

+ [(a) 2, + fav 2'1 ra-y (axi )

6
Table 1. Exchange Coefficient (1) and Source Term

(SO) for Different 0 Variables

Dependent Exchange Source Terms

Variable Coefficients

r

1 0 0 (Continuity)

u efax +ef ax /Jef , rP~f f + a+ Ae ,f f TT-

+ a a av

V neff - - g P - PrI.f + - (p'r Ay

a+. ( a [ a z,] ( o)

W neff au a,

+"•" ." ) +' ' a 

h 4- '''

w /A~eff a 1of 7z

k f- Gk + Gb P4

ok

ar e

Ci Pe-f-f Ri

Ci

7
b ffi -g (5)

The term Gb, in the sources of k and e equations, represents the

buoyancy effect on turbulence [6]. This term is incorporated

into the code, but was not activated in the present development/

demonstration computations, reported in Sections 4 and 5 of this

report.

2.3 Description of Physical Models

The physical models of turbulence, heat transfer, combustion, and

radiation selected and incorporated into the fire code have been

outlined in Section 1. The Salient features of these models are

described below.

Turbulence Model: Two turbulence models are incorporated in the

fire code:

* The k-c model of turbulence [10], with modifications

for buoyancy terms [5].

* The multiple scale turbulence model of Chen and. Kim

[11].

The basis of the k-c model has been well documented [10] and the

model has been verified on a variety of nonreactive and reactive

flow problems. It involves solution of two transport equations

for: kinetic energy of turbulence, k, and turbulence dissipation

rate, e. The source terms for these equations, as shown in Table

1, involve semi-empirical constants which are assigned the

following standard values.

C, = 1.44,

C2 = 1.92,

CA = 0.09,

Ok = 1.0, and

a = 1.3

The major shortcoming of the k-7 model is its inability to

predict the response rate to changes in main field strain rates.

Computational study for flows in homogeneous turbulence in the

uniform shear flow showed that the k-e model correctly predicted

the asymptotic level of turbulence, but the predicted time to

reach that level was an order of magnitude longer than the

experimental one.

Recently, several works on the multiple scale turbulence (MST)

model have been conducted. One of the most advanced ones has

been recently presented by Chen and Kim [11]. The MST model subdivides the kinetic energy of turbulence into two regions: The

large scale (or production) region, kP, and the small (dissipaa
tive) scale region, k . The partitioning, as shown in Figure

2.1, is determined as a part of the solution and depends on local

turbulence intensity, production, transfer, and dissipation rate.

The partition is moved into higher wave numbers when production

is high, and to low wave numbers when production vanishes. The

source terms for the set of four equations are:

Sk = G - pep p

Sk t = p t

S = 2 + Cp 2 G - Cp 3 pep) (6) SCpz CpP1 ep

Se t = P Ct I Cp +t p Ct3 e t

1 6 Et )

Figure 2.1. Partitioning of the Turbulence Energy Spectral

Density in the Multiple Scale Turbulence Model

The total kinetic energy of turbulence is calculated as the sum

of both scales, viz:

and the constitutive equation for turbulent viscosity is:

(7)

Pt = C;,p k 2 /ef (B)

The turbulence constants for the MST model are:

kp'kt. 7 5; p= get .75; Cp = 0.09; Cp = 0.21;

Cp2 = 1.24; Cp 3 = 1.84; Ct, = 0.29; Ct 2 = 1.28; Ct 3 = 1.66.

The MST model has been tested so far only for nonreactive flows.

It is currently being verified for reactive flows in conjuction

9
with the eddy breakup model. It is planned to use the model in

the second stage of the project.

Heat Transfer Model: The basic energy equation is the conservation of stagnation enthalpy h which is defined as:

h = h + h (U 2 + v 2 + w2 ) (9)

The second term on RHS is the kinetic energy and it is generally

important only for high-velocity (compressible) flows, for the

fire code this term is neglected. The static enthalpy, h, is

defined as:

h = Z mj hj (10)

where mj is the j-th species mass fraction, and summation is

taken over all species within the mixture composition. The

static enthalpy of individual species is prescribed as input or

can be computed within the code based on the JANNAF tables, as:

h, =ho +fT CpdT (11)

f.298K 298

An iterative process is required to compute the temperature, T,

from known enthalpy, h, because of the nonlinear relationship

between them.

At present, thermal boundary conditions can be either of prescribed heat flux or prescribed wall temperature. An important

addition, to be made in the Phase II of the project, is to

incorporate a transient heat conduction model for accurate

calculation of wall heat losses through solid walls. This model

will be based on integral profile method and has already been

successfully demonstrated in an earlier NRC/BNL study performed

by the principal investigator of this project.

Combustion Model: The principal equation used in the combustion

model is the conservation equation for chemical species which has

the form:

aaE (pmj) + div(pulmj) = div(r,ff.j grad aj ) + Rj (12)

where mi is the mass fraction of chemical species J, R, is the

mass rate of reaction of species J due to combustion per unit

volume, and r~ffj is the effective exchange coefficient for

species j.

At present, there are two different options considered for the

reaction rate Rj:

- One-step Arrhenius reaction-model; and

- Eddy-break-up model.

10
The Arrhenius reaction rate for fuel can be written as:

Rf= -pAMF (P MFU ( PCX exp -(13)

where A is the reaction rate constant, E is the activation

energy, R is the gas constant, MFu and Mox are molecular weights

of fuel and oxygen, and a and f are reaction order constants.

This model is mainly used for kinetically controlled chemical

reactions.

For most of the turbulent flows in which mixing is the controlling parameter, the eddy-break-up model of Spalding [12] modified

by Magnussen [13) is usually used. The reaction rate for fuel is

written as:

Rf = -CR Pgh e/k (14)

where CR is an empirical constant and g represents the local

mean-square concentration fluctuations, which is expressed

algebraically as follows (13]:

gk = min m U, (15)

where B=4.5 is a constant.

Provisions for more complex reaction models such as two or fourstep reaction models with intermediate species such as CO, C2H4 , etc. are also needed in the fire code, particularly for underventilated and/or small room fires.

Successful modeling of combustion systems depends on a correct

description of the chemical reactions involved. While the

complete chemical kinetics of CHm oxidation are not fully

understood, there are models available which are more sophisticated than those currently used in the present demontration

study. One of these is a four-step kinetic mechanism for

hydrocarbon combustion. The results of Hautman [14] indicated

the viability of representing the C.H2•÷a - aliphatic hydrocarbon

oxidation in a simplified kinetics scheme in the form:

C. H2.+. 2 - C2 H, + H2 (16)

C2 H2 +02 - 2CO + 2H2 (17)

Co + 1⁄2 02 - Co 2 (18)

H2 + 02 - H20 (19)

11
where CZH. represents a class of olefinic intermediates. The

above reactions are controlled by reaction rate constants which

are expressed as:

Ru = -:0L [fu]a [0 2 ]b [C2HA]J exp (-E/RT) (20)

Rc2H4 = -i0 [C 2 H4]a (O2]b [fu]c exp (-E/RT) (21)

RcO = -10x [Cola 0 2 ]b [H 2 0]© exp (-do) exp (-E/RT) (22)

RH2 = -1ox [H2 l] (O]b [CHA]O exp (-E/RT) (23)

where x, a, b, c, and d are constants which are different for

each reaction, o is the equivalence ratio, E is the activation

energy, and R is the universal gas constant.

In the second phase of the project, various advanced reaction

models will be reviewed again and suitable one(s) will be

utilized.

Thermal Radiation Model: The radiation transport equation for an

absorbing-emitting gray medium in local thermodynamic equilibrium

can be written as [15]:

d . AIr = - aI. + a7 -74daE (24)

where d is the unit vector representing direction of radiating

beam, I is the radiation intensity, a is beam absorption, s is

scattering coefficient, E=aT4 where a is the Bolzman Constant.

A flux model is used to solve the above integro-differential

equation where the flux model is the assumed radiation of intensity with direction. Assumption of constant I in a quandrant

centered along each (±) coordinate direction results in six

differential equations for intensities. Therefore, this model is

known as the six-flux radiation model. This model has been used,

successfully in several reactive flow simulations [16, 17].

2.4 Description of Numerical Procedure

A finite-volume approach based on the modified SIMPLEC algorithm

[18, 19, 20] has been used tn the fire code. The code has

provisions for cartesian (x, y, z) and cylindrical (x, y, 0)

coordinates. Grid distributions can be non-uniform to accommodate finer grid spacings in the regions of steep flow gradients. Internal solid blockages due to, say, cabinets, cables,

partitions, etc. are represented by using the "partial porosity"

(or "porosity-resistivity") approach.

12
The fire code uses a staggered grid practice in the so-called

"backward boomerang" mode, shown schematically in Figure 2.2.

The velocity components are calculated at the cell faces while

all other scalar variables are calculated at the grid center.

f vi li+# I vii+1

Asij* I

i-2i Uj-lj Ij Uiii Ui j i+li

Aw A AwA+.j

___________vi-1i WSJ Iv 1i

L XAsij

Uq.... il-I

Figure 2.2 Staggered Grid Notation in the x-y Plane

There are six types of geometrical boundary conditions that can

be specified in the fire code.

- Fixed mass flow rate inlets that can account for forced

or free multiple ventilation inlets; - Fixed pressure inlets or exits; - Solid wall boundaries; - Symmetry planes; - Internal mass sources without momentum; and - Zero gradient boundaries such as adiabatic walls.

Special attention is required for the near wall boundary

conditions. For velocities, classical wall functions, based on

the logarithmic low of wall, are used. Near wall turbulence

parameters, that are based on the assumption of local equilibrium

between total turbulence kinetic energy and second scale dissipa- tion rate, are assumed. For the energy equation, prescribed wall

temperature or heat flux boundary conditions are used. As

mentioned earlier, in Phase II, a transient heat conduction model

for solid walls will be incorporated. This will enable coupled

solution of heat transfer between gas and solid regions.

Finite-Volume Differencing Scheme: The partial differential

equations (1) are integrated over the finite volume, and for

transient calculations, over the finite time interval to obtain

the descretized equations. A hybrid formulation is used for

convection and diffusion terms. This formulation has proven to

be very robust, but does lead to significant numerical diffusion

13
of solutions. Therefore, a more advanced (second-order accurate)

scheme, i.e. central differencing with second and fourth-order

dissipation terms, is also included in the code. The secondorder scheme will be tested in Phase II work.

For transient calculations, a fully implicit formulation is

employed where the values of flow variables are taken to be those

which prevail at the current time step (backward Euler formulation). Integration of the partial differential equaitons results

in linear algebraic equations of the form:

a = a1• + awOw + as$s + SUO + le (25)

where az, aw, ... , etc. are called "link coefficients", SU and SP

are the linearized sources, Ml is the accumulation term, and ap

(diagonal coefficient) is expressed as:

ap = a. + aw + a. + a. - SPO + MO (26)

Solution Procedure: A set of linearized equations, one equation

for each variable at each grid node, needs to be solved for each

time step. These equations are coupled with each other because

of the presence of several dependent variables in each equation.

Also, since the original partial differential equations are nonlinear, the linearized equations need to be solved iteratively

with under-relaxation. The required under-relaxation and number

of iterations strongly depend upon the manner in which the

equations are solved. In the present code, the following

practices have been chosen.

1. Equations are solved in a successive (variable-byvariable) manner.

2. Each variable is solved for over the whole calculation

domain (rather than over a plane or a line) before the

solution of the next variable starts.

3. Provisions are made for the solutions by:

a) an iterative whole-field (3-D) solver, and

b) Jacobi Point-by-Point solver.

The selected whole-field solver [21] is a modified version of

Stone's Implicit' Procedure (22], and is considerably more

efficient especially for ill-conditioned and non-symmetric

matrices. At least, the pressure-correction equations are solved

by using this equation solver. This equation plays a central

part (and consumes major computational effort) in the selected

solution scheme. The solution of other variables, such as

velocities, can be equally efficient (depending upon the various

versions of the SIMPLE algorithm) by the whole-field or point-bypoint solution.

An inertial under-relaxation practice is employed for all

dependent variables except pressure. This practice is achieved

14
by adding an inertia term, viz: I4(q-9") to the finite difference equation in the following manner:

ZaO, + SU + 1 4*

= a. - SP + 10 (27)

where suffix ,n, denotes all cell neighbors, and superscript , , denotes the previous iteration value of 9- The "inertia term"

,Io, is calculated as:

p VOL (28) AtF 0

where p is the fluid density, VOL is the grid cell volume and

AtFý is the "false" time step specified for each dependent variable 0. Typically Att is specified to be of the order of

resident time of fluid in a typical control volume. It should be

noted that in converged solutions 0-0*, and therefore the

addition of the inertia term does not affect the final solution.

In order to monitor the convergence of the iterative solution

procedure, the residual errors are calculated for the equations of all dependent variables and at each control cell, in the

following manner:

cop = EanOn + SU - SPOp, - (Ma.) p*. (29)

where top stands for the residual error in the equation of

variable 0 at point p, and summation is taken over all link

coefficients: N, E, W, S, .. , etc.

The residual errors are used as indicators for the convergence of

the solution. Generally, the sum of absolute values of residuals

for each variable are required to reduce approximately by three

orders of magnitude. In Phase II, after the experience of

several fire test problems, specific recommendations and built-in

convergence criterion will be provided.

3. CODE STRUCTURE

During the process of designing the structure of the fire code,

maximum emphasis was placed on the modularity and the ease of

use. For these purposes, the computer code was divided into

three interconnected modules as previously shown in Figure 1.1.

The pre-processor module provides the means to set up the test

case in either interactive or "batch" mode. It allows specifica- tion of the geometry, grid, and boundary conditions, solution and

printout control, etc. Provisions are also made in the preprocessor for providing physical properties data. The pre-processor

overall structure is shown in Figure 3.1.

15
INPUT

ASCII

FILE

COMMAND

PROCESS ING

MODULE

CONTROL

* Checkout

OLD INPUT

FILE FOR

MAIN CODE

MODULE

INPUT FILE

FOR MAIN CODE

Figure 3.1 Structure of the Pre-Processor

16
The main code module (Equation Solver) includes a large number of

subroutines, each with specific functions such as treatment of

boundary conditions, grid and geometry, initial field, physical

properties, physical models, solution algorithm and output

control. Major parts of physical models are coded in selfcontained, modular subroutines. This module also includes a

subroutine for dumping the data to files for further (restart)

analysis and/or graphical post processing of selected flow

fields.

The post-processor module is prepared as a separate program that

has access to the data files generated by the main code and the

pre-processor modules. Graphical display of grid and geometry, velocity vectors, contours of all variables, and profiles are

available through an interactive mode in the post-processor.

An overall structure of the main code and the post-processor is

shown in Figure 3.2. Further details, such as functions and

detailed flow charts of all major subroutines will be provided in

the final report of Phase II work.

4. CODE VERIFICATION STUDY

4.1 Overview

The selected solution scheme and physical models have been used

successfully in many codes, including several written by the

members of the present project team. However, before a code can

be used for engineering analysis, it must be checked out/verified

systematically for a large number of problems of increasing complexity. This section summarizes the verification problems which have so far been simulated with the fire code.

1. Two-dimensional driven cavity flow with different

Reynolds numbers.

2. Laminar flow over a backward-facing step. 3. Developing turbulent pipe flow.,

4. Two-dimensional axisymmetric turbulent flow over a

backwark-facing step. 5. Two-dimensional natural convection in an enclosed

cavity, with different Rayleigh numbers (10 3

s<Rasl50).

6. Three-dimensional natural convection in an enclosed

cavity (room with hot and cold walls).

7. Shear mixing layer (propane-air non-reacting).

All cases were selected because each of them contains at least

one flow feature which is likely-to be present in fire problems.

For cases 1, 2, 3, 4, 5, and 7 bench-mark data are available for

quantitative comparisons. Case 6 was performed mainly to check

the code capability for three-dimensional buoyancy-dominated flow. This case is designed to have symmetric boundary conditions such that the calculated flow field can be examined for the

symmetry conditions.

17
Figure 3.2 Structure of the Main Code

18
In the following sections, the verification problems are briefly

described and their results are presented.

4.2 Two-Dimensional Driven Cavity Flow

The driven cavity problem is a good test case for two-dimensional

elliptic flows. In this problem, the top wall of the cavity

moves at a constant velocity while the remaining three walls are

stationary (no-slip boundaries). The geometry is illustrated in

Figure 4.1. A sample computational grid with uniformly distributed 80 grid cells in the x-direction and 80 grid cells in

the y-direction is shown in Figure 4.1.

Computations have been performed with several grids in the range

of 10 x 10 to 120 x 120, with uniform and non-uniform distributions, for three Reynolds numbers (i.e. Re=100, 1000, and 10000).

Results for Re=1000 case are presented in Figure 4.2. The

streamline contours show the expected recirculation pattern with

a big central vortex and two secondary weak vortices in the

corners. These results have been compared to the bench-mark

solutions reported by Ghia, et. al. [23], who solved the same

problem with 257 x 257 grid points. The comparison between

Figures 4.2a and 4.2b indicate a good agreement. Further

improvement in accuracy is expected with the use of a secondorder differencing scheme (in Phase II). Similar agreements were

obtained at the two other Reynolds numbers.

4.3 Two-Dimensional Laminar Flow Over a Backward-FacinQ Step

In this case, analysis of laminar flow over a backward-facing

step with 1:2 expansion ratio was performed. The geometry and

computational grids are shown in Figure 4.3. The computational

grid used consisted of 120 grid cells in the x-direction and 40

grid cells in the y-direction. The grid was clustered near the

step and in the expected recirculation zone. The inlet location

was chosen such that, based on a uniform inlet velocity profile,

the flow is fully developed when it reaches the step.

Results were obtained for Re=100 and Re=378. For Re=100, the

predicted reattachment length was found to be x/h=2.7 which is

within five percent of the experimental data reported by Armaly,

et. al. [24]. For Re=398, the fire code overpredicted the

reattachment length by about eight percent. The axial velocity

profiles for both cases (Re=100, 398) are shown in Figure 4.4 and

they indicate good agreement with the data reported in Reference

23.

4.4 Developing Turbulent Pipe Flow

In this verification case, the developing turbulent flow in a

pipe is simulated. The geometry of the pipe is depicted in

Figure 4.5. The length of the pipe, L, was chosen such that the

19
"wall

(a)

(b)

Figure 4.1 Geometry and Computational Grid for Driven

Cavity Problem

20
(a) (b)

Stram function

Contour Contour

Lettcr Value or w number Value or iv

a -1.0 x I0",s 0 1.0 x 10"-

b -. 0x 10" I 1.0X 10"

c -LOx 10"' 2 .OX I0x - d -L.Ox 10"6 3 l.Ox I0"

C -0.0100 4 S.0x tO0'

f -0.0300 5 1.Ox I0"-

£ -0.0500 6 2.5 X 10-4

h -0.0700 7 5.O x I0-,

i -4.0900 9 IOx I0"

j -0.1000 9 I.5 x I0"-

k -0.1100 I0 3.0 X 10"$

I -0.1150

m -0.1175

Figure 4.2 Streamline Contours from, a) Present Fire Code;

b) Ghia, et. al. [23]

21
ge '..- '0

-. 

H I 0=N

(a)

* Note that both graphs are not to scale

(b)

Figure 4.3 Magnified Inlet Region of Backward-Facing Step;

a) Geometry, and b) Computational Grid

22
020 02 

V

O050W 045 01 015s 02

V

Re = 100 Predicted

v Experimental (24]

0014-up 9 /h 4 18

* I~

fuels I I II 4 I II

-02 09 02 04 06 Os

V

Re = 398

Figure 4.4 Comparison of Predicted Axial Velocity Profiles

with Experimental Data of Armaly, et. al. [24]

23
- - - - - - - - - - -- - -

f

I

F- r L I -H

Figure 4.5 Geometry of the Developing Pipe Flow Problem

0.

S.

- Predicted

Measured [25]

0

0

0

Veloo ity

23 30 35

Figure 4.6 Predicted and Measured Axial Velocity Profile

for Fully Developed Turbulent Pipe Flow

24
flow gets fully developed downstream. The standard k-c model was

used to obtain results for Re=Sxl0. The velocity profile in the

fully-developed region is presented and compared to the experimental results of Laufer [25] in Figure 4.6. The comparison

indicates satisfactory agreement.

4.5 Two-Dimensional Axisymmetric Flow Over a Backward-Facing

Two-dimensional non-swirling turbulent flow over an axisymmetric

backward-facing step is considered in this case, as shown in

Figure 4.7. Results were obtained, using the k-c turbulence

model for Re=1.25x10'. The axial velocity profiles in the fully developed region (before the step) and at x/h=2.0 are shown in

Figure 4.8 and compared to the experimental data of Samimy, et.

al. [26]. As shown in both figures, the fire code predictions

have good agreement with the experimental data.

4.6 Two-Dimensional Natural Convection

In this verification case, flow in an enclosed square cavity with

differentially heated walls is considered. The geometry and

boundary conditions of the computational domain are shown in

Figure 4.9. The horizontal (top and bottom) walls are insulated

(DT/ey = 0) while the vertical walls are held at constant

temperatures,

T = Th at x = 0

and

T = T© at x = L.

Results were obtained for Rayleigh numbers of Ra=103 , Ra=10 5 , and

Ra=107. Velocity vectors for Ra=105 are shown in Figure 4.10.

Streamline and temperature contours for three different Rayleigh numbers are shown in Figures 4.11 and 4.12, respectively. For

Ra=103, the streamlines are smooth and no distinct boundary layer

can be determined. As the Rayleigh number increases, secondary

recirculation eddies form as shown in Figure 4.10 (Ra=10 5 , Ra=107) and a significant temperature gradient develops near the

walls (Figures 4.12b and 4.12c).

The predicted behavior of the flow depicted in Figures 4.11 and

4.12 is identical to the benchmark results published in the

literature by G. DeWahl Davis [27] and others.

4.7 Three-Dimensional Natural Convection

The geometry and computational grids for this three-dimensional

problem are shown in Figure 4.13. The grid size employed consists of 30 grid cells in the horizontal x-direction, 30 grid

cells in the vertical y-direction, and 10 in the z-direction.

Boundary conditions are comprised of zero wall velocities and

25
U,

-N - -m mm - -m u ~t - -m

--1⁄2

Figure 4.7 Geometry of the Inlet Section of the 2-D Axisymmetric Backward-Facing Step Problem

Figure 4.8 Predicted and Measured Axial Velocity Profiles

26
8T-= 0

ay

Th Tc

8T

By 0

Figure 4.9 Grid and Boundary Conditions for 2-D Natural

Convection Problem

Figure 4.10 2-D Natural Convection Velocity Vectors for Re=105

27
(a)

(b) (C)

Figure 4.11 Streamlines in a Square Enclosure for a) Ra=103 ;

b) Ra=105 ; c) Ra=10O

28
.( a)

(b) (c)

Figure 4.12 Temperature Contours in a Square Enclosure for

a) Ra=103 ; b) Ra=106; c) Ra=-07

29
I

Cold

(x--o)

Cold

y

z

I

I

Hot

(z=o)

I

x-y Plane

I i

x-z Plane

y-z Plane

Figure 4.13 Grid and Geometry for 3-D Natural Convection

Problem

30
prescribed wall temperatures and zero heat fluxes, specified as

follows:

Cold Walls: T(x=O,y,z) = Te ,

T(x=L,y,z) = T,

Adiabatic Walls: BT(x,y=O,z) 0 0

By

8T(xy=M,Z) 00 

ay

Hot Wall: T(x,y,z=O) = Th , and

Symmetry Plane: 80 (x,y,z=N) = 0

az

As in the previous case, the flow here is dominated by the

buoyancy force. As the temperature near the hot wall (z=o)

increases, the fluid becomes less dense, creating a buoyancy force resulting from density difference. The buoyancy force acts

vertically upward and causes the less dense fluid to rise as

shown in Figure 4.14a. The opposite happens near the two

vertical cold walls because the fluid is more dense, therefore it

falls near both walls creating two circulating eddies near the

middle of the room as shown in Figure 4.14b. Since the flow

rises near the hot vertical plane, it falls near the symmetry vertical plane as shown in Figure 4.14c. Figure 4.15 shows the

velocity vectors in two horizontal planes where the flow is

directed toward the hot wall near the bottom plane and away from

the hot wall near the top wall. Figure 4.16 shows the velocity

vectors in three vertical YZ planes. Near the cold wall, the

flow is falling as expected while some recirculating motion is

observed in the middle of the room.

Temperature contours in three planes are shown in Figure 4.17.

The contours in the mid-vertical plane indicate that the hottest

fluid is near the middle of the room, away from the cold walls.

Figures 4.17b and 4.17c show the highest isotherms near the hot

wall as expected.

4.8 Shear Layer Mixing (Propane-Air Non-Reacting)

In this last verification case, the problem of turbulent nonreacting propane-air jets flow was considered. The geometry and

the computational grids are shown in Figure 4.18. The computational grid consisted of 20 grid cells in the radial direction,

where 5 cells were placed within the inner jet and 15 cells

within the outer jet. The horizontal x-direction was modeled by

31
XY 1

.511 11.5.

• *I I,.

• , t/I •_•..'('\\, , ...

• 5.51 o.oo . . ,...•

(a) Near Hot Plane

Xv 5 Xv to

i

1

i SI

(b). Mid-Vertical Plane (c) Near Symmetry Plane

Figure 4.14 Predicted Vertical (xy) Velocity Distribution

for the 3-D Natural Convection Problem

32
xz I Xz 30

* % % ,' %%%%\\ \ t ,//,//f/,,"/l, - - ,.

* . -. \..\'X\\\\.- \ % .111/1//,/////

* S % %%~%.'s%'%~s%~ I

* ~ I

* ~%\\'%\%WA\V~ I

* ~~%~\\\\\\\\\~ I

* ~~A\\\\\\\\ \~ I

* ~\\\\\\\\\\\~ I

~\\~N\V~ I

* %%~%'%NNNN%~%~\\~ 1

U ~

I ~////W~-,,,.

I lt//////,,~,..

I t/////////,i,.

I //////,~y~~,,

I / / ~

- ,...=== p, . , • % # 0 ' . .

Bottom Plane Top Plane

Figure 4.15 Predicted Horzontal Velocity Distribution

for 3-D Natural Convection Problem

vzI

t

•z %

S.. 'S 8~ ~

1 lb ~. \\\\

I I

I I

II

II

I I

Ii

It

I I

a i

I'

.8 5

£ 1/

il/I

a''''''

a Id///j I a * I I -- - ~' I I a

Yz 15

I . - -. ~. - - -

/ ,' z-..- - ".S

YZ2.Z

, - - • _.• • . S

I . " 

Figure 4.16 Predicted Vertical (yz) Velocity Distribution

for 3-D Natural Convection Problem

33
hot

a)~~~~~~~~ 

er l, pl nlyY t cve 

twx

xY 5

a) vertical (XY) plane, parallel

to hot wall

hot

hot

tL.x

YZ is

b) horizontal (XZ) plane c) vertical (YZ) plane,

parallel to cold walls

Figure 4.17 Calculated Isotherms in Three Mid Planes

34
,/ .10 " e " .41 ,e -.-, .0 ý/ " ý/ " /

SI

D

.1

/I / / / O1

- .1161 . -

Geometry

20 x 20 Grid

Figure 4.18 Geometry and Computational Grid for Shear

Layer Mixing Problem

35
20 grid cells placed uniformly. The inner jet velocity was 53

m/s while the outer jet velocity was maintained at 9 m/s. The

centerline variations in the mixture fraction are shown and

compared to the measured data of Schefer and Dibble [28] in

Figure 4.19. The comparison indicated good agreement upstream

and satisfactory agreement downstream.

5. CODE DEMONSTRATION FOR FIRE PROBLEMS

5.1 Problem Selection

All verification cases, described in the previous section, have

been of non-reactive flow nature. To demonstrate the code's

capability for reactive flows, similar to those in typical fire

problems, the following experimental investigations were considered as the potential source for a demonstration problem.

1. Cable burning experiments conducted by Factory Mutual

Research Corporation (FMRC) under the sponsorship of

the Electric Power Research Institute; EPRI Report NP2660.

2. Tests for Twenty-foot Separation Criteria of NRC,

conducted by Underwriters Laboratory in conjunction

with Sandia National Laboratory; Report NUREG/CR-3192,

SAND83-0306, 1983.

3. Tests for fire in compartments, conducted by the

National Bureau of Standards; Report NBSIR-82-250,

1982.

4. Enclosure Environment Characterization testing for the

Base Line Validation of Computer Fire Simulation Codes,

tests conducted by FMRC, under the guidance of Sandia

National Laboratories for the U.S. Nuclear Regulatory

Commission; Report NUREG/CR-4681, SAND86-1296, 1987.

Of the above, the last series of tests have been found to be the

most suitable because the tests were designed:

a) to simulate fire conditions of nuclear power plant

enclosures; and

b) to provide data for the validation of fire modeling

computer programs.

This test series included 22 tests, however, so far only limited

data (for 3 tests) have been analyzed and reported in Reference

9.

In the following sections, the test cases will be described along

with the numerical setup, and finally, the results will be

presented, analyzed, and compared to the experimental results.

36
I - Predicted

0 Measured [28]

0.

0

20 40

K/D

60

Figure 4.19 Predicted and Measured Mixture Fraction

at the Centerline

37
5.2 Test Case Selected (FNRC/SNL/NRC)

Experimental Setup: The test enclosure has a size of sixty feet

long by forty feet wide by twenty feet high (18 x 12 x 6 meters)

as shown in Figure 5.1. A forced ventilation system with six

inlet ports and one outlet port is used. The inlet ports

extended down through the ceiling of the enclosure to a point

four feet (1.2 meters) below the ceiling.. These inlets were

capped with standard, commercially available four-way air

diffusers. The enclosure surfaces were instrumented at a number

of locations for surface temperature, optical density, surface

heat flux, etc. The fuel/fire sources utilized in the tests

included gas burner fires, liquid pool fires, and simple solid

fuel fires.

A total of eighteen Base Line Validation tests were conducted in

the enclosure with no internal obstructions present, while

another four Base Line Validation tests were conducted in a

control room mock-up configuration of the test enclosure.

Initial test results are available [9] for cases 4 and 5 of the

non-internal obstruction cases, and case 21 which was a mock-up

configuration test. Test case 4 was selected as the base case

for code validation. Cases 5 and 21 have also been considered for

two-dimensional analysis with the fire code. These two-dimensional analyses will be presented later in this section; case 4

and its two-dimensional results are described first.

Case 4 utilized a gas burner fire source located near the middle

of the room. The burner utilized was a 36-inch diameter sand

burner. The gaseous fuel was forced to flow up through a base of

loose sand filling the burner body. The nominal peak value of

the heat release rates was measured to be 516 kW. The location

of the fire along with the location of the instrumentation, for

which results were presented, is shown in Figure 5.2. Figure 5.3

shows the heat release rate profile for the gas burner that was

utilized in case 4. The numerical setup of case 4 is presented

in the following section.

Grid and Boundary Conditions: Due to the symmetry of the

geometry of the enclosure and the location of the fire (middle of

the room), only half of the room was numerically modeled. A

relatively coarse grid of 1980 control cells (9 in the horizontal

x-direction, 11 in the vertical y-direction, and 20 in the other

horizontal z-direction) was used. The grid was constructed in

such a way that the instrumentation locations shown in Figure 5.2

would coincide the grid nodes. The ventilation inlets were

placed in three cells located 1.2 meters below the ceiling. Each

grid cell has an area of 0.46 x 0.46 M2 , which corresponds to

that used in the experiment. The ducts that are used for air

passage inside the room were modeled by blocked cells above the

ventilation inlets. The three planes of the computational grid

are shown in Figure 5.4.

38
.4'

-4:

% %

~ C

(\j ~q

A

Top

View

(

Exhaust Duct Inlet Duct (typ. one of two)

Si de7III_ t

View

t..o CU!

0 ON/M,

60

Figure 5.1 Enclosure Ventilation System of FMRC/SNL/NRC

(1986/87) Tests

39
Vent Inlet Ports -(19"'7)

-(18°)

--(14')

-- +-(lCO)

Exhaust Port

Fire Source

(6')

Instrument

0 Vertical Location

of Instrumentation

Note: a is the instrument location where the temperature

measurements are available [8]

Figure 5.2 Location of Fire Source and Ventilation Inlet and Outlet Ports
1000

800

M600 16k-W PeakHRR

400

200

0

8 2 4 6 8 10 12 14

Time (min)

Figure 5.3 Measured and Simulated Heat Release Rate Profile

41
z Vertical y-z Plane

x I I

•z Horizontal x-z Plane

y L Vertical x-y Plane

Figure 5.4 Computational Grid Chosen for the Demonstration

Calculation

42
The calculations were performed in two steps. First, steadystate (initial t=O, pre-fire) flow distribution is calculated

with the following boundary conditions.

* Ventilation inlet mass flow rate of 0.076 kg/sec.

* Exhaust port zero relative pressure (Poxit = 0).

* Adiabatic walls.

Then the transient analysis is performed with the same boundary

conditions and the fire source modeling. The transient combustion calculations were performed with the time intervals of At =

15 seconds.

Fire Source and Combustion Setup: In the experimental test case

considered, the (non-premixed) gas propylene (C3 H6 ) fuel was

supplied through the sand burner. A diffusion flame resulted

with a 36-inch base (burner diameter) extending upward to the

ceiling of the enclosure. Experimental data for the heat release

rate, HRR, shown in Figure 5.3, was used to compute the fuel mass

flow rate from the following relation:

=HRRr ci a =(30)

QL secj

where QL = 4.58 • 107 j/kg is the heat of combustion of propylene.

In the "growth-mode" heat release test, the resulting flow rate

is proportional to the square of elapsed time from ignition to

reach full intensity at four minutes. The peak value was

maintained for the remaining six minutes, and then full cut off

occured.

In the real diffusion flame environment the heat release rate

will be distributed within the flame zone and will extend up tothe ceiling of the enclosure. The spacial distribution of the

heat release rate was unknown. It was decided, therefore, to use

a simplified, global one-step finite-rate reaction model for the

fuel combustion in the following form:

C3 He + SO2 -- 3C02 + 4H20 (31)

Under real conditions, however, at relatively low temperatures,

the combustion process may result in much slower heat release and

several intermediate reactions, producing: C (soot), CO, C2 H4,, CH2 , H2 , etc. It should be indicated, therefore, that two types

of inaccuracies will result from the assumed combustion model:

1. the overall heat release will be larger in the computational studies than in the experiments. As a result,

the temperature levels may be over-predicted, and

2. single-step reaction will produce faster reaction rates

43
and may result in too large a heat release in the

vicinity of the burner.

For the purpose of the present study, an Arrhenius reaction model

was used to represent the reaction rate as:

RfU= _PA~p u (PMiU) C, [p OX) exp(-.j1 (32)

with rate constants A = 1.1010, and E/R = 18000.

In the second phase of the project, a more adequate multi (two or

four) step reaction model should be used as discussed in Section

2.3.

5.3 Results of 3-D Computations

Steady-State Results: The steady-state solution depicts the flow

distribution under the normal forced ventilation conditions in a

large, empty room with no fuel supply. Figure 5.5 shows velocity

vectors in three planes of the computational domain. Figure 5.5a

shows the velocity vectors at the vertical symmetry plane near

the exit. As expected, the flow is strongest near the exit and

shows an upward motion at the opposite side of the room, resulting from the third ventilation inlet near the end of the room.

Figure 5.5b shows the velocity vectors in the mid-vertical plane

where the ventilation inlets are located. It can be seen that

the air exits from the ports and flows in all directions,

resembling four-way diffuser outlets. The same phenomena is

observed in Figure 5.5c, which shows the velocity vectors in the

horizontal xz plane near the ventilation inlets.

Transient Results: As mentioned earlier, transient calculations

were performed with fifteen second intervals (at=15 sec) for a

total time of six minutes. Presented results from the transient

analysis include velocity vectors and temperature contours at one

minute intervals. Figure 5.6 shows the temperature contours in

the vertical yz plane, where the fire is located. It is observed

that as early as one minute, the hot temperature starts propagating in the room. For the first two minutes, a hot layer is

observed near the ceiling. This hot layer propagates and

descends into the room as time increases to five minutes. The

fire phenomena is observed in all four plots, where the flameshaped contours are observed above the fire. The temperature, as

expected, increases with time due to the addition of heat, as

shown in the profile in Figure '5.3. Figure 5.7 also shows the

temperature contours in the vertical xy plane at the fire

location. The same observations made earlier apply to these

plots. The hot layer is clearly seen in the first three minutes.

The layer descends into the room and the temperature near the

fire increases to above 330"C at the four minute time.

44
(a)

.XZ Plane 1

Vmin z 0.00014 m/z

Vmax z 0.36 m/a

YZ Plane 6

Vmin a 0.0 m/5

Vmax = 0.057 m/s

* S * - - .~ -- -

L{- I ~ - *. * b/%\:jh$V~<

,' ' - - - - - ~ I - I ' - .~ / .* ~

St t - a - - - - - - * - - - - , , ,

(b)

_ - e.~ .= -

.- . .. " • ,. '. ., 

•

-- - : - . . . . ,, I ..

XZ Plane 8

Vmin = 0.0002 m/a

VMax = 0.049 M/8

(c)

Figure 5.5 Calculated Velocity Distribution in the Test Room

for Initial (Steady-State Pre-Fire) Condition

45
Time 1 1.0 Min

113.95 

I

30.65 IC

Time a 3.0 Min

192.15 "C62.25 It '

371.95 ICTime 

a 5 .0 Min

3 "71 .8 5 "C - _

Figure 5.6 Temperature Contours at Near-Fire Vertical

y-z Plane at t = 1, 2, 3, and 5 Minutes

46
22.35 It2

113.95 *C

Time a 2.0 MiJn

331.~o65 IC

192.15 ..- 0

Time a 1.0 Min Time 3.0 Min

Time = 4.0 Min Timm a 5e. MIn Time a 6.0 Kin

Figure 5.7 Temperature Contours at Near-Fire Vertical x-y

Plane at t= 1, 2, 3, 4, 5, and 6 Minutes

47
Figures 5.8, 5.9,, and 5.10 show the velocity vectors in three

vertical yz planes. Figure 5.8 shows the vectors in the first

vertical yz plane near the fire. At one minute, it is observed

that the air is entrained into the fire and accelerated vertically upward toward the ceiling, while the flow has no uniform or

apparent pattern away from the fire. The same phenomena persists

at the two minute time with a more uniform development of the

flow. At four minutes, circulatory (eddy) motion is observed in

different places of the room, with the biggest eddy located near

the exit. Thie same phenomena is observed at five minutes.

Figure 5.9 shows the velocity vecto rs in the mid-vertical yz

plane where the ventilation inlets are located. The same overall

phenomena observed at the first vertical plane is observed here

(viz: the circulatory motion or several flow eddies). The

effect of the ventilation is observed in the first two minutes

where the velocity is directed downward ano away from the inlets

with the circulation eddies being pushed toward the floor of the

room. At four minutes, the effects of the ventilation disappears

as the flow near the inlets is dominated by the circulation

motions created by the fire and no motion from the inlet ports is

observed.

Figure 5.10 shows the velocity vectors in the last vertical yz

plane,, at the end of the room. The circulation motion of the

flow, which is a characteristic of fire flows, is still apparent

and has a more pronounced shape at later times (four and five

minutes).

Figure 5.11 presents a summary plot of the predicted and experimentally measured profiles at a vertical line located at the

symmetry plane of the room close to the burner (the triangle in

Figure 5.2). Results are compared for 2, 3,j 4,, and 5 minute

intervals.

The overall tendency of the temperature pattern is similar in

both computations and experimental measurements. At t=2 minutes

after the ignition,, the temperature profiles are very similar,

with only 2*C over-prediction near the floor (6' level), and 8*C

over-prediction near the ceiling (18' level). As the time

increases, the profile shape remains similar, but the overprediction of the temperature level becomes greater. The maximum

over-prediction of temperature is observed at t=5 minutes, and is

of the order of 306C. Several factors contribute to the inaccuracy. Their relative importance is discussed on the following

pages.

1. Grid Coarseness: With only one grid cell within the

flame region (burner diameter equal to the cell size)

the maximum flame temperature will be smeared within

the cell. The real flame front. temperatures may be as

high as 12000C,, whereas cell averaged combustion zone

temperatures are of the order of 3700C (Figure 5.6).

48
~- 4- - - .~ - - - .. a, * -

* ~ * I ~ - - - ~. -p ~ %

~ , , -

.- / ~p ~

- -p *~ - - %' S..

1~ S., %~

- ~ / - - 'S \

- " / ~ * a '~

- ~ I - - -

-0

'S

S -

- - -p

_ S.

S.'.!'

- .- *!~

A•

,,-,,, "- - • S.

N- ,, 'I

: :.-.'

i S *

I 1J j 

, - .S -- 4. . ~ 5' . .. a. 

a , * 

• t- .--- • Y *5 %t p * - /, / I -. "'

-, S, 9. -

I/ ' I- ,-• 

'• . - l" I i' q l -• " '- l -iVmin z 0.011 m/e

Vmax = 1.772 m/e

Time = 1.0 Min

Vmin w 0.011 m/c

Vmax c 2.470 m/a

Time = 2.0 Min

'. # "''Zki- J!/l s.7. Vrin Vmax c = 0.042 2.268 rn/u m/u

;4 Time - 41.0 Min

-Ia , ,,.,, :: ,. ,,,

• ~~~, " \

- ,,,, -f -. , . -

"" t.. " - - -'v°.,. ,

ILI

"5'%1 'IVain = 0.035 a/u * "Vmax S = 3.183 r/u - 00 Time a=690 Min

,,. .... ,- I . *1 . / I, , -

Figure 5.8 Velocity Contours at the First Vertical y-z

Plane Near Ventilation Exit

49
I

4.. 

a

W 4~~q- ~ ~- 4 ~ ~ p ~- I - __

*% ~ ~9% 

- ~ s.-.

I S - ~

5 '-~ %* \ \ Vrniriz 0.9 rn/a

- //,

-

I;,;I

I -

I ,

/

/

~5/

-- p

-b ~

i=.i

r-N

Vmax 0.463 r/s

0 , Time a 1.0 min

I ' s n ho i

- , , ... --.. -, ,., -,

/-- - 4. - -

..l t, • ..,.• -- .- 9...- ,. -,.• . 4. -. ., .. • •

"- ¶

N1 V a in = 0 . 00 M/ 8

7 V m a x = 0.7 9 4 M/ 2

1i ur 5. Veoity et s at ?Time = 2.0 Min

N a . Vl ( - 6) Ti me 

tma --. 9 /

Nea Vetlto Inlets - (y- 6)

50
I - 4-~ 4--- 46

a- I

A - - ' - a-. S ~

/ , -- % - -

, / - _ - .~ .~ - =

It 1, ~= - - S

/ ,, S ~ \ ~**%.

, I . * - -~

- - - I - -

Vmin 0 0.017 m/a

Vmax = 0.504 M/u

Time = 1.0 Min

Vmin - 0.046 m/6

Vmax x 1.079 m/u

Time = 2.0 Min

Vmin = 0.046 m/u

Vmax 2.483 m/r

Time z 4.0 Min

- ~- -

/4.- uS%. % -. 4..

/4.~

I I

I ' ~ ~ t - PD I -- , I

-.--. ,

b ~ -.

.5

.5

\

%

S a-

-

- - F - -

-

*~ ~ 4-

if- 4- 4~~

"'1

* . /

t I { C~ ~' IVminn .0453 rn/z.

I I • ''/ dVmax i 3.673 m1/

Time = 6.0 Min

Figure 5.10 Velocity Vectors at the Last Vertical

y-z Plane Near the Wall

51
20.

0. 10. 20. 30. 40. 50. 50. 70. 0. 10. 20. 30. 40. 50. 60. 70.

T("C) T('C)

20.

15.

10.

H(ft)

5.

-. /p

I

I

I

I

I

/

U

g

I

a

/

*

I

I

I

I

I

-.---. I

/

I

'I'

I

I

I

S

I

J

4.0 min 5.0 min

0. 10.

m m m m

20. 30. 40. 50. 60. 70. 0. 10. 20. 30. 40. 50. 50. 70.

T("C) T("C)

Figure 5.11 Comparison of Predicted Vertical Temperature

Profiles with Experimental Data

52
The second inaccuracy will result from numerical

diffusion in the region of high property gradients

(flame zone). As a result, high temperatures will be

predicted in the vicinity of the flame as well. This

explains the largest over-prediction of temperatures at

the six foot level, shown in Figure 5.11.

2. Combustion Model: As discussed earlier (section 2.*3),

the one-step reaction model will result in higher

overall energy release and the release will take place

mainly in the burner vicinity. This explains the

higher temperature levels at the six foot elevation,

and higher overall temperature levels. This inaccuracy

can be significantly reduced by employing a multireaction model.

3. Wall Heat Transfer: Thermal radiation is the key

mechanism of quick energy transfer between the high

temperature flame and the solid walls (the flame-gas

exchange is less important). The irradiated walls are

then transferring the energy by convection to the

boundary layer gases. As a result, significant energy

is lost from the flame zone to all walls of the domain.

Not accounting for radiation, in the flame vicinity,

will contribute to higher predicted temperatures at the

six foot level. In addition, the heat transfer between

the irradiated walls and the nearby gas will contribute

to a large vertical temperature gradient. This

phenomena can be observed from experimental data. Lack

of wall heat transfer (by convection and radiation) in

our present demonstration calculations resulted in

smoother temperature profiles, visible from Figure

5.11.

5.4 Investigatory 2-D Computations

To investigate the above-described contributing factors, and to

further examine some of the basic code capabilities, a series of

two-dimensional test cases were considered. The two-dimensional

setup assumes only one plane in the x-direction while it ma intains the same geometry in the vertical y-z plane. The boundary

conditions used were the same as those used in the three-dimensional case. The inlet ports were placed at the same vertical

and horizontal locations as in the three-dimensional problem.

The ventilation rates were reduced by a factor of 6.1 proportional to the size of the room (viz: x-distance = 6.1 meters in 3-D;

x-distance = 1.0 meter in 2-D) and were 0.0125 kg/sec. The fire

source was maintained at the same locations, and the heat release

rate was reduced by the same factor as the ventilation rates.

The reason for choos ing two-dimensional instead of three-dimensional cases was the shorter computer time and fast turn-around

53
of two-dimensional computations. The cases selected include

effects of:

* grid size;

* fire in a room with obstacles;

* sensitivity to ventilation inlet, velocity specif ications, and size of exhaust port; and

* higher ventilation rates.

The following sections describe the cases used for two-dimensional analysis and present the obtained results.

Grid Size Sensitivity Study: In this study, two cases were

considered. In the f irst case, the same grid as in the threedimensional analysis, in the yz plane, was used. In the second

case, the grid was ref ined to 35 x 19 grid cells, as shown in

Figure 5.12b. A steady-state plus six minutes transient

calculations were performed for both cases.

Results from the finer grid case are presented in Figures 5.13

and 5.14. Figure 5.13 shows the velocity vectors at steady-state

and at four 'minutes of the transient. At the steady-state, the

air exits from the inlet ports and flows in all directions with

the overall motion being directed, as expected, toward the exit.

At four minutes into the transient, the ventilation effects are

overwhelmed by the buoyancy dominated flow, and the fire source

is apparent. The phenomena of recirculation is observed as in

the three-dimensional case. Figure .5.14 presents temperature

contours at one and four minutes into the transient. Both plots

show the hot layer near the ceiling. At four minutes, the hot

layer descends toward the bottom of the room.

Results of the coarse grid calculation show the same phenomena.

Figure 5.15 presents a comparison between coarse grid and finer

grid results. For this comparison,, the vertical temperature

profiles, at the same location that was used for comparison with

experimental results in the three-dimensional analysis, are

plotted at 1, 2, 3, 4, 5, and 6 minutes for both cases. it is

observed that the hot layer, near the ceiling, is more pronounced

in the fine grid case. Also, as expected, the temperature at the

six foot level is lower when the f iner grid is used. Both of

those differences make the finer grid results relatively more

agreeable with experimental data. Considerable further refinement is essential. Since -the f low is truly three-dimensional,

further grid-sensitivity studies should be three-dimensional, and

should be conducted after the incorporation of improved physical

models, in Phase II study.

Fire in Room with Obstacles: In this case, the same operational

conditions were maintained along with the finer grid of the

previous case. However, an obstacle was placed above the fire,

as shown in Figure 5.16 by the shaded area.

54
Y• Z (a) 20 x 11 Grid Size

y

t-z

(b) 35 x 19 Grid Size

Figure 5.12 Computational Grids Used for Grid Sensitivity

Study; (2-D Investigatory Problem)

55
ZZ

a I P a

- - - - - ..- ~..e -r - *

- - - - - ~ ~ . .

- ~e.rr~,v[~7.. - .

- - - - i I a I S ~ I S

a * * * * *

- - - - - - a a * a I S S

a - - S * * S

C C - C * I * *

Vmin = 0.0 m/z

Vmax a 0.052 m/a

Steady State

I 

S~ 

~~~~~~~~~~ .- .--

~~ Sz% z 

S- 

--- - Z..* -

IN S - - S * S S - ~ S '

I

- - - - -I - aC- -. %

a.. ---- - .~ '-S -S 5

~ I

I. * , * leg I 5 I. * ~ I.t

'a * - Pd *~ I

I - - - - P * P - - - - C- V .* 'a a 

P -

Vd I' C- C- - - V 11d I V V - - d P - - C

*5.1 PIE - * S

- dv, PUP I S

dPI -

- - - - a a * . .

Vain = .0121 m/s

Vmax = 1.578 m/s

Time = 4.8 Min

Figure 5.13 Velocity Vectors for Fine Grid

,7.35cC Time - 1.o Min

Time - 4.o Min

t. * 132.4 •C

Figure 5.14 Temperature Contours for Fine Grid

56
L. 20 x 11 Grid size

------ 35 x 19 Grid size

20.

18.

Height (ft)

16.

14.

12.

10.

ui 8.0

6.0

4.0

2.0

0.0

0. 5. 10. 15. 20. 25. 30. 35. 40. 45. 50. 55. 60. 65. 70. 75. 80. 85. 90. 95. 100.

Temperature (*C)

Figure 5.15 Vertical Temperature Profiles Near the Fire Source
I I L

Figure 5.16 Fine Grid with Obstacles for Mock-Up

2-D Test Case

A steady-state and three-minute transient calculations were

performed. Figure 5.17 shows the velocity vectors for the

steady-state and at 1, 2, and 3 minutes of the transient. The

velocity vectors in the steady-state are similar to those shown

in the previous case, with a stronger flow above the blockage.

After the fire starts, the flow pattern of the air changes

dramatically, where the ventilation inlet port influence is not

observed any more, even at one minute time. However, two

recirculation regions are formed, one below the blockage and the

other above the blockage, and extend up to the exit. The cold

air is entrained into the fire, then accelerated from under the

blockage. Some of this hot air escapes toward the ceiling,

because of its lower density, creating a big recirculation motion

in the room.

Figure 5.18 shows the temperature contours at the 1, 2, and 3

minute times of the transient. At one minute, the temperature

contours resemble a flame in the middle of the room. However, at

two and three minutes, the flame develops and a hot layer is

formed near the ceiling.

Effects of a Wider Exhaust Outlet: The results of the abovedescribed case showed a large recirculation zone under the exit.

To examine whether or not this recirculation region is influenced

by the size of the exhaust port, it was decided to model the same

case with a wider outlet (exit area almost doubled). Figure 5.19

shows the velocity vectors for this case. At the steady-state,

the same overall pattern of the flow is observed and the effect

of the wider outlet is apparent by four velocity vectors which

are directed outward near the exit. Also, the maximum velocity

is almost halved when the exit is widened. At one and two

minutes, the effect of the size of the exit diminished as the

flow pattern resembles that of the narrower exit (Figure 5.17).

As before, the wide recirculation zone persists.

58
* - S - 0---------a *• . - - -b -e -e -• •• • - - - .

- - - - - - -- a. = a - a. a. - a. a..-a. ~ - . a. - a. ~ % I . S

a. - r P a .

m*.S - . -e- *

* ~ a.

~ - - I' I ~ * * a. --- - - - - , * *a.a.a.r a. *. - - - - . . S *

* - .

____________________________________________ .1

- . - - - - -

• .• . . . . °. a ° . • - - - a.

• - ° • • • -- - - . -

. . .. .... . *.. a iiV " : "'--

• ° " ° " a"". * S " " S :~i - a- - - - -

Vain x 0.0 n/u

Vnax = 0.052 lne

Steady State

Vmin = 0.0 m/s

Vmax = 0.501 M/s

Time = 1.0 Min

Vmin = 0.0 m/6

Vmax = 0.716 n/6

Time = 2.0 Min

'V

I'

*

a

a

a

a

I

I

IS

IS

S.

- -

4-

I'

d

I

I

I

I

S

S

S

I

I

.5

Sb

I 1 --. ON - - . a. ra - --- -- 1.

z z- -- a .- a .a .a .a .a .

- - - -- - - - -0.0

- - - - - - - - - - -

-- a,.. '' . . . . . j .. . . .. . " "" "

y. .. - .. - .. . . . . . . . . . . . a._a .d a .a. -• Ose

I **S S

IS S

IS - S

%S - ~s. - S * - -

IS - - -

= -- ~~--.'--*--~'-' Vmin a 6.0 ma/

- -.... ---- -14. pmVax e.... . Y v..', = 0.905 ,a/R

-------.... . ' " z- '" Time 3.0. Min

. . .. . . . . .. . . .

Figure 5.17 Predicted Velocity Vectors for Fine Grid

with Obstacles

59
Time 1.0 Hin

Time " 2.0 Hin

20.351 u CTime O a 3.0 e

Figure 5.18 Temperature Contours for Flow with Obstacles

60
. . . ..: . ... - . . . . . . . . . . . . ... 

* 9. S t.... -, . --.- . • • .- - ° a - . - . a

9. .~- .-. . . . . . . . - - - -- - -

. . . . - . . . . . . . . . . . , -. *.~% %

. . ~ . . .. ~ . . . . .

Vain = 6.0 a/c

Vmax = 0.022 m/a

Steady State

Vain = 0.0 a/e

Vmax = 0.497 m/a

Time = 1.0 Min

Vain = 0.0 a/c

Vmax = 0.717 a/c

Time = 2.0 Min

--- a

IS *~ I llI,

*1

'A

'I

II'

i.A

I. A

- -S

- I S S - - - --- --- - -- ~ - I - - - - - C - - - - - - - - - - -. - - * . I S S * A S S * * * I S ~ ge - - - - C - p - C - - - Cr - - St -- ~~5 5-5*~~~

~. - n -. - - - - - - ~D ~P

A A A A A - - -

* * - . - - . -

* * * - * * - - uP?

* *......-.- - - - - I I

* C - I I ~ S % '~ ~ I

* - S - - A - * I d * £ - - - - - S 1 1 r~ -, * - - A ~.- CFigure 5.19 Predicted Velocity Distribution for a

Wider Ventilation Outlet

61
Ef fects of Simulatingi Ventilation Air with Finite Momentum: IIn

the above cases, the ventilation air at the inlet ports was

simulated as a mass source, and zero momentum source, i.e. it

distributed equally in all directions. In the present case, the

air was forced to flow vertically downward with finite momentum

(no four-way diffuser effect). A steady-state plus three

minutes of the transient calculation were made. Figure 5.20

shows the velocity vectors for the steady-state, one minute and

two minutes time into the transient. At the steady-state, the

only difference in flow pattern,, from that in Figure 5.17, is

observed near the inlet ports where strong or downward flow

motion is seen. After the fire starts, at one and three minutes,

the effects of the ventilation are insignificant and the overall

pattern resembles that of Figure 5.17.

Effects of Increased Ventilation Flow Rate: To observe the

effect of the increase in the ventilation rate, another twodimensional test case was set up using ten times higher ventilation rate. This increased ventilation rate corresponds to that

used in case 5 of the experiment (9]. The finer grid along with

four-way diffuser effects were used in this case. A steady-state

and three-minute transient simulations were performed. Figure

5.21 shows the velocity vectors for steady-state and at 1, 2, and

3 minutes of the transient. The overall pattern of the flow is

similar to that of the lower inlet rates shown earlier in Figure

5.17. However, the ventilation inlets are still observed at one

and two minutes into the transient. Figure 5.22 shows the

temperature contours at 1, 2, and 3 minutes of the transient. At

two minutes after the fire starts, the level of the contours is,

as expected, somewhat lower when the ventilation rate is higher.

6. CONCLUJSIONS AND RECOMMENDATIONS

6.1 Achievements of Phase I Study

A fire modeling computer code has been developed by using welltested solution techniques. As a result, within the short (six

month) period of the present study, it has been possible to

develop,, verify and demonstrate basic capabilities of the code.

The results presented in Sections 4 and 5 clearly establish the

feasibility of consolidating relevant techniques into a fire

code, which can be widely used by various organizations concerned

with the efficient and safe operation of nuclear power plants.

6.2 Recommendations for Phase II Work

The work completed in Phase I has clearly established the

feasibility of an advanced state-of-the-art fire simulation code.

To convert this code into a reliable and usable code (i.e. fully

validated and well-documented code), the following work items are

recommended for SEIR Phase II study.

62
4

------- - - - ------ - - -

~ '.5. S.. s~* ~ - *~. = - S. 0~*~~~

I _____ I Y~ ~ -. - - -~

\'.~u..,, - -- -~ I - ~ -. - 4 I - - - k ~ I * - - --- - - - r - - --- - -

. - - - - -- .- ~ e~ - - --- ~ - - - * - .. S. a. --- -- - - - - - *~ - -. - - - - - - - - ~-C - - E S - -

. - - - - - - - - - - - - - - S

* . a - . a

- - - -

- - - =

- S.. 5. J Sp.S -5, 5.

I ~

O

Vmin= 0.0 rn/a

Vuax c 0.021 n/s

Steady State

Vnin = 0.0 m//

Vnax a 0.528 n/a

Time z 1.0 Min

Vmin 6 0.0 n/6

Vmax a 0.907 M/u

Time = 3.0 MIN

'S..-

~I I

Ii

I.

'I

'I

Pb

I 'S

I*

'S.

S. *

S..

'~~~ . .. . .... "i"" ' '0"~ 'o 

- - - - - 0• - - - - - W p - - .

................ .v....( I f. "

5. * a - 5 t 05r 

Figure 5.20 Velocity Distribution for Ventilations with

No Diffusers

63
= - - - - - - ~ - - = a ~ ~ ~ - - - - ~ - - - % % I

- - d ~ ~ ' I - . * ~ *l9E~ I ~ ~ - - ' - .- ~ * c~ - -

I - - ~ - -

-~ - - .- ---~ ~.--4- .......-.-z~..1 I ** - - - .-. ~\ .% F ~Is ~ R ~ - r* - --- - - 4 I . * * * , 5 ~ ~ ~ *~ - - - - - - - a I . . * * , - - - - - . S S * * * * 4

--- 4- 4- 4- -

. - - - * , * *

*4~~ - - - -- - - - - C C - 4 - S - S 4 4

F. . . F . S

Vain =z0.0 m//

Vmax : 0.529 m/a

Steady State

Vain a 0.0 m/s

Vmax = 0.542 m/z

Time = 1.0 Min

rr'-.- ---- - -- - -4-4- -

Fr - * - *~ S - r - - - - = - = - - - - ~ ~ ~ - . - - - -

£ * - - - - C C - U C

C4444 ,~ ' S S - - * . - - 4

*1% S. *.--------------------S - *5 5 * S - 0 0 - C C C 4

* . S

- -- -- - -... -- 4- *0~ -- 4- ~ 4~ 9- -

V - - - j 4-.-4-4-~~~ I - £ j % I I C -

- C C I S S~ *~ ~ III*~

* 0 4 F ~ % S S % % % ~ C - 4 5 % '5 *5 ~

- ~ .. * 4 . F - - -

0

@

B

G

41

0

0

0

0

0

e

0

Q

0

I

Q

¢

I

I .. - . .......... , rr . ... . . - S ., . )

,,,,... ............ I ...........

.. - IIs -.......... - - - - - - -

*...............................- ,- •.. .

.... ............. , • . . . .-

a................... •, --.-. • -. .,-.3, , &

Vain a 0.0 a/s

Vmax 0.716 m/u

'Time = 2.0 Min

Vain = 0.0 m/S

Vaax: 0.907 m/u

Time = 3.0 kin

Figure 5.21 Calculated Velocity Vectors for Higher

Ventilation Inlet Flow Rates

64
Time = 1.0 Min

Time z 2.0 Min

40.15 C

Time x 3.0 Min

Figure 5.22 Temperature Contours for Higher Ventilation

Inlet Flow Rates

65
1. Completion and testing of the physical models (described in Section 1.3).

2. Systematic parametric studies to, establish the

sensitivity of results to computational mesh, and

various physical models.

3. Validation of the code against several (at least six)

test cases of the 1986-87 NRC/SNL/FMRC experiments. If

necessary, additional data from earlier tests (published by SNL, NBS, EPRI, FRS (UK), etc.) should also be

used.

4. Finalization of the code, with emphasis on the ease-ofuse, modularity, and portability.

5. Documentation of the code (including all assumptions,

theoretical models, code structure and user's guide).

A detailed work plan with time and cost estimates will be

submitted to NRC in the SBIR Phase II proposal.

6.3 Recommendations for Phase III Work

After the satisfactory completion of Phase II work, to effectively utilize and disseminate the developed code, the following

three activities are recommended to be carried out with commercial and/or non-SBIR government funds.

1. Incorporation of selected data bank, i.e. physical

properties and models for typical fire sources pertinent to nuclear power plants.

2. Enhancement of graphical display of results (improved

color graphics and animation of results) so that the

time evolution of the fire environments can be

presented (and recorded) in the form of movies and/or

video cassettes. Such movie (real-time) presentations

are very helpful in general comprehension of complex

flow developments, and can be very useful in assessing

various safety aspects such as evacuation time.

3. Development of technical training and user support

material, and a provision for code maintenance and

necessary enhancements.

For items 1 and 2 above, necessary support will be sought from

commercial organizations such as: Factory Mutual, Ove Arup &

Partners, and IBM. In the past, these organizations have shown

positive inclinations, however, due to the uncertainty and delay

in NRC's request for 1988 SBIR proposals, no formal agreements

have been established.

For item 3 above, CFDRC plans to provide the necessary support

from its internal funds.

66
7. REFERENCES

1. Cline, D.D., Von Riesemann, W.A., and Chavez, J.M., "Investigation of Twenty-Foot Separation Distance as a Fire

Protection Method as Specified in 10CFR50, Appendix R,"

NUREG/CR-3192, SAND83-0306, October 1983.

2. Singhal, A.K. and Tam, L.T., "Three-Dimensional Transient

Fire Environment Predictions," reported in the PRC/NASA

Report; "Solid Rocket Motor/Vehicle Assembly Building

Inadvertent Ignition Effects Study," NASA KSC-DF-441,

September 1981.

3. Boccio, J.L., Usher, J.L., Singhal, A.K., and Tam, L.T.,

"The Use of a Field Model to Analyze Probable Fire Environments Encountered within the Complex Geometries of Nuclear

Power Plants," paper presented at the 23rd National Heat

Transfer Conference, Denver, Colorado, August 1985,

published in Fire and Combustion Systems, HTD vol. 45, pp.

159-166.

4. Usher, J.L., Boccio, J.L., Singhal, A.K., and Tam, L.T.,

"Fire Environment Determination in the LaSalle NPP Control

Room," paper presented at the International ANS/ENS Topical

Meeting on Operability of Nuclear Power Systems in Normal

and Adverse Environments in Albuquerque, New Mexico,

September 29 - October 3, 1986.

5. Markatos, N.C., Malin, M.R., and Cox, G., "Mathematical

Modeling of Buoyancy-Induced Smoke Flow in Enclosures," Int.

J. Heat Mass Transfer, vol. 25, no. 1, pp. 63-75, 1982.

6. Markatos, N.C. and Cox, G., "Turbulent Buoyant Heat Transfer

in Enclosures Containing a Fire Source," Seventh International Heat Transfer Conference, Munchen, 6, p. 373, 1982.

7. Yang, K.T., Kelleher, M., Yang, H.Q., and Raycraft, Jt,

"Field Modeling of Fire and Smoke Spread in Confined

Spaces,"; Technical Report, Department of Aerospace and

Mechanical Engineering, Unversity of Notre Dame, 1987, paper

under preparation.

8. Boccio, J.L. and Usher, J.L., "The Use of a Field Model to

Assess Fire Behavior in Complex Nuclear Power Plant Enclosures: Present Capabilities and Future Prospects,"

NUREG/CR-4479, BNL-NUREG-51948, December 1985.

9. Nowlen, S.P., "Enclosure Environment Characterization

Testing for the Base Line Validation of Computer Fire

Simulation Codes," prepared by Sandia National Laboratories

for the U.S. Nuclear Regulatory Commission, NUREG/CR-4681,

SAND86-1296, 1987.

10. Launder, B.E. and Spalding, D.B., Mathematical Models of

Turbulence, Academic Press, London, 1972.

11. Kim, S.W. and Chen, C.P. "A Multiple-Time-Scale Turbulence

Model Based on Variable Partitioning of Turbulent Kinetic

Energy Spectrum," AIAA-88-1771, 1988.

12. Spalding, D.B., "Concentraton Fluctuations in a Round

Turbulent Free Jet," Chem. Eng. Sci., vol. 26, pp. 95-107,

1971.

67
13. Magnussen, B.F. and Hjertager, B.M., "On Mathematical

Modeling of Turbulent Combustion with Special Emphasis on

Soot Formation and Combustion," 16th Symposium (International) on Combustion, The Combustion Institute, pp. 719-729,

1976.

14. Hautman, D.J., Dryer, F.L., Shug, K.P., and Glassman, I.,

"Multistep Overall Kinetic Mechanism for the Oxidation of

Hydrocarbons," Comb. Science and Tech., vol. 25, nos. 5, 6,

1981.

15. Chandrasekhar, S., "Radiative Transfer," Dover Publications,

New York, 1960.

16. Lockwood, P.C. and Shah, N.G., "A New Radiation Solution

Method for Incorporation in General Combustion Prediction

Procedures," Eighteenth Symposium (International) on

Combustion, The Combustion Institute, Pittsburgh, Pennsylvania, 1981.

17. Lockwood, F.C. and Spalding, D.B., "Predictions of a

Turbulent Reacting Duct Flow with Significant Radiation,"

Thermodynamic Session, Proc. Colloques d'Evian de la Societe

Francaise de Physique, May 1981.

18. Patankar, S.V. and Spalding, D.B., "A Calculation Procedure

for Heat, Mass and Momentum Transfer in Three-Dimensional

Parabolic Flows," International Journal of Heat and Mass

Transfer, vol. 15, pp. 1787-1806, Pergamon Press, 1972.

19. Patankar, S.V., "Numerical Heat Transfer and Fluid Flow,"

McGraw Hill, 1980.

20. Van Doormaal, J.P. and Raithby, G.D., "Enhancement of the

SIMPLE Method for Predicting Incompressible Fluid Flows,"

Numerical Heat Transfer, vol. 7, pp. 147-163, 1984.

21. Przekwas, A.J., "Whole Field Solution Method Elliptic

Difference Equations with General Boundary Conditions,"

under preparation for submission to Numerical Heat Transfer.

22. Stone, H.L., "Iterative Solution of Implicit Approximation

of Multidimensional Partial Differential Equations," SIAM J.

Num. Annual., vol. 5, 1968.

23. Ghia, U., Ghia, K.N., and Shin, C.T., "High-Re Solutions for

Incompressible Flow Using the Navier-Stokes Equations and a

Multi-Grid Method," J. Comput. Phys., vol. 38, pp. 387-411,

1982.

24. Armaly, B.F., Durst, F., Pereira, J.C.F., and Sch6nung, B.,

"Experimental and Theoretical Investigation of BackwardFacing Step Flow," J. Fluid Mech., vol. 127, pp. 473-496,

1983.

25. Laufer, J., "The Structure of Turbulence in Fully Developed

Pipe Flow," NACA TN-2954, Washington, June 1983.

26. Samimy, M., Nejad, A.S., Langenfeld, C.A., Craig, R.R., and

Vanka, S.V., "Isothermal Swirling Flow in a Dump Combustor,"

AIAA-87-1352, 1987.

27. de Wahl Davis, G., "Natural Convection in a Square CavityA Benchmark Solution," Int. J. Numerical Methods in Fluids,

vol. 3, pp. 249-264, 1983.

28. Schefer, R.W. and Dibble, R.W., "Rayleigh Scattering

68
Measurements of Mixture Fraction in a Turbulent Non-Reacting

Propane Jet," AIAA-86-0278, 1986.

69
NRC FORM as U.L NUCLEARI RjULATORY COMMIgSION I REPORT NUMIERA JieA•pdb TrOC.oWVd .1' No. N ,V 

0•441

NRiCM 1102.

.,= BIBLIOGRAPHIC DATA SHEET NUREG/CR-5233

SEE INSTRUCTIONS ON THE REVERSE

2. TITLE AND SUBTITLE 3 LEAVE BLANK

A Computer Code for Fire Protection and Risk Analysis

of Nuclear Plants 4. DATE REPORT COMPLETED

MONTH YEAR

S. A-THORISI June 1988

S DATE REPORT ISSUED

Ashok K. Singhal, Sami D. Habchi, and MONTH I YEAR

Andrzej J. Przekwas September 1988 7. PERFORMING ORGANIZATION NAME AND MAILING ADDRESS Oft-cw Je C1•*) I PROJECT/TASKAWORK UNIT NUMBER

CFD Research Corporation

3313 Bob Wallace Avenue, Suite 205 ,. FIN OR GRANT NUMBER

Huntsville, Alabama 35805 D204 1

10. SPONSORING ORGANIZATION NAME AND MAILING ADDRESS IHusciu& Z. Cea) I Is TYPE OF REPORT

Division of Engineering Final Report for SBIR

Office of Nuclear Regulatory Research Phase I Study

U.S. Nuclear Regulatory Commission b. PERIOD COVERED .neu.t.w sbad

Washington, DC 20555 Oct. '87 - June '88

12 SUPPLEMENTARY NOTES

13. ABSTRACT 12W wydMs o"r w

A fire modeling computer code has been developed with

sufficient flexibility for accurate representations of geometry,

ventilation and other conditions as may be present in cable

rooms, control rooms, and other enclosures in nuclear power

plants. The computer code is capable of three-dimensional,

transient, turbulent flow and heat transfer calculations with

chemical reaction and radiation. The code has a modular structure, specifically designed for fire problems.

The code employs the latest relevant finite-volume solution

techniques. The code has been applied to a series of benchmark

problems and a recent fire test problem. This has confirmed the

feasibility of the fire code. Considerable further work is

needed to enhance the physical models (to improve the realism of

predicted solutions) and to validate and document the final code.

Specific recommendations are made for Phases II and III of the

Project.

14 DOCUMENT ANALYSIS - i. KEYWORDSIDESCRIPTORS

Fire Safety

Field Model

Fire Modeling

Conputer Simulation

Fire Protection

15 AVAILABILITY

STATEMENT

Unlimited

16 SECUIFTY CLASSIFICATION

ITI,.t pWg)

Unclassifie

I Ine-A~ aq if i Pd 17. NUMBER OF PAGES

1t PRICE

b. IDE'ETIFIERSIOPEN.ENDIO TERMS
UNITED STATES

NUCLEAR REGULATORY COMMISSION

WASHINGTON, D.C. 20555

OFFICIAL BUSINESS

PENALTY FOR PRIVATE USE, $300

120555110976 1 1AN1RP

US NRC-NRR

DIV-ENG & SYST TECH

BRANCH CHIEF

PLANT SYSTEMS BRANCH 4W/F 8-D-1

WASHINGTON DC 20555

SPECIALFOURTH-CLASS RATE IPOSTAGE & FEES PAIDI

1 USNRC 1

PEMI N.G-6