This document has been superseded 

by a later version. For the latest 

version go to the web site: 

http://fire.nist.gov/fds
NIST Special Publication 1018

Fire Dynamics Simulator (Version 4)

Technical Reference Guide

Kevin McGrattan, Editor

NIST Special Publication 1018

Fire Dynamics Simulator (Version 4)

Technical Reference Guide

Kevin McGrattan, Editor

Fire Research Division

Building and Fire Research Laboratory

in cooperation with

VTT Building and Transport, Finland

March 2006

UNITED STATES OF AMERICA

DEPARTMENT OF COMMERCE

U.S. Department of Commerce

Donald L. Evans, Secretary

Technology Administration

Phillip J. Bond, Under Secretary for Technology

National Institute of Standards and Technology

Arden L. Bement, Jr., Director
Certain commercial entities, equipment, or materials may be identified in this

document in order to describe an experimental procedure or concept adequately. Such

identification is not intended to imply recommendation or endorsement by the

National Institute of Standards and Technology, nor is it intended to imply that the

entities, materials, or equipment are necessarily the best available for the purpose.

National Institute of Standards and Technology Special Publication 1018

Natl. Inst. Stand. Technol. Spec. Publ. 1018, 94 pages (March 2006)

CODEN: NSPUE2

U.S. GOVERNMENT PRINTING OFFICE

WASHINGTON: 2004

For sale by the Superintendent of Documents, U.S. Government Printing Office

Internet: bookstore.gpo.gov – Phone: (202) 512-1800 – Fax: (202) 512-2250

Mail: Stop SSOP, Washington, DC 20402-0001
Preface

The use of fire models currently extends beyond the fire research laboratories and into the engineering,

fire service and legal communities. Surveys [1, 2] of available fire models show a significant increase in

number over the last decade. Sufficient evaluation of any model is necessary to ensure that users can judge

the adequacy of its technical basis, appropriateness of its use, and confidence level of its predictions. This

document provides the theoretical basis for the Fire Dynamics Simulator (FDS) and a summary of the work

performed to evaluate the model.

This guide is based in part on the “Standard Guide for Evaluating the Predictive Capability of Deterministic Fire Models,” ASTM E 1355 [3]. ASTM E 1355 defines model evaluation as “the process of

quantifying the accuracy of chosen results from a model when applied for a specific use.” The model evaluation process consists of two main components: verification and validation. Verification is a process to check

the correctness of the solution of the governing equations. Verification does not imply that the governing

equations are appropriate; only that the equations are being solved correctly. Validation is a process to determine the appropriateness of the governing equations as a mathematical model of the physical phenomena

of interest. Typically, validation involves comparing model results with experimental measurement. Differences that cannot be explained in terms of numerical errors in the model or uncertainty in the measurements

are attributed to the assumptions and simplifications of the physical model.

Evaluation is critical to establishing both the acceptable uses and limitations of a model. Throughout

its development, FDS has undergone various forms of evaluation, both at NIST and beyond. This guide

provides a survey of work conducted to date to evaluate FDS. Roughly half of the referenced studies were

aimed primarily at model evaluation, the other half describe limited work to validate FDS for a specific use.

The latter group were performed mostly by practicing engineers who did not have the time or resources to

comprehensively evaluate the model. Collectively, the body of work forms the basis of a model evaluation.

As FDS development continues, the work performed to date will provide a framework for future research.

i
ii
Disclaimer

The US Department of Commerce makes no warranty, expressed or implied, to users of the Fire Dynamics

Simulator (FDS), and accepts no responsibility for its use. Users of FDS assume sole responsibility under

Federal law for determining the appropriateness of its use in any particular application; for any conclusions

drawn from the results of its use; and for any actions taken or not taken as a result of analysis performed

using these tools.

Users are warned that FDS is intended for use only by those competent in the fields of fluid dynamics,

thermodynamics, heat transfer, combustion, and fire science, and is intended only to supplement the informed judgment of the qualified user. The software package is a computer model that may or may not have

predictive capability when applied to a specific set of factual circumstances. Lack of accurate predictions

by the model could lead to erroneous conclusions with regard to fire safety. All results should be evaluated

by an informed user.

Throughout this document, the mention of computer hardware or commercial software does not constitute endorsement by NIST, nor does it indicate that the products are necessarily those best suited for the

intended purpose.

iii
iv
Acknowledgments

The development and maintenance of the Fire Dynamics Simulator has been made possible through a partnership of public and private organizations, both in the United States and abroad. Following is a list of

contributors from the various sectors of the fire research, fire protection engineering and fire services communities:

The Fire Dynamics Simulator has been under development for almost 25 years. At NIST, Howard Baum

and Ronald Rehm laid the theoretical groundwork for the model and devised the basic numerical solvers that

are still used in the code today. Kevin McGrattan added features needed to make the model accessible to

practicing fire protection engineers, and he remains the custodian of the FDS source code. Glenn Forney

developed the visualization tool Smokeview that not only made the public release possible, but it also serves

as the principal diagnostic tool for the continuing development of FDS. Kuldeep Prasad added the multipleblock feature of the model, paving the way for parallel processing. William (Ruddy) Mell has developed

unique applications of the model in areas of microgravity combustion and wildland fire spread. Charles

Bouldin devised the basic framework of the parallel version of the code.

Several postdoctoral fellows at NIST have made significant contributions to the model development.

Jason Floyd implemented the mixture fraction combustion model, Francine Battaglia studied the extension

of the model to look at fire whirls, Javier Trelles extended the large eddy simulation technique to large, windblown fire plumes, Amy Musser considered indoor air quality applications along with her advisor Steve

Emmerich of the Building Environment Division of BFRL, Dike Ezekoye worked on an early combustion

model, Kevin Cassel used early versions of the model to study gravity currents/ceiling jets.

In the labs at NIST, Anthony Hamins has conducted a number of small and large scale experiments to

validate FDS, including a substantial set of experiments in response to the World Trade Center fires and

collapse. Alex Maranghides oversees the Large Fire Laboratory at NIST where these tests were conducted,

and has helped to design the experiments. Dick Gann has overseen the fire part of the WTC effort. Other

participants in the WTC effort are Rodney Bryant, Michelle Donnelly, Eric Johnsson, George Mulholland

and Tom Ohlemiller. Technicians Laurean Delauter, Jay McElroy, Jack Lee, Mike Smith and Bill Twilley

have conducted numerous validation experiments in both large and small test facilities at NIST.

Takashi Kashiwagi, Greg Linteris, and John Yang have worked on small-scale experiments that have

been used to validate FDS. Tony Putorti, through graduate work with Prof. Arvind Atreya at the University

of Michigan, has conducted experiments and analysis for sprinkler applications. William Grosshandler and

Tom Cleary helped to develop the smoke detector activation algorithm that has been implemented in FDS.

Grosshandler is also the developer of RadCal, a library of subroutines that have been incorporated in FDS

to provide the radiative properties of gases and smoke.

Also at NIST, Dan Madrzykowski, Doug Walton, Bob Vettori, Dave Stroup, Steve Kerber, and Nelson

Bryner have applied FDS to the reconstruction of several real fire incidents that led to the death of fire

fighters. Dave Evans and Doug Walton oversaw the development of ALOFT, A Large Outdoor Fire plume

Trajectory model, an early off-shoot of FDS.

At VTT Building and Transport, Finland, Simo Hostikka developed and implemented the radiation

solver and the pyrolysis model for charring materials. Also at VTT, Jukka Hietaniemi, Jukka Vaari and

v
Timo Korhonen have performed validation studies of various sub-models. The VTT Fire Research group

continues to work on model development and validation of the model for various applications.

Support for validation experiments and the preparation of the FDS manuals has been provided by the

Office of Nuclear Regulatory Research of the US Nuclear Regulatory Commission (US NRC). Special

thanks to Mark Salley and Monideep Dey for their efforts. The Office of Nuclear Material Safety and

Safeguards, another branch of the NRC, has supported modeling studies of tunnel fires under the direction

of Chris Bajwa and Allen Hansen.

Within the fire protection engineering community, many individuals have made contributions to the

model, either through development of sub-models or through validation work. Wei Zhang, Vijay d’Souza,

Jason Sutula, Steve Olenick, Doug Carpenter, and Richard Roby of Combustion Science and Engineering

(CSE) implemented a smoke detector activation model within FDS and also validated various parts of the

model. Doug Carpenter and Erin Mack Ashley of CSE conducted the literature review that is now part of

the FDS Technical Reference Guide.

Jason Floyd and Javier Trelles, former NIST post-docs now working at Hughes Associates, continue to

make vital contributions to the model because of their familiarity with the source code.

Chris Wood of ArupFire, Dave Sheppard of the US Bureau of Alcohol, Tobacco and Firearms (ATF),

and Doug Carpenter of Combustion Science and Engineering developed a training course on the use of FDS

and Smokeview on behalf of the Society of Fire Protection Engineers (SFPE). Morgan Hurley of the SFPE

supported this effort. Prof. David McGill of Seneca College, Ontario, Canada has also conducted a remotelearning course on the use of FDS. Prof. Ian Thomas of Victoria University has presented short courses on

the use of FDS in Australia. His students have also performed some validation work on compartment fires.

At the University of Maryland, Professors Fred Mowrer, Jim Quintiere, Jose Torero and Marino di Marzo

have supervised students who have conducted validation work with FDS. In particular, Phil Friday and T. Ma

have performed extensive validation work on the basic hydrodynamic solvers, and Paulo Ruffino has studied

the activation of sprinklers that have been wetted by previously activated sprinklers. Profs. Arnaud Trouve

and Andre Marshall are presently studying the application of FDS to under-ventilated fires.

At Worcester Polytechnic Institute (WPI), Prof. Johnathan Barnett has supervised several students, in

particular Jay Ierardi and Chris Lautenberger, who have experimented with more detailed combustion models for small-scale applications.

At Northwestern University, Prof. Richard Lueptow has advised doctoral students Dave Sheppard and

Jay Schwille whose research involved measurements of water sprays from sprinklers and mist nozzles.

At the University of Illinois, Chicago, Prof. Ishwar Puri along with students Anurag Jhalani and Achintya

Mukhopadhyay validated the model against slot burner measurements. The purpose of their study was to

quantify the effects of flame stretch.

At Virginia Tech, Prof. Uri Vandsburger hosted NIST Post-Doc Jason Floyd during his two-year tenure

at NIST. Doctoral student Chris Wieczorek performed a series of compartment fire experiments used in part

to validate FDS.

At Johns Hopkins University, Baltimore, Maryland, Jason Averill, a member of the NIST staff, has done

graduate work with advisors Profs. Charles Meneveau and Omar Knio. The goal of the study was to evaluate

the potential improvements to the large eddy simulation technique presently implemented in FDS.

At the University of Canterbury, Christchurch, New Zealand, Prof. Charles Fleischmann and his students

have provided an extensive set of material properties of upholstered furnishings, plus performed validation

studies. Doctoral student Jason Clement spent a semester at NIST validating the hydrodynamics of FDS

against salt water data.

At Kingston University in Britain, Prof. Jennifer Wen and students have investigated the Smagorinsky

model within FDS and experimented with various other approaches to model turbulence.

At Factory Mutual Research Corporation, John de Ris has contributed his expertise in developing an

experimental soot growth model in FDS, along with Chris Lautenberger at WPI.

vi
At Tyco Fire Products, Dave LeBlanc and James Golinveaux have contributed information and test data

for various sprinkler designs.

Grantees of the Building and Fire Research Laboratory at NIST who have contributed to FDS development include Prof. Pat Pagni of the University of California, Berkeley, Profs. G. Faeth and Arvind Atreya at

the University of Michigan, and Prof. Jay Gore at Purdue University.

Several guest researchers at NIST have made contributions in various areas, among them Prof. Yuji

Nakamura, Nagoya University, Japan; Ales Jug, University of Iceland; Prof Bakhtier Farouk, Drexel University, Claudia Rexfort, Gerhard-Mercator-University Duisburg, Germany; Prof. Damian Rouson, City

College, New York; Prof. Woe-Chul Park, Pukyong National University, Korea.

At National Fire Protection Association (NFPA), Rick Mulhaupt and Eric Peterson of the Research

Foundation have extended invitations to various researchers to speak at a series of conferences about their

efforts with FDS.

There have been many others who have contributed to the development of FDS, many of whom are

listed in the References of this document.

Kevin McGrattan, September 2005

vii
viii
Contents

Preface i

Disclaimer iii

Acknowledgments v

1 Introduction 1

2 Model and Scenario Definition 5

2.1 Model Documentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.1.1 Name and Version of the Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.1.2 Type of Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.1.3 Model Developers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.1.4 Model Uses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.1.5 Model Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.1.6 Relevant Publications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.1.7 Governing Equations, Assumptions and Numerics . . . . . . . . . . . . . . . . . . 8

2.1.8 Limitations of the Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.1.9 Input Data Required to Run the Model . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.1.10 Property Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.2 Scenarios for Which FDS has been Evaluated . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.2.1 Description of Scenarios or Phenomenon of Interest . . . . . . . . . . . . . . . . . 12

2.2.2 List of Quantities Predicted by FDS upon which Evaluation is Based . . . . . . . . . 13

2.2.3 Degree of Accuracy Required for Each Output Quantity . . . . . . . . . . . . . . . 13

3 Theoretical Basis for the Model 15

3.1 Hydrodynamic Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

3.1.1 The Fundamental Conservation Equations . . . . . . . . . . . . . . . . . . . . . . . 15

3.1.2 Large Eddy Simulation (LES) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3.1.3 Direct Numerical Simulation (DNS) . . . . . . . . . . . . . . . . . . . . . . . . . . 18

3.1.4 The Low Mach Number Assumption . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.1.5 The Momentum Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

3.2 Combustion Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

3.2.1 Mixture Fraction Combustion Model . . . . . . . . . . . . . . . . . . . . . . . . . 22

3.2.2 Enhancements to the Mixture Fraction Model . . . . . . . . . . . . . . . . . . . . . 24

3.2.3 Finite-Rate Reaction (DNS) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

3.3 Thermal Radiation Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

3.4 Thermal Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

ix
3.4.1 Convective Heat Transfer to Walls . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

3.4.2 Thermoplastic Fuels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

3.4.3 Liquid Fuels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

3.4.4 Charring Fuels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

3.5 Sprinklers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

3.5.1 Sprinkler Activation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

3.5.2 Sprinkler Droplet Size Distribution . . . . . . . . . . . . . . . . . . . . . . . . . . 33

3.5.3 Sprinkler Droplet Trajectory in Air . . . . . . . . . . . . . . . . . . . . . . . . . . 35

3.5.4 Sprinkler Droplet Transport on a Surface . . . . . . . . . . . . . . . . . . . . . . . 35

3.5.5 Mass and Energy Transfer from Droplets . . . . . . . . . . . . . . . . . . . . . . . 35

3.5.6 Interaction of Droplets and Radiation . . . . . . . . . . . . . . . . . . . . . . . . . 36

3.5.7 Fire Suppression by Water . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

3.6 Smoke Detectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

3.7 Numerical Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

3.7.1 Simplified Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

3.7.2 Temporal Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

3.7.3 Spatial Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

3.7.4 Large Eddy vs. Direct Numerical Simulation . . . . . . . . . . . . . . . . . . . . . 43

3.7.5 The Mass Transport Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

3.7.6 Discretizing the Divergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

3.7.7 Thermal Radiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

3.7.8 Interaction of Droplets and Radiation . . . . . . . . . . . . . . . . . . . . . . . . . 50

3.7.9 Thermal and Material Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . 51

3.7.10 The Momentum Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

3.7.11 The Pressure Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

3.7.12 Particle Tracking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

3.8 Review of the Theoretical Development of the Model . . . . . . . . . . . . . . . . . . . . . 58

3.8.1 Assessment of the Completeness of Documentation . . . . . . . . . . . . . . . . . . 58

3.8.2 Assessment of Justification of Approaches and Assumptions . . . . . . . . . . . . . 59

3.8.3 Assessment of Constants and Default Values . . . . . . . . . . . . . . . . . . . . . 60

4 Mathematical and Numerical Robustness 61

4.1 Comparison with Analytical Solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

4.2 Code Checking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

4.3 Numerical Tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

5 Model Sensitivity 65

5.1 Grid Sensitivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

5.2 Sensitivity of Large Eddy Simulation Parameters . . . . . . . . . . . . . . . . . . . . . . . 67

5.3 Sensitivity of Radiation Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

5.4 Sensitivity of Thermophysical Properties of Solid Fuels . . . . . . . . . . . . . . . . . . . . 69

5.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

6 Model Validation 71

6.1 Validation Work with Pre-Release Versions of FDS . . . . . . . . . . . . . . . . . . . . . . 71

6.2 Validation of FDS since 2000 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

6.2.1 Comparison with Full-Scale Tests Conducted Specifically for the Chosen Evaluation 73

6.2.2 Comparison with Engineering Correlations . . . . . . . . . . . . . . . . . . . . . . 74

x
6.2.3 Comparisons with Previously Published Full-Scale Test Data . . . . . . . . . . . . . 74

6.2.4 Comparison with Standard Tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

6.2.5 Comparison with Documented Fire Experience . . . . . . . . . . . . . . . . . . . . 81

7 Conclusion 85

Nomenclature 87

References 99

xi
xii
Chapter 1

Introduction

The idea that the dynamics of a fire might be studied numerically dates back to the beginning of the computer age. Indeed, the fundamental conservation equations governing fluid dynamics, heat transfer, and

combustion were first written down over a century ago. Despite this, practical mathematical models of fire

(as distinct from controlled combustion) are relatively recent due to the inherent complexity of the problem.

Indeed, in his brief history of the early days of fire research, Hoyt Hottel noted “A case can be made for fire

being, next to the life processes, the most complex of phenomena to understand” [4].

The difficulties revolve about three issues: First, there are an enormous number of possible fire scenarios

to consider due to their accidental nature. Second, the physical insight and computing power required to

perform all the necessary calculations for most fire scenarios are limited. Any fundamentally based study

of fires must consider at least some aspects of bluff body aerodynamics, multi-phase flow, turbulent mixing

and combustion, radiative transport, and conjugate heat transfer; all of which are active research areas in

their own right. Finally, the “fuel” in most fires was never intended as such. Thus, the mathematical models

and the data needed to characterize the degradation of the condensed phase materials that supply the fuel

may not be available. Indeed, the mathematical modeling of the physical and chemical transformations of

real materials as they burn is still in its infancy.

In order to make progress, the questions that are asked have to be greatly simplified. To begin with,

instead of seeking a methodology that can be applied to all fire problems, we begin by looking at a few

scenarios that seem to be most amenable to analysis. Hopefully, the methods developed to study these “simple” problems can be generalized over time so that more complex scenarios can be analyzed. Second, we

must learn to live with idealized descriptions of fires and approximate solutions to our idealized equations.

Finally, the methods should be capable of systematic improvement. As our physical insight and computing

power grow more powerful, the methods of analysis can grow with them.

To date, three distinct approaches to the simulation of fires have emerged. Each of these treats the

fire as an inherently three dimensional process evolving in time. The first to reach maturity, the “zone”

models, describe compartment fires. Each compartment is divided into two spatially homogeneous volumes,

a hot upper layer and a cooler lower layer. Mass and energy balances are enforced for each layer, with

additional models describing other physical processes appended as differential or algebraic equations as

appropriate. Examples of such phenomena include fire plumes, flows through doors, windows and other

vents, radiative and convective heat transfer, and solid fuel pyrolysis. Descriptions of the physical and

mathematical assumptions behind the zone modeling concept are given in separate papers by Jones [5] and

Quintiere [6], who chronicle developments through 1983. Model development since then has progressed to

the point where documented and supported software implementing these models are widely available [7].

The relative physical and computational simplicity of the zone models has led to their widespread use in

the analysis of fire scenarios. So long as detailed spatial distributions of physical properties are not required,

1
and the two layer description reasonably approximates reality, these models are quite reliable. However, by

their very nature, there is no way to systematically improve them. The rapid growth of computing power

and the corresponding maturing of computational fluid dynamics (CFD), has led to the development of CFD

based “field” models applied to fire research problems. Virtually all this work is based on the conceptual

framework provided by the Reynolds-averaged form of the Navier-Stokes equations (RANS), in particular

the k−ε turbulence model pioneered by Patankar and Spalding [8]. The use of CFD models has allowed the

description of fires in complex geometries, and the incorporation of a wide variety of physical phenomena.

However, these models have a fundamental limitation for fire applications – the averaging procedure at the

root of the model equations.

RANS models were developed as a time-averaged approximation to the conservation equations of fluid

dynamics. While the precise nature of the averaging time is not specified, it is clearly long enough to require

the introduction of large eddy transport coefficients to describe the unresolved fluxes of mass, momentum

and energy. This is the root cause of the smoothed appearance of the results of even the most highly resolved

fire simulations. The smallest resolvable length scales are determined by the product of the local velocity

and the averaging time rather than the spatial resolution of the underlying computational grid. This property

of RANS models is typically exploited in numerical computations by using implicit numerical techniques

to take large time steps.

Unfortunately, the evolution of large eddy structures characteristic of most fire plumes is lost with

such an approach, as is the prediction of local transient events. It is sometimes argued that the averaging

process used to define the equations is an “ensemble average” over many replicates of the same experiment

or postulated scenario. However, this is a moot point in fire research since neither experiments nor real

scenarios are replicated in the sense required by that interpretation of the equations. The application of

“Large Eddy Simulation” (LES) techniques to fire is aimed at extracting greater temporal and spatial fidelity

from simulations of fire performed on the more finely meshed grids allowed by ever faster computers.

The phrase LES refers to the description of turbulent mixing of the gaseous fuel and combustion products with the local atmosphere surrounding the fire. This process, which determines the burning rate in

most fires and controls the spread of smoke and hot gases, is extremely difficult to predict accurately. This

is true not only in fire research but in almost all phenomena involving turbulent fluid motion. The basic

idea behind the LES technique is that the eddies that account for most of the mixing are large enough to be

calculated with reasonable accuracy from the equations of fluid dynamics. The hope (which must ultimately

be justified by comparison to experiments) is that small-scale eddy motion can either be crudely accounted

for or ignored.

The equations describing the transport of mass, momentum, and energy by the fire-induced flows must

be simplified so that they can be efficiently solved for the fire scenarios of interest. The general equations of

fluid dynamics describe a rich variety of physical processes, many of which have nothing to do with fires.

Retaining this generality would lead to an enormously complex computational task that would shed very

little additional insight on fire dynamics. The simplified equations, developed by Rehm and Baum [9], have

been widely adopted by the larger combustion research community, where they are referred to as the “low

Mach number” combustion equations. They describe the low speed motion of a gas driven by chemical

heat release and buoyancy forces. Oran and Boris provide a useful discussion of the technique as applied

to various reactive flow regimes in the chapter entitled “Coupled Continuity Equations for Fast and Slow

Flows” in Ref. [10]. They comment that “There is generally a heavy price for being able to use a single

algorithm for both fast and slow flows, a price that translates into many computer operations per time step

often spent in solving multiple and complicated matrix operations.”

The low Mach number equations are solved numerically by dividing the physical space where the fire

is to be simulated into a large number of rectangular cells. Within each cell the gas velocity, temperature,

etc., are assumed to be uniform; changing only with time. The accuracy with which the fire dynamics can

be simulated depends on the number of cells that can be incorporated into the simulation. This number is

2
ultimately limited by the computing power available. Present day, single processor desktop computers limit

the number of such cells to at most a few million. This means that the ratio of largest to smallest eddy

length scales that can be resolved by the computation (the “dynamic range” of the simulation) is on the

order of 100. Parallel processing can be used to extend this range to some extent, but the range of length

scales that need to be accounted for if all relevant fire processes are to be simulated is roughly 104to 105

because combustion processes take place at length scales of 1 mm or less, while the length scales associated

with building fires are of the order of tens of meters. The form of the numerical equations discussed below

depends on which end of the spectrum one wants to capture directly, and which end is to be ignored or

approximated.

3
4
Chapter 2

Model and Scenario Definition

2.1 Model Documentation

This section provides a short description of the Fire Dynamics Simulator following the framework suggested

by ASTM E 1355 [3]. It is intended to outline the major features of the model, its history, the underlying

physical assumptions, and so on. More detailed information about the algorithm itself can be found in the

next chapter.

2.1.1 Name and Version of the Model

The name of the program is the NIST Fire Dynamics Simulator or FDS. FDS is a Fortran 90 computer

program that solves the governing equations of fluid dynamics, and Smokeview is a companion program

written in C/OpenGL programming language that produces images and animations of the results. Version 1

of FDS was publicly released in February 2000, version 2 in December 2001, and version 3 in November

2002. The present version of FDS is 4, released in July 2004. Changes in the version number correspond

to major changes in the physical model or input parameters. For minor changes and bug fixes, incremental

versions are released, referenced according to fractions of the integer version number. For example, version

3.1 would be an updated release of version 3.0. In addition, day-to-day bug fixes made in response to user

feedback are referenced by a compilation date that is printed at the top of the diagnostic output file.

2.1.2 Type of Model

FDS is a Computational Fluid Dynamics (CFD) model of fire-driven fluid flow. The model solves numerically a form of the Navier-Stokes equations appropriate for low-speed, thermally-driven flow with an

emphasis on smoke and heat transport from fires. The partial derivatives of the conservation equations of

mass, momentum and energy are approximated as finite differences, and the solution is updated in time

on a three-dimensional, rectilinear grid. Thermal radiation is computed using a finite volume technique on

the same grid as the flow solver. Lagrangian particles are used to simulate smoke movement and sprinkler

discharge.

2.1.3 Model Developers

The Fire Dynamics Simulator (FDS) was developed and is currently maintained by the Fire Research Division in the Building and Fire Research Laboratory (BFRL) at the National Institute of Standards and

Technology (NIST). A substantial contribution to the development of the model was made by VTT Building

and Transport in Finland. Additional contributors are cited in the Acknowledgments.

5
2.1.4 Model Uses

Throughout its development, FDS has been aimed at solving practical fire problems in fire protection engineering, while at the same time providing a tool to study fundamental fire dynamics and combustion. FDS

can be used to model the following phenomena:

• Low speed transport of heat and combustion products from fire

• Radiative and convective heat transfer between the gas and solid surfaces

• Pyrolysis

• Flame spread and fire growth

• Sprinkler, heat detector, and smoke detector activation

• Sprinkler sprays and suppression by water

Although FDS was designed specifically for fire simulations, it can be used for other low-speed fluid flow

simulations that do not necessarily include fire or thermal effects. To date, about half of the applications

of the model have been for design of smoke control systems and sprinkler/detector activation studies. The

other half consist of residential and industrial fire reconstructions.

2.1.5 Model Results

FDS computes the temperature, density, pressure, velocity and chemical composition within each numerical

grid cell at each discrete time step. There are typically hundreds of thousands to several million grid cells

and thousands to hundreds of thousands of time steps. In addition, FDS computes at solid surfaces the

temperature, heat flux, mass loss rate, and various other quantities. The user must carefully select what data

to save, much like one would do in designing an actual experiment. Even though only a small fraction of the

computed information can be saved, the output typically consists of fairly large data files. Typical output

quantities for the gas phase include:

• Gas temperature

• Gas velocity

• Gas species concentration (water vapor, CO2, CO, N2)

• Smoke concentration and visibility estimates

• Pressure

• Heat release rate per unit volume

• Mixture fraction (or air/fuel ratio)

• Gas density

• Water droplet mass per unit volume

On solid surfaces, FDS predicts additional quantities associated with the energy balance between gas and

solid phase, including

• Surface and interior temperature

6
• Heat flux, both radiative and convective

• Burning rate

• Water droplet mass per unit area

Global quantities recorded by the program include:

• Total Heat Release Rate (HRR)

• Sprinkler and detector activation times

• Mass and energy fluxes through openings or solids

Time histories of various quantities at a single point in space or global quantities like the fire’s heat release

rate (HRR) are saved in simple, comma-delimited text files that can be plotted using a spreadsheet program.

However, most field or surface data are visualized with a program called Smokeview, a tool specifically

designed to analyze data generated by FDS. FDS and Smokeview are used in concert to model and visualize fire phenomena. Smokeview performs this visualization by presenting animated tracer particle flow,

animated contour slices of computed gas variables and animated surface data. Smokeview also presents

contours and vector plots of static data anywhere within a scene at a fixed time.

A complete list of FDS output quantities and formats is given in Ref. [11]. Details on the use of Smokeview are found in Ref. [12].

2.1.6 Relevant Publications

Each version of FDS and Smokeview are documented by three separate publications – the FDS Technical

Reference Guide, the FDS User’s Guide [11], and the Smokeview User’s Guide [12]. The User’s Guides

only describe the mechanics of using the computer programs. The Technical Reference Guide provides the

theory and algorithm details, plus a description of the verification and validation studies.

There are numerous sources that describe various parts of the model. The basic set of equations solved

in FDS were formulated by Rehm and Baum in the Journal of Research of the National Bureau of Standards [9]. The basic hydrodynamic algorithm evolved at NIST through the 1980s and 1990s, incorporating fairly well-known numerical schemes that are documented in books by Anderson, Tannehill and

Pletcher [13], Peyret and Taylor [14], and Ferziger and Peric [15]. This last book provides a good descrip-  ́

tion of the large eddy simulation technique and provides references to many current publications on the

subject. Numerical techniques appropriate for combustion systems are described by Oran and Boris [10].

The mixture fraction combustion model is described in a review article by Bilger [16]. Basic heat transfer theory is provided by Holman [17] and Incropera [18]. Thermal radiation is described in Siegel and

Howell [19].

Much of the current knowledge of fire science and engineering is found in the SFPE Handbook of

Fire Protection Engineering [20]. Popular textbooks in fire protection engineering include those by Drysdale [21] and Quintiere [22]. On-going research in fire and combustion is documented in several periodicals

and conference proceedings. The International Association of Fire Safety Science (IAFSS) organizes a

conference every two years, the proceedings of which are frequently referenced by fire researchers. Interscience Communications, a London-based publisher of several fire-related journals, hosts a conference

known as Interflam roughly every three years in the United Kingdom. The Combustion Institute hosts an

international symposium on combustion every two years, and in addition to the proceedings of this symposium, the organization publishes its own journal Combustion and Flame. The papers appearing in the

IAFSS conference proceedings, the Combustion Symposium proceedings, and Combustion and Flame are

all peer-reviewed, while those appearing in the Interflam proceedings are selected based on the submission

7
of a short abstract. Both the Society for Fire Protection Engineers (SFPE) and the National Fire Protection

Association (NFPA) publish peer-reviewed technical journals entitled the Journal of Fire Protection Engineering and Fire Technology. Other often-cited, peer-reviewed technical journals include the Fire Safety

Journal, Fire and Materials, Combustion Science and Technology, Combustion Theory and Modeling and

the Journal of Heat Transfer.

Research at NIST is documented in various ways beyond contributions made by staff to external journals

and conferences. NIST publishes several forms of internal reports, special publications, and its own journal

called the Journal of Research of NIST. An internal report, referred to as a NISTIR, is a convenient means

to disseminate information, especially when the quantity of data exceeds what could normally be accepted

by a journal. Often a NISTIR is streamlined and published externally, with the NISTIR serving as a more

complete record of the work performed. Previous versions of the FDS Technical Reference Guide and

User’s Guide were published as NISTIRs. The current FDS and Smokeview manuals are being published

as NIST Special Publications, distinguished from NISTIRs by the fact that they are permanently archived.

Work performed by an outside person or organization working under a NIST grant or contract is published

in the form of a NIST Grant/Contract Report (GCR). All work performed by the staff of the Building and

Fire Research Laboratory at NIST beyond 1993 is permanently stored in electronic form and made freely

available via the Internet and yearly-released compact disks (CDs).

2.1.7 Governing Equations, Assumptions and Numerics

Following is a brief description of the major components of FDS. Detailed information regarding the assumptions and governing equations associated with the model is provided in Section 3.1.

Hydrodynamic Model FDS solves numerically a form of the Navier-Stokes equations appropriate for lowspeed, thermally-driven flow with an emphasis on smoke and heat transport from fires. The core

algorithm is an explicit predictor-corrector scheme, second order accurate in space and time. Turbulence is treated by means of the Smagorinsky form of Large Eddy Simulation (LES). It is possible to

perform a Direct Numerical Simulation (DNS) if the underlying numerical grid is fine enough. LES

is the default mode of operation.

Combustion Model For most applications, FDS uses a mixture fraction combustion model. The mixture

fraction is a conserved scalar quantity that is defined as the fraction of gas at a given point in the flow

field that originated as fuel. The model assumes that combustion is mixing-controlled, and that the

reaction of fuel and oxygen is infinitely fast. The mass fractions of all of the major reactants and

products can be derived from the mixture fraction by means of “state relations,” empirical expressions

arrived at by a combination of simplified analysis and measurement.

Radiation Transport Radiative heat transfer is included in the model via the solution of the radiation transport equation for a non-scattering gray gas. In a limited number of cases, a wide band model can be

used in place of the gray gas model. The radiation equation is solved using a technique similar to a

finite volume method for convective transport, thus the name given to it is the Finite Volume Method

(FVM). Using approximately 100 discrete angles, the finite volume solver requires about 15 % of the

total CPU time of a calculation, a modest cost given the complexity of radiation heat transfer. Water

droplets can absorb thermal radiation. This is important in cases involving mist sprinklers, but also

plays a role in all sprinkler cases. The absorption coefficients are based on Mie theory.

Geometry FDS approximates the governing equations on one or more rectilinear grids. The user prescribes

rectangular obstructions that are forced to conform with the underlying grid.

8
Boundary Conditions All solid surfaces are assigned thermal boundary conditions, plus information about

the burning behavior of the material. Usually, material properties are stored in a database and invoked

by name. Heat and mass transfer to and from solid surfaces is usually handled with empirical correlations, although it is possible to compute directly the heat and mass transfer when performing a Direct

Numerical Simulation (DNS).

Sprinklers and Detectors The activation of sprinklers and heat and smoke detectors are modeled using

fairly simple correlations based on thermal inertia in the case of sprinklers and heat detectors, and the

lag in smoke transport through smoke detectors. Sprinkler sprays are modeled by Lagrangian particles

that represent a sampling of the water droplets ejected from the sprinkler.

2.1.8 Limitations of the Model

Although FDS can address most fire scenarios, there are limitations in all of its various algorithms. Some of

the more prominent limitations of the model are listed here. More specific limitations are discussed as part

of the description of the governing equations in Section 3.1.

Low Speed Flow Assumption The use of FDS is limited to low-speed1 flow with an emphasis on smoke

and heat transport from fires. This assumption rules out using the model for any scenario involving

flow speeds approaching the speed of sound, such as explosions, choke flow at nozzles, and detonations.

Rectilinear Geometry The efficiency of FDS is due to the simplicity of its rectilinear numerical grid and

the use of fast, direct solvers for the pressure field. This can be a limitation in some situations where

certain geometric features do not conform to the rectangular grid, although most building components

do. There are techniques in FDS to lessen the effect of “sawtooth” obstructions used to represent nonrectangular objects, but these cannot be expected to produce good results if, for example, the intent

of the calculation is to study boundary layer effects. For most practical large-scale simulations, the

increased grid resolution afforded by the fast pressure solver offsets the approximation of a curved

boundary by small rectangular grid cells.

Fire Growth and Spread Because the model was originally designed to analyze industrial-scale fires, it

can be used reliably when the heat release rate (HRR) of the fire is specified and the transport of

heat and exhaust products is the principal aim of the simulation. In these cases, the model predicts

flow velocities and temperatures to an accuracy within 5 % to 20 % of experimental measurements,

depending on the resolution of the numerical grid 2. However, for fire scenarios where the heat release

rate is predicted rather than prescribed, the uncertainty of the model is higher. There are several

reasons for this: (1) properties of real materials and real fuels are often unknown or difficult to obtain,

(2) the physical processes of combustion, radiation and solid phase heat transfer are more complicated

than their mathematical representations in FDS, (3) the results of calculations are sensitive to both the

numerical and physical parameters. Current research is aimed at improving this situation, but it is

safe to say that modeling fire growth and spread will always require a higher level of user skill and

judgment than that required for modeling the transport of smoke and heat from prescribed fires.

Combustion For most applications, FDS uses a mixture fraction combustion model. The mixture fraction

is a conserved scalar quantity that is defined as the fraction of gas at a given point in the flow field that

1Mach numbers less than about 0.3

2

It is extremely rare to find measurements of local velocities and/or temperatures from fire experiments that have reported error

estimates that are less than 5 %. Thus, the most accurate calculations using FDS do not introduce significantly greater errors in

these quantities than the vast majority of fire experiments.

9
originated as fuel. The model assumes that combustion is mixing-controlled, and that the reaction

of fuel and oxygen is infinitely fast, regardless of the temperature. For large-scale, well-ventilated

fires, this is a good assumption. However, if a fire is in an under-ventilated compartment, or if a

suppression agent like water mist or CO2 is introduced, fuel and oxygen may mix but may not burn.

Also, a shear layer with high strain rate separating the fuel stream from an oxygen supply can prevent

combustion from taking place. The physical mechanisms underlying these phenomena are complex,

and even simplified models still rely on an accurate prediction of the flame temperature and local

strain rate. Sub-grid scale modeling of gas phase suppression and extinction is still an area of active

research in the combustion community. Until reliable models can be developed for building-scale

fire simulations, simple empirical rules can be used that prevent burning from taking place when the

atmosphere immediately surrounding the fire cannot sustain the combustion. Details are found in

Section 3.2.

Radiation Radiative heat transfer is included in the model via the solution of the radiation transport equation for a non-scattering gray gas, and in some limited cases using a wide band model. The equation

is solved using a technique similar to finite volume methods for convective transport, thus the name

given to it is the Finite Volume Method (FVM). There are several limitations of the model. First, the

absorption coefficient for the smoke-laden gas is a complex function of its composition and temperature. Because of the simplified combustion model, the chemical composition of the smokey gases,

especially the soot content, can effect both the absorption and emission of thermal radiation. Second,

the radiation transport is discretized via approximately 100 solid angles. For targets far away from

a localized source of radiation, like a growing fire, the discretization can lead to a non-uniform distribution of the radiant energy. This can be seen in the visualization of surface temperatures, where

“hot spots” show the effect of the finite number of solid angles. The problem can be lessened by the

inclusion of more solid angles, but at a price of longer computing times. In most cases, the radiative

flux to far-field targets is not as important as those in the near-field, where coverage by the default

number of angles is much better.

2.1.9 Input Data Required to Run the Model

All of the input parameters required by FDS to describe a particular scenario are conveyed via one or two

text files created by the user. These files contain information about the numerical grid, ambient environment,

building geometry, material properties, combustion kinetics, and desired output quantities. The numerical

grid is one or more rectilinear meshes with (usually) uniform cells. All geometric features of the scenario

have to conform to this numerical grid. Objects smaller than a single grid cell are either approximated as a

single cell, or rejected. The building geometry is input as a series of rectangular blocks. Boundary conditions

are applied to solid surfaces as rectangular patches. Materials are defined by their thermal conductivity,

specific heat, density, thickness, and burning behavior. There are various ways that this information is

conveyed, depending on the desired level of detail.

A significant part of the FDS input file directs the code to output various quantities in various ways.

Much like in an actual experiment, the user must decide before the calculation begins what information to

save. There is no way to recover information after the calculation is over if it was not requested at the start.

A complete description of the input parameters required by FDS can be found in the FDS User’s

Guide [11].

2.1.10 Property Data

Any simulation of a real fire scenario involves prescribing material properties for the walls, floor, ceiling, and

furnishings. FDS treats all of these objects as homogeneous solids, thus the physical parameters for many

10
real objects can only be viewed as approximations to the actual properties. Describing these materials in the

input file is the single most challenging task for the user. Thermal properties such as thermal conductivity,

specific heat, density, and thickness can be found in various handbooks, or in manufacturers literature, or

from bench-scale measurements. More difficult is the burning behavior at different heat fluxes. Even though

entire books are devoted to the subject [23], it is still difficult to find information on a particular item.

There is a small file listing thermal properties of various common materials that is distributed with the

FDS software. However, this data is meant to serve only as a starting point. For most practical applications,

the data accompanying the program will not be sufficient. Also, the data is not necessarily validated, and

none of it is to be referenced as NIST data. The users of FDS are responsible for verifying the accuracy

or appropriateness of each input parameter used. The FDS manuals describe the governing equations in

which the material properties are used, but ultimately the user must decide if a given material conforms to

the assumptions underlying the physical model.

11
2.2 Scenarios for Which FDS has been Evaluated

This section provides a description of the scenarios or phenomena of interest that have been included in

past efforts to evaluate FDS. Because FDS is used both for research and for practical applications, there

are some routines within the model that have not been comprehensively validated. Even those parts of the

model for which a substantial amount of validation work has been performed will be pushed beyond these

limits by the need to perform increasingly comlex analyses. Consequently, it is the responsibility of the user

to demonstrate the applicability of the model for scenarios that have not yet been validated.

2.2.1 Description of Scenarios or Phenomenon of Interest

FDS is suited for a wide range of thermally-driven fluid flow scenarios, including fire, both in the open (e.g.

unconfined fire plumes) as well as within the built environment. To date, about half of the applications of

FDS have been for design, and half for forensic reconstruction.

Design applications typically involve an existing building or a building under design. A so-called “design fire” is prescribed either by a regulatory authority or by the engineers performing the analysis. Because

the fire’s heat release rate is known, the role of the model is to predict the transport of heat and combustion

products throughout the room or rooms of interest. Ventilation equipment is often included in the simulation, like fans, blowers, exhaust hoods, HVAC ducts, smoke management systems, etc. Sprinkler and heat

and smoke detector activation are also of interest. The effect of the sprinkler spray on the fire is usually

less of interest since the fire is prescribed rather than predicted. Detailed descriptions of the contents of the

building are usually not necessary because these items are not assumed to be burning, and even if they are,

the burning rate will be fixed, not predicted. Sometimes, it is necessary to predict the heat flux from the

fire to a nearby “target,” and even though the target may heat up to some prescribed ignition temperature,

the subsequent spread of the fire usually goes beyond the scope of the analysis because of the uncertainty

inherent in object to object fire spread.

Forensic reconstructions require the model to simulate an actual fire based on information that is collected after the event, such as eye witness accounts, unburned materials, burn signatures, etc. The purpose

of the simulation is to connect a sequence of discrete observations with a continuous description of the fire

dynamics. Usually, reconstructions involve more gas/solid phase interaction because virtually all objects in

a given room are potentially ignitable, especially when flashover occurs. Thus, there is much more emphasis on such phenomena as heat transfer to surfaces, pyrolysis, flame spread, and suppression. In general,

forensic reconstructions are more challenging simulations to perform because they require more detailed

information about the room contents, and there is much greater uncertainty in the total heat release rate as

the fire spreads from object to object.

Validation studies of FDS to date have focussed more on design applications than reconstructions. The

reason is that design applications usually involve prescribed fires and demand a minimum of thermophysical

properties of real materials. Transport of smoke and heat is the primary focus, and measurements can be

limited to well-placed thermocouples, a few heat flux gauges, gas samplers, etc. Phenomena of importance

in forensic reconstructions, like second item ignition, flame spread, vitiation effects and extinction, are more

difficult to model and more difficult to study with well-controlled experiments. Uncertainties in material

properties and measurements, as well as simplifying assumptions in the model, often force the comparison

between model and measurement to be qualitative at best. Nevertheless, current validation efforts are moving

in the direction of these more difficult issues.

12
2.2.2 List of Quantities Predicted by FDS upon which Evaluation is Based

An FDS simulation predicts the thermal environment of the entire volume of interest by computing thermodynamic quantities in each grid cell of the computational mesh. The number of quantities predicted by the

model is greater than the number of quantities that are typically measured in validation experiments. Following is a list of quantities predicted in each computational cell by the model that are typically compared

to experimental measurement. The order of the items on the list emphasize their relative importance in most

validation studies.

• Gas temperature

• Gas velocity

• Major gas species concentrations (CO2, CO, O2)

• Smoke concentration/obscuration

• Pressure

On solid surfaces, FDS predictions of the following are typically compared with measurements:

• Surface temperature (and sometimes back surface temperature)

• Total heat flux (radiative plus convective)

• Burning (or mass loss) rate

In addition to these point measurements, model and experimental predictions of global quantities, such as

heat release rate and the flux of mass and energy through openings or walls in the compartment, are often

compared. These global quantities can be used to assess the overall energy budget, confirming energy

conservation in the model, and helping to assess systematic uncertainty in the experiments.

2.2.3 Degree of Accuracy Required for Each Output Quantity

The degree of accuracy for each output variable required by the user is highly dependent on the technical

issues associated with the analysis. The user must ask: How accurate does the analysis have to be to answer

the technical question posed? Thus, a generalized definition of the accuracy required for each quantity with

no regard as to the specifics of a particular analysis is not practical and would be limited in its usefulness.

Returning to the earlier definitions of “design” and “reconstruction,” design applications typically are

more accurate because the heat release rate is prescribed rather than predicted, and the initial and boundary

conditions are far better characterized. Mathematically, a design calculation is an example of a “well-posed”

problem in which the solution of the governing equations is advanced in time starting from a known set of

initial conditions and constrained by a known set of boundary conditions. The accuracy of the results is a

function of the fidelity of the numerical solution, which is mainly dependent on the size of the computational

grid. It will be shown in Chapter 6 that validation efforts to date involving well-characterized geometries

and prescribed fires have shown that FDS predictions vary from being within experimental uncertainty to

being about 20 % different than measurements of temperature, heat flux, gas concentration, etc.

A reconstruction is an example of an “ill-posed” problem because the outcome is known whereas the

initial and boundary conditions are not. There is no single, unique solution to the problem, that is, it is

possible to simulate numerous fires that produce the given outcome. There is no right or wrong answer,

but rather a small set of plausible fire scenarios that are consistent with the collected evidence. These

simulations are then used to demonstrate to fire service personnel why the fire behaved as it did based on

13
the current understanding of fire physics incorporated in the model. Most often, the result of the analysis is

only qualitative. If there is any quantification at all, it could be in the time to reach critical events, like a roof

collapse or room flashover.

14
Chapter 3

Theoretical Basis for the Model

This chapter presents the theoretical basis for FDS, followed by a discussion of the model review process.

The physical model is first presented as a set of partial differential equations, with appropriate simplifications

and approximations noted. Then, the finite difference form of the equations are presented. This latter set

constitutes the numerical algorithm.

3.1 Hydrodynamic Model

An approximate form of the Navier-Stokes equations appropriate for low Mach number applications is

used in the model. The approximation involves the filtering out of acoustic waves while allowing for large

variations in temperature and density [9]. This gives the equations an elliptic character, consistent with

low speed, thermal convective processes. The computation can either be treated as a Direct Numerical

Simulation (DNS), in which the dissipative terms are computed directly, or as a Large Eddy Simulation

(LES), in which the large-scale eddies are computed directly and the sub-grid scale dissipative processes

are modeled. The choice of DNS vs. LES depends on the objective of the calculation and the resolution

of the computational grid. If, for example, the problem is to simulate the flow of smoke through a large,

multi-room enclosure, it is not possible to resolve the combustion and transport processes directly. However,

for small-scale combustion experiments, it is possible to compute the transport directly and the combustion

processes to some extent.

3.1.1 The Fundamental Conservation Equations

The conservation equations for mass, momentum and energy for a Newtonian fluid are presented here.

These are the same equations that can be found in almost any textbook on fluid dynamcs or CFD. A particularly useful reference for a description of the equations, the notation used, and the various approximations

employed is Anderson et al. [13].

Conservation of Mass

∂ρ

∂t

+∇· ρu = 0 (3.1)

Conservation of Momentum (Newton’s Second Law)

∂

∂t

(ρu) +∇· ρuu+∇p = ρf+∇· τi j (3.2)

15
Conservation of Energy (First Law of Thermodynamics)

∂

∂t

(ρh) +∇· ρhu =

Dp

Dt

+q ̇

000 −∇·q+Φ (3.3)

Equation of State for a Perfect Gas

p =

ρR T

M

(3.4)

Note that this is a set of partial differential equations consisting of six equations for six unknowns, all

functions of three spatial dimensions and time: the density ρ, the three components of velocity u = (u,v,w),

the temperature T, and the pressure p. The sensible enthalpy h is a function of the temperature: h = R T

T0

cp(T

0

)dT0.

The notation adopted above is deliberately intended to render the equations in as simple and compact

a form as possible. However, when all of the terms are expanded, the equations appear unwieldy. For

the model developers, this is unavoidable, but for the model user, it is sufficient to understand the basic

conservation laws to appreciate what the model does, and does not represent. Following is a brief description

of the terms. Note that boldfaced quantities represent vectors. Boldfaced quantities with the subscripts i j

represent tensors, which can be best understood as 3×3 matrices.

The mass conservation mass equation is often written in terms of the mass fractions of the individual

gaseous species:

∂

∂t

(ρYi) +∇· ρYiu = ∇· ρDi∇Yi +m ̇

000

i

(3.5)

Summing these equations over all species yields the original mass conservation equation. This implies that

the sum of the mass fractions ∑Yi = 1, the sum of the production/loss rates ∑m ̇

000

i = 0, and the sum of the

diffusion terms ∑ρDi∇Yi = 0.

In the momentum equation, the term uu is a tensor formed by multiplying (in the sense of two 3×1

matrices) the vectors u

T

and u. The term ∇ · ρuu is thus a vector formed by applying the vector operator

∇ = ( ∂

∂x

,

∂

∂y

,

∂

∂z

) to the tensor. The force term f in the momentum equation consists of gravity plus other

forces such as the drag exerted by liquid droplets. The stress tensor τi j is defined:

τi j = μ



2 Si j −

2

3

δi j(∇·u)



; δi j =



1 i = j

0 i 6= j

; Si j =

1

2



∂ui

∂xj

+

∂uj

∂xi



i, j = 1,2,3 (3.6)

The term Si j is the strain tensor, written using conventional tensor notation. The symbol μ is the dynamic

viscosity of the fluid.

In the energy equation, note the use of the material derivative, D( )/Dt = ∂( )/∂t +u·∇( ). The term  ̇q

000

is the heat release rate per unit volume from a chemical reaction. The term ∇ · q represents the conductive

and radiative heat fluxes. The term Φ is known as the dissipation function, the rate at which kinetic energy

is transferred to thermal energy due to the viscosity of the fluid:

Φ ≡ τi j ·∇u ≡ μ



2 Si j ·Si j −

2

3

(∇·u)

2



(3.7)

= μ

"

2



∂u

∂x

2

+2



∂v

∂y

2

+2



∂w

∂z

2

+ (3.8)



∂v

∂x

+

∂u

∂y

2

+



∂w

∂y

+

∂v

∂z

2

+



∂u

∂z

+

∂w

∂x

2

−

2

3



∂u

∂x

+

∂v

∂y

+

∂w

∂z

2

#

(3.9)

16
This term is usually neglected because it is very small relative to the heat release rate of the fire. However,

it is included here because it plays a role in the discussion of large eddy simulation.

The equations presented here form the basis for a wide variety of engineering applications, but not

without the further application of simplifying assumptions unique to each field. The only assumptions

made thus far are that the fluid is a perfect gas, that the stress is linearly dependent on the strain, and that

heat is conducted according to Fourier’s Law. Except for extreme applications, like hypersonic flow, these

assumptions are not contested. In the sections to follow, additional assumptions will be imposed on the

governing equations to apply them to fire and other low speed thermal processes. Most important is the

treatment of the diffusive and source terms that differentiates one type of CFD model from another, and

indeed one type of fire model from another.

3.1.2 Large Eddy Simulation (LES)

The most distinguishing feature of any CFD model is its treatment of turbulence. Chapter 1 contains a brief

history of turbulence modeling as it has been applied to the fire problem. Of the three main techniques of

simulating turbulence, FDS contains only Large Eddy Simulation (LES) and Direct Numerical Simulation

(DNS). There is no Reynolds-Averaged Navier-Stokes (RANS) capability in FDS.

LES is a technique used to model the dissipative processes (viscosity, thermal conductivity, material

diffusivity) that occur at length scales smaller than those that are explicitly resolved on the numerical grid.

This means that the parameters μ, k and D in the equations above cannot be used directly in most practical

simulations. They must be replaced by surrogate expressions that “model” their impact on the approximate

form of the governing equations. This section contains a simple explanation of how these terms are modeled

in FDS. Note that this discussion is quite different than what it typically found in the literature, thus the

reader is encouraged to consider other explanations of the technique in the references that are listed in a

review article by Pope [24].

Recall from Section 3.1.1 that there is a small term in the energy equation known as the dissipation

function Φ, the rate at which kinetic energy is converted to thermal energy by viscosity. To understand

where this term originates, form an evolution equation for the kinetic energy of the fluid by taking the dot

product of the momentum equation (3.2) with the velocity vector1:

ρ

Du

Dt

·u ≡ ρ

D



|u|

2/2



Dt

= ρf·u−∇p ·u+∇·(τi j ·u)−Φ (3.10)

As mentioned above Φ is a neglible quantity in the energy equation. However, its functional form is useful

in representing the dissipation of kinetic energy from the resolved flow field. Following the analysis of

Smagorinsky [25], the viscosity μ is modeled

μLES = ρ(Cs ∆)

2



2 S ̄

i j ·S ̄i j −

2

3

(∇·u ̄)

2

1

2

(3.11)

where Csis an empirical constant and ∆ is a length on the order of the size of a grid cell. The bar above

the various quantities denotes that these are the resolved, or filtered, values, meaning that they are computed

on a numerical grid. The other diffusive parameters, the thermal conductivity and material diffusivity, are

related to the turbulent viscosity by

kLES =

μLES cp

Pr

; (ρD)l,LES =

μLES

Sc

(3.12)

The Prandtl number Pr and the Schmidt number Sc are assumed to be constant for a given scenario.

1

In this section it is convenient to work with the Lagrangian form of the conservation equations.

17
The model for the viscosity, μLES, serves two roles: first, it provides a stabilizing effect in the numerical algorithm, damping out numerical instabilities as they arise in the flow field, especially where vorticity

is generated. Second, it has the appropriate mathematical form to describe the dissipation of kinetic energy from the flow. Note the similar mathematical form of μLES and the dissipation function Φ, defined in

Eq. (3.9)2. In the parlance of the turbulence community, the dissipation rate (most often denoted by ε) is

often related to the turbulent kinetic energy (most often denoted by k) by the relation ε ≈ k

3/2/L, where L

is a length scale. Similar scaling arguments show that the dissipation function formed using the modeled

viscosity and the resolvable flow variables, Φ ̄ , has this same three-halves dependence.

There have been numerous refinements of the original Smagorinsky model [26, 27, 28], but it is difficult

to assess the improvements offered by these newer schemes. There are two reasons for this. First, the

structure of the fire plume is so dominated by the large-scale resolvable eddies that even a constant eddy

viscosity gives results almost identical to those obtained using the Smagorinsky model [29]. Second, the lack

of precision in most large-scale fire test data makes it difficult to assess the relative accuracy of each model.

The Smagorinsky model with constant Cs produces satisfactory results for most large-scale applications

where boundary layers are not well-resolved. In fact, experience to date using the simple form of LES

described above has shown that the best results are obtained when the Smagorinsky constant Csis set as low

as possible to maintain numerical stability. In other words, the most realistic flow simulations are obtained

when resolvable eddies are not “damped” by excessive amounts of artificial viscosity. The “best” value of

Cs, a hotly debated topic, has more to do with the differencing scheme of the convective terms than it does

with any physical phenomenon.

3.1.3 Direct Numerical Simulation (DNS)

There are some flow scenarios where it is possible to use the diffusive parameters μ, k and D directly.

Usually, this means that the numerical grid cells are on the order of 1 mm or less, and the simulation is

regarded as a Direct Numerical Simulation (DNS). In FDS, for a DNS, the viscosity, thermal conductivity

and material diffusivity are approximated from kinetic theory because the temperature dependence of each

is important in combustion scenarios. The viscosity of the lth species is given by

μl =

26.69×10−7(Ml T)

1

2

σ

2

l Ωv

kg

m s

(3.13)

where σlis the Lennard-Jones hard-sphere diameter (A ̊) and Ωv is the collision integral, an empirical function of the temperature T. The thermal conductivity of the lth species is given by

kl =

μl cp,l

Pr

W

m K

(3.14)

where the Prandtl number Pr is 0.7. The viscosity and thermal conductivity of a gas mixture are given by

μDNS = ∑

l

Yl μl; kDNS = ∑

l

Yl kl (3.15)

The binary diffusion coefficient of the lth species diffusing into the mth species is given by

Dlm =

2.66×10−7 T

3/2

M

1

2

lm σ

2

lm ΩD

m2

s

(3.16)

2Often in the literature, the Smagorinsky model is written μLES = ρ(Cs ∆)2



S ̄

i j ·S ̄i j 1

2

. Its original formulation in the meteorological literature neglected the divergence as well, but not the factor of 2. This is a minor detail given the arbitrary constant in the

expression, but it is important to maintain a relationship between the dissipation function Φ and the Smagorinky viscosity.

18
where Mlm = 2(1/Ml +1/Mm)

−1

, σlm = (σl +σm)/2, and ΩD is the diffusion collision integral, an empirical

function of the temperature T [30]. It is assumed that nitrogen is the dominant species in any combustion

scenario considered here, thus the diffusion coefficient in the species mass conservation equations is that of

the given species diffusing into nitrogen

(ρD)l,DNS = ρ Dl0 (3.17)

where species 0 is nitrogen.

3.1.4 The Low Mach Number Assumption

Another distinguishing feature of a CFD model is the regime of flow speeds (relative to the speed of sound)

for which it is designed. High speed flow codes involve compressibility effects and shock waves. Low speed

solvers, however, explicitly eliminate compressibity effects that give rise to acoustic (sound) waves.

As written in Section 3.1.1, the Navier-Stokes equations describe the propagation of information at

speeds comparable to that of the fluid flow (for fire, 10-20 m/s), but also at speeds comparable to that of

sound waves (for still air, 300 m/s). Solving a discretized form of the basic N-S equations would require

extremely small time steps in order to account for information traveling at the speed of sound, making

practical simulations impossible.

Following the work of Rehm and Baum [9], an approximation to the Equation of State (3.4) is made

by decomposing the pressure into a “background” component, a hydrostatic component, and a flow-induced

perturbation

p = p0 −ρ∞gz+ p ̃ (3.18)

For most applications, p0 is constant and the other two components are relatively small. Adjustments to this

assumption can be made in the case when the pressure rises due to a fire in a tightly sealed enclosure, or

when the height of the domain is on the order of a kilometer and p0 can no longer be assumed constant and

must be considered a function of the altitude [31].

The purpose of decomposing the pressure is that for low-Mach number flows, it can be assumed that the

temperature and density are inversely proportional, and thus the equation of state can be approximated

p0 = ρTR ∑(Yi/Mi) = ρTR /M (3.19)

The pressure p in the state and energy equations is replaced by the background pressure p0 to filter out

sound waves that travel at speeds that are much faster than typical flow speeds expected in fire applications.

The low Mach number assumption serves two purposes. First, the filtering of acoustic waves means that the

time step in the numerical algorithm is bound only by the flow speed as opposed to the speed of sound, and

second, the modified state equation leads to a reduction in the number of dependent variables in the system

of equations by one. The energy equation (3.3) is never explicitly solved, but its source terms are included

in the expression for the flow divergence, to be derived presently.

The divergence of the flow is obtained by taking the material derivative of the modified equation of

state (3.19), and then substituting terms from the mass and energy conservation equations. First, define

the constant-pressure specific heat of the mixture: cp = ∑l cp,l Yl where cp,lis the temperature-dependent

specific heat of species l. Next, define the enthalpy h = ∑l hl Yl where

hl =

Z T

T0

cp,l(T

0

) dT0(3.20)

Taking the material derivative of Eq. 3.19 and substituting terms from the mass and energy equations, the

divergence of the flow can be expressed in terms of the thermodynamic quantities

∇·u =

1

ρcp T



q ̇

000 −∇·q



+



1

ρcp T

−

1

p0



d p0

dt

19
M

ρ ∑

l

∇· ρDl∇(Yl/Ml)−

1

ρcp T ∑

l

hl ∇· ρDl∇Yl +

1

ρ ∑

l



M

Ml

−

hl

cp T



m ̇

000

l

(3.21)

The pressure rise term on the right hand side of the divergence expression is non-zero only if it assumed that

the enclosure is tightly sealed, in which case the background pressure p0 can no longer be assumed constant

due to the increase (or decrease) in mass and thermal energy within the enclosure. The evolution equation

for the pressure is found by integrating Eq. (3.21) over the entire domain Ω

d p0

dt

=

Z

Ω

1

ρcp T



q ̇

000 −∇·q



dV −

Z

∂Ω

u· dS

,Z

Ω



1

p0

−

1

ρcp T



dV (3.22)

This last equation is rarely needed in practical fire simulations because compartments are rarely sealed, thus

the background pressure is assumed constant.

3.1.5 The Momentum Equation

The momentum equation is simplified to make it easier to solve numerically. First, start with the nonconservative form of the momentum equation introduced above

ρ



∂u

∂t

+ (u·∇)u



+∇p = ρg+f+∇· τi j (3.23)

Next, make the following substitutions:

1. Subtract the hydrostatic pressure gradient, ρ∞g, from both sides. Note that ∇p = ρ∞g+∇p ̃.

2. Apply the vector identity: (u·∇)u = ∇|u|

2/2−u×ω

3. Divide all terms by the density, ρ

4. Decompose the pressure term:

∇p ̃

ρ

=

∇p ̃

ρ∞

+



1

ρ

−

1

ρ∞



∇p ̃

5. Define H = |u|

2/2+ p ̃/ρ∞

Now the momentum equation can be written

∂u

∂t

−u×ω+∇H +



1

ρ

−

1

ρ∞



∇p ̃ =

1

ρ

[(ρ−ρ∞)g+f+∇· τi j] (3.24)

The numerical solution of the pressure equation obtained by taking the divergence of Eq. (3.24) is greatly

simplified by either neglecting the last term on the left hand side, or in cases where it cannot be neglected,

treating it with some care. The decision to either neglect the extra pressure term or to approximate it depends

on its relative contribution to the creation of vorticity. When the momentum equation is integrated around a

closed loop that moves with the fluid, in the absence of any external forces, f, one can readily identify the

sources of vorticity:

dΓ

dt

=

I

1

ρ∞



1−

ρ∞

ρ



∇p ̃· dx+

I

ρ−ρ∞

ρ

g · dx+

I

(∇· τi j)· dx (3.25)

The first term on the right hand side represents the baroclinic torque. The second term is buoyancy-induced

vorticity. The third term represents the vorticity generated by viscosity or sub-grid scale mixing, as in

20
boundary and shear layers. In most practical large scale fire simulations, the fire itself occupies a small part

of the computational domain. Hence, the fire is often not well resolved by the numerical grid, in which case

the vorticity generated in regions where there are large deviations in density is not captured directly. The

mixing of air and combustion products occurs in the plume above the fire where buoyancy is the dominant

source of vorticity. In these calculations, the baroclinic torque can be neglected to simplify the numerical

solution. In simulations where detailed flame dynamics are resolvable, the pressure term responsible for the

baroclinic torque cannot be neglected, but for reasons to be made clear below, must be treated differently

than the other pressure term. In neither of these cases is the Boussinesq approximation invoked. The fluid

is still considered thermally-expandable; the divergence is non-zero; the mass and energy equations are not

modified.

The reason for either neglecting the baroclinic torque or decomposing the pressure term in the momentum equation is to simplify the elliptic partial differential equation obtained by taking the divergence of the

momentum equation

∇

2H = −

∂(∇·u)

∂t

−∇·F ; F = −u×ω+



1

ρ

−

1

ρ∞



∇p ̃ −

1

ρ

((ρ−ρ∞)g+f+∇· τi j) (3.26)

Note that the pressure  ̃p appears on both sides of Eq. (3.26). The pressure on the right hand side is taken from

the previous time step of the overall explicit time-marching scheme. It can be neglected if the baroclinic

torque is not considered important in a given simulation. The pressure on the left hand side (incorporated in

the variable H ) is (always) solved for directly. The reason for the decomposition of the pressure term is so

that the linear algebraic system arising from the discretization of Eq. (3.26) has constant coefficients (i.e. it

is separable) and can be solved to machine accuracy by a fast, direct (i.e. non-iterative) method that utilizes

Fast Fourier Transforms (FFT). No-flux or forced-flow boundary conditions are specified by asserting that

∂H

∂n

= −Fn −

∂un

∂t

(3.27)

where Fn is the normal component of F at the vent or solid wall, and ∂un/∂t is the prescribed rate of change

in the normal component of velocity at a forced vent. Initially, the velocity is zero everywhere.

At open external boundaries the pressure-like term H is prescribed, depending on whether the flow is

outgoing or incoming

H = |u|

2/2 outgoing

H = 0 incoming (3.28)

The outgoing boundary condition assumes that the pressure perturbation  ̃p is zero at an outgoing boundary

and that H is constant along streamlines. The incoming boundary condition assumes that H is zero infinitely

far away. At the boundary between two meshes, the pressure boundary condition is similar to that at an

external open boundary, except that where the flow is incoming, H is taken from the adjacent mesh.

21
3.2 Combustion Model

There are two types of combustion models used in FDS. The choice depends on the resolution of the underlying grid. For a DNS calculation where the diffusion of fuel and oxygen can be modeled directly, a global

one-step, finite-rate chemical reaction is most appropriate. However, in an LES calculation where the grid

is not fine enough to resolve the diffusion of fuel and oxygen, a mixture fraction-based combustion model

is used.

3.2.1 Mixture Fraction Combustion Model

The mixture fraction combustion model is based on the assumption that large-scale convective and radiative

transport phenomena can be simulated directly, but physical processes occurring at small length and time

scales must be represented in an approximate manner. The nature of the approximations employed are

necessarily a function of the spatial and temporal resolution limits of the computation, as well as our current

(often limited) understanding of the phenomena involved.

The actual chemical rate processes that control the combustion energy release are often unknown in fire

scenarios. Even if they were known, the spatial and temporal resolution limits imposed by both present

and foreseeable computer resources places a detailed description of combustion processes beyond reach.

Thus, the model adopted here is based on the assumption that the combustion is mixing-controlled. This

implies that all species of interest can be described in terms of a mixture fraction Z(x,t). The mixture

fraction is a conserved quantity representing the fraction of material at a given point that originated in the

fuel stream. The relations between the mass fraction of each species and the mixture fraction are known

as “state relations”. The state relation for the oxygen mass fraction provides the information needed to

calculate the local oxygen mass consumption rate. The form of the state relation that emerges from classical

laminar diffusion flame theory is a piecewise linear function. This leads to a “flame sheet” model, where

the flame is a two dimensional surface embedded in a three dimensional space. The local heat release rate

is computed from the local oxygen consumption rate at the flame surface, assuming that the heat release

rate is directly proportional to the oxygen consumption rate, independent of the fuel involved. This relation,

originally proposed by Huggett [32], is the basis of oxygen calorimetry.

Start with the most general form of the combustion reaction

νF Fuel+νO O2 → ∑

i

νP,i Products (3.29)

The numbers νi are the stoichiometric coefficients for the overall combustion process that reacts fuel “F”

with oxygen “O” to produce a number of products “P”. The stoichiometric equation (3.29) implies that the

mass consumption rates for fuel and oxidizer are related as follows:

m ̇

000

F

νFMF

=

m ̇

000

O

νOMO

(3.30)

The mixture fraction Z is defined as:

Z =

sYF −(YO −Y

∞

O

)

sYI

F +Y

∞

O

; s =

νOMO

νFMF

(3.31)

By design, it varies from Z = 1 in a region containing only fuel to Z = 0 where the oxygen mass fraction takes

on its un depleted ambient value, Y

∞

O

. Note that Y

I

F

is the fraction of fuel in the fuel stream. The quantities

MF and MO are the fuel and oxygen molecular weights, respectively. The mixture fraction satisfies the

conservation law

ρ

DZ

Dt

= ∇· ρD∇Z (3.32)

22
Figure 3.1: State relations for propane.

obtained from a linear combination of the fuel and oxygen conservation equations. The assumption that the

chemistry is “fast” means that the reactions that consume fuel and oxidizer occur so rapidly that the fuel

and oxidizer cannot co-exist. The requirement that fuel and oxidizer simultaneously vanish defines a flame

surface as:

Z(x,t) = Zf; Zf =

Y

∞

O

sYI

F +Y

∞

O

(3.33)

The assumption that fuel and oxidizer cannot co-exist leads to the “state relation” between the oxygen mass

fraction YO and Z

YO(Z) = 

Y

∞

O

(1−Z/Zf) Z < Zf

0 Z > Zf

(3.34)

State relations for both reactants and products can be derived by considering the following ideal reaction

of a hydrocarbon fuel:

CxHy +η(x+y/4) (O2 +3.76 N2) → max(0,1−η) CxHy+min(1,η) x CO2 +

min(1,η) (y/2) H2O + max(0,η−1) (x+y/4) O2 +η(x+y/4)3.76 N2 (3.35)

Here η is a parameter ranging from 0 (all fuel with no oxygen) to infinity (all oxygen with no fuel). A

correspondence between η and Z is obtained by applying the definition of Z (Eq. 3.31) to the left hand side

of Eq. (3.35). Mass fractions of the products of the infinitely fast reaction (including excess fuel or oxygen)

can be obtained from the right hand side of Eq. (3.35).

An expression for the local heat release rate can be derived from the conservation equations and the state

relation for oxygen. The starting point is Huggett’s [32] relationship for the heat release rate as a function

of the oxygen consumption

q ̇

000 = ∆HO m ̇000

O (3.36)

23
Here, ∆HO is the heat release rate per unit mass of oxygen consumed (about 13,100 kJ/kg for most fuels).

The oxygen mass conservation equation

ρ

DYO

Dt

= ∇· ρD∇YO +m ̇

000

O (3.37)

can be transformed into an expression for the local heat release rate using the conservation equation for the

mixture fraction (3.32) and the state relation for oxygen YO(Z).

−m ̇

000

O = ∇·



ρD

dYO

dZ

∇Z



−

dYO

dZ

∇· ρD∇Z = ρD

d

2YO

dZ2

|∇Z|

2

(3.38)

Neither of these expressions for the local oxygen consumption rate is particularly convenient to apply numerically because of the discontinuity of the derivative of YO(Z) at Z = Zf

. However, an expression for

the oxygen consumption rate per unit area of flame sheet can be derived from Eq. (3.38) by integrating  ̇m

000

O

over a small volume through which the flame sheet cuts. Working with the middle terms of Eq. (3.38) and

noting that dYO/dZ is constant on one side of the flame sheet and zero on the other, the volume integral can

be rewritten as a surface integral over the flame sheet by applying the divergence theorem and seeing that

the two terms cancel at the exterior boundary of the control volume. At this point, it is more convenient to

express the oxygen consumption rate in units of mass per unit time per unit area of flame sheet:

−m ̇

00

O =

dYO

dZ









Z<Zf

ρD ∇Z ·n (3.39)

In the numerical algorithm, the local heat release rate is computed by first locating the flame sheet, then

computing the local heat release rate per unit area, and finally distributing this energy to the grid cells cut

by the flame sheet. In this way, the ideal, infinitely thin flame sheet is smeared out over the width of a grid

cell, consistent with all other gas phase quantities.

3.2.2 Enhancements to the Mixture Fraction Model

The mixture fraction model described in the previous section has several limitations, both numerical and

physical. Its numerical limitations are related to the resolution of the underlying numerical grid. On coarse

grids, the accuracy of the fuel transport and combustion processes is diminished by the high levels of numerical diffusion. The above procedure for determining the local heat release rate works well for calculations

in which the fire is adequately resolved. A measure of how well the fire is resolved is given by the nondimensional expression D

∗/δx, where D∗

is a characteristic fire diameter

D

∗ =



Q ̇

ρ∞ cp T∞

√g

2

5

(3.40)

and δx is the nominal size of a grid cell3. The quantity D

∗/δx can be thought of as the number of computational cells spanning the characteristic (not necessarily the physical) diameter of the fire. The more cells

spanning the fire, the better the resolution of the calculation. For fire scenarios where D

∗

is small relative

to the physical diameter of the fire, and/or the numerical grid is relatively coarse, the stoichiometric surface

Z = Zf will underestimate the observed flame height [33, 34]. It has been found empirically that a good estimate of flame height can be found for crude grids if a different value of Z is used to define the combustion

region

Zf,eff

Zf

= min 1 , C

D

∗

δx



(3.41)

3The characteristic fire diameter is related to the characteristic fire size via the relation Q∗ = (D∗/D)5/2

, where D is the physical

diameter of the fire.

24
Here C is an empirical constant equal to 0.6 to be used for all fire scenarios. As the resolution of the

calculation increases, the Zf,eff approaches the ideal value, Zf, at which point the approximation is no

longer required. The benefit of the expression is that it provides a quantifiable measure of the grid resolution

that takes into account not only the size of the grid cells, but also the size of the fire.

Another consequence of a coarse numerical grid is that a disproportionate amount of the combustion

energy is released near the edges of the burner. From Eq. 3.39, it can be seen that the heat release rate per

unit area of the flame sheet is proportional to the local gradient of the mixture fraction and the local value

of the material diffusivity. The gradient of the mixture fraction is large at the base of the fire because there

a stream of pure fuel meets surrounding air. The diffusivity is large on a coarse grid because it is related to

the Smagorinsky viscosity. To prevent too much of the energy from being released too close to the burner

when a coarse grid is used, there is a maximum bound imposed on the local heat release rate per unit area

of flame sheet. This upper bound is based on a simple analysis in which the fire is assumed to be conical in

shape with surface area, A, and a flame height, H, given by Heskestad’s correlation [35]

H/D = 3.7 Q

∗2/5 −1.02 ; A = πr

p

r

2 +h2

(3.42)

The surface area of a real flame is larger than that of a cone, so the upper bound estimate will prevent too

much energy from being released too close to the fire when a coarse grid is used, but will be high enough

not to interfere with the calculation when the grid is well-resolved. Any energy that is “clipped” off due to

the upper bound is redistributed over the entire flame volume.

The physical limitation of the mixture fraction approach is that it is assumed that fuel and oxygen burn

instantaneously when mixed. For large-scale, well-ventilated fires, this is a good assumption. However, if a

fire is in an under-ventilated compartment, or if a suppression agent like water mist or CO2 is introduced, fuel

and oxygen may mix but may not burn. Also, a shear layer with high strain rate separating the fuel stream

from an oxygen supply can prevent combustion from taking place. The physical mechanisms underlying

these phenomena are complex, and even simplified models still rely on an accurate prediction of the flame

temperature and local strain rate. Sub-grid scale modeling of gas phase suppression and extinction is still an

area of active research in the combustion community. Until reliable models can be developed for buildingscale fire simulations, simple empirical rules can be used that prevent burning from taking place when

the atmosphere immediately surrounding the fire cannot sustain the combustion. Based on the work of

Quintiere [6], Beyler [36] and others, a simple model for flame extinction has been implemented in FDS.

The mixture fraction continues to be used to track the progress of the fuel mixing with the surrounding air,

but now the surrounding volume is assessed to determine if it is more or less likely to support combustion.

Figure 3.2 shows values of temperature and oxygen concentration for which burning can and cannot take

place.

Note that once the gas environment falls in the “No Burn” zone, the state relations (Fig. 3.1) are no

longer valid for values of Z below stoichiometric, since now some fuel may be mixed with the other combustion products. To account for the deviation from the ideal state relations, at least one other scalar quantity

in addition to the mixture fraction would have to be tracked in the calculation. Research is ongoing to add

this second scalar to the model.

3.2.3 Finite-Rate Reaction (DNS)

In a DNS calculation, the diffusion of fuel and oxygen can be modeled directly, thus it is possible to implement a relatively simple one-step chemical reaction. Consider the reaction of oxygen and a hydrocarbon

fuel

νCxHyCxHy +νO2 O2 −→ νCO2CO2 +νH2O H2O (3.43)

25
Figure 3.2: Oxygen-temperature phase space showing where combustion is allowed and not allowed to take place.

The reaction rate is given by the expression

d[CxHy]

dt

= −B[CxHy]

a

[O2]

b

e

−E/RT (3.44)

Suggested values of B, E, a and b for various hydrocarbon fuels are given in Refs. [37, 38]. It should be

understood that the implementation of any of these one-step reaction schemes is still very much a research

exercise because it is not universally accepted that combustion phenomena can be represented by such a simple mechanism. Efforts are currently underway to determine in what cases a one-step reaction mechanism

provides a valid description of the combustion process.

26
3.3 Thermal Radiation Model

The Radiative Transport Equation (RTE) for an absorbing/emitting and scattering medium is

s·∇Iλ(x,s) = −[κ(x,λ) +σs(x,λ)]I(x,s) +B(x,λ) + σs(x,λ)

4π

Z

4π

Φ(s,s

0

) Iλ(x,s

0

) dΩ

0

(3.45)

where Iλ(x,s) is the radiation intensity at wavelength λ, s is the direction vector of the intensity, κ(x,λ) and

σs(x,λ) are the local absorption and scattering coefficients, respectively, and B(x,λ) is the emission source

term. The integral on the right hand side describes the in-scattering from other directions. In the case of a

non-scattering gas the RTE becomes

s·∇Iλ(x,s) = κ(x,λ) [Ib(x)−Iλ(x,s)] (3.46)

where Ib(x) is the source term given by the Planck function. This section describes the radiation transport

in the gas phase.

In practical simulations the spectral dependence can not be solved accurately. Instead, the radiation

spectrum is divided into a relatively small number of bands, and a separate RTE is derived for each band.

The limits of the bands are selected to give an accurate representation of the most important radiation bands

of CO2 and water. The band specific RTE’s are now

s·∇In(x,s) = κn(x) [Ib,n(x)−In(x,s)], n = 1...N (3.47)

where In is the intensity integrated over the band n, and κn is the appropriate mean absorption coefficient

inside the band. The source term can be written as a fraction of the blackbody radiation

Ib,n = Fn(λmin,λmax) σ T

4

/π (3.48)

where σ is the Stefan-Boltzmann constant. The calculation of factors Fn is explained in Ref. [19]. When the

intensities corresponding to the bands are known, the total intensity is calculated by summing over all the

bands

I(x,s) =

N

∑

n=1

In(x,s) (3.49)

From a series of numerical experiments it has been found that six bands are usually enough (N = 6). If the

absorption of the fuel is known to be important, separate bands can be reserved for fuel, and the total number

of bands is increased to ten (N = 10). For simplicity, the fuel is assumed to be CH4. The limits of the bands

are shown in Table 3.1.

Table 3.1: Limits of the spectral bands.

9 Band Model 1 2 3 4 5 6 7 8 9

Soot CO2 CH4 Soot CO2 H2O H2O Soot Soot Major Species

H2O, Soot Soot Soot Soot CH4, Soot

ν (1/cm) 10000 3800 3400 2800 2400 2174 1429 1160 1000 50

λ (μm) 1.00 2.63 2.94 3.57 4.17 4.70 7.00 8.62 10.0 200

6 Band Model 1 2 3 4 5 6

Soot CO2 CH4 CO2 H2O, CH4, Soot Soot Major Species

H2O, Soot Soot Soot

Even with a reasonably small number of bands, the solution of N RTE’s is very time consuming. Fortunately, in most large-scale fire scenarios soot is the most important combustion product controlling the

27
thermal radiation from the fire and hot smoke. As the radiation spectrum of soot is continuous, it is possible

to assume that the gas behaves as a gray medium. The spectral dependence is lumped into one absorption

coefficient (N = 1) and the source term is given by the blackbody radiation intensity

Ib(x) = σT(x)

4

/π (3.50)

In optically thin flames, where the amount of soot is small compared to the amount of CO2 and water, the

gray gas assumption may produce significant overpredictions of the emitted radiation.

For the calculation of the gray or band-mean absorption coefficients, κn, a narrow-band model, RadCal [39], has been implemented in FDS. At the start of a simulation the absorption coefficient(s) are tabulated

as a function of mixture fraction and temperature. During the simulation the local absorption coefficient is

found by table-lookup.

In calculations of limited spatial resolution, the source term, Ib, in the RTE requires special treatment

in the neighborhood of the flame sheet because the temperatures are smeared out over a grid cell and are

thus considerably lower than one would expect in a diffusion flame. Because of its dependence on the

temperature raised to the fourth power, the source term must be modeled in those grid cells cut by the flame

sheet. Elsewhere, there is greater confidence in the computed temperature, and the source term can assume

its ideal value there

κ Ib =



κσT

4/π Outside flame zone

max(χr q ̇

000/4π , κσT4/π) Inside flame zone (3.51)

Here,  ̇q

000 is the chemical heat release rate per unit volume and χr

is an empirical estimate of the local fraction

of that energy emitted as thermal radiation. Near the flame in large scale calculations, neither κ nor T can

be computed reliably, hence the inclusion of the empirical radiation loss term. Note the difference between

the prescription of a local χr and the resulting global equivalent. For a small fire (D < 1 m), the local χris

approximately equal to its global counterpart. However, as the fire increases in size, the global value will

typically decrease due to a net re-absorption of the thermal radiation by the increasing smoke mantle.

The boundary condition for the radiation intensity leaving a gray diffuse wall is given as

Iw(s) = εIbw +

1−ε

π

Z

s

0

·nw<0

Iw(s

0

) |s

0

·nw| dΩ (3.52)

where Iw(s) is the intensity at the wall, ε is the emissivity, and Ibw is the black body intensity at the wall.

The radiative transport equation (3.47) is solved using techniques similar to those for convective transport in finite volume methods for fluid flow [40], thus the name given to it is the Finite Volume Method

(FVM).

To obtain the discretized form of the RTE, the unit sphere is divided into a finite number of solid angles.

In each grid cell a discretized equation is derived by integrating equation (3.47) over the cell i jk and the

control angle δΩl, to obtain

Z

Ωl

Z

Vi jk

s·∇In(x,s)dV dΩ =

Z

Ωl

Z

Vi jk

κn(x) [Ib,n(x)−In(x,s)]dV dΩ (3.53)

The volume integral on the left hand side is replaced by a surface integral over the cell faces using the

divergence theorem. Assuming that the radiation intensity In(x,s) is constant on each of the cell faces, the

surface integral can be approximated by a sum over the cell faces. More detail on the discretization and

solution of the RTE can be found in Section 3.7.7.

The radiant heat flux vector qris defined

qr(x) = Z

4π

sI(x,s) dΩ (3.54)

28
The radiative loss term in the energy equation is

−∇·qr(x) = κ(x) [U(x)−4πIb(x)] ; U(x) = Z

4π

I(x,s)dΩ (3.55)

In words, the net radiant energy gained by a grid cell is the difference between that which is absorbed and

that which is emitted.

29
3.4 Thermal Boundary Conditions

The type of thermal boundary condition applied at any given surface depends on whether that surface is to

heat up and burn, whether the burning rate will simply be prescribed, or whether there is to be any burning

at all.

3.4.1 Convective Heat Transfer to Walls

The heat fluxes to a solid surface consist of gains and losses from convection and radiation. The radiative

flux at the surface is obtained from the boundary condition for the radiation equation, Eq. (3.52).

The calculation of the convective heat flux depends on whether one is performing a Direct Numerical

Simulation (DNS) or a Large Eddy Simulation (LES). In a DNS calculation, the convective heat flux to a

solid surface  ̇q

00

c

is obtained directly from the gas temperature gradient at the boundary

q ̇

00

c = −k

∂T

∂n

(3.56)

where n is the spatial coordinate pointing into the solid. In an LES calculation, the convective heat flux to

the surface is obtained from a combination of natural and forced convection correlations

q ̇

00

c = h ∆T W/m2

; h = max C|∆T|

1

3 ,

k

L

0.037 Re

4

5 Pr

1

3



W/m2/K (3.57)

where ∆T is the difference between the wall and the gas temperature (taken at the center of the grid cell

abutting the wall), C is the coefficient for natural convection (1.52 for a horizontal surface and 1.31 for

a vertical surface) [17], L is a characteristic length related to the size of the physical obstruction, k is the

thermal conductivity of the gas, and the Reynolds Re and Prandtl Pr numbers are based on the gas flowing

past the obstruction. Since the Reynolds number is proportional to the characteristic length, L, the heat

transfer coefficient is weakly related to L. For this reason, L is taken to be 1 m for most calculations.

3.4.2 Thermoplastic Fuels

If it is assumed that the pyrolysis of the solid material occurs at the surface, or if no pyrolysis occurs at all,

the treatment of boundary is as follows.

If the surface material is assumed to be thermally-thick, a one-dimensional heat conduction equation

for the material temperature, Ts(x,t), is applied in the direction x pointing into the solid (the point x = 0

represents the surface)

ρs cs

∂Ts

∂t

=

∂

∂x



ks

∂Ts

∂x



; −ks

∂Ts

∂x

(0,t) = q ̇

00

c +q ̇

00

r −m ̇

00 ∆Hv (3.58)

where ρs, cs and ks are the temperature-dependent density, specific heat and conductivity of the material;  ̇q

00

c

is the convective and  ̇q

00

r

is the (net) radiative heat flux at the surface,  ̇m

00 is the mass loss rate of fuel and ∆Hv

is the heat of vaporization. It is assumed that fuel pyrolysis takes place at the surface, thus the heat required

to vaporize fuel is extracted from the incoming energy flux. The pyrolysis rate is given by an Arrhenius

expression

m ̇

00 = A ρs e−EA/RT (3.59)

R is the universal gas constant. The value of the pre-exponential factor A and the activation energy EA can

be prescribed so that the material burns in the neighborhood of a prescribed temperature. The actual burning

rate is governed by the overall energy balance in the solid. These parameters are often difficult to obtain for

30
Figure 3.3: Fuel pyrolysis rates for two given “ignition” temperatures. See Eq. 3.59.

real fuels; the intent of using the given expression for the mass loss rate is to mimic the behavior of burning

objects when details of their pyrolysis mechanisms are unknown. Figure 3.3 displays a few plots of the mass

flux versus temperature for two different temperatures.

If the surface material is assumed to be thermally-thin, that is, its temperature is assumed uniform across

its width, Ts(t) is affected by gains and losses due to convection, radiation and pyrolysis. The thermal lag of

the material is a function of the product of its density, specific heat and thickness δ

dTs

dt

=

q ̇

00

c +q ̇

00

r −m ̇

00 ∆Hv

ρscsδ

(3.60)

The convective and radiative fluxes are summed over the front and back surfaces of the thin fuel. Unless

otherwise specified, the back surface is assumed to face an ambient temperature void. Note that the individual values of the parameters ρs

, cs and δ are not as important as their product, thus often in the literature

and in the computer program, the three values are lumped together as a product. The pyrolysis rate for a

thermally-thin fuel is the same as for a thermally-thick; see Eq. 3.59.

3.4.3 Liquid Fuels

The rate at which liquid fuel evaporates when burning is a function of the liquid temperature and the concentration of fuel vapor above the pool surface. Equilibrium is reached when the partial pressure of the fuel

vapor above the surface equals the Clausius-Clapeyron pressure

pcc = p0 exp−

hvMf

R



1

Ts

−

1

Tb

 (3.61)

31
where hv is the heat of vaporization, Mf

is the molecular weight, Tsis the surface temperature, and Tb is the

boiling temperature of the fuel [41].

For simplicity, the liquid fuel itself is treated like a thermally-thick solid for the purpose of computing

the heat conduction. There is no computation of the convection of the liquid within the pool.

3.4.4 Charring Fuels

The heat transfer and pyrolysis of charring materials like wood are described using a one-dimensional model,

originally developed by Atreya [42], and further developed by Ritchie et al. [43]. To simplify the original

model, the pyrolysis is assumed to take place over an infinitesimally thin front. The model includes the

conduction of heat inside the material, the evaporation of moisture and the degradation of the virgin material

to gaseous fuel and char. The volatile gases are instantaneously transported to the surface. The governing

equation for energy is

ρc

∂Ts

∂t

=

∂

∂x

ks

∂Ts

∂x

+

∂ρs

∂t

[∆Hpyr −C(T −T0)] + ∂ρm

∂t

[∆Hev −D(T −T0)] (3.62)

where ρsis the total density of the solid and ρm is the moisture density. ∆Hpyr and ∆Hev are the heats of

pyrolysis and water vaporization, respectively. The boundary condition on the surface is

−ks

∂Ts

∂x

(0,t) = q ̇

00

c +q ̇

00

r

(3.63)

Notice that the energy required to pyrolyze the fuel is not extracted from the incoming heat flux, as is done

in the thermoplastic model, but rather the energy is drawn from inside the material at the point where the

pyrolysis is occurring.

The coefficients C and D are defined

C =

ρs0cp,s0 −ρcharcp,char

ρs0 −ρchar

−cp,g ; D = cp,m −cp,g (3.64)

where the subscripts s0 and char refer to the virgin material and char, and the subscripts g and m refer to

the gaseous products of pyrolysis and the moisture. An overbar denotes the average value for the quantity

between the temperature T and ambient temperature T0.

The pyrolysis rate of the material is modeled with a single first order Arrhenius reaction

m ̇

00 = A (ρs0 −ρchar) e−EA/RT (3.65)

The coefficients A and EA are chosen so that the pyrolysis takes place very close to the desired pyrolysis

temperature. The pyrolysis is assumed to take place over a thin front moving inside the material. The

velocity of the front is given by

v =

m ̇

00

ρs0 −ρchar

(3.66)

and the front moves a distance ∆x = v ∆t. The evaporation rate of moisture is modeled in a similar manner.

The following definitions are used to calculate the thermal properties of the material during the drying

and charring processes.

ρc = ρacp,s0 +ρccp,char +ρmcp,m ; ks = ks0

ρa

ρs0

+kchar

ρc

ρchar

(3.67)

ρa = ρs0

ρs −ρchar

ρs0 −ρchar

; ρc = ρs −ρa (3.68)

where the thermal properties of the virgin material and char may be temperature dependent.

The selection of the proper material properties and pyrolysis coefficients is a difficult task. To reduce

the uncertainty related to the coefficients of the pyrolysis rate, a thin pyrolysis front is assumed so that the

burning rate is controlled mainly by the heat of pyrolysis and the heat transfer inside the material.

32
3.5 Sprinklers

Simulating the effects of a sprinkler spray involves a number of pieces: predicting activation, computing the

droplet trajectories and tracking the water as it drips onto the burning commodity.

3.5.1 Sprinkler Activation

The temperature of the sensing element of a given sprinkler is estimated from the differential equation put

forth by Heskestad and Bill [44], with the addition of a term to account for the cooling of the link by water

droplets in the gas stream from previously activated sprinklers

dTl

dt

=

p

|u|

RTI

(Tg −Tl)−

C

RTI

(Tl −Tm)−

C2

RTI

β|u| (3.69)

Here Tlis the link temperature, Tg is the gas temperature in the neighborhood of the link, Tm is the temperature of the sprinkler mount (assumed ambient), and β is the volume fraction of (liquid) water in the gas

stream. The sensitivity of the detector is characterized by the value of RTI. The amount of heat conducted

away from the link by the mount is indicated by the “C-Factor”, C. The constant C2 has been empirically determined by DiMarzo and co-workers [45, 46, 47] to be 6×106 K/(m/s) 1

2 , and its value is relatively constant

for different types of sprinklers.

The algorithm for heat detector activation is exactly the same as for sprinkler activation, except there is

no accounting for conductive losses or droplet cooling. Note that neither sprinklers nor heat detectors are

sensitive to thermal radiation.

3.5.2 Sprinkler Droplet Size Distribution

Once activation is predicted, a sampled set of spherical water droplets is tracked from the sprinkler to either

the floor or the burning commodity. In order to compute the droplet trajectories, the initial size and velocity

of each droplet must be prescribed. This is done in terms of random distributions. The initial droplet

size distribution of the sprinkler spray is expressed in terms of its Cumulative Volume Fraction (CVF),

a function that relates the fraction of the water volume (mass) transported by droplets less than a given

diameter. Researchers at Factory Mutual have suggested that the CVF for an industrial sprinkler may be

represented by a combination of log-normal and Rosin-Rammler distributions [48]

F(d) =







√

1

2π

Z d

0

1

σd

0 e

−

[ln(d

0/dm)]2

2σ2 dd0

(d ≤ dm)

1−e

−0.693(

d

dm

)

γ

(dm < d)

(3.70)

where dm is the median droplet diameter (i.e. half the mass is carried by droplets with diameters of dm

or less), and γ and σ are empirical constants equal to about 2.4 and 0.6, respectively 4. The median drop

diameter is a function of the sprinkler orifice diameter, operating pressure, and geometry. Research at

Factory Mutual has yielded a correlation for the median droplet diameter [49]

dm

D

∝ We−

1

3 (3.71)

where D is the orifice diameter of the sprinkler. The Weber number, the ratio of inertial forces to surface

tension forces, is given by

We =

ρwU

2D

σw

(3.72)

4The Rosin-Rammler and log-normal distributions are smoothly joined if σ = 2/(

√

2π(ln 2) γ) = 1.15/γ .

33
Figure 3.4: Cumulative Volume Fraction and Cumulative Number Fraction functions of the droplet size

distribution from a typical industrial-scale sprinkler. The median diameter dm is 1 mm, σ = 0.6 and γ = 2.4.

where ρw is the density of water, U is the water discharge velocity, and σw is the water surface tension

(72.8×10−3 N/m at 20 ◦C). The discharge velocity can be computed from the mass flow rate, a function of

the sprinkler’s operating pressure and K-Factor. FM reports that the constant of proportionality in Eq. (3.71)

appears to be independent of flow rate and operating pressure. Three different sprinklers were tested in their

study with orifice diameters of 16.3 mm, 13.5 mm, 12.7 mm and the constants were approximately 4.3,

2.9, 2.3, respectively. The strike plates of the two smaller sprinklers were notched, while that of the largest

sprinkler was not [49].

In the numerical algorithm, the size of the sprinkler droplets are chosen to mimic the Rosin-Rammler/lognormal distribution. A Probability Density Function (PDF) for the droplet diameter is defined

f(d) = F

0

(d)

d

3

Z ∞

0

F

0

(d

0

)

d

03

dd0(3.73)

Droplet diameters are randomly selected by equating the Cumulative Number Fraction of the droplet distribution with a uniformly distributed random variable U

U(d) = Z d

0

f(d

0

)dd0(3.74)

Figure 3.4 displays typical Cumulative Volume Fraction and Cumulative Number Fraction functions.

Every droplet from a given sprinkler is not tracked. Instead, a sampled set of the droplets is tracked.

Typically, 1,000 droplets per sprinkler per second are tracked (50 droplets every 0.05 s, depending on user

preference). The procedure for selecting droplet sizes is as follows: Suppose water is leaving the sprinkler at

a mass flow rate of  ̇m. Suppose also that the time interval for droplet insertion into the numerical simulation

is δt, and the number of droplets inserted each time interval is N. Choose N uniformly distributed random

numbers between 0 and 1, call them Ui, obtain N droplet diameters di based on the given droplet size

distribution (Eq. (3.74), and then compute a weighting constant C from the mass balance

m ̇ δt = C

N

∑

i=1

4

3

πρw



di

2

3

(3.75)

34
The mass and heat transferred from each droplet will be multiplied by the weighting factor C.

3.5.3 Sprinkler Droplet Trajectory in Air

For a sprinkler spray, the force term f in Eq. (3.2) represents the momentum transferred from the water

droplets to the gas. It is obtained by summing the force transferred from each droplet in a grid cell and

dividing by the cell volume

f =

1

2

∑ρCDπr

2

d

(ud −u)|ud −u|

δxδyδz

(3.76)

where CD is the drag coefficient, rd is the droplet radius, ud is the velocity of the droplet, u is the velocity of

the gas, ρ is the density of the gas, and δxδyδz is the volume of the grid cell. The trajectory of an individual

droplet is governed by the equation

d

dt

(mdud) = md g−

1

2

ρCD πr

2

d

(ud −u)|ud −u| (3.77)

where md is the mass of the droplet. The drag coefficient is a function of the local Reynolds number

CD =







24/Re Re < 1

241+0.15Re0.687/Re 1 < Re < 1000

0.44 1000 < Re

(3.78)

Re =

ρ|ud −u|2rd

μ

(3.79)

where μ is the dynamic viscosity of air.

3.5.4 Sprinkler Droplet Transport on a Surface

When a water droplet hits a solid horizontal surface, it is assigned a random horizontal direction and moves

at a fixed velocity until it reaches the edge, at which point it drops straight down at the same fixed velocity.

This “dripping” velocity has been measured to be on the order of 0.5 m/s [50, 51]. Penetration of water

into porous materials is handled very crudely by assigning a fraction of the water droplets that strike a solid

horizontal surface to go straight through the solid at a slow velocity. Neither the fraction nor the velocity

has been validated.

3.5.5 Mass and Energy Transfer from Droplets

The evaporation of water droplets is handled semi-empirically. A water droplet suspended in air will evaporate as a function of the droplet equilibrium vapor mass fraction, the local gas phase vapor mass fraction,

the heat transfer to the droplet, and the droplet’s motion relative to the gas. A correlation for the mass loss

rate of a droplet that involves these parameters is given here [52]

dmd

dt

= −2π rd Sh ρD (Yd −Yg) (3.80)

The subscripts d and g refer to the droplet and gas, respectively, md is the droplet mass, D is the diffusion

coefficient for water vapor into air, Y is the water vapor mass fraction, and Sh is the droplet Sherwood

number, given by a correlation involving the Reynolds and Schmidt numbers

Sh = 2+0.6 Re

1

2 Sc

1

3 (3.81)

35
The vapor mass fraction of the gas, Yg, is obtained from the overall set of mass conservation equations and

the vapor mass fraction of the droplet is obtained from the Clausius-Clapeyron equation

Xd = exp

hv Mw

R



1

Tb

−

1

Td

 ; Yd =

Xd

Xd(1−Ma/Mw) +Ma/Mw

(3.82)

where Xd is the droplet water vapor volume fraction, hv is the heat of vaporization, Mw is the molecular

weight of water, Ma is the molecular weight of air, R is the gas constant, Tb is the boiling temperature of

water and Td is the droplet temperature.

In addition to calculating the mass transfer due to evaporation, the transfer of energy must also be

calculated. The droplet heats up due to the convective heat transfer across the surface of the droplet minus

the energy required to evaporate water

md cp,w

dTd

dt

= Ad hd (Tg −Td)−

dmd

dt

hv (3.83)

Here cp,w is the specific heat of water, Ad = 4πr

2

d

is the surface area of the droplet, hd is the heat transfer

coefficient, given by

hd =

Nu k

2rd

; Nu = 2+0.6 Re

1

2 Pr

1

3 (3.84)

Nu is the Nusselt number, k is the thermal conductivity of air, and the Prandtl number, Pr, is about 0.7 for

air. The Sherwood number, Sh, is analogous to the Nusselt number, with the Schmidt number about 0.6

compared to the 0.7 for the Prandtl number.

Finally, the exchange of mass and energy between the droplets and the gas results in an additional term

that must be added to the expression for the divergence, Eq. (3.21)

∇·u = ... +

R

γ p0



ρ ∑(Yi/Mi)

∂T

∂t

+

T

Mw

m ̇

000

w



(3.85)

where  ̇m

000

w

is the water evaporation rate per unit volume. The liquid water droplets are assumed to occupy

no volume, simplifying the analysis.

3.5.6 Interaction of Droplets and Radiation

The attenuation of thermal radiation by water droplets is an important consideration, especially for water

mist systems [53]. Water droplets attenuate thermal radiation through a combination of scattering and absorption [54]. The radiation-droplet interaction must therefore be solved for both the accurate prediction of

the radiation field and for the droplet energy balance.

If the gas phase absorption and emission in Eq. (3.45) are temporarily neglected for simplicity, the

radiative transport equation becomes

s·∇Iλ(x,s) = −[κd(x,λ) +σd(x,λ)]I(x,s) +κd(x,λ) Ib,d(x,λ) + σd(x,λ)

4π

Z

4π

Φ(s,s

0

) Iλ(x,s

0

) dΩ

0

(3.86)

where κd is the droplet absorption coefficient, σd is the droplet scattering coefficient and Ib,d is the emission

term of the droplets. Φ(s,s

0

) is a scattering phase function that gives the scattered intensity from direction

s

0

to s. The local absorption and scattering coefficients are calculated from the local droplet number density

N(x) and mean diameter dm(x) as

κd(x,λ) = N(x)

R ∞

0

f(r,dm(x)) Ca(r,λ) dr

σd(x,λ) = N(x)

R ∞

0

f(r,dm(x)) Cs(r,λ) dr (3.87)

36
where r is the droplet radius and Ca and Cs are absorption and scattering cross sections, respectively, given

by Mie theory. The droplet number density function f(r,dm) is assumed to have the same form as the droplet

size distribution, but a different mean.

An accurate computation of the in-scattering integral on the right hand side of Eq (3.86) would be

extremely time consuming. It is here approximated by dividing the total 4π solid angle to a “forward angle”

δΩland “ambient angle” δΩ∗ = 4π−δΩl. For compatibility with the FVM solver, δΩlis set equal to the

control angle given by the angular discretization. However, it is assumed to be symmetric around the center

of the control angle. Within δΩlthe intensity is Iλ(x,s) and elsewhere it is approximated as

U

∗

(x,λ) = U(x,λ)−δΩl

Iλ(x,s)

δΩ∗

(3.88)

where U(x) is the total integrated intensity. The in-scattering integral can now be written as

σd(x,λ)

4π

Z

4π

Φ(s,s

0

) Iλ(x,s

0

) dΩ

0 = σd(x,λ)[χf Iλ(x,s) + (1−χf)U∗

(x)] (3.89)

where χf = χf(r,λ) is a fraction of the total intensity originally within the solid angle δΩlthat is scattered

into the same angle δΩl. Defining now an effective scattering coefficient section

σ ̄ d(x,λ) = 4πN(x)

4π−δΩl

Z ∞

0

(1−χf) Cs(r,λ) dr (3.90)

the spray RTE becomes

s·∇Iλ(x,s) = −[κd(x,λ) +σ ̄ d(x,λ)]I(x,s) +κd(x,λ) Ib,d(x,λ) + σ ̄ d(x,λ)

4π

U(x,λ) (3.91)

This equation can be integrated over the spectrum to get the band specific RTE’s. The procedure is exactly

the same as what is used for the gas phase RTE. After the band integrations, the spray RTE for band n

becomes

s·∇In(x,s) = −[κd,n(x) +σ ̄ d,n(x)]In(x,s) +κd,n(x) Ib,d,n(x) + σ ̄ d(x,λ)

4π

Un(x) (3.92)

where the source function is based on the average droplet temperature within a cell. The absorption and

scattering cross sections and the scattering phase function are calculated using the MieV code developed

by Wiscombe [55]. Both κd and σ ̄ d are averaged over the possible droplet radii and wavelength before

the actual simulation. A single constant temperature is used in the wavelength averaging. This “radiation

temperature” Trad should be selected to represent a typical radiating flame temperature. A value 1173 K is

used by default. The averaged quantities, being now functions of the droplet mean diameter only, are saved

in one-dimensional arrays. During the simulation, the local properties are calculated as a table lookup using

the local mean droplet diameter [56]. Details of the computation are given in Section 3.7.8.

3.5.7 Fire Suppression by Water

The above two sections describe heat transfer from a droplet of water to a hot gas, a hot solid, or both.

Although there is some uncertainty in the values of the respective heat transfer coefficients, the fundamental

physics are fairly well understood. However, when the water droplets encounter burning surfaces, simple

heat transfer correlations become more difficult to apply. The reason for this is that the water is not only

cooling the surface and the surrounding gas, but it is also changing the pyrolysis rate of the fuel. If the

surface of the fuel is planar, it is possible to characterize the decrease in the pyrolysis rate as a function of

the decrease in the total heat feedback to the surface. Unfortunately, most fuels of interest in fire applications

are multi-component solids with complex geometry at scales unresolvable by the computational grid.

37
To date, most of the work in this area has been performed at Factory Mutual. An important paper on

the subject is by Yu et al. [57]. The authors consider dozens of rack storage commodity fires of different

geometries and water application rates, and characterize the suppression rates in terms of a few global

parameters. Their analysis yields an expression for the total heat release rate from a rack storage fire after

sprinkler activation

Q ̇ = Q ̇

0 e

−k(t−t0)

(3.93)

where Q ̇

0 is the total heat release rate at the time of application t0, and k is a fuel-dependent constant. For

the FMRC Standard Plastic commodity k is given as

k = 0.716  ̇m

00

w −0.0131 s−1

(3.94)

where  ̇m

00

w

is the flow rate of water impinging on the box tops, divided by the area of exposed surface (top

and sides). It is expressed in units of kg/m2/s. For the Class II commodity, k is given as

k = 0.536  ̇m

00

w −0.0040 s−1

(3.95)

Unfortunately, this analysis is based on global water flow and burning rates. Equation (3.93) accounts

for both the cooling of non-burning surfaces as well as the decrease in heat release rate of burning surfaces.

In the FDS model, the cooling of unburned surfaces and the reduction in the heat release rate are computed

locally, thus it is awkward to apply a global suppression rule. However, the exponential nature of suppression

by water is observed both locally and globally, thus it is assumed that the local burning rate of the fuel can

be expressed in the form [50, 51]

m ̇

00

f

(t) = m ̇

00

f,0

(t) e

−

R

k(t)dt (3.96)

Here  ̇m

00

f,0

(t) is the burning rate per unit area of the fuel when no water is applied and k(t) is a linear function

of the local water mass per unit area, m

00

w

, expressed in units of kg/m2,

k(t) = a m00

w

(t) s

−1

(3.97)

Note that a is an empirical constant.

38
3.6 Smoke Detectors

An informative discussion of the issues associated with smoke detection can be found in the SFPE Handbook

chapter “Design of Detection Systems,” by Schifiliti, Meacham and Custer [20]. The authors point out that

the difficulty in modeling smoke detector activation stems from a number of issues: (1) the production and

transport of smoke in the early stage of a fire are not well-understood, (2) detectors often use complex

response algorithms rather than simple threshold or rate-of-change criteria, (3) detectors can be sensitive to

smoke particle number density, size distribution, refractive index, composition, etc., and (4) most computer

models, including FDS, do not provide detailed descriptions of the smoke besides its bulk transport. This

last point is the most important. At best, in its present form, FDS can only provide to an activation algorithm

the velocity and smoke concentration of the ceiling jet flowing past the detector. Regardless of the detailed

mechanism within the device, any activation model included within FDS can only account for the entry

resistance of the smoke due to the geometry of the detector. Issues related to the effectiveness of ionization

or photoelectric type detectors cannot be addressed by FDS.

Consider the simple idealization of a “spot-type” smoke detector. A disk-shaped cover lined with a fine

mesh screen forms the external housing of the device, which is usually mounted to the ceiling. Somewhere

within the device is a relatively small sensing chamber where the smoke is actually detected in some way. A

simple model of this device has been proposed by Heskestad. He suggested that the mass fraction of smoke

in the sensing chamber of the detector Yc lags behind the mass fraction in the external free stream Ye by

a time period δt = L/u, where u is the free stream velocity and L is a length characteristic of the detector

geometry. The change in the mass fraction of smoke in the sensing chamber can be found by solving the

following equation:

dYc

dt

=

Ye(t)−Yc(t)

δt

(3.98)

The detector activates when Yc rises above a detector-specific threshold.

A more detailed model of smoke detection involving two filling times rather than one has also been

proposed. Smoke passing into the sensing chamber must first pass through the exterior housing, then it must

pass through a series of baffles before arriving at the sensing chamber. There is a time lag associated with

the passing of the smoke through the housing and also the entry of the smoke into the sensing chamber. Let

δte be the characteristic filling time of the entire volume enclosed by the external housing. Let δtc be the

characteristic filling time of the sensing chamber. Cleary et al. [58] suggested that each characteristic filling

time is a function of the free-stream velocity u outside the detector

δte = αeu

βe

; δtc = αcu

βc

(3.99)

The α’s and β’s are empirical constants related to the specific detector geometry The change in the mass

fraction of smoke in the sensing chamber Yc can be found by solving the following equation:

dYc

dt

=

Ye(t −δte)−Yc(t)

δtc

(3.100)

Ye is the mass fraction of smoke outside of the detector in the free-stream. A simple interpretation of the

equation is that the concentration of the smoke that enters the sensing chamber at time t is that of the freestream at time t −δte.

An analytical solution for Eq. (3.100) can be found, but it is more convenient to simply integrate it

numerically as is done for sprinklers and heat detectors. The predicted mass fraction of smoke in the sensing

chamber Yc(t) can be expressed as an obscuration per unit length as follows:



1−

I

I0



×100 =



1−e

−κρYcl



×100 (3.101)

39
where κ is the specific extinction coefficient, ρ is the density of the external gases in the ceiling jet, and

l is the preferred unit of length (typically 1 m or 1 ft). For most flaming fuels, a suggested value for κ is

8700 m2/kg ± 1100 m2/kg at a wavelength of 633 nm [59].

The SFPE Handbook has references to various works on smoke detection and suggested values for the

characteristic length L. FDS includes the one parameter Heskestad model as a special case of the four

parameter Cleary model. For the Cleary model, one must set αe, βe, αc, and βc, whereas for the Heskestad

model only L = αc need be prescribed. Eq. (3.100) is still used, with αe = 0 and βe = βc = −1. Proponents

of the four parameter model claim that the two filling times are needed to better capture the behavior of

detectors in a very slow free-stream (u < 0.5 m/s). Rather than declaring one model better than another, the

algorithm included in FDS allows the user to pick these various parameters, and in so doing, pick whichever

model the user feels is appropriate.

40
3.7 Numerical Method

First the equations that are being solved are presented. Each of the conservation equations emphasize the

importance of the velocity divergence and vorticity fields, as well as the close relationship between the

thermally expandable fluid equations [9] and the Boussinesq equations for which the authors have developed

highly efficient solution procedures [60, 61]. All spatial derivatives are approximated by second order finite

differences and the flow variables are updated in time using an explicit second order predictor-corrector

scheme.

3.7.1 Simplified Equations

Regardless of whether one is performing an LES or a DNS calculation, the overall solution algorithm is the

same.

Conservation of Mass

∂ρ

∂t

+u·∇ρ = −ρ∇·u (3.102)

Conservation of Species

∂ρYl

∂t

+u·∇ρYl = −ρYl∇·u+∇· ρD∇Yl +m ̇

000

l

(3.103)

Conservation of Momentum

∂u

∂t

+u×ω+∇H =

1

ρ

((ρ−ρ∞)g+f+∇· τi j) (3.104)

Divergence Constraint

∇·u =

1

ρcp T

 

∇· k∇T +∇·∑

l

R

cp,ldT ρDl∇Yl −∇·qr +q ̇

000!

+



1

ρcp T

−

1

p0



d p0

dt

(3.105)

Equation of State

p0(t) = ρTR ∑

l

Yl/Ml (3.106)

Notice that the source terms from the energy conservation equation have been incorporated into the divergence and ultimately are involved in the mass conservation equation. The temperature is found from the

density and background pressure via the equation of state.

3.7.2 Temporal Discretization

All calculations start with ambient initial conditions. At the beginning of each time step, the quantities ρ

n

,

Y

n

i

, u

n

, H n, and p

n

0

are known. All other quantities can be derived from them. Note that the superscript

(n+1)e refers to an estimate of the value of the quantities at the (n+1)st time step.

1. The thermodynamic quantities ρ, Yi, and p0 are estimated at the next time step with an explicit Euler

step. For example, the density is estimated

ρ

(n+1)e = ρn −δt(un

·∇ρn +ρ

n∇·un

) (3.107)

The divergence (∇ · u)

(n+1)e

is formed from these estimated thermodynamic quantities. The normal

velocity components at boundaries that are needed to form the divergence are assumed known.

41
2. A Poisson equation for the pressure is solved with a direct solver

∇

2Hn = −

(∇·u)

(n+1)e −(∇·u)n

δt

−∇·F

n

(3.108)

Note that the vector F contains the convective, diffusive and force terms of the momentum equation.

These will be described in detail below. Then the velocity is estimated at the next time step

u

(n+1)e = un −δt(Fn +∇Hn

) (3.109)

Note that the divergence of the estimated velocity field is identically equal to the estimated divergence

(∇·u)

(n+1)e

that was derived from the estimated thermodynamic quantities. The time step is checked

at this point to ensure that

δt max

δx

|u|

,

δy

|v|

,

δz

|w|



< 1 ; 2 δt ν



1

δx

2

+

1

δy

2

+

1

δz

2



< 1 (3.110)

If the time step is too large, it is reduced so that it satisfies both constraints and the procedure returns

to the beginning of the time step. If the time step satisfies the stability criteria, the procedure continues

to the corrector step. See Section 3.7.10 for more details on stability.

3. The thermodynamic quantities ρ, Yi, and p0 are corrected at the next time step. For example, the

density is corrected

ρ

n+1 =

1

2



ρ

n +ρ(n+1)e −δt(u(n+1)e

·∇ρ(n+1)e +ρ

(n+1)e∇·u(n+1)e

)



(3.111)

The divergence (∇·u)

(n+1)

is derived from the corrected thermodynamic quantities.

4. The pressure is recomputed using estimated quantities

∇

2H(n+1)e = −

2(∇·u)

n+1 −(∇·u)(n+1)e −(∇·u)n

δt

−∇·F

(n+1)e

(3.112)

The velocity is then corrected

u

n+1 =

1

2

h

u

n +u(n+1)e −δt



F

(n+1)e +∇H(n+1)e

i (3.113)

Note again that the divergence of the corrected velocity field is identically equal to the corrected

divergence.

3.7.3 Spatial Discretization

Spatial derivatives in the governing equations are written as second order accurate finite differences on a

rectilinear grid. The overall domain is a rectangular box that is divided into rectangular grid cells. Each cell

is assigned indices i, j and k representing the position of the cell in the x, y and z directions, respectively.

Scalar quantities are assigned in the center of each grid cell, thus ρ

n

i jk is the density at the nth time step

in the center of the cell whose indices are i, j and k. Vector quantities like velocity are assigned at cell

faces, thus the x component of velocity u is defined at the faces whose normals are parallel to the x-axis,

the y component v is defined at the faces whose normals are parallel to the y-axis, and the z component w is

defined at the faces whose normals are parallel to the z-axis. The quantity u

n

i jk is the x component of velocity

at the forward pointing face of the i jkth cell; u

n

i−1, jk is at the backward pointing face of the i jkth cell.

42
3.7.4 Large Eddy vs. Direct Numerical Simulation

The major difference between an LES and a DNS calculation is the form of the viscosity, and the thermal

and material diffusivities. For a Large Eddy Simulation, the dynamic viscosity is defined at cell centers

μi jk = ρi jk (Cs ∆)

2

|S| (3.114)

where Csis an empirical constant, ∆ = (δxδyδz)

1

3 , and

|S|

2 = 2



∂u

∂x

2

+2



∂v

∂y

2

+2



∂w

∂z

2

+



∂u

∂y

+

∂v

∂x

2

+



∂u

∂z

+

∂w

∂x

2

+



∂v

∂z

+

∂w

∂y

2

−

2

3

(∇·u)

2

(3.115)

The quantity |S| consists of second order spatial differences averaged at cell centers. The thermal conductivity and material diffusivity of the fluid are related to the viscosity by

ki jk =

cp,0 μi jk

Pr

; (ρD)i jk =

μi jk

Sc

(3.116)

where Pr is the Prandtl number and Sc is the Schmidt number, both assumed constant. Note that the specific

heat cp,0 is that of the dominant species of the mixture. Based on simulations of smoke plumes, Csis 0.20,

Pr and Sc are 0.5. There are no rigorous justifications for these choices other than through comparison with

experimental data [62].

The dynamic viscosity, thermal conductivity and diffusion coefficients for a DNS calculation are defined

at cell centers

μi jk = ∑

l

Yl,i jk μl(Ti jk) (3.117)

ki jk = ∑

l

Yl,i jk kl(Ti jk) (3.118)

Dl,i jk = Dl0(Ti jk) (3.119)

where the values for each individual species are approximated from kinetic theory [30]. The term Dl0 is

the binary diffusion coefficient for species l diffusing into the predominant species 0, usually nitrogen. It

is often the case that the numerical grid is too coarse to resolve steep gradients in flow quantities when the

temperature is near ambient. However, as the temperature increases and the diffusion coefficients increase

in value, the situation improves. As a consequence, there is a provision in the numerical algorithm to place

a lower bound on the viscous coefficients to avoid numerical instabilities at temperatures close to ambient.

3.7.5 The Mass Transport Equations

Due to the low Mach number approximation being used in the model, the mass and energy equations are

combined by way of the divergence. The divergence of the flow field contains much of the fire-specific

source terms described above.

Convective and Diffusive Transport

The density at the center of the i jkth cell is updated in time with the following predictor-corrector scheme.

In the predictor step, the density at the (n+1)st time level is estimated based on information at the nth level

ρ

(n+1)e

i jk −ρ

n

i jk

δt

+ (u·∇ρ)

n

i jk = −ρ

n

i jk(∇·u)

n

i jk (3.120)

43
Following the prediction of the velocity and background pressure at the (n+1)st time level, the density is

corrected

ρ

(n+1)

i jk −

1

2



ρ

n

i jk +ρ

(n+1)e

i jk 

1

2

δt

+ (u·∇ρ)

(n+1)e

i jk = −ρ

(n+1)e

i jk (∇·u)

(n+1)e

i jk (3.121)

The species conservation equations are differenced the same way

(ρYl)

(n+1)e

i jk −(ρYl)

n

i jk

δt

+ (u·∇ρYl)

n

i jk = −(ρYl)

n

i jk(∇·u)

n

i jk + (∇· ρD∇Yl)

n

i jk +m ̇

000

i jk (3.122)

at the predictor step, and

(ρYl)

(n+1)

i jk −

1

2



(ρYl)

n

i jk + (ρYl)

(n+1)e

i jk 

1

2

δt

+(u·∇ρYl)

(n+1)e

i jk = −(ρYl)

(n+1)e

i jk (∇·u)

(n+1)e

i jk +(∇·ρD∇Yl)

(n+1)e

i jk +m ̇

000

i jk

(3.123)

at the corrector step.

The convective terms are written as upwind-biased differences in the predictor step and downwindbiased differences in the corrector step. In the expressions to follow, the symbol ± means + in the predictor

step and − in the corrector step. The opposite is true for ∓.

(u·∇ρ)i jk =

1∓εu

2

ui jk

ρi+1, jk −ρi jk

δx

+

1±εu

2

ui−1, jk

ρi jk −ρi−1, jk

δx

+

1∓εv

2

vi jk

ρi, j+1,k −ρi jk

δy

+

1±εv

2

vi, j−1,k

ρi jk −ρi, j−1,k

δy

+

1∓εw

2

wi jk

ρi j,k+1 −ρi jk

δz

+

1±εw

2

wi j,k−1

ρi jk −ρi j,k−1

δz

(3.124)

(u·∇ρYl)i jk =

1∓εu

2

ui jk

(ρYl)i+1, jk −(ρYl)i jk

δx

+

1±εu

2

ui−1, jk

(ρYl)i jk −(ρYl)i−1, jk

δx

+

1∓εv

2

vi jk

(ρYl)i, j+1,k −(ρYl)i jk

δy

+

1±εv

2

vi, j−1,k

(ρYl)i jk −(ρYl)i, j−1,k

δy

+

1∓εw

2

wi jk

(ρYl)i j,k+1 −(ρYl)i jk

δz

+

1±εw

2

wi j,k−1

(ρYl)i jk −(ρYl)i j,k−1

δz

(3.125)

Note that without the inclusion of the ε’s, these are simple central difference approximations. The ε’s are

local Courant numbers, εu = uδt/δx, εv = vδt/δy, and εw = wδt/δz, where the velocity components are those

that immediately follow. Their role is to bias the differencing upwind. Where the local Courant number is

near unity, the difference becomes nearly fully upwinded. Where the local Courant number is much less

than unity, the differencing is more centralized [63].

3.7.6 Discretizing the Divergence

The divergence in both the predictor and corrector step is discretized

(∇·u)i jk =

1

ρT cp

 

q ̇

000

i jk + (∇· k∇T)i jk +∑

l

(∇· hlρD∇Yl)i jk!+



1

ρT cp

−

1

p0



d p0

dt

(3.126)

The thermal and material diffusion terms are pure central differences, with no upwind or downwind bias,

thus they are differenced the same way in both the predictor and corrector steps. For example, the thermal

44
conduction term is differenced as follows:

(∇· k∇T)i jk =

1

δx



ki+

1

2

, jk

Ti+1, jk −Ti jk

δx

−ki−

1

2

, jk

Ti jk −Ti−1, jk

δx



+

1

δy



ki, j+

1

2

,k

Ti, j+1,k −Ti jk

δy

−ki, j−

1

2

,k

Ti jk −Ti, j−1,k

δy



+

1

δz



ki j,k+

1

2

Ti j,k+1 −Ti jk

δz

−ki j,k−

1

2

Ti jk −Ti j,k−1

δz



(3.127)

The temperature is extracted from the density via the equation of state

Ti jk =

p0

ρi jkR ∑

N

l=0

(Yl,i jk/Ml)

(3.128)

Because only species 1 through N are explicitly computed, the summation is rewritten

N

∑

l=0

Yl,i jk

Ml

=

1

M0

+

N

∑

l=1



1

Ml

−

1

M0



Yl (3.129)

In isothermal calculations involving multiple species, the density can be extracted from the average molecular weight

ρi jk =

p0

T∞R ∑

N

l=0

Yl,i jk/Ml

(3.130)

Again, because only species 1 through N are explicitly computed, this expression can be written

ρi jk =

M0 p0

T∞R

+

N

∑

l=1



1−

M0

Ml



(ρYl)i jk (3.131)

Heat Release Rate (Mixture Fraction)

When using the mixture fraction-based combustion model, we must extract the local heat release rate per

unit volume from the computed mixture fraction field. Ideally, the reaction of fuel and oxygen is assumed

to occur at an infinitesimally thin sheet where Z = Zf, the stoichiometric value of the mixture fraction. In

the finite difference calculation, the location of the sheet must be determined, and then the energy that is

released must be distributed to the appropriate grid cells. To motivate the numerical procedure, we will work

with the exact (i.e. non-discrete) form of the relevant expressions.

The heat release rate per unit area of flame sheet is

q ̇

00 = ∆HO

dYO

dZ









Z<Zf

ρD ∇Z ·n (3.132)

where n is the outward facing unit normal. Note that both dYO/dZ and ∇Z · n are negative. The finite

difference approximation for this expression must be accurate to ensure that the integrated heat release rate

over the flame sheet is consistent with the prescribed fuel flow rate. For the sake of simplicity, take a

simple gas burner with a constant fuel flow rate  ̇mb. At any instant in time, the mass flow of fuel (plus fuel

by-products) through the flame sheet above the burner is given by

m ̇ f(t) = Z

f

ρZu·n dS−

Z

f

ρD ∇Z ·n dS (3.133)

45
where the surface integrals are over the flame sheet whose normal direction points away from the fuel side.

Averaged over time, the mass of fuel and its by-products flowing through the flame will equal the mass of

fuel flowing out of the burner. At the flame sheet (Z = Zf) the convective term can be simplified

Z

f

ρZu·n dS = Zf

Z

f

ρu·n dS ; Zf =

Y

∞

O

sYI

F +Y

∞

O

; s =

νOMO

νFMF

(3.134)

Averaged over time, both the integral R

f

ρu · n dS and the quantity  ̇mf(t) will equal the mass flow of fuel

from the burner  ̇mb, yielding

m ̇ b (1−Zf) = −

Z

f

ρD∇Z ·n dS (3.135)

Using the identities

dYO

dZ









Z<Zf

= −

Y

∞

O

Zf

= −(sYI

F +Y

∞

O ) ; 1−Zf =

sYI

F

sYI

F +Y

∞

O

(3.136)

we arrive at

m ̇ b Y

I

F

s ∆HO = ∆HO

dYO

dZ









Z<Zf

Z

f

ρD ∇Z ·n dS (3.137)

The quantity  ̇mb Y

I

F

is the mass flow of fuel (allowing for the possibility that the fuel is diluted), and s ∆HO

is the heat of combustion. This exercise confirms that if we integrate Eq. (3.132) over the flame sheet and

time-average the result, we obtain the total heat release rate of the fire.

Suppose now we repeat the exercise with the finite difference form of the equations. Ultimately, we seek

a discretization of Eq. (3.132) so that when we numerically integrate it over the flame sheet, we obtain the

desired heat release rate. The flame sheet is now approximated by a stair-stepped surface formed by grid cell

interfaces such that cells inside the surface have mixture fraction values greater than stoichiometric (Z > Zf)

and cells outside the surface have values that are less (Z < Zf). Remember that the mixture fraction and all

other scalar quantities are defined at cell centers. The mass flux of mixture fraction at the interface between

cell i jk and cell i+1, jk is written in conservative form (excluding the subscript jk to reduce clutter)

m ̇

00

f = ui(ρZ)

i+

1

2

−(ρD)

i+

1

2

Zi+1 −Zi

δx

(3.138)

Note that the horizontal component of velocity uiis defined at the interface between the cell i and cell (i+1),

and that the subscript (i +

1

2

) denotes the average of the values of a given quantity in cell i and (i + 1). It

is important to remember that the value of the mixture fraction at the cell interface is close, but not exactly

equal to, Zf. To account for this, we rewrite the convective component of the mass flux as follows:

ui(ρZ)

i+

1

2

= ui ρi+

1

2

Zf +ui ρi+

1

2

(Z −Zf) ; Z =

(ρZ)

i+

1

2

ρi+

1

2

=

ρi+1Zi+1 +ρiZi

ρi+1 +ρi

(3.139)

The term ui ρi+

1

2

Zf when summed over the stair-stepped flame surface will equal exactly (via telescoping

sums) Zf multiplied by the mass flow through the surface. When averaged over time, the mass flow through

the flame surface will be equal to the mass flow of fuel through the burner. This reasoning is completely

analogous to that above which resulted in Eq. (3.135). The term ui ρi+

1

2

(Z −Zf) does not appear in the analytical derivation above. It represents the error associated with approximating the flame sheet at the interface

between two grid cells. It must be accounted for in the finite difference approximation to Eq. (3.132):

q ̇

00 = ∆HO

dYO

dZ









Z<Zf

(

(ρD)

i+

1

2

(Zi+1 −Zi)/δx−ui ρi+

1

2

(Z −Zf) Zi > Zf > Zi+1

(ρD)

i+

1

2

(Zi −Zi+1)/δx+ui ρi+

1

2

(Z −Zf) Zi < Zf < Zi+1

(3.140)

The energy released at the cell interfaces that make up the approximate flame surface is divided between the

two cells on either side of the interface.

46
Heat Release Rate (Finite-Rate Reaction)

In a DNS calculation (usually), a one-step, finite-rate reaction of a hydrocarbon fuel is assumed

νCxHyCxHy +νO2 O2 −→ νCO2CO2 +νH2O H2O (3.141)

The single step finite rate reaction model provides a chemical “switch” that can turn the idealized combustion process on or off without arbitrarily introducing rules governing combustible oxygen concentrations or

ignition temperatures. In any cell where the reaction is “on,” the chemical reaction time scale is much shorter

than any convective or diffusive transport time scale. Thus, it makes sense to calculate the consequences of

the reaction assuming all other processes are frozen in a state corresponding to the beginning of the time

step. For each grid cell, at the start of a time step where t = t

n

and Y

n

CxHy,i jk ≡ YF(t

n

) and Y

n

O2,i jk ≡ YO(t

n

),

the following ODE is solved numerically with a 2nd order Runge-Kutta scheme

dYF

dt

= −

Bρ

a+b−1

i jk

Mb

O Ma−1

F

YF(t)

aYO(t)b

e

−E/RTi jk (3.142)

dYO

dt

= −

νO MO

νF MF

dYF

dt

(3.143)

The temperature Ti jk and density ρi jk are fixed at their values at time t

n

and the ODE is iterated from t

n

to

t

n+1

in about 100 time steps. The pre-exponential factor B, the activation energy E, and the exponents a and

b are input parameters. The average heat release rate over the entire time step is given by

q ̇

000n

i jk = ∆H ρ

n

i jk

YF(t

n

)−YF(t

n+1

)

δt

(3.144)

where δt = t

n+1 −tn

. The species mass fractions are adjusted at this point in the calculation (before the

convection and diffusion update)

Y

n

l,i jk = Yl(t

n

)−

νl Ml

νF MF

(YF(t

n

)−YF(t

n+1

)) (3.145)

3.7.7 Thermal Radiation

The discretized RTE is derived by integrating equation (3.46) over the grid cell i jk and the control angle δΩl

Z

Ωl

Z

Vi jk

s·∇I(x,s)dV dΩ =

Z

Ωl

Z

Vi jk

κ(x) [Ib(x)−I(x,s)]dV dΩ (3.146)

The volume integral on the left hand side is replaced by a surface integral over the cell faces using the

divergence theorem. Assuming that the radiation intensity I(x,s) is constant on each of the cell faces, the

surface integral can be approximated by a sum over the cell faces. Assuming further that I(x,s) is constant

within the volume Vi jk and over the angle δΩl we obtain

6

∑

m=1

Am I

l

m

Z

Ωl

(s·nm)dΩ = κi jk hIb,i jk −I

l

i jki

Vi jk δΩl(3.147)

where

I

l

i jk radiant intensity in direction l

I

l

m radiant intensity at cell face m

Ib,i jk radiant blackbody Intensity in cell

47
δΩl

solid angle corresponding to direction l

Vi jk volume of cell i jk

Am area of cell face m

nm unit normal vector of the cell face m

It must be noticed, that while the intensity is assumed constant within the angle δΩl, its direction covers the

angle δΩlexactly.

In Cartesian coordinates5, the normal vectors nm are the base vectors of the coordinate system and the

integrals over the solid angle do not depend on the physical coordinate, but the direction only. The intensities

on the cell boundaries, I

l

m, are calculated using a first order upwind scheme. If the physical space is swept in

the direction s

l

, the intensity I

l

i jk can be directly solved from an algebraic equation. This makes the numerical

solution of the FVM very fast. Iterations are needed only to account for the reflective boundaries. However,

this is seldom necessary in practice, because of the small time step set by the flow solver.

The spatial discretization for the RTE solver is the same as for the fluid solver. The coordinate system

used to discretize the solid angle is shown in Figure 3.5. The discretization of the solid angle is done by

φ

s

x

z

y

θ

Figure 3.5: Coordinate system of the angular discretization.

dividing first the polar angle, θ, into Nθ bands, where Nθ is an even integer. Each θ-band is then divided

into Nφ(θ) parts in the azimuthal (φ) direction. Nφ(θ) must be divisible by 4. The numbers Nθ and Nφ(θ)

are chosen to give the total number of angles NΩ as close to the value defined by the user as possible. NΩ is

calculated as

NΩ =

Nθ

∑

i=1

Nφ(θi) (3.148)

The distribution of the angles is based on empirical rules that try to produce equal solid angles δΩl = 4π/NΩ.

The number of θ-bands is

Nθ = 1.17 N

1/2.26

Ω

(3.149)

rounded to the nearest even integer. The number of φ-angles on each band is

Nφ(θ) = max4,0.5NΩ



cos(θ

−)−cos(θ+)

	 (3.150)

rounded to the nearest integer that is divisible by 4. θ

− and θ+ are the lower and upper bounds of the

θ-band, respectively. The discretization is symmetric with respect to the planes x = 0, y = 0, and z = 0.

5

In the axisymmetric case equation (3.147) becomes a little bit more complicated, as the cell face normal vectors nm are not

always constant. However, the computational efficiency can still be retained.

48
This symmetry has three important benefits: First, it avoids the problems caused by the fact that first order

upwind scheme, used to calculate intensities on the cell boundaries, is more diffusive in non-axial directions

than axial. Second, the treatment of the mirror boundaries becomes very simple, as will be shown later.

Third, it avoids so called “overhang” situations, where s·i, s·j or s·k changes sign inside the control angle.

These “overhangs” would make the resulting system of linear equations more complicated.

In the axially symmetric case these “overhangs” can not be avoided, and a special treatment, developed

by Murthy and Mathur [64], is applied. In these cases Nφ(θi) is kept constant, and the total number of angles

is NΩ = Nθ ×Nφ. In addition, the angle of the vertical slice of the cylinder is chosen to be same as δφ.

The cell face intensities, I

l

m appearing on the left hand side of (3.147) are calculated using a first order

upwind scheme. Consider, for example, a control angle having a direction vector s. If the radiation is

traveling in the positive x-direction, i.e. s· i ≥ 0, the intensity on the upwind side, I

l

xu is assumed to be the

intensity in the neighboring cell, I

l

i−1 jk, and the intensity on the downwind side is the intensity in the cell

itself I

l

i jk.

On a rectilinear grid, the normal vectors nm are the base vectors of the coordinate system and the

integrals over the solid angle can be calculated analytically. Equation (3.147) can be simplified

a

l

i jkI

l

i jk = a

l

x

I

l

xu +a

l

y

I

l

yu +a

l

z

I

l

zu +b

l

i jk (3.151)

where

a

l

i jk = Ax|D

l

x

|+Ay|D

l

y

|+Az|D

l

z

|+κi jk Vi jkδΩl(3.152)

a

l

x = Ax|D

l

x

| (3.153)

a

l

y = Ay|D

l

y

| (3.154)

a

l

z = Az

|D

l

z

| (3.155)

b

l

i jk = κi jk Ib,i jk Vi jk δΩl

(3.156)

δΩl =

Z

Ωl

dΩ =

Z

δφ

Z

δθ

sinθ dθ dφ (3.157)

D

l

x =

Z

Ωl

(s

l

·i)dΩ (3.158)

=

Z

δφ

Z

δθ

(s

l

·i)sinθ dθdφ

=

Z

δφ

Z

δθ

cosφsinθsinθ dθdφ

=

1

2



sinφ

+ −sinφ−

 ∆θ−

cosθ

+ sinθ+ −cosθ− sinθ−



D

l

y =

Z

Ωl

(s

l

·j)dΩ (3.159)

=

Z

δφ

Z

δθ

sinφsinθsinθ dθdφ

=

1

2



cosφ

− −cosφ+

 ∆θ−

cosθ

+ sinθ+ −cosθ− sinθ−



D

l

z =

Z

Ωl

(s

l

·k)dΩ (3.160)

49
=

Z

δφ

Z

δθ

cosθsinθ dθdφ

=

1

2

∆φhsinθ

+

2 −

sinθ

−

2

i

Here i, j and k are the base vectors of the Cartesian coordinate system. θ

+, θ−, φ+ and φ+ are the upper

and lower boundaries of the control angle in the polar and azimuthal directions, respectively, and ∆θ =

θ

+ −θ− and ∆φ = φ+ −φ−. The solution method of (3.151) is based on an explicit marching sequence [65].

The marching direction depends on the propagation direction of the radiation intensity. As the marching is

done in the “downwind” direction, the “upwind” intensities in all three spatial directions are known, and the

intensity I

l

i jk can be solved directly. Iterations may be needed only with the reflective walls and optically

thick situations. Currently, no iterations are made.

The boundary condition on a solid wall is given as

I

l

w = ε

σT

4

w

π

+

1−ε

π ∑

Dl

0

w<0

I

l

0

w

|D

l

0

w

| (3.161)

where D

l

0

w =

R

Ωl

0(s·nw)dΩ. The constraint D

l

0

w < 0 means that only the “incoming” directions are taken into

account when calculating the reflection. The net radiative heat flux on the wall is

qw =

NΩ

∑

l=1

I

l

w

Z

δΩl

(s·nw)dΩ =

NΩ

∑

l=1

I

l

wD

l

n

(3.162)

where the coefficients D

l

n

are equal to ±D

l

x

, ±D

l

y or ±D

l

z

, and can be calculated for each wall element at the

start of the calculation.

The open boundaries are treated as black walls, where the incoming intensity is the black body intensity

of the ambient temperature. On mirror boundaries the intensities leaving the wall are calculated from the

incoming intensities using a predefined connection matrix.

I

l

w,i jk = I

l

0

(3.163)

Computationally intensive integration over all the incoming directions is avoided by keeping the solid angle

discretization symmetric x, y and z planes. The connection matrix associates one incoming direction l

0

to

each mirrored direction on each wall cell.

The local incident radiation intensity is

Ui jk =

NΩ

∑

l=1

I

l

i jkδΩl

(3.164)

3.7.8 Interaction of Droplets and Radiation

The computation of χf for a similar but simpler situation has been derived in Ref. [66]. It can be shown that

here χf becomes

χf =

1

δΩl

Z μ

l

0

Z μ

l

0

Z μd,π

μd,0

P0(θd)

(1−μ2)(1−μ

02)−(μd −μμ0)2

dμd dμdμ0(3.165)

where μd is a cosine of the scattering angle θd and P0(θd) is a single droplet scattering phase function

P0(θd) = λ

2



|S1(θd)|

2 +|S2(θd)|2



2Cs(r,λ)

(3.166)

50
S1(θd) and S2(θd) are the scattering amplitudes, given by Mie-theory. The integration limit μ

l

is a cosine of

the polar angle defining the boundary of the symmetric control angle δΩl

μ

l = cos(θl

) = 1−

2

NΩ

(3.167)

The limits of the innermost integral are

μd,0 = μμ0 +

p

1−μ

2

p

1−μ

02

; μd,π = μμ0 −

p

1−μ

2

p

1−μ

02

(3.168)

When χfis integrated over the droplet size distribution to get an averaged value, it is multiplied by Cs(r,λ).

It is therefore |S1|

2 +|S2|2

, not P0(θd), that is integrated. Physically, this means that intensities are added,

not probabilities [55].

3.7.9 Thermal and Material Boundary Conditions

Four types of thermal boundary conditions are applied at solid surfaces. The first, and simplest, is an

adiabatic boundary condition that states that there is no temperature gradient normal to the surface. It is

implemented by assigning to the grid cell that is embedded in the solid (the ghost cell) the same temperature

as the first cell in the gas(the gas cell).

The second type of boundary condition is where the solid surface has a prescribed temperature (usually

this prescribed temperature is a function of time).

The third type of boundary condition assumes the solid to be thermally-thin. The surface temperature is

updated in time according to

T

n+1

w = T

n

w +δts

q ̇

00

c +q ̇

00

r

ρscsδ

(3.169)

where Tw is the wall temperature, δtsis the time step used when updating the thermal boundary conditions

(usually greater than the hydrodynamic time step δt), and ρs, cs, δ are the input density, specific heat and

thickness of the wall. In a DNS calculation where the boundary layer is resolved, the convective flux to the

wall is given by

q ̇

00

c = −k

Tgas −Tw

δn/2

(3.170)

where δn is the size of a grid cell in the normal direction to the wall. In an LES calculation where the

boundary layer is not resolved,

q ̇

00

c = C|Tgas −Tw|

1

3 (Tgas −Tw) W/m2

(3.171)

whereC is an empirical coefficient (1.31 for vertical surface; 1.52 for horizontal), and Tgas is the temperature

of the gas in the cell bordering the wall.

The fourth type of thermal boundary condition is for a thermally-thick solid. In this case, a one dimensional heat transfer calculation is performed at each boundary cell designated as thermally-thick. The width

of the solid δ is partitioned into N cells, clustered near the front face. The cell boundaries are located at

points xi

xi = f(ξi) = δ

e

sξi/δ −1

e

s −1

(3.172)

where 0 ≤ i ≤ N, ξi = iδξ, δξ = δ/N, and 0 < s ≤ 1 is a measure of the degree of clustering of the cells at the

front face. The width of each cell is δxi = f

0

(ξi−

1

2

)δξ, 1 ≤ i ≤ N where ξi−

1

2

= (i−

1

2

)δξ. The temperature at

the center of the ith cell is denoted Ts,i. The (temperature-dependent) thermal conductivity of the solid at the

51
center of the ith cell is denoted ks,i

. The temperatures are updated in time using an implicit Crank-Nicholson

scheme

T

n+1

s,i −T

n

s,i

δt

=

1

2(ρscs)iδxi

 

ks,i+

1

2

T

n

s,i+1 −T

n

s,i

δxs,i+

1

2

− ks,i−

1

2

T

n

s,i −T

n

s,i−1

δxs,i−

1

2

+

ks,i+

1

2

T

n+1

s,i+1 −T

n+1

s,i

δxi+

1

2

− ks,i−

1

2

T

n+1

s,i −T

n+1

s,i−1

δxi−

1

2

!

(3.173)

for 1 ≤ i ≤ N. The boundary condition is discretized

−ks,1

T

n+1

s,1 −T

n+1

s,0

δx 1

2

= q ̇

00

c +q ̇

00

r −4 ε σ T

n

3

s,

1

2



T

n+1

s,

1

2

−T

n

s,

1

2



(3.174)

where Ts,

1

2

= (Ts,1 +Ts,0)/2 is the temperature at the front face. Notice that the radiative emission term has

been linearized

T

(n+1)

4

s,

1

2

−T

n

4

s,

1

2

≈ 4 T

n

3

s,

1

2



T

n+1

s,

1

2

−T

n

s,

1

2



(3.175)

The wall temperature is defined Tw ≡ Ts,

1

2

= (Ts,0 +Ts,1)/2.

Regardless of how the wall temperature is determined, there are two ways of coupling the wall temperature with the fluid calculation. Gas phase temperatures are defined at cell centers; the wall is defined at

the boundary of the bordering gas phase cell and a “ghost” cell inside the wall. As far as the gas phase calculation is concerned, the normal temperature gradient at the wall is expressed in terms of the temperature

difference between the “gas” cell and the “ghost” cell. The wall temperature affects the gas phase calculation through the prescription of the ghost cell temperature. This ghost cell temperature has no physical

meaning on its own. Only the difference between ghost and gas cell temperatures matters, for this defines

the heat transfer to the wall.

In a DNS calculation, the wall temperature is assumed to be an average of the ghost cell temperature

and the temperature of the first cell in the gas, thus the ghost cell temperature is defined

Tghost = 2Tw −Tgas (3.176)

For an LES calculation, the heat lost to the boundary is equated with an empirical expression

k

Tgas −Tghost

δn

= C|Tgas −Tw|

1

3 (Tgas −Tw) (3.177)

where δn is the distance between the center of the ghost cell and the center of the gas cell. This equation is

solved for Tghost, so that when the conservation equations are updated, the amount of heat lost to the wall is

equivalent to the empirical expression on the right hand side. Note that Tghost is purely a numerical construct.

It does not represent the temperature within the wall, but rather establishes a temperature gradient at the wall

consistent with the empirical correlation.

At solid walls there is no transfer of mass, thus the boundary condition for the lth species at a wall is

simply

Yl,ghost = Yl,gas (3.178)

where the subscripts “ghost” and “gas” are the same as above since the mass fraction, like temperature, is

defined at cell centers. At forced flow boundaries either the mass fraction Yl,w or the mass flux  ̇m

00

l

of species

l may be prescribed. Then the ghost cell mass fraction can be derived because, as with temperature, the

52
normal gradient of mass fraction is needed in the gas phase calculation. For cases where the mass fraction

is prescribed

Yl,ghost = 2Yl,w −Yl,gas (3.179)

For cases where the mass flux is prescribed, the following equation must be solved iteratively

m ̇

00

l = un

ρghostYl,ghost +ρgasYl,gas

2

−ρD

Yl,gas −Yl,ghost

δn

∓

δt u2

n

2

ρgasYl,gas −ρghostYl,ghost

δn

(3.180)

where  ̇m

00

l

is the mass flux of species l per unit area, un is the normal component of velocity at the wall

pointing into the flow domain, and δn is the distance between the center of the ghost cell and the center of

the gas cell. Notice that the last term on the right hand side is subtracted at the predictor step and added at

the corrector step, consistent with the biased unwinding introduced earlier.

Once the temperature and species mass fractions have been defined in the ghost cell, the density in the

ghost cell is computed from the equation of state

ρghost =

p0

R Tghost ∑l(Yl,ghost/Ml)

(3.181)

3.7.10 The Momentum Equation

The three components of the momentum equation are

∂u

∂t

+Fx +

∂H

∂x

= 0 ; Fx = wωy −vωz −

1

ρ



fx +

∂τxx

∂x

+

∂τxy

∂y

+

∂τxz

∂z



(3.182)

∂v

∂t

+Fy +

∂H

∂y

= 0 ; Fy = uωz −wωx −

1

ρ



fy +

∂τyx

∂x

+

∂τyy

∂y

+

∂τyz

∂z



(3.183)

∂w

∂t

+Fz +

∂H

∂z

= 0 ; Fz = vωx −uωy −

1

ρ



fz +

∂τzx

∂x

+

∂τzy

∂y

+

∂τzz

∂z



(3.184)

The spatial discretization of the momentum equations take the form

∂u

∂t

+Fx,i jk +

Hi+1, jk −Hi jk

δx

= 0 (3.185)

∂v

∂t

+Fy,i jk +

Hi, j+1,k −Hi jk

δy

= 0 (3.186)

∂w

∂t

+Fz,i jk +

Hi j,k+1 −Hi jk

δz

= 0 (3.187)

where Hi jk is taken at center of cell i jk, ui jk and Fx,i jk are taken at the side of the cell facing in the forward x

direction, vi jk and Fy,i jk at the side facing in the forward y direction, and wi jk and Fz,i jk at the side facing in

the forward z (vertical) direction. In the definitions to follow, the components of the vorticity (ωx,ωy,ωz) are

located at cell edges pointing in the x, y and z directions, respectively. The same is true for the off-diagonal

terms of the viscous stress tensor: τzy = τyz, τxz = τzx, and τxy = τyx. The diagonal components of the stress

tensor τxx, τxx, and τxx; the external force components (fx, fy, fz); and the unwinding bias terms εu, εv, and

εw are located at the respective cell faces.

Fx,i jk =



1∓εw

2

wi+

1

2

, jk ωy,i jk +

1±εw

2

wi+

1

2

, j,k−1 ωy,i j,k−1



−



1∓εv

2

vi+

1

2

, jk ωz,i jk +

1±εv

2

vi+

1

2

, j−1,k ωz,i, j−1,k



53
−

1

ρi+

1

2

, jk



fx,i jk +

τxx,i+1, jk −τxx,i jk

δx

+

τxy,i jk −τxy,i, j−1,k

δy

+

τxz,i jk −τxz,i, j,k−1

δz



(3.188)

Fy,i jk =



1∓εu

2

ui, j+

1

2

,k ωz,i jk +

1±εu

2

ui−1, j+

1

2

,k ωz,i−1, jk

−



1∓εw

2

wi, j+

1

2

,k ωx,i jk +

1±εw

2

wi, j+

1

2

,k−1 ωx,i j,k−1



−

1

ρi, j+

1

2

,k



fy,i jk +

τyx,i jk −τyx,i−1, jk

δx

+

τyy,i, j+1,k −τyy,i jk

δy

+

τyz,i jk −τyz,i, j,k−1

δz



(3.189)

Fz,i jk =



1∓εv

2

vi j,k+

1

2

ωx,i jk +

1±εv

2

vi, j−1,k+

1

2

ωx,i, j−1,k



−



1∓εu

2

ui j,k+

1

2

ωy,i jk +

1±εu

2

ui−1, j,k+

1

2

ωy,i−1, jk

−

1

ρi j,k+

1

2



fz,i jk +

τzx,i jk −τzx,i−1, jk

δx

+

τzy,i jk −τzy,i, j−1,k

δy

+

τzz,i j,k+1 −τzz,i jk

δz



(3.190)

ωx,i jk =

wi, j+1,k −wi jk

δy

−

vi j,k+1 −vi jk

δz

(3.191)

ωy,i jk =

ui j,k+1 −ui jk

δz

−

wi+1, jk −wi jk

δx

(3.192)

ωz,i jk =

vi+1, jk −vi jk

δx

−

ui, j+1,k −ui jk

δy

(3.193)

τxx,i jk = μi jk 2

ui jk −ui−1, jk

δx

−

2

3

(∇·u)i jk≡ μi jk 

4

3

(∇·u)i jk −2

vi jk −vi, j−1,k

δy

−2

wi jk −wi j,k−1

δz



(3.194)

τyy,i jk = μi jk 2

vi jk −vi, j−1,k

δy

−

2

3

(∇·u)i jk≡ μi jk 

4

3

(∇·u)i jk −2

ui jk −ui−1, jk

δx

−2

wi jk −wi j,k−1

δz



(3.195)

τzz,i jk = μi jk 2

wi jk −wi j,k−1

δz

−

2

3

(∇·u)i jk≡ μi jk 

4

3

(∇·u)i jk −2

ui jk −ui−1, jk

δx

−2

vi jk −vi j−1,k

δy



(3.196)

τxy,i jk = τyx,i jk = μi+

1

2

, j+

1

2

,k



ui, j+1,k −ui jk

δy

+

vi+1, jk −vi jk

δx



(3.197)

τxz,i jk = τzx,i jk = μi+

1

2

, j,k+

1

2



ui j,k+1 −ui jk

δz

+

wi+1, jk −wi jk

δx



(3.198)

τyz,i jk = τzy,i jk = μi, j+

1

2

,k+

1

2



vi j,k+1 −vi jk

δz

+

wi, j+1,k −wi jk

δy



(3.199)

εu =

uδt

δx

(3.200)

εv =

vδt

δy

(3.201)

εw =

wδt

δz

(3.202)

The variables εu, εv and εw are local Courant numbers evaluated at the same locations as the velocity component immediately following them, and serve to bias the differencing of the convective terms in the upwind

direction. The subscript i+

1

2

indicates that a variable is an average of its values at the ith and the (i+1)th

cell. The divergence defined in Eq. (3.126) is identically equal to the divergence defined by

(∇·u)i jk =

ui jk −ui−1, jk

δx

+

vi jk −vi, j−1,k

δy

+

wi jk −wi j,k−1

δz

(3.203)

54
The equivalence of the two definitions of the divergence is a result of the form of the discretized equations,

the time-stepping scheme, and the direct solution of the Poisson equation for the pressure.

Force Terms

The external force term components, in addition to including the effects of buoyancy, may also include the

drag force from sprinkler droplets.

fx,i jk =

1

2

∑ρCDπr

2

d

(ud −ui jk)|ud −u|

δxδyδz

−(ρi+

1

2

, jk −ρ∞)gx (3.204)

fy,i jk =

1

2

∑ρCDπr

2

d

(vd −vi jk)|ud −u|

δxδyδz

−(ρi, j+

1

2

,k −ρ∞)gy (3.205)

fz,i jk =

1

2

∑ρCDπr

2

d

(wd −wi jk)|ud −u|

δxδyδz

−(ρi j,k+

1

2

−ρ∞)gz (3.206)

where g = (gx,gy,gz) is the gravity vector, rd is the radius of a droplet, u = (ud,vd,wd) the velocity of a

droplet, CD the drag coefficient, and δxδyδz the volume of the i jkth cell. The summations represent all

droplets within a grid cell centered about the x, y and z faces of a grid cell respectively.

Time Step

The time step is constrained by the convective and diffusive transport speeds via two conditions. The first is

known as the Courant-Friedrichs-Lewy (CFL) condition:

δt max

|ui jk|

δx

,

|vi jk|

δy

,

|wi jk|

δz



< 1 (3.207)

The estimated velocities u

(n+1)e

, v

(n+1)e and w(n+1)e are tested at each time step to ensure that the CFL

condition is satisfied. If it is not, then the time step is set to 0.8 of its allowed maximum value and the

estimated velocities are recomputed (and checked again). The CFL condition asserts that the solution of the

equations cannot be updated with a time step larger than that allowing a parcel of fluid to cross a grid cell.

For most large-scale calculations where convective transport dominates diffusive, the CFL condition is most

important.

However, in small, finely-gridded domains, a second condition often dominates:

2 maxν,D,

k

ρcp



δt



1

δx

2

+

1

dy2

+

1

δz

2



< 1 (3.208)

Note that this constraint applied to the momentum, mass and energy equations via the relevant diffusion

parameter – viscosity, material diffusivity or thermal conductivity. This constraint on the time step, often

referred to as the Von Neumann criterion, is typical of any explicit, second-order numerical scheme for

solving a parabolic partial differential equation. The Von Neumann criterion is only invoked for DNS

calculations or for LES calculations with grid cells smaller than 5 mm.

3.7.11 The Pressure Equation

The divergence of the momentum equation yields a Poisson equation for the pressure

Hi+1, jk −2Hi jk +Hi−1, jk

δx

2

+

Hi, j+1,k −2Hi jk +Hi, j−1,k

δy

2

+

Hi j,k+1 −2Hi jk +Hi j,k−1

δz

2

= −

Fx,i jk −Fx,i−1, jk

δx

−

Fy,i jk −Fy,i, j−1,k

δy

−

Fz,i jk −Fz,i j,k−1

δz

−

∂

∂t

(∇·u)i jk (3.209)

55
The lack of a superscript implies that all quantities are to be evaluated at the same time level. This

elliptic partial differential equation is solved using a direct (non-iterative) FFT-based solver [67] that is part

of a library of routines for solving elliptic PDEs called CRAYFISHPAK6. To ensure that the divergence of

the fluid is consistent with the definition given in Eq. (3.21), the time derivative of the divergence is defined

∂

∂t

(∇·u)i jk =

(∇·u)

(n+1)e

i jk −(∇·u)

n

i jk

δt

(3.210)

at the predictor step, and then

∂

∂t

(∇·u)i jk =

2(∇·u)

n+1

i jk −(∇·u)

(n+1)e

i jk −(∇·u)

n

i jk

δt

(3.211)

at the corrector step. The discretization of the divergence was given in Eq. (3.126).

Direct Poisson solvers are most efficient if the domain is a rectangular region, although other geometries

such as cylinders and spheres can be handled almost as easily. For these solvers, the no-flux condition

(3.212) is simple to prescribe at external boundaries. For example, at the floor, z = 0, the Poisson solver is

supplied with the Neumann boundary condition

Hi j,1 −Hi j,0

δz

= −Fz,i j,0 (3.212)

However, many practical problems involve more complicated geometries. For building fires, doors and

windows within multi-room enclosures are very important features of the simulations. These elements may

be included in the overall domain as masked grid cells, but the no-flux condition (3.212) cannot be directly

prescribed at the boundaries of these blocked cells. Fortunately, it is possible to exploit the relatively small

changes in the pressure from one time step to the next to enforce the no-flux condition. At the start of a time

step, the components of the convection/diffusion term F are computed at all cell faces that do not correspond

to walls. At those cell faces that do correspond to solid walls, prescribe

Fn = −

∂H

∂n

∗

+

un

δt

(3.213)

where Fn is the normal component of F at the wall.

The asterisk indicates the most recent value of the pressure. Obviously, the pressure at this particular

time step is not known until the Poisson equation is solved. Equation (3.213) asserts that following the

solution of the Poisson equation for the pressure, the normal component of velocity un will be driven nearly

(but not exactly) to zero. This is approximate because the true value of the velocity time derivative depends on the solution of the pressure equation, but since the most recent estimate of pressure is used, the

approximation is very good. Also, even though there are small errors in normal velocity at solid surfaces,

the divergence of each blocked cell remains exactly zero for the duration of the calculation. In other words,

the total flux into a given obstruction is always identically zero, and the error in normal velocity is usually

at least several orders of magnitude smaller than the characteristic flow velocity. When implemented as part

of a predictor-corrector updating scheme, the no-flux condition at solid surfaces is maintained remarkably

well.

At open boundaries (say i = I), H is prescribed depending on whether the flow is incoming or outgoing

HI+

1

2

, jk = (u

2

I, jk +v

2

I, j−

1

2

,k

+w

2

I, j−

1

2

,k

)/2 uI, jk > 0

HI+

1

2

, jk = 0 uI, jk < 0

(3.214)

6CRAYFISHPAK, a vectorized form of the elliptic equation solver FISHPAK, was originally developed at the National Center

for Atmospheric Research (NCAR) in Boulder, Colorado.

56
where I is the index of the last gas phase cell in the x direction and uI, jk is the x component of velocity at the

boundary. The value of H in the ghost cell is

HI+1, jk = 2HI+

1

2

, jk −HI, jk (3.215)

3.7.12 Particle Tracking

Thermal elements are introduced into the flow field as a means of introducing heat and as a way to visualize

the flow. The position xp of each thermal element is governed by the equations

dxp

dt

= u (3.216)

The thermal element positions are updated according to the same predictor-corrector scheme that is applied

to the other flow quantities. Briefly, the position xp of a given thermal element is updated according to the

two step scheme

x

(n+1)e

p = x

n

p +δt u

n

(3.217)

x

n+1

p =

1

2



x

n

p +x

(n+1)e

p +δt u

(n+1)e



(3.218)

where the bar over the velocity vector indicates that the velocity of the fluid is interpolated at the element’s

position.

57
3.8 Review of the Theoretical Development of the Model

ASTM E 1355 [3] requires that the theoretical basis of the model be reviewed by one or more recognized

experts fully conversant with the chemistry and physics of fire phenomenon, but not involved with the

production of the model. FDS is reviewed in two ways, internally and externally. All documents issued

by the National Institute of Standards and Technology are formally reviewed internally by members of the

staff. The theoretical basis of FDS is laid out in the present document, and is subject to internal review by

staff members who are not active participants in the development of the model, but who are members of the

Fire Research Division and are considered experts in the fields of fire and combustion. Externally, papers

detailing various parts of FDS are regularly published in peer-reviewed journals and conference proceedings.

In addition, FDS is used world-wide by fire protection engineering firms who review the technical details

of the model related to their particular application. Some of these firms also publish in the open literature

reports documenting internal efforts to validate the model for a particular use. Many of these studies are

referenced in the present document.

3.8.1 Assessment of the Completeness of Documentation

As discussed in Section 2.1.6, there are three relevant documents associated with FDS and Smokeview: the

FDS Technical Reference Guide, the FDS User’s Guide [11], and the Smokeview User’s Guide [12]. The

Technical Reference Guide provides a description of the equations, assumptions, and approximations of

the various sub-models. The User’s Guides provide instruction on how to actually run the model, plus a

description of the physical and numerical parameters.

The FDS Technical Reference Guide contains a description of the governing equations and how they

are solved numerically. However, it does not contain a detailed description of the derivation of the equations or the numerical techniques. For example, the basic conservation equations are solved numerically

using Large Eddy Simulation (LES), a technique common to the CFD community. The documentation of

FDS does not provide a detailed description of this technique, but rather makes appropriate references to

supporting papers and reports. There are numerous sources that describe various parts of the model. The

basic set of equations solved in FDS were formulated by Rehm and Baum in the Journal of Research of the

National Bureau of Standards [9]. The basic hydrodynamic algorithm evolved at NIST through the 1980s

and 1990s, incorporating fairly well-known numerical schemes that are documented in books by Anderson,

Tannehill and Pletcher [13], Peyret and Taylor [14], and Ferziger and Peric [15]. This last book provides a  ́

good description of the large eddy simulation technique and provides references to many current publications on the subject. Numerical techniques appropriate for combustion systems are described by Oran and

Boris [10]. The mixture fraction combustion model is described in a review article by Bilger [16]. Basic

heat transfer theory is provided by Holman [17] and Incropera [18]. Thermal radiation is described in Siegel

and Howell [19].

Much of the current knowledge of fire science and engineering is found in the SFPE Handbook of

Fire Protection Engineering [20]. Popular textbooks in fire protection engineering include those by Drysdale [21] and Quintiere [22]. On-going research in fire and combustion is documented in several periodicals

and conference proceedings. The International Association of Fire Safety Science (IAFSS) organizes a

conference every two years, the proceedings of which are frequently referenced by fire researchers. Interscience Communications, a London-based publisher of several fire-related journals, hosts a conference

known as Interflam roughly every three years in the United Kingdom. The Combustion Institute hosts an

international symposium on combustion every two years, and in addition to the proceedings of this symposium, the organization publishes its own journal Combustion and Flame. The papers appearing in the

IAFSS conference proceedings, the Combustion Symposium proceedings, and Combustion and Flame are

all peer-reviewed, while those appearing in the Interflam proceedings are selected based on the submission

58
of a short abstract. Both the Society for Fire Protection Engineers (SFPE) and the National Fire Protection

Association (NFPA) publish peer-reviewed technical journals entitled the Journal of Fire Protection Engineering and Fire Technology. Other often-cited, peer-reviewed technical journals include the Fire Safety

Journal, Fire and Materials, Combustion Science and Technology, Combustion Theory and Modeling and

the Journal of Heat Transfer.

Research at NIST is documented in various ways beyond contributions made by staff to external journals

and conferences. NIST publishes several forms of internal reports, special publications, and its own journal

called the Journal of Research of NIST. An internal report, referred to as a NISTIR, is a convenient means

to disseminate information, especially when the quantity of data exceeds what could normally be accepted

by a journal. Often a NISTIR is streamlined and published externally, with the NISTIR serving as a more

complete record of the work performed. Previous versions of the FDS Technical Reference Guide and

User’s Guide were published as NISTIRs. The current FDS and Smokeview manuals are being published

as NIST Special Publications, distinguished from NISTIRs by the fact that they are permanently archived.

Work performed by an outside person or organization working under a NIST grant or contract is published

in the form of a NIST Grant/Contract Report (GCR). All work performed by the staff of the Building and

Fire Research Laboratory at NIST beyond 1993 is permanently stored in electronic form and made freely

available via the Internet and yearly-released compact disks (CDs).

Finally, the quality of the FDS and Smokeview user’s manuals is checked implicitly by the fact that

the majority of model users have not taken a formal training course in the actual use of the model, but are

able to read the supporting documents, perform a few sample simulations, and then systematically build

up a level of expertise appropriate for their applications. The developers receive daily feedback from users

on the clarity of the documentation and add clarifications when needed. Before new versions of the model

are released, there is a several month “beta test” period in which users test the new version using the updated documentation. This process is similar, although less formal, to that which most computer software

programs undergo.

3.8.2 Assessment of Justification of Approaches and Assumptions

The technical approach and assumptions of the model have been presented in the peer-reviewed scientific

literature and at technical conferences cited in the previous section. The major assumptions of the model, for

example the large eddy simulation technique and the mixture fraction combustion model, have undergone

a roughly 40 year development and are now documented in popular introductory text books. More specific

sub-models, like the sprinkler spray routine or the various pyrolysis models, have yet to be developed to

this extent. As a consequence, all documents produced by NIST staff are required to go through an internal

editorial review and approval process. This process is designed to ensure compliance with the technical

requirements, policy, and editorial quality required by NIST. The technical review includes a critical evaluation of the technical content and methodology, statistical treatment of data, uncertainty analysis, use of

appropriate reference data and units, and bibliographic references. The FDS and Smokeview manuals are

first reviewed by a member of the Fire Research Division, then by the immediate supervisor of the author

of the document, then by the chief of the Fire Research Division, and finally by a reader from outside the

division. Both the immediate supervisor and the division chief are technical experts in the field. Once the

document has been reviewed, it is then brought before the Editorial Review Board (ERB), a body of representatives from all the NIST laboratories. At least one reader is designated by the Board for each document

that it accepts for review. This last reader is selected based on technical competence and impartiality. The

reader is usually from outside the division producing the document and is responsible for checking that

the document conforms with NIST policy on units, uncertainty and scope. He/she does not need to be a

technical expert in fire or combustion.

Besides formal internal and peer review, FDS is subjected to continuous scrutiny because it is available

59
free of charge to the general public and is used internationally by those involved in fire safety design and

post-fire reconstruction. The source code for FDS is also released publicly, and has been used at various

universities world-wide, both in the classroom as a teaching tool as well as for research. As a result, flaws

in the theoretical development and the computer program itself have been identified and corrected. As FDS

continues to evolve, the user base will continue to serve as a means to evaluate the model, as important to

its development as the formal internal and external peer-review processes.

3.8.3 Assessment of Constants and Default Values

A comprehensive assessment of the numerical and physical parameters used in FDS is not available in one

document. Instead, specific parameters have been exercised in various verification and validation studies

performed at NIST and elsewhere. Numerical parameters are described in the Technical Reference Guide

and are subject to the internal review process at NIST, but many physical parameters are extracted from the

open literature and do not undergo a formal review. The model user is expected to assess the appropriateness

of default values provided by FDS and make changes to the default values if need be.

60
Chapter 4

Mathematical and Numerical Robustness

ASTM E 1355 [3] outlines methods to evaluate the mathematical and numerical robustness of deterministic

fire models. This process, often referred to as model verification, ensures the accuracy of the numerical solution of the governing equations. The methods include comparison with analytical solutions, code checking,

and numerical tests.

4.1 Comparison with Analytical Solutions

Most complex combustion processes, including fire, are turbulent and time-dependent. There are no closedform mathematical solutions for the fully-turbulent, time-dependent Navier-Stokes equations. CFD provides

an approximate solution for the non-linear partial differential equations by replacing them with discretized

algebraic equations that can be solved using a powerful computer. While there is no general analytical

solution for fully-turbulent flows, certain sub-models address phenomenon that do have analytical solutions,

for example, one-dimensional heat conduction through a solid. These analytical solutions can be used to

test sub-models within a complex code such as FDS. The developers of FDS routinely use such practices

to verify the correctness of the coding of the model [68, 69]. Such verification efforts are relatively simple

and routine and the results may not always be published nor included in the documentation. Examples of

routine analytical testing include:

• The radiation solver has been verified with scenarios where simple objects, like cubes or flat plates,

are positioned in simple, sealed compartments. All convective motion is turned off, the object is given

a fixed surface temperature and emissivity of one (making it a black body radiator). The heat flux

to the cold surrounding walls is recorded and compared to analytical solutions. These studies help

determine the appropriate number of solid angles to be set as the default.

• Solid objects are heated with a fixed heat flux, and the interior and surface temperatures as a function

of time are compared to analytical solutions of the one-dimensional heat transfer equation. These

studies help determine the number of nodes to use in the solid phase heat transfer model. Similar

studies are performed to check the pyrolysis models for thermoplastic and charring solids.

• Early in its development, the hydrodynamic solver that evolved to form the core of FDS was checked

against analytical solutions of simplified fluid flow phenomena. These studies were conducted at

the National Bureau of Standards (NBS)1 by Rehm, Baum and co-workers [70, 71, 72, 73]. The

emphasis of this early work was to test the stability and consistency of the basic hydrodynamic solver,

especially the velocity-pressure coupling that is vitally important in low Mach number applications.

1The National Institute of Standards and Technology (NIST) was formerly known as the National Bureau of Standards.

61
Many numerical algorithms developed up to that point in time were intended for use in high-speed flow

applications, like aerospace. Many of the techniques adopted by FDS were originally developed for

meteorological models, and as such needed to be tested to assess whether they would be appropriate

to describe relatively low-speed flow within enclosures.

• A fundamental decision made by Rehm and Baum early in the FDS development was to use a direct

(rather than iterative) solver for the pressure. In the low Mach number formulation of the NavierStokes equations, an elliptic partial differential equation for the pressure emerges, often referred to as

the Poisson equation. Most CFD methods use iterative techniques to solve the governing conservation

equations to avoid the necessity of directly solving the Poisson equation. The reason for this is that

the equation is time-consuming to solve numerically on anything but a rectilinear grid. Because FDS

is designed specifically for rectilinear grids, it can exploit fast, direct solvers of the Poisson equation,

obtaining the pressure field with one pass through the solver to machine accuracy. FDS employs

double-precision (8 byte) arithmetic, meaning that the relative difference between the computed and

the exact solution of the discretized Poisson equation is on the order of 10−12. The fidelity of the

numerical solution of the entire system of equations is tied to the pressure/velocity coupling because

often simulations can involve hundreds of thousands of time steps, with each time step consisting of

two solutions of the Poisson equation to preserve second-order accuracy. Without the use of the direct

Poisson solver, build-up of numerical error over the course of a simulation could produce spurious

results. Indeed, an attempt to use single-precision (4 byte) arithmetic to conserve machine memory

led to spurious results simply because the error per time step built up to an intolerable level.

4.2 Code Checking

An examination of the structure of the computer program can be used to detect potential errors in the numerical solution of the governing equations. The coding can be verified by a third party either manually or

automatically with profiling programs to detect irregularities and inconsistencies [3].

At NIST and elsewhere, FDS has been compiled and run on computers manufactured by IBM, HewlettPackard, Sun Microsystems, Digital Equipment Corporation, Apple, Silicon Graphics, Dell, Compaq, and

various other personal computer vendors. The operating systems on these platforms include Unix, Linux,

Microsoft Windows, and Mac OSX. Compilers used include Lahey Fortran, Digital Visual Fortran, Intel

Fortran, IBM XL Fortran, HPUX Fortran, Forte Fortran for SunOS, the Portland Group Fortran, and several

others. Each combination of hardware, operating system and compiler involves a slightly different set of

compiler and run-time options and a rigorous evaluation of the source code to test its compliance with

the Fortran 90 ISO/ANSI standard [74]. Through this process, out-dated and potentially harmful code is

updated or eliminated, and often the code is streamlined to improve its optimization on the various machines.

However, simply because the FDS source code can be compiled and run on a wide variety of platforms does

not guarantee that the numerics are correct. It is only the starting point in the process because it at least rules

out the possibility that erratic or spurious results are due to the platform on which the code is running.

Beyond hardware issues, there are several useful techniques for checking the FDS source code that

have been developed over the years. One of the most best ways is to exploit symmetry. FDS is filled with

thousands of lines of code in which the partial derivatives in the conservation equations are approximated

as finite differences. It is very easy in this process to make a mistake. Consider, for example, the finite

difference approximation of the thermal diffusion term in the i jkth cell of the three-dimensional grid:

(∇· k∇T)i jk ≈

1

δx



ki+

1

2

, jk

Ti+1, jk −Ti jk

δx

−ki−

1

2

, jk

Ti jk −Ti−1, jk

δx



+

62
1

δy



ki, j+

1

2

,k

Ti, j+1,k −Ti jk

δy

−ki, j−

1

2

,k

Ti jk −Ti, j−1,k

δy



+

1

δz



ki j,k+

1

2

Ti j,k+1 −Ti jk

δz

−ki j,k−

1

2

Ti jk −Ti j,k−1

δz



which is written as follows in the Fortran source code:

DTDX = (TMP(I+1,J,K)-TMP(I,J,K))*RDXN(I)

KDTDX(I,J,K) = .5*(KP(I+1,J,K)+KP(I,J,K))*DTDX

DTDY = (TMP(I,J+1,K)-TMP(I,J,K))*RDYN(J)

KDTDY(I,J,K) = .5*(KP(I,J+1,K)+KP(I,J,K))*DTDY

DTDZ = (TMP(I,J,K+1)-TMP(I,J,K))*RDZN(K)

KDTDZ(I,J,K) = .5*(KP(I,J,K+1)+KP(I,J,K))*DTDZ

DELKDELT = (KDTDX(I,J,K)-KDTDX(I-1,J,K))*RDX(I) +

. (KDTDY(I,J,K)-KDTDY(I,J-1,K))*RDY(J) +

. (KDTDZ(I,J,K)-KDTDZ(I,J,K-1))*RDZ(K)

This is one of the simpler constructs because the pattern that emerges within the lines of code make it fairly

easy to check. However, a mis-typing of an I or a J, a plus or a minus sign, or any of a hundred different

mistakes can cause the code to fail, or worse produce the wrong answer. A simple way to eliminate many of

these mistakes is to run simple scenarios that have perfectly symmetric initial and boundary conditions. For

example, put a hot cube in the exact center of a larger cold compartment, turn off gravity, and watch the heat

diffuse from the hot cube into the cold gas. Any simple error in the coding of the energy equation will show

up almost immediately. Then, turn on gravity, and in the absence of any coding error, a perfectly symmetric

plume will rise from the hot cube. This checks both the coding of the energy and the momentum equations.

Similar checks can be made for all of the three dimensional finite difference routines. So extensive are these

types of checks that the release version of FDS has a routine that generates a tiny amount of random noise

in the initial flow field so as to eliminate any false symmetries that might arise in the numerical solution.

The process of adding new routines to FDS is as follows: typically the routine is written by one person

(not necessarily a NIST staffer) who takes the latest version of the source code, adds the new routine, and

writes a theoretical and numerical description for the FDS Technical Reference Guide, plus a description

of the input parameters for the FDS User’s Guide. The new version of FDS is then tested at NIST with

a number of benchmark scenarios that exercise the range of the new parameters. Provisional acceptance

of the new routine is based on several factors: (1) it produces more accurate results when compared to

experimental measurement, (2) the theoretical description is sound, and (3) any empirical parameters are

obtainable from the open literature or standard bench-scale apparatus. If the new routine is accepted, it is

added to a test version of the software and evaluated by external users and/or NIST grantees whose research

is related to the subject. Assuming that there are no intractable issues that arise during the testing period,

the new routine eventually becomes part of the release version of FDS.

Even with all the code checking performed at NIST, it is still possible for errors to go unnoticed. One

remedy is the fact that the source code for FDS is publicly released. Although it consists of on the order

of 10,000 lines of Fortran statements, various researchers outside of NIST have been able to work with

it, add enhancements needed for very specific applications or for research purposes, and report back to the

developers bugs that have been detected. The source code is organized into 14 separate files, each containing

subroutines related to a particular feature of the model, like the mass, momentum, and energy conservation

equations, sprinkler activation and sprays, the pressure solver, etc. The lengthiest routines are devoted to

input, output and initialization. Most of those working with the source code do not concern themselves with

these lengthy routines but rather focus on the finite-difference algorithm contained in a few of the more

63
important files. Most serious errors are found in these files, for they contain the core of the algorithm. The

external researchers provide feedback on the organization of the code and its internal documentation, that is,

comments within the source code itself. Plus, they must compile the code on their own computers, adding to

its portability. Some of the work performed by researchers who have modified the source code is discussed

in Chapter 6. However, most of the routine error reports are via electronic mail and are undocumented.

Most of the current error reports involve routines that are not frequently used by the FDS developers. For

example, the opening of compartment doors or the breaking of windows, especially upon activation of a

heat detector, is a feature of FDS commonly used in the fire protection engineering community but less so

at NIST. As a result, numerous reports have been made over the years in which a complicated sequence of

events prescribed by the user is not carried out by the program. The errors are easy to fix, but the number

of possible permutations of events make it difficult to check them all. Another problem reported by users

are scenarios that extend the parameter range beyond which the model was originally conceived. Walls

made of foam, fires in refrigerators, gas leaks, fuel spills, etc., are just some of the phenomena that users

have attempted to model but have run into difficulty because the model parameters either have never been

exercised (e.g. very low thermal conductivities) or are not allowed (e.g. temperatures below ambient). These

reports by the users help to improve and extend the use of the model.

4.3 Numerical Tests

Numerical techniques used to solve the governing equations within a model can be a source of error in

the predicted results. The hydrodynamic model within FDS is second-order accurate in space and time.

This means that the error terms associated with the approximation of the spatial partial derivatives by finite

differences is of the order of the square of the grid cell size, and likewise the error in the approximation of

the temporal derivatives is of the order of the square of the time step. As the numerical grid is refined, the

“discretization error” decreases, and a more faithful rendering of the flow field emerges. The issue of grid

sensitivity is extremely important to the proper use of the model and will be taken up in the next chapter.

A common technique of testing flow solvers is to systematically refine the numerical grid until the

computed solution does not change, at which point the calculation is referred to as a Direct Numerical

Solution (DNS) of the governing equations. For most practical fire scenarios, DNS is not possible on

conventional computers. However, FDS does have the option of running in DNS mode, where the NavierStokes equations are solved without the use of sub-grid scale turbulence models of any kind. Because the

basic numerical method is the same for LES and DNS, DNS calculations are a very effective way to test the

basic solver, especially in cases where the solution is steady-state. Throughout its development, FDS has

been used in DNS mode for special applications. For example, FDS (or its core algorithms) have been used at

a grid resolution of roughly 1 mm to look at flames spreading over paper in a microgravity environment [75,

76, 77, 78, 79, 80], as well as ”g-jitter” effects aboard spacecraft [81]. Simulations have been compared

to experiments performed aboard the US Space Shuttle. The flames are laminar and relatively simple in

structure, and the comparisons are a qualitative assessment of the model solution. Similar studies have been

performed comparing DNS simulations of a simple burner flame to laboratory experiments [82]. Another

study compared FDS simulations of a counterflow diffusion flames to experimental measurements and the

results of a one-dimensional multi-step kinetics model [83].

Early work with the hydrodynamic solver compared two-dimensional simulations of gravity currents

with salt-water experiments [60]. In these tests, the numerical grid was systematically refined until almost

perfect agreement with experiment was obtained. Such convergence would not be possible if there were a

fundamental flaw in the hydrodynamic solver.

64
Chapter 5

Model Sensitivity

A sensitivity analysis considers the extent to which uncertainty in model inputs influences model output.

Model parameters can be the physical properties of solids and gases, boundary conditions, initial conditions,

etc. The parameters can also be purely numerical, like the size of the numerical grid. FDS typically requires

the user to provide several dozen different types of input parameters that describe the geometry, materials,

combustion phenomena, etc. By design, the user is not expected to provide numerical parameters besides the

grid size, although the optional numerical parameters are described in both the Technical Reference Guide

and the User’s Guide.

FDS does not limit the range of most of the input parameters because applications often push beyond

the range for which the model has been validated. FDS is still used for research at NIST and elsewhere,

and the developers do not presume to know in all cases what the acceptable range of any parameter is. Plus,

FDS solves the fundamental conservation equations and is much less susceptible to errors resulting from

input parameters that stray beyond the limits of simpler empirical models. However, the user is warned that

he/she is responsible for the prescription of all parameters. The FDS manuals can only provide guidance.

The grid size is the most important numerical parameter in the model, as it dictates the spatial and temporal accuracy of the discretized partial differential equations. The heat release rate is the most important

physical parameter, as it is the source term in the energy equation. Property data, like the thermal conductivity, density, heat of vaporization, heat capacity, etc., ought to be assessed in terms of their influence on

the heat release rate. Validation studies (see Chapter 6) have shown that FDS predicts well the transport of

heat and smoke when the HRR is prescribed. In such cases, minor changes in the properties of bounding

surfaces do not have a significant impact on the results. However, when the HRR is not prescribed, but

rather predicted by the model using the thermophysical properties of the fuels, the model output is sensitive

to even minor changes in these properties.

The sensitivity analyses described in this chapter are all performed in basically the same way. For a given

scenario, best estimates of all the relevant physical and numerical parameters are made, and a “baseline”

simulation is performed. Then, one by one, parameters are varied by a given percentage, and the changes

in predicted results are recorded. This is the simplest form of sensitivity analysis. More sophisticated

techniques that involve the simultaneous variation of several parameters are impractical with a CFD model

because the computation time is too long and the number of parameters too large to perform the necessary

number of calculations to generate decent statistics.

5.1 Grid Sensitivity

The most important decision made by a model user is the size of the numerical grid. In general, the finer the

numerical grid, the better the numerical solution of the equations. FDS is second-order accurate in space and

65
time, meaning that halving the grid cell size will decrease the discretization error in the governing equations

by a factor of 4. Because of the non-linearity of the equations, the decrease in discretization error does not

necessarily translate into a comparable decrease in the error of a given FDS output quantity. To find out

what effect a finer grid has on the solution, model users usually perform some form of grid sensitivity study

in which the numerical grid is systematically refined until the output quantities do not change appreciably

with each refinement. Of course, with each halving of the grid cell size, the time required for the simulation

increases by a factor of 24 = 16 (a factor of two for each spatial coordinate, plus time). In the end, a

compromise is struck between model accuracy and computer capacity.

Some grid sensitivity studies have been documented and published. Since FDS was first publicly released in 2000, significant changes in the combustion and radiation routines have been incorporated into the

model. However, the basic transport algorithm is the same, as is the critical importance of grid sensitivity. In

compiling sensitivity studies, only those that examined the sensitivity of routines no longer used have been

excluded.

As part of a project to evaluate the use of FDS version 1 for large scale mechanically ventilated enclosures, Friday [84] performed a sensitivity analysis to find the approximate calculation time based on varying

grid sizes. A propylene fire with a nominal heat release rate was modeled in FDS. There was no mechanical

ventilation and the fire was assumed to grow as a function of the time from ignition squared. The compartment was a 3 m by 3 m by 6.1 m space. Temperatures were sampled 12 cm below the ceiling. Four grid

sizes were chosen for the analysis: 30 cm, 15 cm, 10 cm, 7.5 cm. Temperature estimates were not found to

change dramatically with different grid dimensions.

Using FDS version 1, Bounagui et al. [85] studied the effect of grid size on simulation results to determine the nominal grid size for future work. A propane burner 0.1 m by 0.1 m was modeled with a heat

release rate of 1500 kW. The grid sizes used in the modeling are given in Table 5.1. The results of the model

were compared against Heskestad’s plume correlation. Table 5.1 shows the plume centerline temperature

comparisons. Obviously, the smaller grids provide better agreements among predictive methods. A similar

analysis was performed using Alpert’s ceiling jet correlation [86] that also showed better predictions and

reasonable agreement with smaller grid sizes. In a related study, Bounagui et al. [87] used FDS to evaluate

the emergency ventilation strategies in the Louis-Hippolyte-La Fontaine Tunnel in Montreal, Canada.

Table 5.1: Grid Sizes for Computational Domain used in Ref. [85]

Cases Grid Size (m)

Case 1 0.20 x 0.20 x 0.20

Case 2 0.14 x 0.14 x 0.14

Case 3 0.10 x 0.10 x 0.10

Case 4 0.08 x 0.08 x 0.08

Table 5.2: Plume Centerline Temperature Comparisons [85]

Heskestad Case 1 Case 2 Case 3 Case 4

839 C 479 C 593 C 962 C 967 C

Xin [88] used FDS to model a methane fueled square burner (1 m by 1 m) in the open. Engineering

correlations for plume centerline temperature and velocity profiles were compared with model predictions to

assess the influence of the numerical grid and the size of the computational domain. The results showed that

66
FDS is sensitive to grid size effects, especially in the region near the fuel surface, and domain size effects

when the domain width is less than twice the plume width. FDS uses a constant pressure assumption at open

boundaries. This assumption will affect the plume behavior if the boundary of the computational domain is

too close to the plume.

Ierardi and Barnett [89] used FDS version 3 to model a 0.3 m square methane diffusion burner with heat

release rate values in the range of 14.4 kW to 57.5 kW. The physical domain used was 0.6 m by 0.6 m with

uniform grid spacings of 15, 10, 7.5, 5, 3, 1.5 cm for all three coordinate directions. For both fire sizes, a

grid spacing of 1.5 cm was found to provide the best agreement when compared to McCaffrey’s centerline

plume temperature and velocity correlations [35]. Two similar scenarios that form the basis for Alpert’s

ceiling jet correlation were also modeled with FDS. The first scenario was a 1 m by 1 m, 670 kW ethanol

fire under a 7 m high unconfined ceiling. The planar dimensions of the computational domain were 14 m by

14 m. Four uniform grid spacings of 50, 33.3, 25, and 20 cm were used in the modeling. The best agreement

for maximum ceiling jet temperature was with the 33.3 cm grid spacing. The best agreement for maximum

ceiling jet velocity was for the 50 cm grid spacing. The second scenario was a 0.6 m by 0.6 m 1000 kW

ethanol fire under a 7.2 m high unconfined ceiling. The planar dimensions of the computational domain

were 14.4 m by 14.4 m. Three uniform grid spacings of 60, 30, and 20 cm were used in the modeling. The

results show that the 60 cm grid spacing exhibits the best agreement with the correlations for both maximum

ceiling jet temperature and velocity on a qualitative basis.

Petterson [90] also completed work assessing the optimal grid size for FDS version 2. The FDS model

predictions of varying grid sizes were compared to two separate fire experiments: The University of Canterbury McLeans Island Tests and the US Navy Hangar Tests in Hawaii. The first set of tests utilized a room

with approximate dimensions of 2.4 m by 3.6 m by 2.4 m and fire sizes of 55 kW and 110 kW. The Navy

Hangar tests were performed in a hangar measuring 98 m by 74 m by 15 m in height and had fires in the

range of 5.5 MW to 6.6 MW. The results of this study indicate that FDS simulations with grids of 0.15 m

had temperature predictions as accurate as models with grids as small as 0.10 m. Each of these grid sizes

produced results within 15 % of the University of Canterbury temperature measurements. The 0.30 m grid

produced less accurate results. For the comparison of the Navy Hangar tests, grid sizes ranging from 0.60 m

to 1.80 m yielded results of comparable accuracy.

Musser et al. [91] investigated the use of FDS for course grid modeling of non-fire and fire scenarios.

Determining the appropriate grid size was found to be especially important with respect to heat transfer at

heated surfaces. The convective heat transfer from the heated surfaces was most accurate when the near

surface grid cells were smaller than the depth of the thermal boundary layer. However, a finer grid size

produced better results at the expense of computational time. Accurate contaminant dispersal modeling required a significantly finer grid. The results of her study indicate that non-fire simulations can be completed

more quickly than fire simulations because the time step is not limited by the large flow speeds in a fire

plume.

McGrattan et al. [92] showed that flame height is very sensitive to the grid size in simulations using

coarse grids. To compensate, they implemented an alternative flame surface to account for the decrease in

flame height with increasing grid cell size. This approximate flame surface is not needed when the grid is

sufficiently refined. Details can be found in Section 3.2.2.

5.2 Sensitivity of Large Eddy Simulation Parameters

FDS uses the Smagorinsky form of the Large Eddy Simulation (LES) technique. This means that instead of

using the actual fluid viscosity, the model uses a viscosity of the form

μLES = ρ(Cs ∆)

2

|S|

1

2 (5.1)

67
where Cs

is an empirical constant, ∆ is a length on the order of the size of a grid cell, and the deformation

term |S| is related to the Dissipation Function, given by Eq. (3.9). Related to the “turbulent viscosity” are

comparable expressions for the thermal conductivity and material diffusivity:

kLES =

μLES cp

Pr

; (ρD)LES =

μLES

Sc

(5.2)

The Prandtl number Pr and the Schmidt number Sc are likewise considered to be “turbulent” values. Thus,

Cs, Pr and Sc are a set of empirical constants. Most FDS users simply use the default values of (0.2,0.5,0.5),

but some have explored their effect on the solution of the equations.

In an effort to validate FDS with some simple room temperature data, Zhang et al. [93] tried different

combinations of the Smagorinsky parameters, and suggested the current default values. Of the three parameters, the Smagorinsky constant Cs

is the most sensitive. Smagorinsky [25] originally proposed a value

of 0.23, but researchers over the past three decades have used values ranging from 0.1 to 0.23. There are

also refinements of the original Smagorinsky model [26, 27, 28] that do not require the user to prescribe the

constants, but rather generate them automatically as part of the numerical scheme.

5.3 Sensitivity of Radiation Parameters

Radiative heat transfer is included in FDS via the solution of the radiation transport equation for a nonscattering gray gas, and in some limited cases using a wide band model. The equation is solved using a

technique similar to finite volume methods for convective transport, thus the name given to it is the Finite

Volume Method (FVM). There are several limitations of the model. First, the absorption coefficient for

the smoke-laden gas is a complex function of its composition and temperature. Because of the simplified

combustion model, the chemical composition of the smokey gases, especially the soot content, can effect

both the absorption and emission of thermal radiation. Second, the radiation transport is discretized via

approximately 100 solid angles. For targets far away from a localized source of radiation, like a growing

fire, the discretization can lead to a non-uniform distribution of the radiant energy. This can be seen in the

visualization of surface temperatures, where “hot spots” show the effect of the finite number of solid angles.

The problem can be lessened by the inclusion of more solid angles, but at a price of longer computing

times. In most cases, the radiative flux to far-field targets is not as important as those in the near-field, where

coverage by the default number of angles is much better.

Hostikka et al. examined the sensitivity of the radiation solver to changes in the assumed soot production, number of spectral bands, number of control angles, and flame temperature. Some of the more

interesting findings were:

• Changing the soot yield from 1 % to 2 % increased the radiative flux from a simulated methane burner

about 15 %

• Lowering the soot yield to zero decreased the radiative flux about 20 %.

• Increasing the number of control angles by a factor of 3 was necessary to ensure the accuracy of the

model at the discrete measurement locations.

• Changing the number of spectral bands from 6 to 10 did not have a strong effect on the results.

• Errors of 100 % in heat flux were caused by errors of 20 % in absolute temperature.

The sensitivity to flame temperature and soot composition are consistent with combustion theory, which

states that the source term of the radiative transport equation is a function of the absorption coefficient multiplied by the absolute temperature raised to the fourth power. The number of control angles and spectral

68
bands are user-controlled numerical parameters whose sensitivities ought to be checked for each new scenario. The default values in FDS are appropriate for most large scale fire scenarios, but may need to be

refined for more detailed simulations such as a low-sooting methane burner.

5.4 Sensitivity of Thermophysical Properties of Solid Fuels

An extensive amount of verification and validation work with FDS version 4 has been performed by Hietaniemi, Hostikka, and Vaari at VTT, Finland [94]. The case studies are comprised of fire experiments

ranging in scale from the cone calorimeter (ISO 5660-1) to full-scale fire tests such as the room corner test

(ISO 9705). Comparisons are also made between FDS results and data obtained in the SBI (Single Burning

Item) Euro-classification test apparatus (EN 13823) as well as data obtained in two ad hoc experimental

configurations: one is similar to the room corner test but has only partial linings and the other is a space to

study fires in building cavities.

All of the case studies involve real materials whose properties must be prescribed so as to conform to

the assumption in FDS that solids are of uniform composition backed by a material that is either cold or

totally insulating. Sensitivity of the various physical properties and the boundary conditions were tested.

Some of the findings were:

• The measured burning rates of various materials often fell between two FDS predictions in which cold

or insulated backings were assumed for the solid surfaces. FDS lacks a multi-layer solid model.

• The ignition time of upholstery is sensitive to the thermal properties of the fabric covering, but the

steady burning rate is sensitive to the properties of the underlying foam.

• Moisture content of wooden fuels is very important and difficult to measure.

• Flame spread over complicated objects, like cables laid out in trays, can be modeled if the surface

area of the simplified object is comparable to that of the real object. This suggests sensitivity not

only to physical properties, but also geometry. It is difficult to quantify the extent of the geometrical

sensitivity.

There is little quantification of the observed sensitivities in the study. Fire growth curves can be linear to

exponential in form, and small changes in fuel properties can lead to order of magnitude changes in heat

release rate for unconfined fires. The subject will be revisited in Chapter 6 where it is noted in many of the

studies that predicting fire growth is difficult.

5.5 Summary

The basis of large eddy simulation is that accuracy increases as the numerical mesh is refined. For fire

applications, the grid sensitivity studies have shown that the accuracy of the model is a function of the

characteristic fire diameter D

∗ divided by the grid cell size. It is not enough to describe the resolution of the

calculation solely in terms of the grid cell size, but rather the grid cell size relative to the heat release rate.

For non-fire applications, there are no simple means to evaluate “good resolution.”

As a rule of thumb, in simulations of limited resolution FDS predictions are more reliable in the farfield because the substantial numerical diffusion mimics the unresolved sub-grid scale mixing. This is hard

to quantify other than through comparisons with experiment. In some of the sensitivity studies discussed

above, the authors conclude that the model works best with a cell size of a given value, and often this cell

is not the smallest one tested. In these cases, the authors have found a flow scenario where the unresolved

convective mixing is almost exactly offset by numerical diffusion. This is fortuitous, but the conclusion does

69
not necessarily extend to other scenarios. The disadvantage of any turbulence model, large eddy simulation

included, is that good results are not guaranteed on grids of limited resolution. The advantage of LES over

other turbulence models is that the solution of the actual governing equations, not a temporal or spatial

average, is obtained as the mesh is refined.

The same can be said for phenomena closer in to the fire. However, grid resolution is more critical

for near-field phenomena because numerical diffusion near the fire on coarse grids does not have the same

fortuitous effect as it does on far-field results. In general, coarse resolution will decrease temperatures and

velocities by smearing the values over the large grid cells. This can affect the radiative flux, convection to

surrounding solids, and ultimately flame spread and fire growth.

70
Chapter 6

Model Validation

Model validation is a process where an one simulates a controlled experiment and assesses the degree of

difference between measured and predicted quantities. The assessment involves the model’s input parameters, its mathematical formulation of the physical phenomena, and its interpretation of the measurements. It

is often difficult to sort out the major source of error, especially if the experiment is complicated. It is best

if the experiment is designed purely to validate a model, in which case one can eliminate error sources by

simplifying the geometry and materials. However, most fire experiments are performed for reasons other

than simply to validate a model, and they often involve real materials burning in complicated spaces.

In this chapter, a survey of FDS validation work will be presented. Some of the work has been performed

at NIST, some by its grantees and some by engineering firms using the model. Because each organization

has its own reasons for validating the model, the referenced papers and reports do not follow any particular

guidelines. Some of the works only provide a qualitative assessment of the model, concluding that the

model agreement with a particular experiment is “good” or “reasonable.” Sometimes, the conclusion is that

the model works well in certain cases, not as well in others. These studies are included in the survey because

the references are useful to other model users who may have a similar application and are interested in even

qualitative assessment. It is important to note that some of the papers point out flaws in early releases of

FDS that have been corrected or improved in more recent releases. Some of the issues raised, however, are

still subjects of active research. The research agenda for FDS is greatly influenced by the feedback provided

by users, often through publication of validation efforts.

6.1 Validation Work with Pre-Release Versions of FDS

FDS was officially released in 2000. However, for two decades various CFD codes using the basic FDS

hydrodynamic framework were developed at NIST for different applications and for research. In the mid

1990s, many of these different codes were consolidated into what eventually became FDS. Before FDS, the

various models were referred to as LES, NIST-LES, LES3D, IFS (Industrial Fire Simulator), and ALOFT

(A Large Outdoor Fire Plume Trajectory).

The NIST LES model describes the transport of smoke and hot gases during a fire in an enclosure

using the Boussinesq approximation, where it is assumed that the density and temperature variations in

the flow are relatively small [9, 70, 73, 95]. Such an approximation can be applied to a fire plume away

from the fire itself. Much of the early work with this form of the model was devoted to the formulation

of the low Mach number form of the Navier-Stokes equations and the development of the basic numerical

algorithm. Early validation efforts compared the model with salt water experiments [61, 60, 96], and fire

plumes [97, 98, 99, 29]. Clement validated the hydrodynamic model in FDS by measuring salt water flows

using Laser Induced dye Fluorescence (LIF) [100]. An interesting finding of this work was that the transition

71
from a laminar to a turbulent plume is very difficult to predict with any technique other than DNS.

Eventually, the Boussinesq approximation was dropped and simulations began to include more firespecific phenomena. Simulations of enclosure fires were compared to experiments performed by Steckler [69]. Mell et al. [68] studied small helium plumes, with particular attention to the relative roles of

baroclinic torque and buoyancy as sources of vorticity. Cleary et al. [101] used the LES model to simulate

the environment seen by multi-sensor fire detectors and performed some simple validation work to check

the model before using it. Large fire experiments were performed by NIST at the FRI test facility in Japan,

and at US Naval aircraft hangars in Hawaii and Iceland [102]. Room airflow applications were considered

by Emmerich and McGrattan [103, 104].

These early validation efforts were encouraging, but still pointed out the need to improve the hydrodynamic model by introducing the Smagorinsky form of large eddy simulation. This addition improved the

stability of the model because of the relatively simple relation between the local strain rate and the turbulent

viscosity. There is both a physical and numerical benefit to the Smagorinsky model. Physically, the viscous

term used in the model has the right functional form to describe sub-grid mixing processes. Numerically,

local oscillations in the computed flow quantities are damped if they become large enough to threaten the

stability of the entire calculation.

During the 1980s and 1990s, the Building and Fire Research Laboratory at NIST studied the burning

of crude oil under the sponsorship of the US Minerals Management Service. The aim of the work was to

assess the feasibility of using burning as a means to remove spilled oil from the sea surface. As part of the

effort, Rehm and Baum developed a special application of the LES model called ALOFT. The model was

a spin-off of the two-dimensional LES enclosure model, in which a three-dimensional steady-state plume

was computed as a two-dimensional evolution of the lateral wind field generated by a large fire blown in

a steady wind. The ALOFT model is based on large eddy simulation in that it attempts to resolve the

relevant scales of a large, bent-over plume. Validation work was performed by simulating the plumes from

several large experimental burns of crude oil in which aerial and ground sampling of smoke particulate

was performed [105]. Yamada [106] performed a validation of the ALOFT model for 10 m oil tank fire.

The results indicate that the prediction of the plume cross section 500 m from the fire agree well with the

experimental observations.

6.2 Validation of FDS since 2000

There is an on-going effort at NIST and elsewhere to validate FDS as new capabilities are added. To date,

most of the validation work has evaluated the model’s ability to predict the transport of heat and exhaust

products from a fire through an enclosure. In these studies, the heat release rate is usually prescribed, along

with the production rates of various products of combustion. More recently, validation efforts have moved

beyond just transport issues to consider fire growth, flame spread, suppression, sprinkler/detector activation,

and other fire-specific phenomena.

The validation work discussed below can be organized into several categories: Comparisons with fullscale tests conducted especially for the chosen evaluation, comparisons with previously published full-scale

test data, comparisons with standard tests, comparisons with documented fire experience, and comparisons

with engineering correlations. There is no single method by which the predictions and measurements are

compared. Formal, rigorous validation exercises are time-consuming and expensive. Most validation exercises are done simply to assess if the model can be used for a very specific purpose. While not comprehensive

on their own, these studies collectively constitute a valuable assessment of the model.

72
6.2.1 Comparison with Full-Scale Tests Conducted Specifically for the Chosen Evaluation

As part of the NIST investigation of the World Trade Center fires and collapse, a series of large scale fire

experiments were performed specifically to validate FDS [107]. The tests were performed in a rectangular

compartment 7.2 m long by 3.6 m wide by 3.8 m tall. The fires were fueled by heptane for some tests

and a heptane/toluene mixture for the others. The fuel was sprayed from a nozzle into a steel pan. The

compartment was heavily instrumented.

FDS simulations were performed before testing began to guide the design of the compartment and also

to provide a baseline set of ”blind” predictions1. A uniform numerical grid whose cells are 10 cm on a side

was chosen based on a grid resolution study. The heat release rate (HRR) of the simulated burner was set to

that which was measured in the experiments. No attempt was made to model the spray burner.

To quantify the accuracy of the predictions, the measured and predicted gas temperatures of the upper

layer were compared for each test. Because the uncertainty of the measured HRR is often the major source

of disagreement between model and experiment, an analysis was performed to assess how sensitive the

measured temperatures were to changes in the HRR. According to an empirical correlation by McCaffrey,

Quintiere and Harkleroad [108], the rise in the upper layer gas temperature ∆Tg in a compartment is related

to the overall HRR by the relation

∆Tg = 6.85 

Q ̇ 2

A0

√

H0 hk AT

1

3

(6.1)

where Q ̇ is the HRR (kW), A0 is the area of the opening (m2), H0 is the height of the opening (m), hk is

the thermal conductivity of the walls (kW/m/K), and AT is the total area of the compartment walls (m2).

What is of importance here is the fact that the temperature rise is proportional to the HRR to the 2/3 power.

The reported uncertainty in the heat release rate measurement was 5 % (one standard deviation). The 5 %

uncertainty in the HRR corresponds to a 2/3 x 5 % = 3.3 % uncertainty in the temperature rise. For upper

layer temperatures of approximately 600 ◦C, this translates to roughly 20 ◦C. The difference in measured

and predicted upper layer temperatures ranged from 5 ◦C to 20 ◦C, meaning that the predictions were

within the uncertainty range of the HRR. Even though there are uncertainties in the measurement of the

temperature itself, the discrepancy between measurement and prediction can be explained solely in terms of

the uncertainty in the HRR measurement.

To extend this analysis further, it was observed that the heat flux onto surfaces in the upper layer was

very nearly given by σT

4 where σ = 5.67 × 10−11 kW/m2

/K4and T is the temperature in degrees K. In

other words, the emissivity of the upper layer gases is nearly 1, not surprising given the high level of soot.

Given the uncertainty in the upper layer temperature rise of 3.3 %, the uncertainty in absolute temperature

at 600 ◦C is 2.3 %, leading to an estimate for the uncertainty in heat flux of 4 × 2.3 % = 9.2 %. In most

instances, the difference between measurement and prediction was within 10 %, confirming that the model

is within experimental uncertainty.

In the discussion above, it was shown that the FDS predicted upper layer temperatures and heat fluxes

were within the uncertainty range of the experiment. For this exercise it was assumed that the uncertainty

in the experiment was based solely on the uncertainty of the heat release rate measurement. Thus, it was

shown how sensitive the upper layer temperature and heat flux measurements were to the heat release rate.

In the numerical simulations there are dozens of input parameters prescribed by the model user. Often there

is no way to assess the sensitivity of these parameters except by numerical experiment; that is, running the

model with small changes to the base parameters to see what effect these have on the predictions. To this

1A blind prediction simply means that the calculation is performed before the experiment is conducted. At NIST, the results of

a blind prediction are given to the engineers conducting the experiment so that they can be compared with the measurements immediately after the experiment is completed. This is an effective test the model, and also a good way to ensure that the measurements

make sense.

73
end, the data set was exploited to assess the sensitivity of FDS to grid cell size, wall thermal properties, fire

soot yield, radiation fraction, and various other quantities. Details can be found in Reference [107].

A second set of experiments to validate FDS for use in the World Trade Center investigation is documented in Ref. [109]. The intent of these tests was to evaluate the ability of the model to simulate the growth

of a fire burning 3 office workstations within a compartment of dimensions 11 m by 7 m by 4 m, open at

one end to mimic the ventilation of 5 windows similar to those in WTC 1 and 2. Six tests were performed

with various initial conditions exploring the effect of jet fuel spray and ceiling tiles covering the surface of

the desks and carpet. Measurements were made of the heat release rate and compartment gas temperatures

at four locations using vertical thermocouple arrays. Six different material samples were tested in the NIST

cone calorimeter: desk, chair, paper, computer case, privacy panel, and carpet. Data for the carpet, desk and

privacy panel were input directly into FDS, with the other 3 materials lumped together to form an idealized

fuel type. Open burns of single workstations were used to calibrate the simplified fuel package. Then FDS

was used to make blind predictions of the 3 workstation fires within the compartment. Peak heat release

rates and temperatures were predicted to within 20 % for all tests.

6.2.2 Comparison with Engineering Correlations

There are several examples of fire flows that have been extensively studied, so much so that a set of engineering correlations combining the results of many experiments have been developed. These correlations

are useful to modelers because of their simplicity. The most studied phenomena include fire plumes, ceiling

jets, and flame heights.

Although much of the early validation work before FDS was released involved fire plumes, it remains

an active area of interest. One study by Chow and Yin [110] surveys the performance of various models in

predicting plume temperatures and entrainment. They compare various correlations, a RANS (ReynoldsAveraged Navier-Stokes) model, and FDS. Simulations were carried out that replicated a 470 kW fire with

a diameter of 1 m and an unbounded ceiling. A numerical grid size of 96 by 96 by 96 cells was used in

the FDS calculation and provided results that agreed well with those predicted by the RANS model and the

various correlations.

Battaglia et al. [111] used FDS to simulate fire whirls. First, the model was shown to reproduce the

McCaffrey correlation of a fire plume, then it was shown to reproduce qualitatively certain features of fire

whirls. At the time, FDS used Lagrangian elements to introduce heat from the fire (no longer used), and this

combustion model could not replicate the extreme stretching of the core of the flame zone.

Quintiere and Ma [33, 34] compared predicted flame heights and plume centerline temperatures to

empirical correlations. For plume temperature, the Heskestad correlation [35] was chosen. Favorable agreement was found in the plume region, but the results near the flame region were found to be grid-dependent,

especially for low Q

∗ fires. At this same time, researchers at NIST were reaching similar conclusions, and

it was noticed by both teams that a critical parameter for the model is D

∗/δx, where D∗

is the characteristic

fire diameter and δx is the grid cell size. If this parameter is sufficiently large, the fire can be considered

well-resolved and agreement with various flame height correlations was found. If the parameter is not large

enough, the fire is not well-resolved and adjustments must be made to the combustion routine to account for

it. See Section 3.2.2 for more details.

6.2.3 Comparisons with Previously Published Full-Scale Test Data

Experiments conducted solely for model validation are somewhat rare. More common are validation studies

that use data from past experiments. This section contains brief descriptions of work published comparing

FDS with past experiments or correlations of experimental data.

74
Pool Fires

Xin et al. [112] used FDS to model a 1 m diameter methane pool fire. The computational domain was 2 m

by 2 m by 4 m with a uniform grid size of 2.5 cm. The predicted results were compared to experimental

data and found to qualitatively and quantitatively reproduce the velocity field. The same authors performed

a similar study of a 7.1 cm methane burner [113] and a helium plume [114].

Hostikka et al. [115] modeled small pool fires of methane and methanol to test the FDS radiation solver

for low-sooting fires. They conclude that the predicted radiative fluxes for both fuels are higher than measured values, especially at small heat release rates, due to an over-prediction of the gas temperature.

Hietaniemi, Hostikka and Vaari [94] consider heptane pool fires of various diameters. Predictions of the

burning rate as a function of diameter follow the trend observed in a number of experimental studies. Their

results show an improvement in the model over the earlier work with methanol fires, due to improvements

in the radiation routine and the fact that heptane is more sooty than methanol, simplifying the treatment of

radiation. The authors point out that reliable predictions of the burning rate of liquid fuels require roughly

twice as fine a grid spanning the burner than would be necessary to predict plume velocities and temperatures. The reason for this is the prediction of the heat feedback to the burning surface necessary to predict

rather than to prescribe the burning rate.

Airflows in Non-Fire Compartments

The low Mach number assumption in FDS is appropriate not only to fire, but to most building ventilation

scenarios. An example of how the model can be used to assess indoor air quality is presented by Musser et

al. [91]. The test compartment was a displacement ventilation test room that contained computers, furniture, and lighting fixtures as well as heated rectangular boxes intended to represent occupants. A detailed

description of the test configuration is given by Yuan et al. [116]. The room is ventilated with cool supply

air introduced via a diffuser that is mounted on a side wall near the floor. The air rises as it is warmed by

heat sources and exits through a return duct located in the upper portion of the room. The flow pattern is intended to remove contaminants by sweeping them upward at the source and removing them from the room.

Sulphur hexafluoride, SF6, was introduced into the compartment during the experiment as a tracer gas near

the breathing zone of the occupants. Temperature, tracer concentration, and velocity were measured during

the experiments. For temperature, the two finest grids (50 by 36 by 24 and 64 by 45 by 30) produced results

in which the agreement between the measurement and prediction was considered acceptable. The agreement

for the tracer concentrations were not as good. It was suggested that the difference could be related to the

way the source of the tracer gas was modeled. The comparison of velocity data was deemed reasonable,

given the limitations of the velocity probes at low velocities.

In another study, Musser and Tan [117] used FDS to assess the design of ventilation systems for facilities

in which train locomotives operate. Although there is only a limited amount of validation, the study is useful

in demonstrating a practical use of FDS for a non-fire scenario.

Mniszewski [118] used FDS to model the release of flammable gases in simple enclosures and open

areas. In this work, the gases were not ignited.

Kerber and Walton provided a comparison between FDS version 1 and experiments on positive pressure

ventilation in a full-scale enclosure without a fire. The model predictions of velocity were within 10 % to

20 % of the experimental values [119].

Wind Engineering

Most applications of FDS involve fires within buildings. However, it can be used to model thermal plumes

in the open and wind impinging on the exterior of a building. Rehm, McGrattan, Baum and Simiu [120]

use the LES solver to estimate surface pressures on simple rectangular blocks in a crosswind, and compare

75
these estimates to experimental measurements. In a subsequent paper [121], they consider the qualitative

effects of multiple buildings and trees on a wind field.

A different approach to wind is taken by Wang and Joulain [122]. They consider a small fire in a

wind tunnel 0.4 m wide and 0.7 m tall with flow speeds of 0.5 m/s to 2.5 m/s. Much of the comparison

with experiment is qualitative, including flame shape, lean, length. They also use the model to determine

the predominant modes of heat transfer for different operating conditions. To assess the combustion, they

implement an “Eddy Break-up” combustion model [123] and compare it to the mixture fraction approach

used by FDS. The two models perform better or worse, depending on the operating conditions. Some of

the weaknesses of the mixture fraction model as implemented in FDS version 2 are addressed in subsequent

versions. The Eddy Break-up approach has not been implemented in the official version of FDS.

Growing Fires

Vettori [124] modeled two different fire growth rates in an obstructed ceiling geometry. The rectangular

compartment was 9.2 m by 5.6 m by 2.4 m with a hollow steel door to the outside that remained closed

during the tests. An open wooden stairway led to an upper floor with the same dimensions as the fire

compartment below. Wooden joists measuring 0.038 m by 0.24 m were spaced at 0.41 m intervals across

the ceiling and were supported by a single steel beam that spanned the width of the room. A rectangular

methane gas burner measuring 0.7 m by 1.0 m by 0.31 m was placed in the corner of the chamber. Slow

and fast burning fires that reached 1055 kW in 600 s and 150 s, respectively, were monitored. Four vertical

arrays of Type K thermocouples were used to measure temperatures during the tests. The FDS model used

four grid refinements and piecewise linear grid spacing for each fire growth rate (slow and fast). For the fast

growing fire, the predicted temperatures were within 20 % of the measured values and within 10 % for the

slow growing fire. In general, finer grids produced better agreement.

In a follow-up report, Vettori [125] extended his study to include sloped ceilings, with and without

obstructions. He found that the difference between predicted and measured sprinkler activation times varied

between 4 % and 26 % for all cases studied. He also noted that FDS was able to predict the first activation

of a sprinkler twice as far from the fire as another; caused presumably by the re-direction of smoke by the

beams on the ceiling.

Floyd [126, 127] validated FDS by comparing the modeling results with measurements from fire tests at

the Heiss-Dampf Reaktor (HDR) facility. The structure was originally the containment building for a nuclear

power reactor in Germany. The cylindrical structure was 20 m in diameter and 50 m in height topped by

a hemispherical dome 10 m in radius. The building was divided into eight levels. The total volume of the

building was approximately 11,000 m3. From 1984 to 1991, four fire test series were performed within

the HDR facility. The T51 test series consisted of eleven propane gas tests and three wood crib tests. To

avoid permanently damaging the test facility, a special set of test rooms were constructed, consisting of a

fire room with a narrow door, a long corridor wrapping around the reactor vessel shield wall, and a curtained

area centered beneath a maintenance hatch. The fire room walls were lined with fire brick. The doorway and

corridor walls had the same construction as the test chamber. Six gas burners were mounted in the fire room.

The fuel source was propane gas mixed with 10 % air fed at a constant rate to one of the six burners. For

comparison with the FDS model, only the fire room, hallway and curtained region was input into the model,

for a total of 450,000 grid cells. The burners were defined within FDS as separate vents with a constant inlet

velocity. Two sets of burners were created, the first set at the physical location of the burners as the source

of fuel and second set directly above the first set as a source for ambient air. The data was presented at

fifteen minutes into the fire. The results comparing the measured data and the predicted data are presented

in Table 6.1. The FDS model predicted the layer height and temperature of the space to within 10 % of the

experimental values [126].

FDS predictions of fire growth and smoke movement in large spaces were presented by Kashef [128].

76
Table 6.1: T51 Test Results [126]

Quantity Experimental Data FDS

Upper Temperature ◦C 730 808

Upper Velocity (m/s) 4.5 2.8

Lower Velocity (m/s) 2.3 1.1

Layer Height (m) 1.0 0.9

O2 Concentration (%) 11.6 7.6

The experiments were conducted at the National Research Council Canada. The tests were performed in

a compartment with dimensions of 9 m by 6 m by 5.5 m with 32 exhaust inlets and a single supply fan.

A burner generated fires ranging in size from 15 kW to 1000 kW. FDS produced good predictions of the

experimental layer temperatures and interface heights, but there was some disagreement in the shape of the

temperature profiles.

Flame Spread

FDS was evaluated to predict the heat transfer to the wall from an adjacent pool fire. The experimental results

were based on the work by Beck et al. The predicted heat flux was in agreement with the experimental

results. The temperatures are within 30 % of the measured values near the base of the wall but decrease

more rapidly than the experimental measurements. The difference between the experimental and predicted

values can be attributed to the combustion model within FDS.

The flame spread calculations from FDS were compared to the vertical flame spread over a 5 m slab of

PMMA performed by Factory Mutual Research Corporation (FMRC). The predicted flame spread rate was

within 0.3 m/s for any point in time during the analysis. The comparison at the quasi-steady burning rate

once the full slab was burning shows that FDS over-estimated the burning rate [33, 34].

A charring model was implemented in FDS by Hostikka and McGrattan [129]. The model is a simplification of work done at NIST by Ritchie et al. [43]. The charring model was first used to predict the burning

rate of a small wooden sample in the cone calorimeter. The results were more favorable for higher imposed

heat fluxes. For low imposed fluxes, the heat transfer at the edge of the sample was more pronounced, and

more difficult to model accurately. Full-scale room tests with wood paneling were modeled, but the results

were judged to be grid-dependent. This was likely a consequence of the gas phase spatial resolution, rather

than the solid phase. The authors concluded that it is difficult to predict the growth rate of a fire in a woodlined room without “tuning” the pyrolysis rate coefficients. For real wood products, it is unlikely that all

of the necessary properties can be obtained easily. Thus, grid sensitivity and uncertain material properties

make blind predictions of fire growth on real materials beyond the reach of the current version of the model.

However, the model can still be used for a qualitative assessment of fire behavior as long as the uncertainty

in the flame spread rate is recognized.

Response of Active and Passive Fire Protection

A significant validation effort for sprinkler activation and suppression was a project entitled the International Fire Sprinkler, Smoke and Heat Vent, Draft Curtain Fire Test Project organized by the National Fire

Protection Research Foundation [130]. Thirty-nine large scale fire tests were conducted at Underwriters

Laboratories in Northbrook, IL. The tests were aimed at evaluating the performance of various fire protection systems in large buildings with flat ceilings, like warehouses and “big box” retail stores. All the tests

77
were conducted under a 30 m by 30 m adjustable-height platform in a 37 m by 37 m by 15 m high test bay.

At the time, FDS had not been publicly released and was referred to as the Industrial Fire Simulator (IFS),

but it was essentially the same as FDS version 1.

For model validation of sprinkler activation, the most valuable experiments performed were a series of

heptane spray fires. With the spray burner in different locations, with and without draft curtains, with and

without vertical vents, the model made predictions of sprinkler activation and upper layer temperatures. For

all tests, the first ring of sprinklers surrounding the fire activated within 15 % of the experimental times;

within 25 % for the second ring. The gas temperatures near the ceiling were predicted to within about 15 %

of the measured values.

Most of the full-scale experiments performed during the project used a heptane spray burner to generate

controlled fires of 1 MW to 10 MW. However, 5 experiments were performed with 6 m high racks containing

the Factory Mutual Standard Plastic Commodity, or Group A Plastic. To model these fires, bench scale

experiments were performed to characterize the burning behavior of the commodity, and larger test fires

provided validation data with which to test the model predictions of the burning rate and flame spread

behavior [50, 51]. Two to four tier configurations were evaluated. For the period of time prior to application

of water, the simulated heat release rate was within 20 % of the experimental heat release rates. It should

be noted that the model was very sensitive to the thermal parameters and the numerical grid when used to

model the fire growth in the piled commodity tests.

High rack storage fires of pool chemicals were modeled by Olenick et al. [131] to determine the validity

of sprinkler activation predictions of FDS. The model was compared to full-scale fires conducted in January,

2000 at Southwest Research Institute in San Antonio, Texas. The results indicated that the model accurately

predicted sprinkler activation and the over-pressurization of the test compartment.

FDS has been used to study the behavior of a fire undergoing suppression by a water mist system. Kim

and Ryou [132, 133] compared FDS predictions to results of compartment fire tests with and without the

application of a water mist. The cooling and oxygen dilution were predicted to within about 10 % of the

measurements, but the simulations failed to predict the complete extinguishment of a hexane pool fire. The

authors suggest that this is a result of the combustion model rather than the spray or droplet model.

Another study of water mist suppression using FDS was conducted by Hume at the University of Canterbury, Christchurch, New Zealand [134]. Full-scale experiments were performed in which a fine water

mist was combined with a displacement ventilation system to protect occupants and electrical equipment in

the event of a fire. Simulations of these experiments with FDS showed qualitative agreement, but the version

of the model used in the study (version 3) was not able to predict accurately the decrease in heat release rate

of the fire.

Airflows in Fire Compartments

Cochard [135] used FDS to study the ventilation within a tunnel. He compared the model results with

a full-scale tunnel fire experiment conducted as part of the Massachusetts Highway Department Memorial

Tunnel Fire Ventilation Test Program. The test consisted of a single point supply of fresh air through a 28 m2

opening in a 135 m tunnel. The ventilation was started 2 min after the ignition of a 40 MW fire. Fifteen

temperature measurement trees were placed within the tunnel and replicated within the model. Depending

on location, the difference between predicted and measured temperature rise ranged from 10 % to 20 %.

McGrattan and Hamins [136] also applied FDS to simulate two of the Memorial Tunnel Fire Tests

as validation for the use of the model in studying an actual fire in the Howard Street Tunnel, Baltimore,

Maryland, July 2001. The experiments chosen for the comparison were unventilated. One experiment was a

20 MW fire; the other a 50 MW fire. FDS predictions of peak near-ceiling temperatures were within 50 ◦C

of the measured peak temperatures, which were 600 ◦C and 800 ◦C, respectively.

Friday studied the use of FDS in large scale mechanically ventilated spaces. The ventilated enclosure

78
was provided with air injection rates of 1 to 12 air changes per hour and a fire with heat release rates ranging

from 0.5 MW to 2 MW. The test measurements and model output were compared to assess the accuracy of

FDS [84].

Zhang et al. [93] utilized the FDS model to predict turbulence characteristics of the flow and temperature

fields due to fire in a compartment. The experimental data was acquired through tests that replicated a halfscale ISO Room Fire Test. Two cases were explored: the heat source in the center of the room and the

heat source adjacent to a wall. The heat source was a heating element with an output of 12 kW/m2and

was assumed stable after 300 s. For the first case, the predicted average velocity and temperature profiles

were found to “agree reasonably well.” Near the ceiling, the model under-predicted temperature and overpredicted velocity. The predicted intensity of the temperature fluctuation “agree[d] very well” at all points

except those directly adjacent to the burner. The turbulent heat flux was found to be larger in the region

above the heat source.

The second case also used a burner with a 12 kW/m2 heat source located at the wall. As with the first

case, the predicted mean velocities agree with the experimental results except near the ceiling. The temperatures near the ceiling were found to be over-predicted by FDS. The intensity of the velocity fluctuation

was found to “agree well” with the experimental data except near the ceiling. The predicted intensity of the

temperature fluctuation agrees “very well” with the experimental data except in the region near the middle

of the room. This might be due to the influence of the door sill. Overall, in both cases, the predicted values

agreed well with the experimental values in all regions except near the ceiling.

The ability of version 1 of FDS to accurately predict smoke detector activation was studied by D’Souza [137].

The smoke transport model within FDS was tested and compared with UL 217 test data. The second step

in this research was to further validate the model with full-scale multi-compartment fire tests. The results

indicated that FDS is capable of predicting smoke detector activation when used with smoke detector lag

correlations that correct for the time delay associated with smoke having to penetrate the detector housing.

Another study of smoke detector activation was carried out by Brammer at the University of Canterbury,

New Zealand [138]. Two fire tests from a series performed in a two-story residence were simulated, and

smoke detector activation times were predicted using three different methods. The methods consisted of

either a temperature correlation, a time-lagged function of the optical density, or a thermal device much like

a heat detector. The purpose was to identify ways to reliably predict smoke detector activation using typical

model output like temperature and smoke concentration. It was remarked that simulating the early stage of

the fire is critical to reliable prediction.

Cleary [139] also provided a comparison between FDS computed gas velocity, temperature and concentrations at various detector locations. The research concluded that multi-room fire simulations with the FDS

model can accurately predict the conditions that a sensor might experience during a real fire event. The FDS

model was able to predict the smoke and gas concentrations, heat, and flow velocities at various detector

locations to within 15 % of measurements.

Piergoirgio et al. [140] provided a qualitative analysis of FDS applied to a truck fire within a tunnel.

The goal of their analysis was to describe the spread of the toxic gases within the tunnels, to determine the

places not involved in the spreading of combustion products and to quantify the oxygen, carbon monoxide

and hydrochloric acid concentrations during the fire.

Edwards et al. [141, 142] used FDS to determine the critical air velocity for smoke reversal in a tunnel

as a function of the fire intensity, and his results compared favorably with experimental results. Such information is necessary for planning and implementation of ventilation changes during mine fire fighting and

rescue operations.

79
Table 6.2: Comparison of FDS and measured quantities in a half-scale ISO 9705 compartment [143]

Parameter VTFRL FDS

Room TC 1 191 ◦C 134 ◦C

Door TC 5 348 ◦C 365 ◦C

Room TC 5 585 ◦C 529 ◦C

Door TC 8 227 ◦C 409 ◦C

Room TC 8 625 ◦C 559 ◦C

Door Vel 1 4.4 m/s 3.3 m/s

Door O2 1 14.6 % 13.9 %

Door Vel 5 -0.54 m/s -0.75 m/s

Door O2 5 21.0 % 21.0 %

Room HF1 29 kW/m2 26 kW/m2

Door CO21 3.5 % 3.8 %

Door HF2 3.6 kW/m2 4.1 kW/m2

Door CO24 0.2 % 2.2 %

Combustion Model

Floyd et al. [143, 127] compared the radiation model of FDS version 2 with full-scale data from the Virginia

Tech Fire Research Laboratory (VTFRL). The test compartment was outfitted with equipment capable of

taking temperature, air velocity, gas concentrations, unburned hydrocarbon and heat flux measurements. The

test facility consisted of a single compartment geometrically similar to the ISO 9705 standard compartment

with dimensions of 1.2 m by 1.8 m by 1.2 m in height. The ceiling and walls were constructed of fiberboard

over a steel shell with a floor of concrete. Three baseline experiments were completed with fires ranging in

size from 90 kW to 440 kW. A sample of the test results are presented in Table 6.2.

Overall, FDS predicted the temperatures to within 15 % of the measured temperatures. The FDS velocity

measurements followed the trend of the test data but did not replicate it. The outgoing velocities were underpredicted by 30 % to 40 % and the incoming velocities were over-predicted by 40 %. FDS predicted the

heat flux gauge response to within 10 % of the measured values. The radiation model in FDS predicted the

measured fluxes to within 15 %. The radiation to and from the compartment wall was estimated well. The

mixture fraction model was also successful. For well ventilated tests, the predictions were excellent. The

quality of the FDS prediction decreased with the under-ventilated cases. The main source of error with the

model predictions in the under-ventilated cases was the over-prediction of compartment gas temperatures

and the size of the upper layer [143, 127].

Xin and Gore [144] compared FDS predictions and measurements of the spectral radiation intensities of

small fires. The fuel flow rates for methane and ethylene burners were selected so that the Froude numbers

matched that of liquid toluene pool fires. The heat release rate was 4.2 kW for the methane flame and

3.4 kW for the ethylene flame. Line of sight spectral radiation intensities were measured at six downstream

locations. The spectral radiation intensity calculations were performed by post-processing the transient

scalar distributions provided by FDS. The calculated and measured spectral radiation intensities were found

to be in “excellent” agreement for the gas radiation bands.

Zhang et al. [62] compared the experimental results of a circular methane gas burner to predictions computed by FDS. The compartment was 2.8 m by 2.8 m by 2.2 m high with natural ventilation from a standard

door. Good agreement was found for the temperature prediction at the doorway where the radiation model

was used. The FDS model predicted the temperatures at the corner of the room better than other models

80
compared by the group. It was found that, overall, FDS predicted temperatures well but the prescribed

turbulent Prandtl and Schmidt numbers play an important role in determining the accuracy of the model.

Dillon and Hamins [145] studied the use of FDS in providing data and correlations for fire investigators

to support their investigations. A paraffin wax candle was placed within a 0.61 m by 0.61 m by 0.76 m

plexi-glass enclosure. The chamber was raised 20 mm off the surface to reveal 44 uniformly spaced 6 mm

diameter holes. The holes provided oxygen to the flame without subjecting the flame to a draft. A 150 mm

hole was provided at the top of the enclosure to allow for the heat and combustion products to exit the space.

The heat flux from the candle flame was modeled with FDS. The model provides a prediction of the heat

flux of the candle at a height of 56 mm above the base of the flame with an accuracy of 5 %. The flux is

under predicted by 16 % at 76 mm above the base of the flame. The remainder of the predictions show flux

measurements were under-predicted by 15 % to 40 % of the measured values.

6.2.4 Comparison with Standard Tests

Standard fire tests are performed at various testing laboratories and universities around the world. While

most were not designed as validation tools, they nevertheless can be used as relatively simple, well characterized fire experiments.

An extensive amount of validation work with FDS version 4 has been performed by Hietaniemi, Hostikka,

and Vaari at VTT, Finland [94]. The case studies are comprised of fire experiments ranging in scale from the

cone calorimeter (ISO 5660-1, 2002) to full-scale fire tests such as the room corner test (ISO 9705, 1993).

Comparisons are also made between FDS 4 results and data obtained in the SBI (Single Burning Item)

Euro-classification test apparatus (EN 13823, 2002) as well as data obtained in two ad hoc experimental

configurations: one is similar to the room corner test but has only partial linings and the other is a space

to study fires in building cavities. In the study of upholstered furniture, the experimental configurations are

the cone and furniture calorimeters, and the ISO room. For liquid pool fires, comparison is made to data

obtained by numerous researchers. The burning materials include spruce timber, MDF (Medium Density

Fiber) board, PVC wall carpet, upholstered furniture, cables with plastic sheathing, and heptane. The cases

studied are summarized in Table 6.3.

The scope of the VTT work is considerable. Assessing the accuracy of the model must be done on

a case by case basis. In some cases, predictions of the burning rate of the material were based solely

on its fundamental properties, as in the heptane pool fire simulations. In other cases, some properties of

the material are unknown, as in the spruce timber simulations. Thus, some of the simulations are true

predictions, some are calibrations. The intent of the authors was to provide guidance to engineers using

the model as to appropriate grid sizes and material properties. In many cases, the numerical grid was made

fairly coarse to account for the fact that in practice, FDS is used to model large spaces of which the fuel may

only comprise a small fraction.

6.2.5 Comparison with Documented Fire Experience

Documented fire experience includes known behavior of materials in fires, eyewitness accounts of real fires,

observed post fire conditions, and other means. To date, several actual fires have been reconstructed using

FDS. One case study performed by NIST is documented in Ref. [146]. Two fire fighters were killed and one

severely injured in a townhouse fire in Washington, D.C. during the evening of May 30, 1999. Questions

arose about the injuries the fire fighters had sustained, the lack of thermal damage in the living room where

a fallen fire fighter was found and why the fire fighters never opened their hose lines to protect themselves

or to extinguish the fire.

To answer some of the questions, a rectangular volume of 10 m by 6 m by 5.1 m was divided into 76,500

cells in the FDS model. The FDS results that best replicated the observed fire behavior indicated that the

81
Table 6.3: Summary of materials and test methods considered in VTT study [94]

Burning material Experimental Configuration

10 mm thick spruce timber Cone calorimeter

SBI test

Room corner test

22 mm thick spruce timber Modified room corner test

6 m by 1.1 m by 0.5 m cavity

12 mm thick MDF board Cone calorimeter

SBI test

Room corner test

PVC wall carpet on Cone calorimeter

gypsum plasterboard SBI test

Room corner test

Upholstered furniture: Cone calorimeter

a chair with PU padding Furniture calorimeter

and PP fabric ISO room test

Cables with plastic 6 m by 1.2 m by 0.6 m cavity

sheathing lined with non-combustible board

Heptane Pool fires of various sizes

opening of the basement sliding glass door provided oxygen to a pre-heated, under-ventilated fire. Flashover

was estimated to occur in less than 60 s following the entry of fire fighters into the basement. The resulting

fire gases flowed up the basement stairs and moved across the living room ceiling towards the back wall of

the townhouse. These hot gases came in direct contact with the fire fighters who were killed. The hot gases

traversed the townhouse in less than 2 s, giving the fire fighters little time to respond. The model showed

that the oxygen level was too low to support flaming and, therefore, the fire fighters did not have a visual

cue of the thermal conditions until it was too late. Results of the FDS study were shared with the D.C. fire

department and have been made available via a multi-media CD-ROM to other fire departments across the

country.

Another case study performed at NIST involved a fire in a Houston restaurant [147]. On the morning

of February 14, 2000, a fire started in the office area of a fast food restaurant. Two fire fighters died when

the roof collapsed. The FDS model was used to simulate the fire. The fuel was assumed to be the contents

of a typical office, and the fire was assumed to have a slowly growing heat release rate peaking at 6 MW.

Multiple vents were modeled and the time at which they opened replicated the fire fighters’ actions after

arrival. The model provided a visual representation of the fire during the initial phases until the collapse of

the roof.

NIST also performed a case study on a fire that killed three children and three fire fighters on the morning

of December 22, 1999 [148]. The fire started on top of a stove in a two-story residence. FDS was used to

simulate the fire. The fuel packages consisted of several furniture items in the kitchen and living room with

heat release rates reaching 5.2 MW. The model results indicated the critical event in the fire was flashover

of the kitchen. The fire became a multi-room event after flashover with temperatures increasing to over

600 ◦C. The hot gases spread quickly from the living room to the stairway on the second floor trapping the

fire fighters.

Outside of NIST, FDS has been used to investigate many actual fires, but very few of these studies

are documented in the literature. Exceptions include a study by Rein et al. [149] looking at several fire

82
events using an analytical fire growth model, the NIST zone model CFAST, and FDS. A similar study was

performed several years earlier by Spearpoint et al. [150] as a class exercise at the University of Maryland.

During the SFPE Professional Development Week in the fall of 2001, a workshop was held in which several

engineers related their experiences using FDS as a forensic tool [151]. The role of carbon monoxide in the

deaths of three fire fighters was studied by Christensen and Icove [152]. There is little quantitative validation

of the model afforded by these studies. However, the degree to which the model is able to reproduce observed

behavior can be used as an indicator of the model’s strengths and weaknesses.

83
84
Chapter 7

Conclusion

The equations and numerical algorithm described in this document form the core of an evolving fire model.

As research into specific fire-related phenomena continues, the relevant parts of the model can be improved.

Because the model was originally designed to predict the transport of heat and exhaust products from fires,

it can be used reliably when the fire is prescribed and the numerical grid is sufficiently resolved to capture

enough of the flow structure for the application at hand. It is the job of the user to determine what level

of accuracy is needed. For simulations that simply involve mass and heat transport, the model predicts

flow velocities and temperatures to an accuracy of 5 % to 20% compared to experimental measurements,

depending on the fidelity of the underlying grid. Current research efforts focus on improving both the gas

phase and solid phase descriptions of combustion so that phenomena such as flame spread, fire growth,

under-ventilation, and suppression can be handled more reliably.

Any user of the numerical model must be aware of the assumptions and approximations being employed.

There are two issues for any potential user to consider before embarking on calculations. First, for both real

and simulated fires, the growth of the fire is very sensitive to the thermal properties (conductivity, specific

heat, density, burning rate, etc.) of the surrounding materials. Second, even if all the material properties are

known, the physical phenomena of interest may not be simulated due to limitations in the model algorithms

or numerical grid. Except for those few materials that have been studied to date at NIST, the user must supply

the thermal properties of the materials, and then validate the performance of the model with experiments to

ensure that the model has the necessary physics included. Only then can the model be expected to predict

the outcome of fire scenarios that are similar to those that have actually been tested.

85
86
Nomenclature

As water droplet surface area

B pre-exponential factor for Arrhenius reaction

C Sprinkler C-Factor

CD drag coefficient

Cs Smagorinsky constant (LES)

cp constant pressure specific heat

D diffusion coefficient

D

∗

characteristic fire diameter

dm median volumetric droplet diameter

E activation energy

f external force vector (excluding gravity)

g acceleration of gravity

H total pressure divided by the density

h enthalpy; heat transfer coefficient

hi enthalpy of ith species

h

0

i

heat of formation of ith species

I radiation intensity

Ib radiation blackbody intensity

k thermal conductivity; suppression decay factor

M molecular weight of the gas mixture

Mi molecular weight of ith gas species

m ̇

00

f

fuel mass flux

m ̇

000

i mass production rate of ith species per unit volume

m ̇

00

w water mass flux

m

00

w water mass per unit area

m ̇

00

O

oxygen consumption rate per unit area

Nu Nusselt number

Pr Prandtl number

p pressure

p0 background pressure

p ̃ pressure perturbation

qr radiative heat flux vector

q ̇

000 heat release rate per unit volume

q ̇

00

r

radiative flux to a solid surface

q ̇

00

c

convective flux to a solid surface

Q ̇ total heat release rate

Q

∗

characteristic fire size

R universal gas constant

87
Re Reynolds number

rd water droplet radius

RTI Response Time Index of sprinkler

s unit vector in direction of radiation intensity

Sc Schmidt number

Sh Sherwood number

T temperature

t time

tb thermal element burn-out time (LES)

U integrated radiant intensity

u = (u,v,w) velocity vector

m ̇

000

i

production rate of ith species per unit volume

We Weber number

x = (x,y,z) position vector

Xi volume fraction of ith species

Yi mass fraction of ith species

Y

∞

O mass fraction of oxygen in the ambient

Y

I

F mass fraction of fuel in the fuel stream

ys soot yield

Z mixture fraction

γ ratio of specific heats; Rosin-Rammler exponent

∆H heat of combustion

∆HO energy released per unit mass oxygen consumed

δ wall thickness

κ absorption coefficient

μ dynamic viscosity

νi stoichiometric coefficient, species i

Φ dissipation function

ρ density

τi j viscous stress tensor

χr radiative loss fraction

σ Stefan-Boltzmann constant

σd droplet scattering coefficient

σs scattering coefficient

ω = (ωx,ωy,ωz) vorticity vector

88
Bibliography

[1] S. Olenick and D. Carpenter. An Updated International Survey of Computer Models for Fire and

Smoke. Journal of Fire Protection Engineering, 13:87–110, May 2003.

[2] R. Friedman. Survey of Computer Models for Fire and Smoke. Technical report, Factory Mutual

Research Corporation, Norwood, Massachusetts, March 1990.

[3] American Society for Testing and Materials, West Conshohocken, Pennsylvania. ASTM E 1355-04,

Standard Guide for Evaluating the Predictive Capabilities of Deterministic Fire Models, 2004.

[4] H.C. Hottel. Stimulation of Fire Research in the United States After 1940 (A Historical Account).

Combustion Science and Technology, 39:1–10, 1984.

[5] W.W. Jones. A Review of Compartment Fire Models. NBSIR 83-2684, National Bureau of Standards

(now NIST), Gaithersburg, Maryland, 1983.

[6] J. Quintiere. A Perspective on Compartment Fire Growth. Combustion Science and Technology,

39:11–54, 1984.

[7] G.P. Forney and W.F. Moss. Analyzing and Exploiting Numerical Characteristics of Zone Fire Models. Fire Science and Technology, 14:49–60, 1994.

[8] S.V. Patankar. Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing, New York, 1980.

[9] R.G. Rehm and H.R. Baum. The Equations of Motion for Thermally Driven, Buoyant Flows. Journal

of Research of the NBS, 83:297–308, 1978.

[10] E.S. Oran and J.P. Boris. Numerical Simulation of Reactive Flow. Elsevier Science Publishing

Company, New York, 1987.

[11] K.B. McGrattan and G.P. Forney. Fire Dynamics Simulator (Version 4), User’s Guide. NIST Special Publication 1019, National Institute of Standards and Technology, Gaithersburg, Maryland, July

2004.

[12] G.P. Forney and K.B. McGrattan. User’s Guide for Smokeview Version 4. NIST Special Publication

1017, National Institute of Standards and Technology, Gaithersburg, Maryland, July 2004.

[13] D.A. Anderson, J.C. Tannehill, and R.H. Pletcher. Computational Fluid Mechanics and Heat Transfer. Hemisphere Publishing Corporation, Philadelphia, Pennsylvania, 1984.

[14] R. Peyret and T.D. Taylor. Computational Methods for Fluid Flow. Springer-Verlag, New York,

1983.

89
[15] J.H. Ferziger and M. Peric. ́ Computational Methods for Fluid Dynamics. Springer-Verlag, Berlin,

2nd edition, 1999.

[16] R.W. Bilger. Turbulent Diffusion Flames. Annual Review of Fluid Mechanics, 21:101–135, 1989.

[17] J.P. Holman. Heat Transfer. McGraw-Hill, New York, 5th edition, 1989.

[18] F.P. Incropera and D.P. De Witt. Fundamentals of Heat and Mass Transfer. John Wiley and Sons,

New York, 4th edition, 1996.

[19] R. Siegel and J. R. Howell. Thermal Radiation Heat Transfer. Taylor & Francis, New York, 4th

edition, 2002.

[20] P.J. DiNenno, editor. SFPE Handbook of Fire Protection Engineering. National Fire Protection

Association, Quincy, Massachusetts, 3rd edition, 2002.

[21] D. Drysdale. An Introduction to Fire Dynamics. John Wiley and Sons, New York, 2nd edition, 2002.

[22] J.G. Quintiere. Principles of Fire Behavior. Delmar Publishers, Albany, New York, 1998.

[23] V. Babrauskas. Ignition Handbook. Fire Science Publishers, Issaquah, Washington USA, 1st edition,

2003. Co-published by the Society of Fire Protection Engineers.

[24] S.B. Pope. Ten questions concerning the large-eddy simulation of turbulent flows. New Journal of

Physics, 6:1–24, 2004.

[25] J. Smagorinsky. General Circulation Experiments with the Primitive Equations. I. The Basic Experiment. Monthly Weather Review, 91(3):99–164, March 1963.

[26] J.W. Deardorff. Numerical Investigation of Neutral and Unstable Planetary Boundary Layers. Journal

of Atmospheric Sciences, 29:91–115, 1972.

[27] M. Germano, U. Piomelli, P. Moin, and W.H. Cabot. A Dynamic Subgrid-Scale Eddy Viscosity

Model. Physics of Fluids A, 3(7):1760–1765, 1991.

[28] D.K. Lilly. A Proposed Modification of the Germano Subgrid-Scale Closure Method. Physics of

Fluids A, 4(3):633–635, 1992.

[29] H.R. Baum, K.B. McGrattan, and R.G. Rehm. Three Dimensional Simulations of Fire Plume Dynamics. Journal of the Heat Transfer Society of Japan, 35(139):45–52, 1997.

[30] B.E. Poling, J.M. Prausnitz, and J.P. O’Connell. The Properties of Gases and Liquids. McGraw-Hill,

New York, 5th edition, 2000.

[31] H.R. Baum and K.B. McGrattan. Simulation of Large Industrial Outdoor Fires. In Fire Safety Science – Proceedings of the Sixth International Symposium. International Association for Fire Safety

Science, 2000.

[32] C. Huggett. Estimation of Rate of Heat Release by Means of Oxygen Consumption Measurements.

Fire and Materials, 4(2):61–65, 1980.

[33] K.M. Liang, T. Ma, J.G. Quintiere, and D. Rouson. Application of CFD Modeling to Room Fire

Growth on Walls. NIST GCR 03-849, National Institute of Standards and Technology, Gaithersburg,

Maryland, April 2003.

90
[34] T. Ma and J. Quintiere. Numerical Simulation of Axi-Symmetric Fire Plumes: Accuracy and Limitations. Fire Safety Journal, 38:467–492, 2003.

[35] G. Heskestad. SFPE Handbook of Fire Protection Engineering, chapter Fire Plumes, Flame Height

and Air Entrainment. National Fire Protection Association, Quincy, Massachusetts, 3rd edition, 2002.

[36] C. Beyler. SFPE Handbook of Fire Protection Engineering, chapter Flammability Limits of Premixed

and Diffusion Flames. National Fire Protection Association, Quincy, Massachusetts, 3rd edition,

2002.

[37] I.K. Puri and K. Seshadri. Extinction of Diffusion Flames Burning Diluted Methane and Diluted

Propane in Diluted Air. Combustion and Flame, 65:137–150, 1986.

[38] C.K. Westbrook and F.L. Dryer. Simplified Reaction Mechanisms for the Oxidation of Hydrocarbon

Fuels in Flames. Combustion Science and Technology, 27:31–43, 1981.

[39] W. Grosshandler. RadCal: A Narrow Band Model for Radiation Calculations in a Combustion Environment. NIST Technical Note TN 1402, National Institute of Standards and Technology, Gaithersburg, Maryland, 1993.

[40] G.D. Raithby and E.H. Chui. A Finite-Volume Method for Predicting Radiant Heat Transfer in Enclosures with Participating Media. Journal of Heat Transfer, 112(2):415–423, 1990.

[41] K. Prasad, C. Li, K. Kailasanath, C. Ndubizu, R. Ananth, and P.A. Tatem. Numerical modelling of

methanol liquid pool fires. Combustion Theory and Modelling, 3:743–768, 1999.

[42] A. Atreya. Pyrolysis, Ignition and Fire Spread on Horizontal Surfaces of Wood. NBS GCR 83-449,

National Bureau of Standards (now NIST), Gaithersburg, Maryland, 1983.

[43] S.J. Ritchie, K.D. Steckler, A. Hamins, T.G. Cleary, J.C. Yang, and T. Kashiwagi. The Effect of

Sample Size on the Heat Release Rate of Charring Materials. In Fire Safety Science - Proceedings of

the 5th International Symposium, pages 177–188. International Association For Fire Safety Science,

1997.

[44] G. Heskestad and R.G. Bill. Quantification of Thermal Responsiveness of Automatic Sprinklers

Including Conduction Effects. Fire Safety Journal, 14:113–125, 1988.

[45] P. Ruffino and M. di Marzo. Temperature and Volumetric Fraction Measurements in a Hot Gas Laden

with Water Droplets. Journal of Heat Transfer, 125(2):356–364, April 2003.

[46] P. Ruffino and M. di Marzo. The Effect of Evaporative Cooling on the Activation Time of Fire

Sprinklers. In Fire Safety Science – Proceedings of the Seventh International Symposium, pages

481–492. International Association for Fire Safety Science, 2002.

[47] F. Gavelli, P. Ruffino, G. Anderson, and M. di Marzo. Effect of Minute Water Droplets on a Simulated Sprinkler Link Thermal Response. NIST GCR 99-776, National Institute of Standards and

Technology, Gaithersburg, Maryland, July 1999.

[48] T.S. Chan. Measurements of Water Density and Droplet Size Distributions of Selected ESFR Sprinklers. Journal of Fire Protection Engineering, 6(2):79–87, 1994.

[49] H.Z. Yu. Investigation of Spray Patterns of Selected Sprinklers with the FMRC Drop Size Measuring

System. In Fire Safety Science – Proceedings of the First International Symposium, pages 1165–1176.

International Association For Fire Safety Science, 1986.

91
[50] A. Hamins and K.B. McGrattan. Reduced-Scale Experiments to Characterize the Suppression of

Rack Storage Commodity Fires. NISTIR 6439, National Institute of Standards and Technology,

Gaithersburg, Maryland, 1999.

[51] A. Hamins and K.B. McGrattan. Reduced-Scale Experiments on the Water Suppression of a RackStorage Commodity Fire for Calibration of a CFD Fire Model. In Fire Safety Science – Proceedings

of the Seventh International Symposium, pages 457–468. International Association for Fire Safety

Science, 2002.

[52] N. Cheremisinoff. Encyclopedia of Fluid Mechanics, Volume 3: Gas-Liquid Flows. Gulf Publishing

Company, Houston, Texas, 1986.

[53] T. Ravigururajan and M. Beltran. A Model for Attenuation of Fire Radiation Through Water Droplets.

Fire Safety Journal, 15:171–181, 1989.

[54] A. Tuntomo, C. Tien, and S. Park. Optical Constants of Liquid Hydrocarbon Fuels. Combustion

Science and Technology, 84:133–140, 1992.

[55] W.J. Wiscombe. Improved Mie Scattering Algorithms. Applied Optics, 19(9):1505–1509, 1980.

[56] G.M. Hale and M.R. Querry. Optical constants of water in the 200 nm to 200 mum wavelength region.

Applied Optics, 12(3):555–563, 1973.

[57] H.Z. Yu, J.L. Lee, and H.C. Kung. Suppression of Rack-Storage Fires by Water. In Fire Safety Science

– Proceedings of the Fourth International Symposium, pages 901–912. International Association For

Fire Safety Science, 1994.

[58] T. Cleary, A. Chernovsky, W. Grosshandler, and M. Anderson. Particulate Entry Lag in Spot-Type

Smoke Detectors. In Fire Safety Science – Proceedings of the Sixth International Symposium, pages

779–790. International Association for Fire Safety Science, 1999.

[59] G.W. Mulholland and C. Croarkin. Specific Extinction Coefficient of Flame Generated Smoke. Fire

and Materials, 24:227–230, 2000.

[60] K.B. McGrattan, R.G. Rehm, and H.R. Baum. Fire-Driven Flows in Enclosures. Journal of Computational Physics, 110(2):285–291, 1994.

[61] H.R. Baum, O.A. Ezekoye, K.B. McGrattan, and R.G. Rehm. Mathematical Modeling and Computer

Simulation of Fire Phenomenon. Theoretical and Computational Fluid Dynamics, 6:125–139, 1994.

[62] W. Zhang, N. Ryder, R.J. Roby, and D. Carpenter. Modeling of the Combustion in Compartment

Fires Using Large Eddy Simulation Approach. In Proceedings of the 2001 Fall Technical Meeting,

Eastern States Section. Combustion Institute, Pittsburgh, Pennsylvania, December 2001.

[63] G. Continillo, F.M Denaro, and F.S. Marra. Accuracy and Stability Analysis of Time-Integrated

Schemes for Advection-Diffusion-Reaction Equations. In Seventh International Conference on Numerical Combustion, page 99, 1998.

[64] J.Y. Murthy and S.R. Mathur. Radiative Heat Transfer in Axisymmetric Geometries Using an Unstructured Finite-Volume Method. Numerical Heat Transfer, Part B, 33:397–416, 1998.

[65] S.H. Kim and K.Y. Huh. Assessment of The Finite-Volume Method and The Discrete Ordinate

Method for Radiative Heat Transfer In a Three-Dimensional Rectangular Enclosure. Numerical Heat

Transfer, Part B, 35:85–112, 1999.

92
[66] W. Yang, T. Parker, H.D. Ladouceur, and R.J. Kee. The Interaction of Thermal Radiation and Water

Mist in Fire Suppression. Fire Safety Journal, 39:41–66, 2004.

[67] R.A. Sweet. Direct Methods for the Solution of Poisson’s Equation on a Staggered Grid. Journal of

Computational Physics, 12:422–428, 1973.

[68] W. Mell, K.B. McGrattan, and H. Baum. Numerical Simulation of Combustion in Fire Plumes. In

Twenty-Sixth Symposium (International) on Combustion, pages 1523–1530. Combustion Institute,

Pittsburgh, Pennsylvania, 1996.

[69] K.B. McGrattan, H.R. Baum, and R.G. Rehm. Large Eddy Simulations of Smoke Movement. Fire

Safety Journal, 30:161–178, 1998.

[70] H.R. Baum, R.G. Rehm, P.D. Barnett, and D.M. Corley. Finite Difference Calculations of Buoyant

Convection in an Enclosure, Part I: The Basic Algorithm. SIAM Journal of Scientific and Statistical

Computing, 4(1):117–135, March 1983.

[71] H.R. Baum and R.G. Rehm. Finite Difference Solutions for Internal Waves in Enclosures. SIAM

Journal of Scientific and Statistical Computing, 5(4):958–977, December 1984.

[72] H.R. Baum and R.G. Rehm. Calculations of Three Dimensional Buoyant Plumes in Enclosures.

Combustion Science and Technology, 40:55–77, 1984.

[73] R.G. Rehm, P.D. Barnett, H.R. Baum, and D.M. Corley. Finite Difference Calculations of Buoyant

Convection in an Enclosure: Verification of the Nonlinear Algorithm. Applied Numerical Mathematics, 1:515–529, 1985.

[74] J.C. Adams, W.S. Brainerd, J.T. Martin, B.T. Smith, and J.L. Wagener. Fortran 95 Handbook: Complete ISO/ANSI Reference. MIT Press, Cambridge, Massachusetts, 1997.

[75] K.B. McGrattan, T. Kashiwagi, H.R. Baum, and S.L. Olson. Effects of Ignition and Wind on the

Transition to Flame Spread in a Microgravity Environment. Combustion and Flame, 106:377–391,

1996.

[76] T. Kashiwagi, K.B. McGrattan, S.L. Olson, O. Fujita, M. Kikuchi, and K. Ito. Effects of Slow

Wind on Localized Radiative Ignition and Transition to Flame Spread in Microgravity. In TwentySixth Symposium (International) on Combustion, pages 1345–1352. Combustion Institute, Pittsburgh,

Pennsylvania, 1996.

[77] W. Mell and T. Kashiwagi. Dimensional Effects on the Transition from Ignition to Flame Spread

in Microgravity. In Twenty-Seventh Symposium (International) on Combustion, pages 2635–2641.

Combustion Institute, Pittsburgh, Pennsylvania, 1998.

[78] W. Mell, S.L. Olson, and T. Kashiwagi. Flame Spread Along Free Edges of Thermally-Thin Samples

in Microgravity. In Twenty-Eighth Symposium (International) on Combustion, pages 2843–2849.

Combustion Institute, Pittsburgh, Pennsylvania, 2000.

[79] K. Prasad, Y. Nakamura, S.L. Olson, O. Fujita, K. Nishizawa, K. Ito, and T. Kashiwagi. Effect of

Wind Velocity on Flame Spread in Microgravity. In Twenty-Ninth Symposium (International) on

Combustion, pages 2553–2560. Combustion Institute, Pittsburgh, Pennsylvania, 2002.

[80] Y. Nakamura, T. Kashiwagi, K.B. McGrattan, and H.R. Baum. Enclosure Effects on Flame Spread

over Solid Fuels in Microgravity. Combustion and Flame, 130:307–321, 2002.

93
[81] W.E. Mell, K.B. McGrattan, and H.R. Baum. g-Jitter Effects on Spherical Diffusion Flames. Microgravity Science and Technology, 15(4):12–30, 2004.

[82] A. Mukhopadhyay and I.K. Puri. An Assessment of Stretch Effects on Flame Tip Using the Thin

Flame and Thick Formulations. Combustion and Flame, 133:499–502, 2003.

[83] A. Hamins, M. Bundy, I.K. Puri, K.B. McGrattan, and W.C. Park. Suppression of Low Strain Rate

Non-Premixed Flames by an Agent. In Proceedings of the 6th International Microgravity Combustion

Workshop, NASA/CP-2001-210826, pages 101–104. National Aeronautics and Space Administration,

Lewis Research Center, Cleveland, Ohio, May 2001.

[84] P. Friday and F. W. Mowrer. Comparison of FDS Model Predictions with FM/SNL Fire Test Data.

NIST GCR 01-810, National Institute of Standards and Technology, Gaithersburg, Maryland, April

2001.

[85] A. Bounagui, N. Benichou, C. McCartney, and A. Kashef. Optimizing the Grid Size Used in CFD

Simulations to Evaluate Fire Safety in Houses. In 3rd NRC Symposium on Computational Fluid

Dynamics, High Performance Computing and Virtual Reality, pages 1–8, Ottawa, Ontario, Canada,

December 2003. National Research Council, Canada.

[86] R.L. Alpert. SFPE Handbook of Fire Protection Engineering, chapter Ceiling Jet Flows. National

Fire Protection Association, Quincy, Massachusetts, 3rd edition, 2003.

[87] A. Bounagui, A. Kashef, and N. Benichou. Simulation of the Dynamics of the Fire for a Section of

the L.H.-La Fontaine Tunnel. IRC-RR- 140, National Research Council Canada, Ottawa, Canada,

K1A0R, September 2003.

[88] Y. Xin. Assessment of Fire Dynamics Simulation for Engineering Applications: Grid and Domain

Size Effects. In Proceedings of the Fire Suppression and Detection Research Application Symposium,

Orlando, Florida. National Fire Protection Association, Quincy, Massachusetts, 2004.

[89] J.A. Ierardi and J.R. Barnett. A Quantititive Method for Calibrating CFD Model Calculations. In

Proceedings of the CIB-CTBUH International Conference on Tall Buildings, pages 507–514. International Council for Research and Innovation in Building and Construction (CIB), 2003.

[90] N.M. Petterson. Assessing the feasibility of reducing the grid resolution in fds field modeling. Fire

Engineering Research Report 2002/6, University of Canterbury, Christchurch, New Zealand, March

2002.

[91] A. Musser, K. B. McGrattan, and J. Palmer. Evaluation of a Fast, Simplified Computational Fluid

Dynamics Model for Solving Room Airflow Problems. NISTIR 6760, National Institute of Standards

and Technology, Gaithersburg, Maryland, June 2001.

[92] K.B. McGrattan, J.E. Floyd, G.P. Forney, H.R. Baum, and S. Hostikka. Improved Radiation and

Combustion Routines for a Large Eddy Simulation Fire Model. In Fire Safety Science – Proceedings

of the Seventh International Symposium, pages 827–838. International Association for Fire Safety

Science, 2002.

[93] W. Zhang, A. Hamer, M. Klassen, D. Carpenter, and R. Roby. Turbulence Statistics in a Fire Room

Model by Large Eddy Simulation. Fire Safety Journal, 37:721–752, 2002.

94
[94] J. Hietaniemi, S. Hostikka, and J. Vaari. FDS Simulation of Fire Spread – Comparison of Model Results with Experimental Data. VTT Working Paper 4, VTT Building and Transport, Espoo, Finland,

2004.

[95] R.G. Rehm, H.R. Baum, D.W. Lozier, H. Tang, and J. Sims. Buoyant Convection in an Inclined

Enclosure. In Fire Safety Science – Proceedings of the Third International Symposium, pages 313–

323. International Association for Fire Safety Science, 1991.

[96] R.G. Rehm, K.B. McGrattan, H.R. Baum, and K.W. Cassel. Transport by Gravity Currents in Building Fires. In Fire Safety Science – Proceedings of the Fifth International Symposium, pages 391–402.

International Association for Fire Safety Science, 1997.

[97] H.R. Baum, K.B. McGrattan, and R.G. Rehm. Three Dimensional Simulations of Fire Plume Dynamics. In Fire Safety Science – Proceedings of the Fifth International Symposium, pages 511–522.

International Association for Fire Safety Science, 1997.

[98] H.R. Baum, K.B. McGrattan, and R.G. Rehm. Large Eddy Simulations of Smoke Movement in

Three Dimensions. In Proceedings of the Seventh International Interflam Conference, pages 189–

198. Interscience Communications, London, 1996.

[99] H.R. Baum and B.J. McCaffrey. Fire Induced Flow Field – Theory and Experiment. In Fire Safety

Science – Proceedings of the Second International Symposium, pages 129–148. International Association for Fire Safety Science, 1989.

[100] J.M. Clement and C.M. Fleischmann. Experimental Verification of the Fire Dynamics Simulator

Hydrodynamic Model. In Fire Safety Science – Proceedings of the Seventh International Symposium,

pages 839–862. International Association for Fire Safety Science, 2002.

[101] T. Cleary, M. Anderson, J. Averill, and W. Grosshandler. Evaluating Multi-Sensor Fire Detectors in

the Fire Emulator/Detector Evaluator. In Proceedings of the Eighth International Interflam Conference, pages 453–464. Interscience Communications, June 1999.

[102] W. Davis, K Notarianni, and K. McGrattan. Comparison of Fire Model Predictions with Experiments

Conducted in a Hangar with a 15 Meter Ceiling. NISTIR 5927, National Institute of Standards and

Technology, Gaithersburg, Maryland, December 1996.

[103] S.J. Emmerich and K.B. McGrattan. Application of a Large Eddy Simulation Model to Study Room

Airflow. ASHRAE Transactions, 104(1):1–9, 1998.

[104] S.J. Emmerich. Use of Computational Fluid Dynamics to Analyze Indoor Air Quality Issues. NISTIR

5997, National Institute of Standards and Technology, Gaithersburg, Maryland, April 1997.

[105] K.B. McGrattan, H.R. Baum, and R.G. Rehm. Numerical Simulation of Smoke Plumes from Large

Oil Fires. Atmospheric Environment, 30(24):4125–4136, 1996.

[106] T. Yamada. Smoke Plume Trajectory from In-Situ Burning of Crude Oil in Tomakomai. Technical

report, National Research Institute of Fire and Disaster, Japan, November 1998.

[107] A. Hamins, A. Maranghides, K.B. McGrattan, E. Johnsson, T. Ohlemiller, M. Donnelly, J. Yang,

G. Mulholland, K. Prasad, S. Kukuck, R. Anleitner, and T. McAllister. Report on Experiments to

Validate Fire Dynamic and Thermal-Structural Models for Use in the World Trade Center Investigation. NIST Special Publication 1000-B, National Institute of Standards and Technology, Gaithersburg,

Maryland, September 2004.

95
[108] W.D. Walton and P.H. Thomas. SFPE Handbook of Fire Protection Engineering, chapter Estimating

Temperatures in Compartment Fires. National Fire Protection Association, Quincy, Massachusetts,

3rd edition, 2002.

[109] A. Hamins, A. Maranghides, K.B. McGrattan, Ohlemiller, and R. Anletiner. Experiments to Validate

Models of Fire Growth and Spread for use in the World Trade Center Investigation. NIST Special Publication 1000-10E, National Institute of Standards and Technology, Gaithersburg, Maryland,

September 2004.

[110] W.K. Chow and R. Yin. Discussion on Two Plume Formulae with Computational Fluid Dynamics.

Journal of Fire Sciences, 20:179–201, May 2002.

[111] F. Battaglia, K. McGrattan, R. Rehm, and H. Baum. Simulating Fire Whirls. Combustion Theory and

Modeling, 4:123–138, 2000.

[112] Y. Xin. Baroclinic Effects on Fire Flow Field. In Proceedings of the Fourth Joint Meeting of the U.S.

Sections of the Combustion Institute. Combustion Institute, Pittsburgh, Pennsylvania, March 2005.

[113] Y. Xin, J.P. Gore, K.B. McGrattan, R.G. Rehm, and H.R. Baum. Fire dynamics simulation of a

turbulent buoyant flame using a mixture-fraction-based combustion model. Combustion and Flame,

141:329–335, 2005.

[114] Y. Xin, J.P. Gore, K.B. McGrattan, R.G. Rehm, and H.R. Baum. Large Eddy Simulation of Buoyant

Turbulent Pool Fires. In Twenty-Ninth Symposium (International) on Combustion, pages 259–266.

Combustion Institute, Pittsburgh, Pennsylvania, 2002.

[115] S. Hostikka, K.B. McGrattan, and A. Hamins. Numerical Modeling of Pool Fires using Large Eddy

Simulation and Finite Volume Method for Radiation. In Fire Safety Science – Proceedings of the

Seventh International Symposium, pages 383–394. International Association for Fire Safety Science,

2002.

[116] X. Yuan, Q. Chen, L.R. Glicksman, Y. Hu, and X. Yang. Measurements and Computations of Room

Airflow with Displacement Ventilation. Technical Report RP-949, American Society of Heating,

Refrigerating, Air-Conditioning Engineers, 1999.

[117] A. Musser and L. Tan. Control of Diesel Exhaust Fumes in Enclosed Locomotive Facilities. Technical

Report RP-1191, American Society of Heating, Refrigerating, Air-Conditioning Engineers, 2003.

[118] K. Mniszewski. The Use of FDS for Estimation of Flammable Gas/Vapor Concentrations. In Proceedings of the 3rd Technical Symposium on Computer Applications in Fire Protection Engineering,

pages 143–155. Society for Fire Protection Engineers, Bethesda, Maryland, September 2001.

[119] S. Kerber and W. Walton. Characterizing Positive Pressure Ventilation using Computational Fluid

Dynamics. NISTIR 7065, National Institute of Standards and Technology, Gaithersburg, Maryland,

February 2003.

[120] R. Rehm, K. McGrattan, H. Baum, and E. Simiu. An Efficient Large Eddy Simulation Algorithm

for Computational Wind Engineering: Application to Surface Pressure Computations on a Single

Building. NISTIR 6371, National Institute of Standards and Technology, August 1999.

[121] R.G. Rehm, K.B. McGrattan, and H.R. Baum. Large Eddy Simulation of Flow over a Wooded

Building Complex. Wind & Structures, 5(2):291–300, 2002.

96
[122] H.Y. Wang and P. Joulain. Numerical Simulation of Wind-Aided Turbulent Fires in a Ventilated

Model Tunnel. In Fire Safety Science – Proceedings of the Seventh International Symposium, pages

161–172. International Association for Fire Safety Science, 2002.

[123] F. Magnussen and B.H. Hjertager. On Mathematical Modeling of Turbulent Combustion with Special

Emphasis on Soot Formation and Combustion. In Proceedings of the Sixteenth Symposium (International) on Combustion, pages 719–729. Combustion Institute, Pittsburgh, Pennsylvania, 1977.

[124] R. Vettori. Effect of an Obstructed Ceiling on the Activation Time of a Residential Sprinkler. NISTIR

6253, National Institute of Standards and Technology, Gaithersburg, Maryland, November 1998.

[125] R. Vettori. Effect of a Beamed, Sloped, and Sloped Beamed Ceilings on the Activation Time of a

Residential Sprinkler. NISTIR 7079, National Institute of Standards and Technology, Gaithersburg,

Maryland, December 2003.

[126] J. Floyd. Comparison of CFAST and FDS for Fire Simulation with the HDR T51 and T52 Tests.

NISTIR 6866, National Institute of Standards and Technology, Gaithersburg, Maryland, March 2002.

[127] J.E. Floyd, K.B. McGrattan, S. Hostikka, and H.R. Baum. CFD Fire Simulation Using Mixture Fraction Combustion and Finite Volume Radiative Heat Transfer. Journal of Fire Protection Engineering,

13:11–26, February 2003.

[128] A. Kashef, N. Benichou, G.D. Lougheed, and C. McCartney. A Computational and Experimental

Study of Fire Growth and Smoke Movement in Large Spaces. Technical Report NRCC-45201, National Research Council Canada, 2002.

[129] S. Hostikka and K.B. McGrattan. Large Eddy Simulations of Wood Combustion. In Proceedings

of the Ninth International Interflam Conference, pages 755–762. Interscience Communications, London, 2001.

[130] K.B. McGrattan, A. Hamins, and D. Stroup. Sprinkler, Smoke & Heat Vent, Draft Curtain Interaction — Large Scale Experiments and Model Development. NISTIR 6196-1, National Institute of

Standards and Technology, Gaithersburg, Maryland, September 1998.

[131] S. Olenick, M. Klassen, and R.J. Roby. Validation Study for FDS for a High Rack Storage Fire

Involving Pool Chemicals. In Proceedings of the 3rd Technical Symposium on Computer Applications

in Fire Protection Engineering. Society of Fire Protection Engineers, Bethesda, Maryland, September

2001.

[132] S.C. Kim and H.S. Ryou. An Experimental and Numerical Study on Fire Suppression using a Water

Mist in an Enclosure. Building and Environment, 38:1309–1316, 2003.

[133] S.C. Kim and H.S. Ryou. The Effects of Water Mist on the Compartment Fire. International Journal

of Air-Conditioning and Refrigeration, 12(1):30–36, 2004.

[134] B.P. Hume. Water Mist Suppression in Conjunction with Displacement Ventilation. Master’s thesis,

University of Canterbury, Christchurch, New Zealand, February 2003. Fire Engineering Research

Report 03/4.

[135] S. Cochard. Validation of Fire Dynamics Simulator (Version 2.0) Freeware. Tunnel Management

International Journal, 6(4), December 2003.

97
[136] K. B. McGrattan and A. Hamins. Numerical Simulation of the Howard Street Tunnel Fire, Baltimore,

Maryland, July 2001. NISTIR 6902, National Institute of Standards and Technology, Gaithersburg,

Maryland, January 2003. Joint Publication of NIST and the US Nuclear Regulatory Commission

(NUREG/CR-6793).

[137] V. D’Souza, J.A. Sutula, S.M. Olenick, W. Zhang, and R.J. Roby. Use of Fire Dynamics Simulator

to Predict Smoke Detector Activation. In Proceedings of the 2001 Fall Technical Meeting, Eastern

States Section, pages 175–178. Combustion Institute, Pittsburgh, Pennsylvania, December 2001.

[138] D.R. Brammer. A Comparison between Predicted and Actual Behaviour of Domestic Smoke Detectors in a Realistic House Fire. Master’s thesis, University of Canterbury, Christchurch, New Zealand,

2002.

[139] T. Cleary, M. Donnelly, G. Mulholland, and B. Farouk. Fire Detector Performance Predictions in

a Simulated Multi-Room Configuration. In Proceedings of the 12th International Conference on

Automatic Fire Detection (AUBE ’01). National Institute of Standards and Technology, Gaithersburg,

Maryland, March 2001. NIST SP 965.

[140] A. Piergiorgio, D. Giuseppe, F. Dino, G. Zappellini, and A. Ferrari. CFD Simulations of a Truck Fire

in the Underground Gran Sasso National Laboratory. In Proceedings of the 5th Italian Conference

on Chemical and Process Engineering, volume 5 of AIDIC Conference Series. Associazione Italiana

Di Ingegneria Chimica (AIDIC), Elsevier, May 2002. Papers presented at ICeaP-5, Florence, Italy,

May 20-23, 2001.

[141] J.C. Edwards, R.A. Franks, G.F. Friel, and L. Yuan. Experimental and Modeling Investigation of

the Effect of Ventilation on Smoke Rollback in a Mine Entry. In Proceedings of the SME Annual

Meeting, pages 1–6. Society for Mining, Metallurgy and Exploration, 2005.

[142] C.C. Hwang and J.C. Edwards. The critical ventilation velocity in tunnel fires – a computer simulation. Fire Safety Journal, 40:213–244, 2005.

[143] J.E. Floyd, C. Wieczorek, and U. Vandsburger. Simulations of the Virginia Tech Fire Research Laboratory Using Large Eddy Simulation with Mixture Fraction Chemistry and Finite Volume Radiative

Heat Transfer. In Proceedings of the Ninth International Interflam Conference. Interscience Communications, London, 2001.

[144] Y. Xin and J.P. Gore. Measurements and Calculations of Spectral Radiation Intensities in Buoyant

Turbulent Flames. In Proceedings of the Third Joint Meeting of the U.S. Sections of the Combustion

Institute. Combustion Institute, Pittsburgh, Pennsylvania, 2003.

[145] S. Dillon and A. Hamins. Ignition Propensity and Heat Flux Profiles of Candle Flames for Fire

Investigation. In Vytenis Babrauskas, editor, Fire Science Applications to Fire Investigations, pages

363–376. Interscience Communications, London, 2003.

[146] D. Madrzykowski and R.L. Vettori. Simulation of the Dynamics of the Fire at 3146 Cherry Road

NE Washington D.C., May 30, 1999. NISTIR 6510, National Institute of Standards and Technology,

Gaithersburg, Maryland, April 2000.

[147] R.L. Vettori, D. Madrzykowski, and W.D. Walton. Simulation of the Dynamics of a Fire in a OneStory Restaurant – Texas, February 14, 2000. NISTIR 6923, National Institute of Standards and

Technology, Gaithersburg, Maryland, October 2002.

98
[148] D. Madrzykowski, G.P. Forney, and W.D. Walton. Simulation of the Dynamics of a Fire in a TwoStory Duplex – Iowa, December 22, 1999. NISTIR 6854, National Institute of Standards and Technology, Gaithersburg, Maryland, January 2002.

[149] G. Rein, A. Bar-Ilan, N. Alvares, and A.C. Fernandez-Pello. Estimating the Performance of Enclosure

Fire Models by Correlating Forensic Evidence of Accidental Fires. In Proceedings of the Tenth

International Interflam Conference, pages 1183–1194. Interscience Communications, London, 2004.

[150] M. Spearpoint, F.W. Mowrer, and K. McGrattan. Simulation of a Compartment Flashover Fire Using Hand Calculations, Zone Models and a Field Model. In Proceedings of the Third International

Conference on Fire Research and Engineering (ICFRE3), pages 3–14. Society of Fire Protection

Engineers, Bethesda, Maryland, 1999.

[151] D.J. Carpenter and C.B. Wood, editors. Proceedings of the 3rd Technical Symposium on Computer

Applications in Fire Protection Engineering. Society for Fire Protection Engineers, Bethesda, Maryland, September 2001.

[152] A.M. Christensen and D.J. Icove. The Application of NIST’s Fire Dynamics Simulator to the Investigation of Carbon Monoxide Exposure in the Deaths of Three Pittsburgh Fire Fighters. Journal of

Forensic Sciences, 49(1):1–4, 2004.

99