NIST Special Publication 1018-5

Fire Dynamics Simulator (Version 5)

Technical Reference Guide

Kevin McGrattan

Simo Hostikka

Jason Floyd

Howard Baum

Ronald Rehm

In cooperation with:

VTT Technical Research Centre of Finland

NIST Special Publication 1018-5

Fire Dynamics Simulator (Version 5)

Technical Reference Guide

Kevin McGrattan

Howard Baum

Ronald Rehm

Fire Research Division

Building and Fire Research Laboratory

Simo Hostikka

VTT Technical Research Centre of Finland

Espoo, Finland

Jason Floyd

Hughes Associates, Inc.

Baltimore, Maryland, USA

October 1, 2007

SVNRepository Revision : 727

UNITED STATES OF AMERICA

DEPARTMENT OF COMMERCE

U.S. Department of Commerce

Carlos M. Gutierrez, Secretary

National Institute of Standards and Technology

James M. Turner, Acting 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-5

Natl. Inst. Stand. Technol. Spec. Publ. 1018-5, 86 pages (October 2007)

CODEN: NSPUE2

U.S. GOVERNMENT PRINTING OFFICE

WASHINGTON: 2007

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), following the general

framework set forth in the “Standard Guide for Evaluating the Predictive Capability of Deterministic Fire

Models,” ASTM E 1355 [3].

A separate document, Fire Dynamics Simulator, User’s Guide [4] describes how the FDS software is

actually used.

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
About the Authors

Kevin McGrattan is a mathematician in the Building and Fire Research Laboratory of NIST. He received

a bachelors of science degree from the School of Engineering and Applied Science of Columbia

University in 1987 and a doctorate at the Courant Institute of New York University in 1991. He joined

the NIST staff in 1992 and has since worked on the development of fire models, most notably the Fire

Dynamics Simulator.

Simo Hostikka is a Senior Research Scientist at VTT Technical Research Centre of Finland. He is the

principal developer of the radiation and solid phase sub-models within FDS

Jason Floyd is a Senior Engineer at Hughes Associates, Inc., in Baltimore, Maryland. He received a bachelors of science and Ph.D. in the Nuclear Engineering Program of the University of Maryland. After

graduating, he won a National Research Council Post-Doctoral Fellowship at the Building and Fire

Research Laboratory of NIST, where he developed the combustion algorithm within FDS. He is currently funded by NIST under grant 60NANB5D1205 from the Fire Research Grants Program (15 USC

278f). He is the principal developer of the multi-parameter mixture fraction combustion model and

control logic within FDS.

Howard Baum is a NIST Fellow emeritus. He received a bachelors of science degree from Brooklyn

Polytechnic Institute in 1957 and a doctorate from Harvard University in 1964. After joining the staff

of the National Bureau of Standards (now NIST), he, along with Ronald Rehm, developed the low

Mach number formulation of the Navier-Stokes equations that is still the core transport algorithm

within FDS.

Ronald Rehm is a former NIST fellow, currently retired but still active in a number of projects at NIST.

He received a bachelors of science degree from Purdue University in 1960 and a doctorate from the

Massachusetts Institute of Technology in 1965. Along with Howard Baum, he formulated the original

numerical scheme that evolved over 25 years into the current FDS solver.

v
vi
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:

FDS is supported financially via internal funding at both NIST and VTT, Finland. In addition, support

is provided by other agencies of the US Federal Government, most notably the Nuclear Regulatory Agency

Office of Research. The US NRC Office of Research has funded key validation experiments, the preparation

of the FDS manuals, and the development of various sub-models that are of importance in the area of nuclear

power plant safety. Special thanks to Mark Salley and Jason Dreisbach for their efforts and support. 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.

Another source of support for FDS development has been the Microgravity Combustion Program of the

National Aeronautics and Space Administration (NASA).

Originally, the basic hydrodynamic solver was designed by Ronald Rehm and Howard Baum with programming help from Darcy Barnett, Dan Lozier and Hai Tang of the Computing and Applied Mathematics

Laboratory (CAML) at NIST, and Dan Corley of the Building and Fire Research Laboratory (BFRL). Jim

Sims of CAML developed the original visualization software. The direct pressure solver was written by

Roland Sweet of the National Center for Atmospheric Research (NCAR), Boulder, Colorado. Kuldeep

Prasad added the multiple-mesh data structures, paving the way for parallel processing. William (Ruddy)

Mell has added special routines to extend the model into areas such as microgravity combustion and wildland

fire spread. Charles Bouldin devised the basic framework of the parallel version of the code.

At NIST, 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.

William Grosshandler and Tom Cleary, both currently at NIST, developed an enhancement to the smoke

detector activation algorithm, originally conceived by Gunnar Heskestad of Factory Mutual. Steve Olenick

of Combustion Science and Engineering (CSE) implemented the smoke detector model into FDS.

William 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.

Prof. Nick Dembsey of Worcester Polytechnic Institute (WPI) has provided valuable feedback about the

pyrolysis model used within FDS.

Professor Fred Mowrer of the University of Maryland provided a simple of model of gas phase extinction

to FDS.

Chris Lautenburger of the University of California, Berkeley, and Jose Torero and Guillermo Rein of

the University of Edinburg. have provided valuable insight in the development of the solid phase model.

Finally, thanks to Randall McDermott for a very thorough review of the manuscript and suggestions for

improvements to the various routines.

vii
viii
Contents

Preface i

Disclaimer iii

About the Authors v

Acknowledgments vii

1 Introduction 1

2 Model Overview 5

2.1 Basic Description of FDS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.1.1 Type of Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.1.2 Version History . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.1.3 Model Developers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.1.4 Development Process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.1.5 Intended Uses of FDS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.1.6 Required Input for FDS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.1.7 FDS Output Quantities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

2.1.8 FDS Governing Equations, Assumptions and Numerics . . . . . . . . . . . . . . . . 9

2.1.9 Limitations of FDS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.2 Peer Review Process . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.2.1 Survey of the Relevant Fire and Combustion Literature . . . . . . . . . . . . . . . . 12

2.2.2 Review of the Theoretical Basis of the Model . . . . . . . . . . . . . . . . . . . . . 13

3 Governing Equations 15

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

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

3.1.2 The Low Mach Number Assumption and the Equation of State . . . . . . . . . . . . 17

3.1.3 Combination of the Mass and Energy Equations via the Divergence . . . . . . . . . 18

3.1.4 The Momentum Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.1.5 The Equation for Pressure (Poisson’s Equation) . . . . . . . . . . . . . . . . . . . . 20

3.1.6 Large Eddy Simulation (LES) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

3.1.7 Direct Numerical Simulation (DNS) . . . . . . . . . . . . . . . . . . . . . . . . . . 22

3.2 Combustion Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

3.2.1 Mixture Fraction Combustion Model . . . . . . . . . . . . . . . . . . . . . . . . . 23

3.2.2 Finite-Rate Reaction (DNS) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

3.3 Thermal Radiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

ix
3.4 Solid Phase Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

3.4.1 The Heat Conduction Equation for a Solid . . . . . . . . . . . . . . . . . . . . . . . 33

3.4.2 Component-Averaged Thermal Properties . . . . . . . . . . . . . . . . . . . . . . . 33

3.4.3 Pyrolysis Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

3.4.4 Radiation Heat Transfer to Solids . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

3.4.5 Convective Heat Transfer to Solids . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

3.4.6 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

3.5 Fire Detection Devices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

3.5.1 Sprinkler Activation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

3.5.2 Heat Detectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

3.5.3 Smoke Detectors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

3.6 Liquid Sprays . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

3.6.1 Droplet Size Distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

3.6.2 Droplet Transport in the Gas Phase . . . . . . . . . . . . . . . . . . . . . . . . . . 42

3.6.3 Droplet Transport on a Surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

3.6.4 Absorption and Scattering of Thermal Radiation by Droplets . . . . . . . . . . . . . 43

3.6.5 Heating and Evaporation of Liquid Droplets . . . . . . . . . . . . . . . . . . . . . . 44

3.6.6 Fire Suppression by Water . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

4 Numerical Method 47

4.1 Hydrodynamic Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

4.1.1 Simplified Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

4.1.2 The Mass Transport Equations and the Divergence . . . . . . . . . . . . . . . . . . 50

4.1.3 The Momentum Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

4.1.4 The Pressure Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

4.2 Combustion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

4.2.1 Heat Release Rate (Mixture Fraction) . . . . . . . . . . . . . . . . . . . . . . . . . 60

4.2.2 Heat Release Rate (Finite-Rate Reaction) . . . . . . . . . . . . . . . . . . . . . . . 61

4.3 Thermal Radiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

4.3.1 Discretization of the Radiative Heat Transfer Equation . . . . . . . . . . . . . . . . 62

4.3.2 Radiation Heat Transfer to Solid Objects . . . . . . . . . . . . . . . . . . . . . . . 64

4.4 Solid Phase Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

4.4.1 Solid Surface and Interior Temperatures . . . . . . . . . . . . . . . . . . . . . . . . 66

4.4.2 Convective Heat Transfer to Solid Objects . . . . . . . . . . . . . . . . . . . . . . . 66

4.4.3 Coupling the Gas and Solid Phase . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

4.5 Liquid Sprays . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

4.5.1 Interaction of Droplets and Radiation . . . . . . . . . . . . . . . . . . . . . . . . . 69

4.5.2 Heating and Evaporation of Liquid Droplets . . . . . . . . . . . . . . . . . . . . . . 69

5 Conclusion 71

References 76

A Nomenclature 77

B Derivation of the Velocity Divergence 79

C A Simple Model of Flame Extinction 85

x
Chapter 1

Introduction

Howard Baum, NIST Fellow Emeritus

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” [5].

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 [6] and

Quintiere [7], 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 [8].

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

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

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 [9]. 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 [10],

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. [11]. 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

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

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 Overview

This chapter presents general information about the Fire Dynamics Simulator, following the basic framework set forth in ASTM E 1355 [3]. This chapter also provides information about how the FDS software

is developed and maintained, often referred to as a Configuration Management Plan. Formal guidelines for

such a plan are detailed in IEEE Standard 828-1998.

2.1 Basic Description of FDS

2.1.1 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, sprinkler

discharge, and fuel sprays.

Smokeview is a companion program to FDS that produces images and animations of the results. In

recent years, its developer, Glenn Forney, has added to Smokeview the ability to visualize fire and smoke in

a fairly realistic way. In a sense, Smokeview now is, via its three-dimensional renderings, an integral part of

the physical model, as it allows one to assess the visibility within a fire compartment in ways that ordinary

scientific visualization software cannot.

Although not part of the FDS/Smokeview suite maintained at NIST, there are several third-party and

proprietary “add-ons” to FDS either available commercially or privately maintained by individual users.

Most notably, there are several Graphical User Interfaces (GUIs) that can be used to create the input file

containing all the necessary information needed to perform a simulation.

2.1.2 Version History

Version 1 of FDS was publicly released in February 2000, version 2 in December 2001, version 3 in November 2002, and version 4 in July 2004. The present version of FDS is 5, first released in October, 2007.

Starting with FDS 5, a formal revision management system has been implemented to track changes

to the FDS source code. The open-source program development tools are provided by an Internet-based

organization known as Google Code (code.google.com).

The version number for FDS has three parts. For example, FDS 5.2.12 indicates that this is FDS 5, the

fifth major release. The 2 indicates a significant upgrade, but still within the framework of FDS 5. The 12

5
indicates the twelveth minor upgrade of 5.2, mostly bug fixes and minor user requests.

2.1.3 Model Developers

Currently, FDS is maintained by the Building and Fire Research Laboratory (BFRL) of National Institute

of Standards and Technology. The developers at NIST have formed a loose collaboration of interested

stakeholders, including:

• VTT Technical Research Centre of Finland, a research and testing laboratory similar to NIST

• The Society of Fire Protection Engineers (SFPE) who conduct training classes on the use of FDS

• Fire protection engineering firms that use the software

• Engineering departments at various universities with a particular emphasis on fire

BFRL awards grants on a competitive basis to external organizations who conduct research in fire science

and engineering. Some of these grants have been used to assist the development of FDS. The role of the

grantee in supporting day to day development varies. Not all of the developers outside of NIST are grantees.

Starting with Version 5, the FDS development team uses an Internet-based development environment

called GoogleCode, a free service of the search engine company, Google. GoogleCode is a widely used

service designed to assist open source software development by providing a repository for source code,

revision control, program distribution, bug tracking, and various other very useful services.

Each member of the FDS development team has an account and password access to the FDS repository.

In addition, anonymous access is available to all interested users, who can receive the latest versions of the

source code, manuals, and other items. Anonymous users simply do not have the power to commit changes

to any of these items. The power to commit changes to FDS or its manuals can be granted to anyone on a

case by case basis.

The FDS manuals are typeset using LATEX, specifically, PDF LATEX. The LATEXfiles are essentially text

files that are under SVN (Subversion) control. The figures are either in the form of PDF or jpeg files,

depending on whether they are vector or raster format. There are a variety of LATEXpackages available,

including MiKTeX. The FDS developers edit the manuals as part of the day to day upkeep of the model.

Different editions of the manuals are distinguished by date.

2.1.4 Development Process

Changes are made to the FDS source code daily, and tracked via revision control software. However, these

daily changes do not constitute a change to the version number. After the developers determine that enough

changes have been made to the source, they release a new minor upgrade, 5.2.12 to 5.2.13, for example. This

happens every few weeks. A change from 5.2 to 5.3 might happen only a few times a year, when significant

improvements have been made to the model physics.

There is no formal process by which FDS is updated. Each developer works on various routines, and

makes changes as warranted. Minor bugs are fixed without any communication (the developers are in

different locations), but more significant changes are discussed via email or telephone calls. A suite of

simple verification calculations (included in this document) are routinely run to ensure that the daily bug

fixes have not altered any of the important algorithms. A suite of validation calculations (also included here)

are run with each significant upgrade. Significant changes to FDS are made based on the following criteria,

in no particular order:

6
Better Physics: The goal of any model is to be predictive, but it also must be reliable. FDS is a blend of

empirical and deterministic sub-models, chosen based on their robustness, consistency, and reliability.

Any new sub-model must demonstrate that it is of comparable or superior accuracy to its empirical

counterpart.

Modest CPU Increase: If a proposed algorithm doubles the calculation time but yields only a marginal

improvement in accuracy, it is likely to be rejected. Also, the various routines in FDS are expected

to consume CPU time in proportion to their overall importance. For example, the radiation transport

algorithm consumes about 25 % of the CPU time, consistent with the fact that about one-fourth to

one-third of the fire’s energy is emitted as thermal radiation.

Simpler Algorithm: If a new algorithm does what the old one did using less lines of code, it is almost

always accepted, so long as it does not decrease functionality.

Increased or Comparable Accuracy: The validation experiments that are part of this guide serve as the

metric for new routines. It is not enough for a new algorithm to perform well in a few cases. It must

show clear improvement across the suite of experiments. If the accuracy is only comparable to the

previous version, then some other criteria must be satisfied.

Acceptance by the Fire Protection Community: Especially in regard to fire-specific devices, like sprinklers and smoke detectors, the algorithms in FDS often are based on their acceptance among the

practicing engineers.

2.1.5 Intended Uses of FDS

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.6 Required Input for FDS

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

text file created by the user. The file contains information about the numerical grid, ambient environment,

building geometry, material properties, combustion kinetics, and desired output quantities. The numerical grid consists of one or more rectilinear meshes with (usually) uniform cells. All geometric features

of the scenario must 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.

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

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

and furnishings. FDS treats all of these objects as multi-layered solids, thus the physical parameters for

many 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 conductivity,

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

or from bench-scale measurements. The burning behavior of materials at different heat fluxes is more

difficult to describe, and the properties more difficult to obtain. Even though entire books are devoted to the

subject [12], it is still difficult to find information on a particular item.

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 [4].

2.1.7 FDS Output Quantities

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 millions of 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

• Heat flux, both radiative and convective

• Burning rate

• Water droplet mass per unit area

Global quantities recorded by the program include:

8
• 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. [4]. Details on the use of Smokeview are found in Ref. [13].

2.1.8 FDS 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 that is 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 combustion model based on the mixture fraction

concept. 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 originates as fuel. Unlike versions of FDS prior to 5, the reaction of

fuel and oxygen is not necessarily instantaneous and complete, and there are several optional schemes

that are designed to predict the extent of combustion in under-ventilated spaces. The mass fractions

of all of the major reactants and products can be derived from the mixture fraction by means of “state

relations,” 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 gray gas. In a limited number of cases, a wide band model can be used in place

of the gray gas model to provide a better spectral accuracy. 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 20 % of the total CPU time of a calculation, a modest cost given the complexity of

radiation heat transfer. Water droplets can absorb and scatter thermal radiation. This is important in

cases involving mist sprinklers, but also plays a role in all sprinkler cases. The absorption and scattering coefficients are based on Mie theory. The scattering from the gaseous species and soot is not

included in the model.

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.

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

9
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 is modeled using fairly

simple correlations of thermal inertia for sprinklers and heat detectors, and transport lag for smoke

detectors. Sprinkler sprays are modeled by Lagrangian particles that represent a sampling of the water

droplets ejected from the sprinkler.

2.1.9 Limitations of FDS

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 a fast, direct solver 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 10 % 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-based 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. In its simplest form, the model assumes that combustion is mixingcontrolled, and that the reaction of fuel and oxygen is infinitely fast, regardless of the temperature. For

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 10 %. Thus, the most accurate calculations using FDS do not introduce significantly greater errors in

these quantities than the vast majority of fire experiments.

10
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 are

allowed to mix and not burn, according to a few empirically-based criteria. The physical mechanisms

underlying these phenomena are complex, and are tied closely to the flame temperature and local strain

rate, neither of which are readily-available in a large scale fire simulation. Subgrid-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 (RTE) for a gray gas, and in some limited cases using a wide band model. The RTE 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, although the user may choose to use more

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 error is called “ray effect” and 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.

11
2.2 Peer Review Process

FDS is reviewed both 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 FDS

Validation Guide [?].

2.2.1 Survey of the Relevant Fire and Combustion Literature

FDS has three separate manuals – the FDS Technical Reference Guide, the FDS Validation Guide [?], and

the FDS User’s Guide [4]. Smokeview has its own User’s Guide [13]. The User Guides only describe the

mechanics of using the computer programs. The Technical Reference Guide provides the theory, algorithm

details, and verification work. The Validation Guide is new with FDS 5. It describes validation studies of

FDS, both those performed by NIST and by others.

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

in FDS was formulated by Rehm and Baum in the Journal of Research of the National Bureau of Standards [10]. 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 [14], Peyret and Taylor [15], and Ferziger and Peric [16]. 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 [11].

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

Howell [20].

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

Fire Protection Engineering [21]. Popular textbooks in fire protection engineering include those by Drysdale [22] and Quintiere [23]. 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

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

12
called the Journal of Research of NIST. An internal report, referred to as a NISTIR (NIST Inter-agency

Report), is a convenient means to disseminate information, especially when the quantity of data exceeds

what could normally be accepted by a journal. Often parts of a NISTIR are published externally, with the

NISTIR itself serving as the 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) or other electronic

media.

2.2.2 Review of the Theoretical Basis of the Model

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.

Recently, the US Nuclear Regulatory Commission (US NRC) published a seven-volume report on its

own verification and validation study of five different fire models used for nuclear power plant applications [24]. Two of the models are essentially a set of empirically-based correlations in the form of engineering “spread sheets.” Two of the models are classic two-zone fire models, one of which is the NIST

developed CFAST. FDS is the sole CFD model in the study. More on the study and its results can be found

in the FDS Validation Guide [?].

Besides formal internal and peer review, FDS is subjected to continuous scrutiny because it is available

free of charge to the general public and is used internationally by those involved in fire safety design and

post-fire reconstruction. The quality of the FDS and Smokeview User Guides 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

13
software programs undergo. Also, the source code for FDS is 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. We consider this

process as important to the development of FDS as the formal internal and external peer-review processes.

14
Chapter 3

Governing Equations

This chapter presents the governing equations of FDS. The numerical algorithm is presented in the next

chapter. The governing equations are presented as a set of partial differential equations, with appropriate

simplifications and approximations noted. The numerical method essentially consists of a finite difference

approximation of the governing equations and a procedure for updating these equations in time.

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 [10]. 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 subgrid-scale dissipative processes

are modeled. The numerical algorithm is designed so that LES becomes DNS as the grid is refined. Most

applications of FDS are LES. For example, in simulating the flow of smoke through a large, multi-room

enclosure, it is not possible to resolve the combustion and transport processes directly. However, for smallscale combustion experiments, it is possible to compute the transport and combustion processes directly.

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. [14].

Conservation of Mass

∂ρ

∂t

+∇· ρu = m ̇

000

b

(3.1)

Conservation of Momentum (Newton’s Second Law)

∂

∂t

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

15
Conservation of Energy (First Law of Thermodynamics)

∂

∂t

(ρh) +∇· ρhu =

Dp

Dt

+q ̇

000 −q ̇000

b −∇·q ̇

00 +ε (3.3)

Equation of State for a Perfect Gas

p =

ρR T

W

(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 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. 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 equation is often written in terms of the mass fractions of the individual gaseous

species, Yα:

∂

∂t

(ρYα) +∇· ρYαu = ∇· ρDα∇Yα +m ̇

000

α +m ̇

000

b,α

(3.5)

Summing these equations over all species yields the original mass conservation equation because ∑Yα = 1

and ∑m ̇

000

α = 0 and ∑m ̇

000

b,α = m ̇

000

b

, by definition, and because it is assumed that ∑ρDα∇Yα = 0. This last

assertion is not true, in general. However, transport equations are solved for total mass and all but one of

the species, implying that the diffusion coefficient of the implicit species is chosen so that the sum of all the

diffusive fluxes is zero.

In the momentum equation, the term uu is a diadic 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 fb in the momentum equation represents external 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

000

b

is the energy transferred to

the evaporating droplets. The term q ̇

00 represents the conductive and radiative heat fluxes:

q ̇

00 = −k∇T −∑α

hαρDα∇Yα +q ̇

00

r

(3.7)

where k is the thermal conductivity.

The term ε in the energy equation is known as the dissipation rate. It is 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.8)

16
= μ

"

2



∂u

∂x

2

+2



∂v

∂y

2

+2



∂w

∂z

2

+ (3.9)



∂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.10)

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, that heat is

conducted according to Fourier’s Law, and that gas species diffuse according to Fick’s Law. 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. The most important of these assumptions involves the treatment of the

diffusion and source terms that differentiates one type of CFD model from another.

3.1.2 The Low Mach Number Assumption and the Equation of State

A 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 these equations would require extremely small

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

difficult.

Following the work of Rehm and Baum [10], an approximation to the equation of state (3.4) is made by

decomposing the pressure into a “background” component and a perturbation. The original version of FDS

assumed that the background component of the pressure applied to the entire computational domain, most

often a single compartment. Starting in FDS version 5, it is now assumed that the background component

of pressure can differ from compartment to compartment. If a volume within the computational domain

is isolated from other volumes, except via leak paths or ventilation ducts, it is referred to as a “pressure

zone” and assigned its own background pressure. The pressure within the mth zone, for example, is a linear

combination of its background component and the flow-induced perturbation:

p(x,t) = pm(z,t) + p ̃(x,t) (3.11)

Note that the background pressure is a function of z, the vertical spatial coordinate, and time. For most

compartment fire applications, pm changes very little with height or time. However, for situations where the

pressure increases due to a fire in a tightly sealed enclosure, or when the height of the domain is significant,

pm takes these effects into account [25]. The ambient pressure field is denoted p0(z). Note that the subscript

0 denotes the exterior of the computational domain, not time 0. This is the assumed atmospheric pressure

stratification that serves as both the initial and boundary condition for the governing equations.

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 (in the mth pressure zone)

can be approximated

pm = ρTR ∑(Yα/Wα) = ρTR /W (3.12)

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

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

17
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 stratification of the atmosphere is derived from the relation

d p0

dz

= −ρ0(z)g (3.13)

where ρ0 is the background density. Using Eq. (3.12), the background pressure can be written as a function

of the background temperature, T0(z)

p0(z) = p∞ exp −

Z z

z∞

W g

R T0

dz(3.14)

where the subscript infinity generally refers to the ground. Usually, the temperature stratification of the

atmosphere is specified by the user, in which case p0and ρ0 are derived from Eqs. (3.14) and (3.12),

respectively.

3.1.3 Combination of the Mass and Energy Equations via the Divergence

Because of the low Mach number assumption, the divergence of the flow, ∇ · u, plays a very important

role in the overall solution scheme. The divergence is obtained by taking the material (total) derivative of

the modified Equation of State (3.12), and then substituting terms from the mass and energy conservation

equations. A few assumptions are needed. First, define the constant-pressure specific heat of the mixture:

cp = ∑α cp,αYα where cp,α is the temperature-dependent specific heat of species α. Next, define the enthalpy

h = ∑α hαYα where

hα =

Z T

Tref

cp,α(T

0

) dT0(3.15)

Finally, note that the material derivative of the background pressure, pm(z,t), is written:

Dpm

Dt

=

∂pm

∂t

+w

∂pm

∂z

=

∂pm

∂t

−wρmg (3.16)

The divergence can now be written:

∇·u = D +P

∂pm

∂t

(3.17)

where

D =

m ̇

000

b

ρ

W

Wb

+

R

W cp pm



q ̇

000 −q ̇000

b −∇·q ̇

00

−P wρmg +

W

ρ ∑α

∇· ρDα∇(Yα/Wα)−

R

W cp pm

∑α

hα∇· ρDα∇Yα +

1

ρ ∑α



W

Wα

−

hα

cpT



m ̇

000

α

P =

1

pm



R

W cp

−1



Contributions to the divergence of the flow include the heat release rate of the fire,  ̇q

000, heat losses to evaporating droplets,  ̇m

000

b

, the net heat flux from thermal conduction and radiation, ∇ · q ̇

00, updrafts of air over

considerable heights of the atmosphere, the net mass flux from gas species diffusion and production, and

18
global pressure changes. The change in the background pressure with time, ∂pm/∂t, is non-zero only if it

assumed that the compartment is tightly sealed, in which case the background pressure, pm, can no longer

be assumed constant due to the increase (or decrease) in mass and thermal energy within the enclosure. The

time derivative of the background pressure of the mth pressure zone, Ωm, is found by integrating Eq. (3.17)

over the zone volume:

∂pm

∂t

=

Z

Ωm

D dV −

Z

∂Ωm

u· dS

.Z

Ωm

P dV (3.18)

Equation (3.18) is essentially a consistency condition, ensuring that blowing air or starting a fire within a

sealed compartment leads to an appropriate decrease in the divergence within the volume.

3.1.4 The Momentum Equation

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

ρ



∂u

∂t

+ (u·∇)u



+∇p = ρg+fb +∇· τi j (3.19)

Next, we make the following substitutions:

1. Subtract the hydrostatic pressure gradient of the nth pressure zone, ρn(z,t)g, from both sides. Note

that ∇p = ρng+∇p ̃.

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

2/2−u×ω

3. Divide all terms by the density, ρ

4. Decompose the pressure term, noting that ρ∞ is the constant density at the ground:

∇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

ρ

h

(ρ−ρn)g+fb +∇· τi ji(3.20)

The numerical solution of the pressure equation obtained by taking the divergence of Eq. (3.20) 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. An evolution equation for the circulation, Γ =

H

u·dx,

reveals the sources of vorticity in the absence of any external forces, fb:

dΓ

dt

=

I 

1

ρ∞

−

1

ρ



∇p ̃· dx+

I

ρ−ρn

ρ

g · dx+

I

∇· τi j

ρ

· dx (3.21)

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 molecular and subgrid-scale stresses, as in

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

19
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; and the mass and energy equations are

not modified.

3.1.5 The Equation for Pressure (Poisson’s Equation)

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

ρ



(ρ−ρ0)g+fb +∇· τi j(3.22)

Note that the pressure  ̃p appears on both sides of Eq. (3.22). 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 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.22) 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.23)

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.24)

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.

3.1.6 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

20
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 [26].

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

ε, 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

= ρfb ·u−∇p ·u+∇·(τi j ·u)−ε (3.25)

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 [27], the viscosity μ is modeled

μLES = ρ(Cs ∆)

2



2 Si j ·Si j −

2

3

(∇·u)

2

1

2

(3.26)

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

Prt

; (ρD)l,LES =

μLES

Sct

(3.27)

The turbulent Prandtl number Prt and the turbulent Schmidt number Sct are assumed to be constant for a

given scenario.

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 rate, ε, defined in Eq. (3.10).

In the parlance of the turbulence community, the dissipation rate is related to the turbulent kinetic energy

(most often denoted by k) by the relation ε ≈ k

3/2/L, where L is a length scale.

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

to assess the improvements offered by these newer schemes for fires. 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 comparable to those obtained using the Smagorinsky model [31]. 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.

1

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

21
3.1.7 Direct Numerical Simulation (DNS)

There are some flow scenarios where it is possible to use the molecular properties μ, 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). 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 species α is given by

μα =

26.69×10−7(Wα T)

1

2

σ2

α Ωv

kg

m s

(3.28)

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

kα =

μα cp,α

Pr

W

m K

(3.29)

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

μDNS = ∑αYα μα ; kDNS = ∑αYα kα (3.30)

The binary diffusion coefficient of species α diffusing into species β is given by

Dαβ =

2.66×10−7 T

3/2

W

1

2

αβ σ

2

αβ ΩD

m2

s

(3.31)

where Wαβ = 2(1/Wα+1/Wβ)

−1

, σαβ = (σα+σβ)/2, and ΩD is the diffusion collision integral, an empirical

function of the temperature T [32]. 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)α,DNS = ρ Dα0 (3.32)

where species 0 is nitrogen.

22
3.2 Combustion Model

There are two types of combustion models used in FDS. The default model makes use of the mixture fraction,

a quantity representing the fuel and the products of combustion. For the second model, individual gas species

react according to specified Arrhenius reaction parameters. This latter model is most often used in a direct

numerical simulation (DNS) where the diffusion of fuel and oxygen can be modeled directly. However,

most often for large eddy simulations (LES), where the grid is not fine enough to resolve the diffusion of

fuel and oxygen, the mixture fraction-based combustion model is assumed.

3.2.1 Mixture Fraction Combustion Model

Given a volume containing a mixture of gas species, a mixture fraction can be defined that is the ratio of

the mass of a subset of the species to the total mass present in the volume. In combustion, the mixture

fraction is a conserved quantity traditionally defined as the (mass) fraction of the gas mixture that originates

in the fuel stream. Thus, at a burner surface the mixture fraction is 1 and in fresh air it is 0. In a region

where combustion has occurred this fraction will be comprised of any unburned fuel and that portion of

the combustion products that came from the fuel. The mixture fraction is a function of space and time,

commonly denoted Z(x,t). If it can be assumed that, upon mixing, the reaction of fuel and oxygen occurs

rapidly and completely, the combustion process is referred to as “mixing-controlled.” This implies that all

species of interest can be described in terms of the mixture fraction alone. The correspondence between the

mass fraction of an individual species and the mixture fraction is called its “state relation.” FDS versions

2 through 4 assumed that the gas mixture could be uniquely determined by the mixture fraction alone, an

assumption that implies that fuel and oxygen react instantaneously upon mixing.

For many applications, “mixed is burned” is a reasonable assumption. However, for fire scenarios where

it cannot be assumed that fuel and oxygen react completely upon mixing, for example in under-ventilated

compartments. The mixture fraction itself remains a valid quantity, but it can no longer be assumed that it

completely defines the composition of the gas mixture. If fuel and oxygen are to mix and not burn, at least

two scalar variables are needed to describe the extent to which the fuel and oxygen react. The strategy for

moving beyond the “mixed is burned” model is as follows. Instead of solving a single transport equation for

the mixture fraction Z, multiple transport equations are solved for components of the mixture fraction Zα.

Fuel mass is still conserved, since ∑Zα = Z. For example, if Z1 represents the (unburned) fuel mass fraction,

YF, and we define Z2 = Z −Z1, then Z2 is the the mass fraction of burned fuel and is the component of Z that

originates from the combustion products. With this approach it is possible to account for the mixing of fuel

and oxygen without burning. In the sections to follow, various multi-step reaction mechanisms are discussed

using this system of accounting. First, however, the single-step mixture fraction model is described.

A Single-Step, Instantaneous Reaction

Consider a simple, one-step reaction of fuel and oxygen:

CxHyOzNaMb +νO2 O2 → νCO2 CO2 +νH2O H2O+νCO CO+νS S+νN2 N2 +νMM (3.33)

Note that the nitrogen in the fuel molecule is assumed to form N2 only. Addition product species can be

specified as some number of moles of an average molecular weight species M. These products are presumed

to not consume oxygen during their formation. Soot is assumed to be a mixture of carbon and hydrogen

with the hydrogen atomic fraction given by XH. The stoichiometric coefficient, νS, represents the amount of

fuel that is converted to soot. It is related to the soot yield, yS, via the relation:

νS =

WF

WS

yS ; WS = XHWH + (1−XH)WC (3.34)

23
Figure 3.1: State relations for methane.

Likewise, the stoichiometric coefficient of CO, νCO, is related to the CO yield, yCO, via:

νCO =

WF

WCO

yCO (3.35)

The yields of soot and CO are based on “well-ventilated” or “post-flame” measurements. The increased

production of CO and soot in an under-ventilated compartment will be addressed in the following sections.

The mixture fraction, Z, can be defined in terms of the mass fraction of fuel and the carbon-carrying

products of combustion:

Z = YF +

WF

xWCO2

YCO2 +

WF

xWCO

YCO +

WF

xWS

YS (3.36)

The mixture fraction satisfies the conservation equation

ρ

DZ

Dt

= ∇· ρD∇Z (3.37)

obtained by taking a linear combination of the transport equations for the fuel and the carbon carrying

products. If it is assumed that combustion occurs so rapidly that the fuel and oxygen cannot co-exist, then

both simultaneously vanish at a flame surface:

Z(x,t) = Zf; Zf =

Y

∞

O2

sYI

F +Y

∞

O2

(3.38)

and all species can be related to Z via the “state relations” shown in Fig. 3.1. In versions of FDS prior to

5, this one-step, instantaneous reaction of fuel and oxygen was assumed. However, starting in version 5, a

more generalized formulation has been implemented and is described next.

24
A Single-Step Reaction, but with Local Extinction

The physical limitation of the single-step reaction model described in the previous section is that it assumes

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, or if the shear layer between fuel and oxidizing streams has a sufficiently large local

strain rate, fuel and oxygen may mix but may not burn. 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. Subgrid-scale modeling of gas phase suppression and extinction is still an area of active

research in the combustion community.

Simple empirical rules can be used to predict local extinction based on the oxygen concentration and

temperature of the gases in the vicinity of the flame sheet. Figure 3.2 shows values of temperature and

oxygen concentration for which burning can and cannot take place. A derivation of the model, based on the

critical flame temperature concept, is given in Appendix C. Note that once the local state of the gases falls

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

into 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. Essentially, there are now

two reactions to consider – the “null” reaction, where fuel and oxygen simply mix and do not burn; and the

“complete” reaction, where fuel and oxygen react and form products according to Eq. (3.33). Note that the

term “complete” does not imply that no soot or CO is formed, but rather that their respective production

rates are proportional to the fuel consumption rate.

With the definition of the mixture fraction, Eq. (3.36), in mind, consider a partitioning of Z into the

following components:

Z1 = YF (3.39)

Z2 =

WF

xWCO2

YCO2 +

WF

xWCO

YCO +

WF

xWS

YS (3.40)

such that Z = Z1 + Z2. Transport equations are required for both Z1 and Z2. At the burner surface, Z1 is

25
assigned the mass flux of fuel, while the mass flux for Z2 is zero. In other words, no combustion products

are generated at the fuel source. Where fuel and oxygen co-exist, a reaction occurs if conditions are favorable

in the sense shown by Fig. 3.2. If a reaction occurs, Z1 is converted to Z2 representing the conversion of fuel

to products. The heat release rate of the fire is obtained by multiplying the fuel consumption rate by the heat

of combustion.

We expand upon this by allowing for the fuel stream to be diluted with an inert gas (nitrogen) and we

note that since we are using our extra parameter to determine extinction, that the yields of CO and soot are

still fixed. Thus we can restate Eq. (3.39) and Eq. (3.40) as:

Z1 =

YF

Y

I

F

(3.41)

Z2 =

WF



x−νCO −(1−XH)νS



WCO2

YCO2

Y

I

F

(3.42)

where Y

I

F

is the fuel mass fraction at the burner surface (Y

I

N2 would be the diluent mass fraction at the burner

surface).

Since the mixture fraction variables result from linear combinations of the species transport equations,

the converse is also true – that species mass fractions are linear combinations of the mixture fraction variables. The mass fractions of the species in the mixture, Yα(Z1,Z2), are found via:

YF = Y

I

F Z1

YN2 = (1−Z) Y

∞

N2 +Y

I

N2

Z1 +

νN2WN2

WF

Y

I

F Z2

YO2 = (1−Z) Y

∞

O2 −

νO2WO2

WF

Y

I

F Z2

YCO2 =

νCO2WCO2

WF

Y

I

F Z2

YH2O =

νH2OWH2O

WF

Y

I

F Z2 (3.43)

YCO =

νCOWCO

WF

Y

I

F Z2 (3.44)

YS =

νSWS

WF

Y

I

F Z2 (3.45)

YM =

νMWM

WF

Y

I

F Z2 (3.46)

The stoichiometric coefficients are defined:

νN2 =

a

2

νO2 = νCO2 +

νCO +νH2O −z

2

νCO2 = x−νCO −(1−XH)νS

νM = b

νH2O =

y

2

−XHνS (3.47)

νCO =

WF

WCO

yCO (3.48)

νS =

WF

WS

yS (3.49)

It is important to note that the definitions of Z1 and Z2, unlike in the single parameter model, do not imply

anything regarding the rate of combustion, only that the combustion occurs in a single step.

CO Production (Two-Step Reaction with Extinction)

The previous section describes the “complete” reaction as the conversion of fuel to products such that the

production rate of each product species is proportional to the fuel consumption rate. This means that for

each fuel molecule, fixed amounts of CO2, H2O, CO, and soot are formed and these products persist in the

26
plume indefinitely with no further reaction. This is not an unreasonable assumption if the purpose of the fire

simulation is to assess the impact of the fire on the larger space. However, in under-ventilated fires, soot and

CO are produced at higher rates, and exist within the fuel-rich flame envelope at higher concentrations, than

would otherwise be predicted with a single set of fixed yields that are based on post-flame measurements.

To account for the production of CO and its eventual oxidation at the flame envelope or within a hot upper

layer, an additional reaction is now needed:

CxHyOzNaMb +ν

0

O2O2 → νH2OH2O+ (ν

0

CO +νCO)CO+νS S+νN2 N2 +νMM (3.50)

ν

0

CO h

CO+

1

2

O2 → CO2

i

(3.51)

The brackets around the second reaction are there merely to emphasize that the sum of the two reactions

equal Eq. (3.33). There are two stoichiometric coefficients for CO – the first, ν

0

CO = x − (1 − XH)νS − νCO,

represents CO that is produced in the first step of the reaction that can potentially be converted to CO2

assuming the conditions are favorable. ν

0

CO is equivalent to νCO2

in Eq. (3.33). The second coefficient,

νCO, is the so-called “well-ventilated,” or “post-flame,” value that was introduced in the previous section.

The proposed model of CO production still does not contain the necessary kinetic mechanism to predict the

“post-flame” concentration of CO without the prescription of the measured value of the post-flame CO yield.

Rather, the proposed model includes the production of large amounts of CO in the first step of a two-step

reaction, followed by a partial conversion to CO2 if there is a sufficient amount of oxygen present.

To describe the composition of the gas species, the mixture fraction, Z, must now be decomposed into

three components:

Z1 =

YF

Y

I

F

(3.52)

Z2 =

WF



x−(1−XH)νS



WCO

YCO

Y

I

F

(3.53)

Z3 =

WF



x−(1−XH)νS



WCO2

YCO2

Y

I

F

(3.54)

Here, x−(1−XH)νS represents the number of carbon atoms in the fuel molecule that are oxidized; that is, the

carbon atoms that are not converted to soot. The decomposition of Z into three components is numerically

convenient. However, the recovery of the individual species mass fractions requires some care. The mass

fraction of any species in the mixture, Yα(Z1,Z2,Z3), is still found via linear combinations of the mixture

fraction variables:

YF = Z1 Y

I

F

YN2 = (1−Z) Y

∞

N2 +Y

I

N2

Z1 +

νN2WN2

WF

Y

I

F

(Z2 +Z3)

YO2 = (1−Z) Y

∞

O2 −

WO2 Y

I

F

WF

(ν

0

O2

Z2 +νO2 Z3)

YCO2 =

νCO2WCO2

WF

Y

I

F Z3

YH2O =

νH2OWH2O

WF

Y

I

F

(Z2 +Z3) (3.55)

YCO =

(νCO +ν

0

CO)WCO2

WF

Y

I

F Z2 (3.56)

YS =

νSWS

WF

Y

I

F

(Z2 +Z3) (3.57)

YM =

νMWM

WF

Y

I

F

(Z2 +Z3) (3.58)

The stoichiometric coefficients are defined:

27
νN2 =

a

2

ν

0

O2 =

ν

0

CO +νH2O

2

νO2 = νCO2 +

νCO +νH2O −z

2

νCO2 = x−νCO −(1−XH)νS

νM = b

νH2O =

y

2

−XHνS (3.59)

ν

0

CO = x−νCO −(1−XH)νS (3.60)

νCO =

WF

WCO

yCO (3.61)

νS =

WF

WS

yS (3.62)

Although these formulae appear complicated, most are determined directly from the composition of the fuel

molecule. The only information expected of the modeler are the fuel composition, the soot and CO yields,

and the atomic fraction of hydrogen in the soot.

Heat Release Rate

The discussion of the various multi-step reactions above is essentially book-keeping, the accounting of the

gas molecules formed in the combustion process. But what of the heat released?

When tracking species with two mixture fraction parameters, there is a single step reaction, the conversion of fuel into a fixed, predefined set of combustion products. Fuel and oxygen react within a grid cell and

release energy according to the rate of fuel consumption:

q ̇

000 = m ̇000

F ∆HF (3.63)

The amount of fuel in the grid cell is obtained directly from the definition of Z1 whereas the amount of

O2 is obtained using the state relationships discussed in the prior section. Combustion is either allowed

or disallowed2 using the relation shown in Fig. 3.2. If combustion is allowed to occur in a grid cell, the

single step combustion is assumed to be infinitely fast. Combustion will consume all of either the fuel or the

oxygen in the grid cell.

In the case of the two-step, three parameter mixture fraction model, the first step converts the fuel to CO

and other combustion products, and the second step oxidizes the CO into CO2. The first step is determined

as it is for the single step reaction. The second step, however, is assumed to be rate dependent. Determining

if the second step can take place in a grid cell is done in one of two ways. For most LES simulations, flame

temperatures will not be resolved on the grid because the heat released is smeared out over a grid cell much

larger than the flame sheet thickness. Since a grid cell with heat released from the first step has a flame,

CO conversion is always allowed in these grid cells. In all other grid cells, the local temperature is used to

perform a finite-rate computation that determines the rate of CO conversion.

On-Going Research

One can envision extending the three parameter mixture fraction model to a fourth parameter for soot production:

Z4 =

WF

xWS

YS (3.64)

This extension would allow for combustion models that predict soot formation along with extinction and

CO production.

2Note that the user has control over the parameters associated with local gas phase extinction.

28
Additional partitionings of the mixture fraction could be added as well; however, there is a point of

diminishing returns. The present advantage of using a mixture fraction approach is to track the movement

of a number of gas species using a smaller number of scalar variables at the expense of having to use state

relationships to obtain species mass fractions. As more mixture fraction variables are added, the net savings

in fewer species transport equations vs. the expense of state relationships will diminish.

29
3.2.2 Finite-Rate Reaction (DNS)

In a DNS calculation, the fine grid resolution enables the direct modeling of the diffusion of chemical species

(fuel, oxygen, and combustion products). Since the flame is being resolved in a DNS calculation, the local

gas temperatures can be used to determine the reaction kinetics. Thus, it is possible to implement a relatively

simple set of one or more chemical reactions to model the combustion. Consider the reaction of oxygen and

a hydrocarbon fuel

CxHy +νO2 O2 −→ νCO2 CO2 +νH2O H2O (3.65)

If this were modeled as a single-step reaction, the reaction rate would be given by the expression

d[CxHy]

dt

= −B[CxHy]

a

[O2]

b

e

−E/RT (3.66)

Suggested values of B, E, a and b for various hydrocarbon fuels are given in Refs. [33, 34]. 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. Improved predictions of the heat release rate may be possible by considering a multistep set of reactions. However, each additional gas species defined in the computation incurs a roughly 5 %

increase in the CPU time.

30
3.3 Thermal Radiation

Energy transport consists of convection, conduction and radiation. Convection of heat is accomplished via

the solution of the basic conservation equations. Gains and losses of heat via conduction and radiation are

represented by the divergence of the heat flux vector in the energy equation, ∇ · q ̇

00. This section describes

the equations associated with the radiative part, q ̇

00

r

.

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

s·∇Iλ(x,s) = −

h

κ(x,λ) +σs(x,λ)

i

Iλ(x,s) +B(x,λ) + σs(x,λ)

4π

Z

4π

Φ(s,s

0

) Iλ(x,s

0

) ds

0

(3.67)

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,λ)

h

Ib(x)−Iλ(x,s)

i

(3.68)

where Ib(x) is the source term given by the Planck function (see below).

In practical simulations the spectral (λ) dependence cannot 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 band specific RTE is

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

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.70)

where σ is the Stefan-Boltzmann constant. The calculation of factors Fn is explained in Ref. [20]. 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.71)

Even with a reasonably small number of bands, solving multiple RTEs is very time consuming. Fortunately, in most large-scale fire scenarios soot is the most important combustion product controlling the

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 then lumped into one

absorption coefficient (N = 1) and the source term is given by the blackbody radiation intensity

Ib(x) = σT(x)

4

/π (3.72)

This is the default mode of FDS and appropriate for most problems of fire engineering. 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. From a series of numerical experiments it

has been found that six bands (N = 6) are usually enough to improve the accuracy in these cases. The limits

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

and water. If the absorption of the fuel is known to be important, separate bands can be reserved for fuel,

31
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

and the total number of bands is increased to nine (N = 9). For simplicity, the fuel is assumed to be CH4.

The limits of the bands are shown in Table 3.1.

For the calculation of the gray or band-mean absorption coefficients, κn, a narrow-band model, RadCal [35], 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 fourth-power dependence

on the temperature, 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 be computed directly

κ Ib =



κσT

4/π Outside flame zone

max(χr q ̇

000/4π , κσT4/π) Inside flame zone (3.73)

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.3 Near the flame in large scale calculations, neither κ nor T can

be computed reliably, hence the inclusion of the empirical radiation loss term which is designed to partition

the fire’s heat release rate in accordance with measured values.

The boundary condition for the radiation intensity leaving a gray diffuse wall is given as

Iw(s) = εσT

4

w

π

+

1−ε

π

Z

s

0

·nw<0

Iw(s

0

) |s

0

·nw| ds

0

(3.74)

where Iw(s) is the intensity at the wall, ε is the emissivity, and Tw is the wall surface temperature. The radiant

heat flux vector q ̇

00

r

is defined

q ̇

00

r

(x) = Z

4π

s

0

I(x,s

0

) ds

0

(3.75)

The gas phase contribution to the radiative loss term in the energy equation is

−∇·q ̇

00

r

(x)(gas) = κ(x) [U(x)−4πIb(x)] ; U(x) = Z

4π

I(x,s

0

)ds

0

(3.76)

In words, the net radiant energy gained by a grid cell is the difference between that which is absorbed and

that which is emitted.

3The radiative fraction, χr, is a useful quantity in fire science. Usually, it is understood to be the fraction of the total heat release

rate that takes the form of thermal radiation. For most combustibles, χr is between 0.3 and 0.4. However, in Eq. (3.73), χr is

interpreted as the fraction of energy radiated from the combustion region. For a small fire (D < 1 m), the local χr is approximately

equal to its global counterpart. However, as the fire increases in size, the global value will typically decrease due to a net reabsorption of the thermal radiation by the increasing smoke mantle.

32
3.4 Solid Phase Model

FDS assumes that solid obstructions consist of multiple layers, with each layer composed of multiple material components that can undergo multiple thermal degradation reactions. Each reaction forms a combination of solid residue (i.e. another material component), water vapor, and/or fuel vapor. Heat conduction is

assumed only in the direction normal to the surface. This section describes the single mass and energy conservation equation for solid materials, plus the various coefficients, source terms, and boundary conditions.

3.4.1 The Heat Conduction Equation for a Solid

A one-dimensional heat conduction equation for the solid phase temperature Ts(x,t) is applied in the direction x pointing into the solid (the point x = 0 represents the surface)4

ρscs

∂Ts

∂t

=

∂

∂x

ks

∂Ts

∂x

+q ̇

000

s

(3.78)

Section 3.4.2 describes the component-averaged material properties, ks and ρscs. The source term,  ̇q

000

s

,

consists of chemical reactions and radiative absorption:

q ̇

000

s = q ̇

000

s,c +q ̇

000

s,r

(3.79)

Section 3.4.3 describes the term  ̇q

000

s,c

, which is essentially the heat production (loss) rate given by the pyrolysis models for different types of solid and liquid fuels. Section 3.4.4 describes the term  ̇q

000

s,r

, the radiative

absorption and emission in depth. Section 3.4.5 describes the convective heat transfer to the solid surface.

Finally, Section 3.4.6 describes how the radiative and convective heat fluxes are applied as boundary conditions for Eq. (3.78).

3.4.2 Component-Averaged Thermal Properties

The conductivity and volumetric heat capacity of the solid are defined

ks =

Nm

∑

α=1

Xα ks,α ; ρscs =

Nm

∑

α=1

ρs,α cs,α (3.80)

Nm is the number of material components forming the solid. ρs,α is the component density

ρs,α = ρsYα (3.81)

where ρsis the density of the composite material and Yα is the mass fraction of material component α. The

solid density is the sum of the component densities

ρs =

Nm

∑

α=1

ρs,α (3.82)

4

In cylindrical and spherical coordinates, the heat conduction equation is written

ρscs

∂Ts

∂t

=

1

r

∂

∂r



rks

∂Ts

∂r



+q ̇

000

s

; ρscs

∂Ts

∂t

=

1

r

2

∂

∂r



r

2

ks

∂Ts

∂r



+q ̇

000

s

(3.77)

FDS offers the user these options, with the assumption that the obstruction is not actually recti-linear, but rather cylindrical or

spherical in shape. This option is useful in describing the behavior of small, complicated “targets” like cables or heat detection

devices.

33
Xα is the volume fraction of component α

Xα =

ρs,α

ρα

, Nm

∑

α0=1

ρs,α0

ρα0

(3.83)

where ρα is the density of material α in its pure form. Multi-component solids are defined by specifying the

mass fractions, Yα, and densities, ρα, of the individual components of the composite.

3.4.3 Pyrolysis Models

This section describes how solid phase reactions and the chemical source term in the solid phase heat conduction equation,  ̇q

000

s,c

, are modeled. This is what is commonly referred to as the “pyrolysis model,” but

it actually can represent any number of reactive processes, including evaporation, charring, and internal

heating.

Specified Heat Release Rate

Often the intent of a fire simulation is merely to predict the transport of smoke and heat from a specified fire.

In other words, the heat release rate is a user input, not something the model predicts. In these instances, the

desired HRR is translated into a mass flux for fuel at a given solid surface, which can be thought of as the

surface of a burner:

m ̇

00

f =

f(t) q ̇

00

user

∆H

(3.84)

Usually, the user specifies a desired heat release rate per unit area,  ̇q

00

user, plus a time ramp, f(t), and the

mass loss rate is computed accordingly.

Solid Fuels

Solids can undergo simultaneous reactions under the following assumptions:

• instantaneous release of volatiles from solid to the gas phase,

• local thermal equilibrium between the solid and the volatiles,

• no condensation of gaseous products, and

• no porosity effects5

Each material component may undergo several competing reactions, and each of these reactions may produce some other solid component (residue), gaseous fuel, and/or water vapor according to the yield coefficients νs

, νf and νw, respectively. These coefficients should usually satisfy νs +νf +νw = 1, but smaller

yields may also be used to take into account the gaseous products that are not explicitly included in the

simulation.

Consider material component α that undergoes Nr,α separate reactions. We will use the index β to

represent one of these reactions:

Materialα → νs,αβ Residueαβ +νw,αβ H2O+νf,αβ HC (3.85)

5Although porosity effects are not explicitly included in the model, it is possible to account for it because the volume fractions

defined by Eq. (3.83) need not sum to unity, in which case the thermal conductivity and absorption coefficient are effectively

reduced.

34
The local density of material component α evolves in time according to the solid phase species conservation

equation

∂

∂t



ρs,α

ρs0



= −

Nr,α

∑

β=1

rαβ +Sα (3.86)

which says that the mass of conponent α is consumed by the solid phase reactions rαβ and produced by

other reactions. rαβ is the rate of reaction β in units (1/s) and ρs0 is the initial density of the material layer.

Sα is the production rate of material component α as a result of the reactions of the other components. The

reaction rates are functions of local mass concentration and temperature, and calculated as a combination of

Arrhenius and power functions:

rαβ =



ρs,α

ρs0

ns,αβ

Aαβ exp−

Eαβ

RTs



max0,Ts −Tthr,αβnt,αβ (3.87)

where Tthr,αβ is a threshold temperature that can be used to dictate that the reaction must not occur below a

user-specified temperature. By default, the term is deactivated (Tthr,αβ = 0 K).

The production term Sα is the sum over all the reactions where the solid residue is material α

Sα =

Nm

∑

α0=1

Nr,α0

∑

β=1

νs,α0β rα0β (where Residueα0β = Materialα) (3.88)

The volumetric production rate of fuel gases and water vapor are

m ̇

000

f = ρs0

Nm

∑

α=1

Nr,α

∑

β=1

νf,αβ rαβ (3.89)

m ̇

000

w = ρs0

Nm

∑

α=1

Nr,α

∑

β=1

νw,αβ rαβ (3.90)

It is assumed that the fuel gases and water vapor are transported instantaneously to the surface, where the

mass fluxes are given by:

m ̇

00

f =

Z L

0

m ̇

000

f

(x)dx ;  ̇m

00

w =

Z L

0

m ̇

000

w

(x)dx (3.91)

where L is the thickness of the surface. The chemical source term of the heat conduction equation consists

of the heats of reaction and the differences in latent heats between the original material and the products

q ̇

000

s,c

(x) = ρs0

Nm

∑

α=1

Nr,α

∑

β=1

rαβ(x)



Hr,αβ −

Z T(x)

T0

∆c(θ)dθ



(3.92)

It is possible to define a reaction where the sum of the residue, fuel and water yields is less than one. FDS

then assumes that the solid is converted to a gas that is not tracked by the gas phase calculation. In the

computation of ∆c(T), the specific heat of this dummy species is assumed to be equal to the specific heat of

the original material cs,α

∆c(T) = cs,α(T)−νs,αβ cs,αβ(T)−νf,αβ cf(T)−νw,αβ cw(T)−(1− ∑

γ=s, f,w

νγ,αβ) cs,α(T) (3.93)

cs,αβ is the specific heat of the residue of the βth reaction of material α. The temperature integrals of the

latent heat contents are computed in advance for all the reactions of each material and the values of  ̇q

000

s,c

are tabulated between T0 and T0 + 2000 K. During the pyrolysis calculation, the values are found from a

“look-up” table.

35
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. According to the Clausius-Clapeyron relation, the volume

fraction of the fuel vapor above the surface is a function of the liquid boiling temperature

Xf = exp−

hvWf

R



1

Ts

−

1

Tb

 (3.94)

where hv is the heat of vaporization, Wfis the molecular weight, Tsis the surface temperature, and Tb is the

boiling temperature of the fuel [36].

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 Radiation Heat Transfer to Solids

If it is assumed that the thermal radiation from the surrounding gases is absorbed within an infinitely thin

layer at the surface of the solid obstruction, then the net radiative heat flux is the sum of incoming and

outgoing components,  ̇q

00

r = q ̇

00

r,in −q ̇

00

r,out, where:

q ̇

00

r,in = ε

Z

s

0

·nw<0

Iw(s

0

) |s

0

·nw| dΩ (3.95)

q ̇

00

r,out = εσT

4

w

(3.96)

However, many common materials do not have infinite optical thickness. Rather, the radiation penetrates

the material to some finite depth. The radiative transport within the solid (or liquid) can be described as

a source term in Eq. (3.78). A “two-flux” model based on the Schuster-Schwarzschild approximation [20]

assumes the radiative intensity is constant inside the “forward” and “backward” hemispheres. The transport

equation for the intensity in the “forward” direction is

1

2

dI+(x)

dx

= κs



Ib −I

+(x)



(3.97)

where x is the distance from the material surface and κsis the absorption coefficient

κs =

Nm

∑

α=1

Xα κs,α (3.98)

A corresponding formula can be given for the “backward” direction. Multiplying Eq. 3.97 by π gives us the

“forward” radiative heat flux into the solid

1

2

dq ̇

+

r

(x)

dx

= κs



σT

4

s −q ̇

+

r

(x)



(3.99)

The radiative source term in the heat conduction equation is a sum of the “forward” and “backward” flux

gradients

q ̇

000

s,r

(x) = dq ̇

+

r

(x)

dx

+

dq ̇

−

r

(x)

dx

(3.100)

The boundary condition for Eq. 3.99 at the solid (or liquid) surface is given by

q ̇

+

r

(0) = q ̇

00

r,in + (1−ε)q ̇

−

r

(0) (3.101)

where  ̇q

−

r

(0) is the “backward” radiative heat flux at the surface. In this formulation, the surface emissivity

and the internal absorption are assumed to be independent properties of the material.

36
3.4.5 Convective Heat Transfer to Solids

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.102)

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.103)

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) [18], 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.6 Boundary Conditions

The boundary condition on the front surface is

−ks

∂Ts

∂x

(0,t) = q ̇

00

c +q ̇

00

r

(3.104)

If the internal radiation is solved for a solid, the radiation boundary condition  ̇q

00

r

is not used.

On the back surface, two possible boundary condition types may be prescribed by the user. (1) If the

back surface is assumed to be open either to an ambient void or to another part of the computational domain,

the back side boundary condition is similar to that of the front side. (2) If the back side is assumed to be

perfectly insulated, an adiabatic boundary condition is used

ks

∂Ts

∂x

= 0 (3.105)

37
3.5 Fire Detection Devices

FDS predicts the thermal environment resulting from a fire, but it relies on various empirical models that

describe the activation of various fire detection devices. These models are described in this section.

3.5.1 Sprinkler Activation

The temperature of the sensing element (or “link”) of an automatic fire sprinkler is estimated from the

differential equation put forth by Heskestad and Bill [37], 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.106)

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 RTI and C-Factor are determined

experimentally. The constant C2 has been empirically determined by DiMarzo and co-workers [38, 39, 40]

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 the sprinkler nor heat detector

models account for thermal radiation.

3.5.2 Heat Detectors

As far as FDS is concerned, a heat detector is just a sprinkler with no water spray. In other words, the

activation of a heat detector is governed by Eq. (3.106), but with just the first term on the right hand side:

dTl

dt

=

p

|u|

RTI

(Tg −Tl) (3.107)

Both the RTI and activation temperature are determined empirically.

3.5.3 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 [21]. 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 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

38
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 [21]. 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)

L/u

(3.108)

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. [41] 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.109)

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.110)

where 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

free-stream at time t −δte.

An analytical solution for Eq. (3.110) can be found, but it is more convenient to simply integrate it

numerically as is done for sprinklers and heat detectors. Then, the predicted mass fraction of smoke in the

sensing chamber, Yc(t), can be converted into an expression for the percent obscuration per unit length by

computing:



1−e

−κρYcl



×100 (3.111)

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 [42].

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 needs to be prescribed. Eq. (3.110) 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 [43].

Additionally, FDS can model the behavior of beam and aspiration smoke detectors. For a beam detector

the user specifies the emitter and receiver positions and the total obscuration at which the detector will alarm.

FDS will then integrate the obscuration over the path length using the predicted soot concentration in each

grid cell along the path. For an aspiration detector the user specifies the sampling locations, the flow rate at

39
each location, the transport time from each sampling point to the detector, the flow rate of any bypass flow,

and the total obscuration at which the detector will alarm. FDS will compute that soot concentration at the

detector by weighting the predicted soot concentrations at the sampling locations with their flow rates after

applying the appropriate time delay.

40
3.6 Liquid Sprays

There are several uses for liquid sprays within FDS, most notably for sprinklers, but also for liquid fuel

sprays. The basic physics is the same in each case, but the liquid properties are different. Also, simulating

the effects of a sprinkler spray involves a number of elements beyond just activation: computing the droplet

trajectories and tracking the water as it drips onto the burning surface.

3.6.1 Droplet Size Distribution

A spray consists of a sampled set of spherical droplets. The size distribution is expressed in terms of its Cumulative Volume Fraction (CVF), a function that relates the fraction of the liquid volume (mass) transported

by droplets less than a given diameter. Researchers at Factory Mutual have suggested that the CVF for a

sprinkler may be represented by a combination of log-normal and Rosin-Rammler distributions [44]

F(d) =







√

1

2π

Z d

0

1

σd

0 e

−

[ln(d

0/dm)]2

2σ2 dd

0

(d ≤ dm)

1−e

−0.693(

d

dm

)

γ

(dm < d)

(3.112)

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.6 The median droplet

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 [45]

dm

D

∝ We−

1

3 (3.113)

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 =

ρdu

2

dD

σd

(3.114)

where ρd is the density of liquid, ud is the discharge velocity, and σd is the liquid surface tension (72.8×10−3

N/m at 20 ◦C for water). The discharge velocity can be computed from the mass flow rate, a function of the

sprinkler’s operating pressure and orifice coefficient known as the “K-Factor.” FM reports that the constant

of proportionality in Eq. (3.113) 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 [45].

In the numerical algorithm, the size of the droplets are chosen to mimic the Rosin-Rammler/log-normal

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

dd

0

; F

0 ≡

dF

dd

(3.115)

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

)dd

0

(3.116)

Figure 3.3 displays typical Cumulative Volume Fraction and Cumulative Number Fraction functions.

6The Rosin-Rammler and log-normal distributions are smoothly joined if σ = 2/(

√

2π(ln 2) γ) = 1.15/γ .

41
Figure 3.3: 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.

3.6.2 Droplet Transport in the Gas Phase

Every droplet from a given sprinkler or nozzle is not tracked. Instead, a sampled set of the droplets is

tracked. The procedure for selecting droplet sizes is as follows: Suppose the mass flow rate of the liquid

is  ̇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.116), and then compute a weighting constant C from the mass balance

m ̇ δt = C

N

∑

i=1

4

3

πρw



di

2

3

(3.117)

The mass and heat transferred from each droplet will be multiplied by the weighting factor C.

For a spray, the force term fb in Eq. (3.2) represents the momentum transferred from the 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

fb =

1

2

∑ρCDπr

2

d

(ud −u)|ud −u|

δxδyδz

(3.118)

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.119)

where md is the mass of the droplet. The drag coefficient is a function of the local Reynolds number (based

on droplet diameter)

CD =







24/ReD ReD < 1

241+0.15Re0.687

D



/ReD 1 < ReD < 1000

0.44 1000 < ReD

(3.120)

42
ReD =

ρ|ud −u|2rd

μ(T)

(3.121)

where μ(T) is the dynamic viscosity of air at temperature T.

3.6.3 Droplet Transport on a Surface

When a liquid 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 for water to be on the order of 0.5 m/s [46, 47]. While attached

to a surface, the “droplet” is assumed to form a thin film of liquid that transfers heat to the solid, and heat

and mass to the gas. Details are included in the next section.

3.6.4 Absorption and Scattering of Thermal Radiation by Droplets

The attenuation of thermal radiation by liquid droplets is an important consideration, especially for water

mist systems [48]. Liquid droplets attenuate thermal radiation through a combination of scattering and

absorption [49]. 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.67) 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

) ds

0

(3.122)

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.123)

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 initial

droplet size distribution, but a mean diameter depending on the location x. For the numerical implementation, the above equations are written in the form

κd(x,λ) = Ad(x)

R ∞

0

f(r,dm(x)) Ca(r,λ)

π(dm(x)/2)

2 dr

σd(x,λ) = Ad(x)

R ∞

0

f(r,dm(x)) Cs(r,λ)

π(dm(x)/2)

2 dr

(3.124)

where Ad is the total cross sectional area per unit volume of the droplets. Ad is approximated as

Ad ≈

ρd(x)

2ρwdm(x)/3

(3.125)

An accurate computation of the in-scattering integral on the right hand side of Eq. (3.122) 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

43
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.126)

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.127)

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 an effective scattering coefficient section

σd(x,λ) = 4πN(x)

4π−δΩl

Z ∞

0

(1−χf) Cs(r,λ) dr (3.128)

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.129)

This equation can be integrated over the spectrum to get the band specific RTE’s. The procedure is exactly

the same as that 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.130)

where the source function is based on the average droplet temperature within a cell. The droplet contribution

to the radiative loss term in the energy equation is

−∇·q ̇

00

r

(x)(droplets) = κd(x) [U(x)−4πIb,d(x)] (3.131)

For each individual droplet, the radiative heating/cooling power is computed as

q ̇r =

md

ρd(x)

κd(x) [U(x)−4πIb,d(x)] (3.132)

where md is the mass of the droplet and ρd(x) is the total density of droplets in the cell.

The absorption and scattering cross sections and the scattering phase function are calculated using the

MieV code developed by Wiscombe [50]. Both κd and σd are averaged over the possible droplet radii and

wavelength before the actual simulation. A constant “radiation” temperature, Trad, is used in the wavelength

averaging. Trad should be selected to represent a typical radiating flame temperature. A value of 1173 K

is used by default. The averaged quantities, 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 look-up using

the local mean droplet diameter [51]. Details of the computation are given in Section 4.5.1.

3.6.5 Heating and Evaporation of Liquid Droplets

Liquid “droplets” are represented either as discrete airborne spheres propelled through the gas, or as rectangular blocks that collectively form a thin liquid film on solid objects. These “droplets” are still individually

tracked as lagrangian particles, but the heat and mass transfer coefficients are different. In the discussion to

follow, the term “droplets” will be used to describe either form.

44
Over the course of a time step of the gas phase solver, the droplets in a given grid cell evaporate as a

function of the liquid equilibrium vapor mass fraction, Yl, the local gas phase vapor mass fraction, Yg, the

(assumed uniform) liquid temperature, Tl, and the local gas temperature, Tg. If the droplet is attached to a

surface, Tsis the solid temperature. The mass and energy transfer between the gas and the liquid can be

described by the following set of equations [52]

dml

dt

= −A hm ρ(Yl −Yg) (3.133)

ml cl

dTl

dt

= A h(Tg −Tl) +A hs (Ts −Tl) +q ̇r +

dml

dt

hv (3.134)

Here, mlis the mass of the liquid droplet (or that fraction of the surface film associated with the formerly

airborne droplet), A is the surface area of the liquid droplet (or that fraction of the film exposed to the gas

and to the wall), hm is the mass transfer coefficient to be discussed below, ρ is the gas density, clis the

liquid specific heat, h is the heat transfer coefficient between the liquid and the gas, hsis the heat transfer

coefficient between the liquid and the solid surface,  ̇qris the rate of radiative heating of the droplet, and hv

is the latent heat of vaporization of the liquid. The vapor mass fraction of the gas, Yg, is obtained from the

gas phase mass conservation equations, and the liquid equilibrium vapor mass fraction is obtained from the

Clausius-Clapeyron equation

Xl = exp

hvWl

R



1

Tb

−

1

Tl

 ; Yl =

Xl

Xl(1−Wa/Wl) +Wa/Wl

(3.135)

where Xd is the equilibrium vapor volume fraction, Wlis the molecular weight of the evaporated liquid, Wa is

the molecular weight of air, R is the universal gas constant, and Tb is the boiling temperature of the liquid.

Mass and heat transfer between liquid and gas are described with analogous empirical correlations. The

mass transfer coefficient, hm, is described by the empirical relationships [19]:

hm =

Sh Dlg

L

; Sh =







2+0.6 Re

1

2

D

Sc

1

3 droplet

0.037 Re

4

5

L

Sc

1

3 film

(3.136)

Sh is the Sherwood number, Dlg is the binary diffusion coefficient between the liquid vapor and the surrounding gas (usually assumed air), L is a length scale equal to either the droplet diameter or 1 m for a

surface film, ReD is the Reynolds number of the droplet (based on the diameter, D, and the relative airdroplet velocity), ReL is the Reynolds number based on the length scale L, and Sc is the Schmidt number

(ν/Dlg, assumed 0.6 for all cases).

An analogous relationship exists for the heat transfer coefficient:

h =

Nu k

L

; Nu =







2+0.6 Re

1

2

D

Pr

1

3 droplet

0.037 Re

4

5

L

Pr

1

3 film

(3.137)

Nu is the Nusselt number, k is the thermal conductivity of the gas, and Pr is the Prandtl number (assumed

0.7 for all cases).

3.6.6 Fire Suppression by Water

The previous 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

45
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.

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. [53]. 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.138)

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.139)

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.140)

Unfortunately, this analysis is based on global water flow and burning rates. Equation (3.138) 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 [46, 47]

m ̇

00

f

(t) = m ̇

00

f,0

(t) e

−

R

k(t)dt (3.141)

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.142)

Note that a is an empirical constant.

46
Chapter 4

Numerical Method

This chapter presents a description of the numerical methods used to solve the equations of the previous

chapter, with sections in roughly the same order.

4.1 Hydrodynamic Model

Each of the conservation equations emphasizes the importance of the velocity divergence and vorticity fields,

as well as the close relationship between the thermally expandable fluid equations [10] and the Boussinesq

equations for which the authors have developed highly efficient solution procedures [54, 55]. 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.

4.1.1 Simplified Equations

Conservation of Mass

∂ρ

∂t

+u·∇ρ = −ρ∇·u (4.1)

Conservation of Species

∂ρYα

∂t

+u·∇ρYα = −ρYα∇·u+∇· ρDα∇Yα +m ̇

000

α

(4.2)

Conservation of Momentum

∂u

∂t

+u×ω+∇H =

1

ρ

((ρ−ρ0)g+fb +∇· τi j) (4.3)

Pressure Equation

∇

2H = −

∂(∇·u)

∂t

−∇·F ; F = u×ω−

1

ρ

((ρ−ρm)g−fb −∇· τi j) (4.4)

Equation of State

pm(z,t) = ρTR ∑αYα/Wα (4.5)

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.

47
Temporal Discretization

All calculations start with ambient initial conditions. At the beginning of each time step, the density, ρ

n

,

species mass fractions, Y

n

α

, velocity, u

n

, modified pressure, H n, and background pressure in zone m, p

n

m

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 ρ, Yα, and pm are estimated at the next time step with an explicit Euler

step. For example, the density is estimated by

ρ

(n+1)e = ρn −δt(un

·∇ρn +ρ

n∇·un

) (4.6)

2. 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.

3. A Poisson equation for the pressure is solved with a direct solver

∇

2Hn = −

(∇·u)

(n+1)e −(∇·u)n

δt

−∇·F

n

(4.7)

Note that the vector F contains the convective, diffusive and force terms of the momentum equation.

These will be described in detail below.

4. The velocity is estimated at the next time step

u

(n+1)e = un −δt(Fn +∇Hn

) (4.8)

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

|u|

δx

,

|v|

δy

,

|w|

δz



< 1 ; 2 δt ν



1

δx

2

+

1

δy

2

+

1

δz

2



< 1 (4.9)

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 4.1.3 for more details on stability.

5. The thermodynamic quantities ρ, Yα, and pm 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

)



(4.10)

6. The corrected divergence (∇·u)

(n+1)

is computed from the corrected thermodynamic quantities.

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

(4.11)

8. The velocity is corrected

u

n+1 =

1

2

h

u

n +u(n+1)e −δt



F

(n+1)e +∇H(n+1)e

i (4.12)

Note again that the divergence of the corrected velocity field is identically equal to the corrected

divergence.

48
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 their

appropriate cell faces. For example, u

n

i jk is the x-component of velocity at the positive-oriented face of the

i jkth cell; u

n

i−1, jk is defined at the negative-oriented face of the same cell.

49
4.1.2 The Mass Transport Equations and the Divergence

Due to the use of the low Mach number approximation, the mass and energy equations are combined through

the divergence. The divergence of the flow field contains many of the fire-specific source terms described

above.

Discretizing the Convective and Diffusive Transport Terms

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 (4.13)

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 (4.14)

The species conservation equations are differenced the same way, with the addition of diffusion and source

terms:

(ρYα)

(n+1)e

i jk −(ρYα)

n

i jk

δt

+··· = ···+ (∇· ρDα∇Yα)

n

i jk +m ̇

000

α,i jk (4.15)

at the predictor step, and

(ρYα)

(n+1)

i jk −

1

2



(ρYα)

n

i jk + (ρYα)

(n+1)e

i jk 

1

2

δt

+··· = ···+ (∇· ρDα∇Yα)

(n+1)e

i jk +m ̇

000

α,i jk (4.16)

at the corrector step.

The convective terms are written as upwind-biased differences in the predictor step and downwindbiased differences in the corrector step [56]. In the definition 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

(4.17)

The convective term in the species transport equation, (u·∇ρYα)i jk is differenced the exact same way. 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 at the predictor step. 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.

50
Flux Correction

The second-order finite-differencing scheme used in FDS cannot fully resolve sharp gradients on a relatively

coarse grid. Instead, steep gradients cause local over-shoots and under-shoots of quantities like temperature,

density and species mass fraction. For mass fraction in particular, this can result in a solution where the mass

fraction exceeds its permissible limits (i.e. the numerical method transports into or out of a grid cell more

mass than is physically possible). The overall numerical scheme is still mass conserving, but non-physical in

regions of high gradients. This problem can be reduced, but not eliminated, by using higher order numerical

methods. These methods are, however, more expensive. Another solution is to perform a flux transport

correction. This involves examining the solution and locating regions where a non-physical solution exists

and then redistributing mass to correct it. Typically this results in some increased numerical diffusion;

however, this is partially mitigated since at any one time step, the correction is applied to a small number of

cells.

The flux correction scheme is performed in both the predictor and corrector steps after updating the

species mass fractions, Yα. For each species, two loops are performed over each computational mesh. The

first loop searches for and corrects under-shoots, and the second loop searches for and corrects over-shoots.

An under-shoot occurs if the mass fraction of a species is less than its permissible minimum, typically 0, or

if there was an outflux of that species in the prior time step, ∇ ρYα,i jk < 0 , and the new species mass fraction

is less than all of its surrounding cells. An over-shoot occurs if the mass fraction of a species is greater than

its permissible maximum, or if there was an influx of that species in the prior time step, ∇ ρYα,i jk > 0 , and

the new species mass fraction is greater than all of its surrounding cells. In each loop a temporary array is

used to store the corrected values which are then applied globally at the end of each loop. Using a temporary

array rather than a cell by cell immediate correction ensures that cells requiring correction are not bypassed

due to the sweep direction of the loop. Under-shoots are corrected first because in a typical mixture fraction

computation it is more likely to have an under-shoot with a value less than the absolute minimum for the

species than it is for an over-shoot to exceed the absolute maximum for the species.

Discretizing the Divergence

The divergence (see Eq. (3.17)) in the mth pressure zone in both the predictor and corrector step is discretized

(∇·u)i jk =

R

W cp pm



q ̇

000

i jk + (∇· k∇T)i jk +...

+

1

pn



R

W cp

−1

∂pm

∂t

−wi jkρ0,kg



(4.18)

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

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



(4.19)

The temperature is extracted from the density via the equation of state

Ti jk =

pm

ρi jkR ∑

Ns

l=0

(Yα,i jk/Wα)

(4.20)

51
Because only species 1 through Ns are explicitly computed, the summation is rewritten

W ≡

Ns

∑

α=0

Yα,i jk

Wα

=

1

W0

+

Ns

∑

α=1



1

Wα

−

1

W0



Yα (4.21)

In isothermal calculations involving multiple species, the density can be extracted from the average molecular weight

ρi jk =

pm

T∞R W

=

W0 pm

T∞R

+

Ns

∑

α=1



1−

W0

Wα



(ρYα)i jk (4.22)

52
4.1.3 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



(4.23)

∂v

∂t

+Fy +

∂H

∂y

= 0 ; Fy = uωz −wωx −

1

ρ



fy +

∂τyx

∂x

+

∂τyy

∂y

+

∂τyz

∂z



(4.24)

∂w

∂t

+Fz +

∂H

∂z

= 0 ; Fz = vωx −uωy −

1

ρ



fz +

∂τzx

∂x

+

∂τzy

∂y

+

∂τzz

∂z



(4.25)

The spatial discretization of the momentum equations takes the form

∂u

∂t

+Fx,i jk +

Hi+1, jk −Hi jk

δx

= 0 (4.26)

∂v

∂t

+Fy,i jk +

Hi, j+1,k −Hi jk

δy

= 0 (4.27)

∂w

∂t

+Fz,i jk +

Hi j,k+1 −Hi jk

δz

= 0 (4.28)

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, τyy, and τzz; the external force components (fx, fy, fz); and the Courant numbers εu, εv, and εw are

located at their 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



−

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



(4.29)

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



(4.30)

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



(4.31)

ωx,i jk =

wi, j+1,k −wi jk

δy

−

vi j,k+1 −vi jk

δz

(4.32)

53
ωy,i jk =

ui j,k+1 −ui jk

δz

−

wi+1, jk −wi jk

δx

(4.33)

ωz,i jk =

vi+1, jk −vi jk

δx

−

ui, j+1,k −ui jk

δy

(4.34)

τxx,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



(4.35)

τyy,i jk = μi jk 

4

3

(∇·u)i jk −2

ui jk −ui−1, jk

δx

−2

wi jk −wi j,k−1

δz



(4.36)

τzz,i jk = μi jk 

4

3

(∇·u)i jk −2

ui jk −ui−1, jk

δx

−2

vi jk −vi, j−1,k

δy



(4.37)

τ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



(4.38)

τ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



(4.39)

τ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



(4.40)

εu =

uδt

δx

(4.41)

εv =

vδt

δy

(4.42)

εw =

wδt

δz

(4.43)

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. By construction, the divergence defined in Eq. (4.18) 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

(4.44)

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.

Viscous Terms (LES)

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| (4.45)

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

(4.46)

The quantity |S| consists of second order spatial differences averaged at cell centers. For example

∂u

∂x

≈

ui jk −ui−1, jk

δxi

(4.47)

∂u

∂y

≈

1

2

 

ui, j+1,k −ui jk

δyj+

1

2

+

ui jk −ui, j−1,k

δyj−

1

2

!

(4.48)

54
The divergence is described in Section 4.1.2.

The thermal conductivity and material diffusivity of the fluid are related to the viscosity by

ki jk =

cp,0 μi jk

Prt

; (ρD)i jk =

μi jk

Sct

(4.49)

where Prtis the turbulent Prandtl number and Sctis the turbulent 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, Prt and Sct are 0.5. There are no rigorous justifications for these choices other

than through comparison with experimental data [57].

Viscous Terms (DNS)

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) (4.50)

ki jk = ∑

l

Yl,i jk kl(Ti jk) (4.51)

Dl,i jk = Dl0(Ti jk) (4.52)

where the values for each individual species are approximated from kinetic theory [32]. 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.

Velocity Boundary Conditions at a Solid Surface

At solid boundaries, the tangential components of velocity need to be set inside the solid in order to establish

a velocity gradient for the numerical velocity field. As an example, consider a solid surface forming the

lower boundary of the computational domain (z = 0). The horizontal velocity components ui j,1 and vi j,1 are

computed at the vertical faces of the first layer of gas phase grid cells above the solid surface. What should

the values of ui j,0 and vi j,0 be set to such that the flow near the wall is impeded by the same amount as a real

turbulent flow? Assuming that the numerical grid is too coarse to resolve the boundary layer (of nominal

thickness d), and assuming a turbulent velocity profile

U(z) = 

U0(z/d)

1/7

z < d

U0 z > d

(4.53)

a simple slip boundary condition can be worked out such that

ui j,0 = ui j,1



1−

d

δz



; vi j,0 = vi j,1



1−

d

δz



(4.54)

where δz is the height of the grid cell. For coarse grids, the boundary condition is very similar to a “free-slip”

condition because d << δz.

Past and present versions of FDS have employed a default boundary condition for LES calculations

where ui j,0 and vi j,0 are set to one-half the values of ui j,1 and vi j,1, respectively. Although from the discussion

55
above this implies that the thickness of the boundary layer d is assumed to be half that of the grid cell δz, in

practice, it really implies that real walls are covered with small obstructions that are half a cell thick or less

in size, in which case they do not appear explicitly in the simulation. The “half-slip” condition is an attempt

to account for this wall clutter.

In a DNS calculation, this discussion is moot. The inner-wall “ghost” cell velocities are set so that the

velocity at the wall surface is zero.

Force Terms

The external force term components are defined at cell faces, just like the velocity components. For example,

the drag force components from sprinkler droplets are:

fx,i jk =

1

2

∑ρCDπr

2

d

(ud −ui jk)|ud −u|

δxδyδz

(4.55)

fy,i jk =

1

2

∑ρCDπr

2

d

(vd −vi jk)|ud −u|

δxδyδz

(4.56)

fz,i jk =

1

2

∑ρCDπr

2

d

(wd −wi jk)|ud −u|

δxδyδz

(4.57)

where 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 (4.58)

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

restricts the time step.

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 (4.59)

Note that this constraint is 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. To save CPU time, the Von Neumann criterion is only

invoked for DNS calculations or for LES calculations with grid cells smaller than 5 mm.

56
4.1.4 The Pressure Equation

The momentum equation is

∂u

∂t

+F+∇H = 0 ; H =

|u|

2

2

+

p ̃

ρ

(4.60)

The divergence operator applied to the three components of the discretized momentum equation (4.26–

4.28) yields a single elliptic partial differential equation for the modified pressure, H , known as the Poisson

equation:

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 (4.61)

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 [58] that is part of a

library of routines for solving elliptic PDEs called CRAYFISHPAK1. To ensure that the divergence of the

fluid is consistent with the definition given in Eq. (3.17), the time derivative of the divergence is defined

∂

∂t

(∇·u)i jk =

(∇·u)

(n+1)e

i jk −(∇·u)

n

i jk

δt

(4.62)

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

(4.63)

at the corrector step. The discretization of the divergence was given in Eq. (4.18).

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, a no-flux condition 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 (4.64)

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 (4.64) 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, we prescribe

Fn = −

∂H

∂n

∗

+

un

δt

(4.65)

where Fn is the normal component of F at the wall.

The asterisk indicates the most recent value of the pressure. The pressure at this particular time step

is not known until the Poisson equation is solved. Equation (4.65) asserts that following the solution of

1CRAYFISHPAK, a vectorized form of the elliptic equation solver FISHPAK, was originally developed at the National Center

for Atmospheric Research (NCAR) in Boulder, Colorado.

57
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

(4.66)

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 (4.67)

Multiple Mesh Considerations (On-Going Research)

Solving the Poisson equation on muliple meshes is considerably more difficult than on a single mesh. On a

single mesh, the linear system of equations involves a simple block tri-diagonal matrix for which efficient

solvers have been developed. For multiple mesh calculations, the same solver is used for each mesh as

described above, with the following (Dirichlet) boundary conditions applied at mesh interfaces (in this case,

the interface x = xmax = x

+):

HI+

1

2

, jk = HI, jk −

δx

2



u

+

0, jk−uI, jk

δt +Fx,I jk

uI, jk > 0

HI+

1

2

, jk = H

+

−

1

2

, jk uI, jk < 0

(4.68)

Here, the term u

+ indicates that the value is from the other mesh. In short, if the flow is incoming, then H

is taken from the other mesh, if the flow is outgoing, then H is chosen so that the value of u at the next time

step is driven towards the value of its counterpart (u

+) in the other mesh. For abutting meshes of the same

cell size, uI, jk and u

+

0, jk overlap each other, and the goal is to keep their values as close as possible.

This procedure is fairly stable, and provides a reasonable exchange of momentum mesh to mesh. However, there is no guarantee that the individual solutions of the Poisson equation on each mesh lead to consistent flows into and out of neighboring meshes. Using similar boundary conditions as described above

for open boundaries, it is possible to ensure a continuous solution of the pressure term, H , across mesh

boundaries, but it is not possible to ensure a continuous gradient. Thus, it is not possible to ensure that the

flow out of one mesh is the same as the flow into the next.

The strategy of approximating the solution of the Poisson equation on multiple meshes is as follows.

First, decompose the pressure term H = H 0 +H 0. Solve

∇

2H0 = −∇·F−

∂

∂t

(∇·u) (4.69)

on each mesh, the same way as described above. At this point, each mesh has a pressure field, but when

the velocity is updated in each mesh, the normal components will not match at mesh interfaces. At the very

least, we need ensure that the volume flux, V ̇ =

R

u · dS, at each mesh interface is consistent mesh to mesh.

What should the volume flux be? There is no way to get it locally. Instead, re-discretize the original Poisson

58
equation over a coarse mesh spanning the entire domain, and apply the divergence theorem to the coarse

grid cells:

Z

∂Ωm

∇H · dS = −

Z

∂Ωm

F· dS−

Z

Ωm

∂(∇·u)

∂t

dV (4.70)

The subscript m refers to individual cells of the coarse mesh. Introduce a single scalar H m for each coarse

grid cell that satisfies:

H m+ −H m

δx

Ax+ −

H m −H m−

δx

Ax− + ··· = −

Z

x+

Fx dy dz+

Z

x−

Fx dy dz+ ··· −Z

Ωm

∂(∇·u)

∂t

dV (4.71)

The subscripts m

+ and m− refer to coarse neighboring cells, and Ax+ and Ax− and so on refer to the area of

the coarse cell interfaces. At a mesh interface, say an x interface, there are two values of Fx. Denote the

average of the two values Fx, and use this in the coarse Poisson solve.

Solve the coarse grid Poisson equation with whatever sparse method that is best. At the interface between two neighboring meshes, we now have an estimate of the rate of change of the volume flux. For

example, at an x interface:

dV ̇

dt









x+

= −

Z

x+

Fx dy dz−

H m+ −H m

δx

Ax+ (4.72)

We want the fine grid solution of the Poisson equation on the individual abutting meshes to satisfy this

constraint at their interface. To force this, solve

∇

2H0 = 0 (4.73)

with Neumann boundary conditions on each of the abutting meshes such that at the interface the gradient of

H 0is constant. For example, at the x

+ interface, the gradient of H 0 on the mesh to the left must satisfy:

dV ̇

dt









x+

=

Z

x+

 

Fx,I jk +

H 0

I+1, jk −H 0I jk

δx

!

dy dz+

Z

x+

∂H 0

∂x

dy dz (4.74)

A similar constraint must be satisfied by the mesh to the right. In this way, the resulting volume fluxes at the

next time step V ̇ =

R

u·dS into and out of each mesh match. It does not guarantee, however, that the normal

component of velocity at each cell matches its counterpart in the other mesh.

The question now is – what values should be assigned to ∇H 0at each cell of each mesh interface? The

only constraint is that the integral R∇H 0· dS at each face satisfy an equation like (4.74) above. Various

strategies have been tried, but as yet no stable, robust strategy has been found. The research on this problem

is on-going.

59
4.2 Combustion

4.2.1 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. The mixture fraction, Z, is partitioned into at least

2 components, Zα, such that the sum of the components equals the mixture fraction. Each component is

tracked via a transport equation, and the conversion of mass from one component to another represents a

reaction step and an associated release of energy.

Two-Parameter Mixture Fraction

When the mixture fraction is divided into two components, Z1 and Z2, there is one chemical reaction that

converts Z1 to Z2. Recall from Section 3.2.1 that this represents single-step combustion with the possibility

of local extinction. From the mixture fraction variables one can determine the amount of fuel and oxygen

present in a grid cell. Ideally, we could use these values in a finite-rate computation to determine the heat

release rate; however, for most computations, the grid resolution is too coarse to resolve the flame. Thus,

flame temperatures will not be realized and a finite-rate computation will not succeed.

Instead, if any grid cell surrounding one containing both fuel and oxygen satisfies the “Burn” criteria

depicted in Fig. 3.2, combustion is assumed to occur. The heat release rate is computed as

q ̇

000 = min

max(ρYF,sρYO2

)

δt

∆H , q ̇

000

max

; s =

WF

νO2WO2

(4.75)

where  ̇q

000

max provides an upper bound of the volumetric heat release rate. It is an empirical parameter that is

tied to the resolution of the underlying grid:

q ̇

000

max =

q ̇

00

max

(δxδyδz)

1

3

;  ̇q

00

max = 200 kW/m2

(4.76)

The 200 kW/m2 value is a user-controlled parameter. It is an estimate of the maximum heat release rate per

unit area of flame sheet.

Once the heat release rate is computed, the mixture fraction variables are updated:

Z1

n+1 = Z1n −

q ̇

000∆t

ρ∆H

; Z2

n+1 = Z2n +

q ̇

000∆t

ρ∆H

(4.77)

Note that the total fuel mass is conserved in this process; Z = Z1 +Z2 is still a conserved quantity.

Three-Parameter Mixture Fraction

When the mixture fraction is divided into three components, Z1, Z2, and Z3, there are two chemical reactions

that convert Z1 to Z2 and Z2 to Z3. Recall from Section 3.2.1 that this represents two-step combustion

(fuel to CO and CO to CO2). The first step occurs as it does for the two-parameter mixture fraction with

a modified heat of combustion that accounts for the conversion of fuel to CO rather than CO2. The second

step is performed for all grid cells that contain CO and O2. If  ̇q

000 6= 0 in a grid cell after the first step, then

additional heat is released according to

q ̇

000

CO = min

max(ρZ2,sρYO2)

δt

∆HCO , q ̇

000

max −q ̇

000

(4.78)

60
If  ̇q

000 = 0 after the first step, then it is presumed that the cell is out of the combustion region (say in the upper

layer of smoke-filled compartment), and a finite-rate reaction computation is performed to convert CO to

CO2 (see the next section for a discussion of the algorithm for computing a finite-rate reaction). The  ̇q

000

CO

computed using the finite-rate reaction is still limited by  ̇q

000

max. Once  ̇q

000

CO is computed the mixture fraction

variables are updated:

Z2

n+1 = Z2

n −

q ̇

000

CO∆t

ρ∆HCO

; Z3

n+1 = Z3

n +

q ̇

000

CO∆t

ρ∆HCO

(4.79)

4.2.2 Heat Release Rate (Finite-Rate Reaction)

For simulations that do not employ a mixture fraction-based combustion model, multi-step, finite-rate reactions are assumed, taking the form:

νCxHy CxHy +νO2 O2 −→ νCO2 CO2 +νH2O H2O (4.80)

It is assumed that 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/ρi jkWF ≡ XF(t

n

) and Y

n

O2,i jk/ρi jkWO2 ≡ XO2

(t

n

), the following set of ODEs is

solved numerically with a second-order Runge-Kutta scheme

dXF

dt

= −B XF(t)

a XO2

(t)

b

e

−E/RTi jk (4.81)

dXO2

dt

= −

νO2

νF

dXF

dt

(4.82)

The temperature Ti jk and density ρi jk are fixed at their values at time t

n

and the ODEs are iterated from t

n

to

t

n+1

in about 20 time steps. The pre-exponential factor B, the activation energy E, and the exponents a and

b are input parameters which are in typically assigned the values of νF and νO2

. At the end of each sub-time

step the values of XF and XO2are updated.

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

(4.83)

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

α,i jk = Yα(t

n

)−

ναWα

νFWF



YF(t

n

)−YF(t

n+1

)



(4.84)

If multiple chemical reactions have been specified, equations 4.81 and 4.82 are evaluated for each reaction

during each of the 100 time steps. The reactions are evaluated in the order that they are entered in the input

file.

61
4.3 Thermal Radiation

This section describes how ∇·q ̇

00

r

(the radiative loss term) is computed at all gas-phase cells. Plus, it describes

how the radiative and convective heat fluxes,  ̇q

00

r

and  ̇q

00

c

, are computed at solid boundaries.

4.3.1 Discretization of the Radiative Heat Transfer Equation

The radiative transport equation (3.69) is solved using techniques similar to those for convective transport

in finite volume methods for fluid flow [59], 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 Eq. (3.69) over the volume of cell i jk and the

control angle δΩl, to obtain

Z

δΩl

Z

Vi jk

s

0

·∇I(x

0

,s

0

)dx

0

ds

0 =

Z

δΩl

Z

Vi jk

κ(x

0

)



Ib(x

0

)−I(x

0

,s

0

)



dx

0

ds

0

(4.85)

The volume integral on the left hand side is replaced by a surface integral over the cell faces using the

divergence theorem. Note that the procedure outlined below is appropriate for each band of a wide band

model, thus the subscript n has been removed for clarity.

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

0

·nm)ds

0 = κi jk h

Ib,i jk −I

l

i jki

Vi jk δΩl(4.86)

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

δΩlsolid 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

Note that while the intensity is assumed constant within the angle δΩl, its direction covers the angle δΩl

exactly.

In Cartesian coordinates2, 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 obtained 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 the time step set by the flow solver is small.

The grid used 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 4.1. The discretization of the solid angle is done by dividing

first the polar angle, θ, into Nθ bands, where Nθ is an even integer. Each θ-band is then divided into Nφ(θ)

2

In the axisymmetric case equation (4.86) 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.

62
φ

s

x

z

y

θ

Figure 4.1: Coordinate system of the angular discretization.

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) (4.87)

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

Ω

(4.88)

rounded to the nearest even integer. The number of φ-angles on each band is

Nφ(θ) = max4,0.5NΩ



cos(θ

−)−cos(θ+)

	 (4.89)

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. This

symmetry has three important benefits: First, it avoids the problems caused by the fact that the 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 [60], 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 (4.86) 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 (4.86) 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 (4.90)

63
where

a

l

i jk = Ax|D

l

x

|+Ay|D

l

y

|+Az|D

l

z

|+κi jk Vi jkδΩl(4.91)

a

l

x = Ax|D

l

x

| (4.92)

a

l

y = Ay|D

l

y

| (4.93)

a

l

z = Az

|D

l

z

| (4.94)

b

l

i jk = κi jk Ib,i jk Vi jk δΩl

(4.95)

δΩl =

Z

Ωl

dΩ =

Z

δφ

Z

δθ

sinθ dθ dφ (4.96)

D

l

x =

Z

Ωl

(s

l

·i)dΩ (4.97)

=

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Ω (4.98)

=

Z

δφ

Z

δθ

sinφsinθsinθ dθdφ

=

1

2



cosφ

− −cosφ+

 ∆θ−

cosθ

+ sinθ+ −cosθ− sinθ−



D

l

z =

Z

Ωl

(s

l

·k)dΩ (4.99)

=

Z

δφ

Z

δθ

cosθsinθ dθdφ

=

1

2

∆φh

sinθ

+

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 (4.90) is based on an explicit marching sequence [61]. 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.

4.3.2 Radiation Heat Transfer to Solid Objects

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

| (4.100)

64
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

q ̇

00

r =

NΩ

∑

l=1

I

l

w

Z

δΩl

(s

0

·nw)ds

0 =

NΩ

∑

l=1

I

l

wD

l

n

(4.101)

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

(4.102)

Computationally intensive integration over all the incoming directions is avoided by keeping the solid angle

discretization symmetric on the 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

(4.103)

65
4.4 Solid Phase Model

4.4.1 Solid Surface and Interior Temperatures

A one dimensional heat transfer calculation is performed at each solid boundary cell for which the user has

prescribed thermal properties. The solid can consist of multiple layers of materials. Each layer is partitioned

into non-uniform cells, clustered near the front and back faces. The smallest cells are chosen based on the

criteria

δx < Ss

s

ks

ρscs

(4.104)

where Ssis a cell size factor defined by the user. By default, Fsis 1.0. Interior cells increase in size by a

user-defined stretch factor when moving inwards from the surfaces. By default, the stretch factor is 2.0. The

cell boundaries are located at points xi. The temperature at the center of the ith cell is denoted Ts,i. The

(temperature-dependent) thermal conductivity of the solid at the center of the ith cell is denoted ks,i. The

temperatures are updated in time using an implicit Crank-Nicolson 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

!

(4.105)

for 1 ≤ i ≤ N. The width of each cell is δxi. The distance from the center of cell i to the center of cell i+1

is δxi+

1

2

.

The boundary condition is discretized

−ks,1

T

n+1

s,1 −T

n+1

s,0

δx 1

2

= q ̇

00

c

(n+1) +q ̇00

r

(n+1)

(4.106)

The convective and radiative fluxes at the next time step are approximated

q ̇

00

c

(n+1) ≈ q ̇00

c

n −h



T

n+1

s,

1

2

−T

n

s,

1

2



(4.107)

q ̇

00

r

(n+1) ≈ q ̇00

r

n −4 ε σ Tn

3

s,

1

2



T

n+1

s,

1

2

−T

n

s,

1

2



(4.108)

The wall temperature is defined Tw ≡ Ts,

1

2

= (Ts,0 +Ts,1)/2.

4.4.2 Convective Heat Transfer to Solid Objects

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

(4.109)

where Tgas is the gas temperature in the center of the first gas-phase grid cell, Tw is the surface or “wall”

temperature, δn is the size of a grid cell in the normal direction to the wall, and k is the temperaturedependent thermal conductivity of the gas mixture.

In an LES calculation, where the boundary layer is not resolved, an empirical heat transfer coefficient,

h, is used:

q ̇

00

c = h(Tgas −Tw) W/m2

(4.110)

66
The heat transfer coefficient, h, is the maximum of its natural and forced empirical values:

h = maxC|Tgas −Tw|

1

3 ,

kair

L

0.037Re4/5

L

Pr1/3



W/m2/K (4.111)

C is an empirical coefficient (1.31 for vertical surface; 1.52 for horizontal), kair is the thermal conductivity of

air (not the “turbulent” value used in the LES solver), L is a length scale (assumed to be 1 m), the Reynolds

number is based on the density and velocity of the gases in the middle of the first grid cell and the length

scale L, and the Prandtl number is assumed to be 0.7.

4.4.3 Coupling the Gas and Solid Phase

Gas phase temperatures are defined at cell centers; solid surfaces lie at the interface of the bordering gas

phase cell and a “ghost” cell inside the solid. As far as the gas phase calculation is concerned, the normal

temperature gradient at the surface is expressed in terms of the temperature difference between the “gas”

cell and the “ghost” cell. The solid surface temperature is not used directly in the gas phase calculation.

Rather, the ghost cell temperature is used to couple the gas and solid phases. The 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 (4.112)

For an LES calculation, the numerical expression for the heat lost to the boundary is equated with the

empirical convective heat transfer

kLES

Tgas −Tghost

δn

−cpunρ

Tgas +Tghost

2

= h (Tgas −Ts)−cpunρTs (4.113)

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 (4.114)

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

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 (4.115)

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

(4.116)

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

67
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 upwinding 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/Wl)

(4.117)

68
4.5 Liquid Sprays

Lagrangian particles are sometimes introduced into the flow field as a means of visualization or simulation

of droplets, aerosols, etc.. The position xp of each particle is governed by the equations

dxp

dt

= u (4.118)

The particle 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 particle is updated according to the two step scheme

x

(n+1)e

p = x

n

p +δt u

n

(4.119)

x

n+1

p =

1

2



x

n

p +x

(n+1)e

p +δt u

(n+1)e



(4.120)

where the bar over the velocity vector indicates that the velocity of the fluid is linearly interpolated at the

particle’s position.

4.5.1 Interaction of Droplets and Radiation

The computation of χf for a similar but simpler situation has been derived in Ref. [62]. 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(4.121)

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,λ)

(4.122)

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Ω

(4.123)

The limits of the innermost integral are

μd,0 = μμ0 +

p

1−μ

2

p

1−μ

02

; μd,π = μμ0 −

p

1−μ

2

p

1−μ

02

(4.124)

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 [50].

4.5.2 Heating and Evaporation of Liquid Droplets

The exchange of mass and energy between liquid droplets and the surrounding gases (or solid surfaces) is

computed droplet by droplet. After the temperature of each droplet is computed, the appropriate amount of

vaporized liquid is added to the given mesh cell, and the cell gas temperature is reduced slightly based on

the energy lost to the droplet.

69
Equation (3.134) is solved semi-implicitly over the course of a gas phase time step as follows. Note that

a few terms have been left out to make the algorithm more clear.

T

n+1

l −T

n

l

δt

=

1

ml cl

"

A h Tg −

T

n+1

l +T

n

l

2

!

−A hm ρ

 

Y

n+1

l +Y

n

l

2

−Yg

!

hv

#

(4.125)

The equilibrium vapor mass fraction, Y

n

l

, is a function of T

n

l

via Eq. (3.135), and its value at the next time

step is approximated via

Y

n+1

l ≈ Y

n

l +



dYl

dTl

n 

T

n+1

l −T

n

l



(4.126)

where the derivative of Yl with respect to temperature is obtained via the chain rule:

dYl

dTl

=

dYl

dXl

dXl

dTl

=

Wa/Wl

(Xl(1−Wa/Wl) +Wa/Wl)

2

hvWl

R T

2

l

exp

hvWl

R



1

Tb

−

1

Tl

 (4.127)

The amount of evaporated liquid is given by

δml = δt A hm ρ



Y

n

l +

1

2



dYl

dTl

n 

T

n+1

l −T

n

l





(4.128)

The amount of heat extracted from the gas is

δq = δt A h Tg −

T

n

l +T

n+1

l

2

!

(4.129)

70
Chapter 5

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.

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.

71
72
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] K.B. McGrattan, B.W. Klein, S. Hostikka, and J.E. Floyd. Fire Dynamics Simulator (Version 5),

User’s Guide. NIST Special Publication 1019-5, National Institute of Standards and Technology,

Gaithersburg, Maryland, October 2007.

[5] H.C. Hottel. Stimulation of Fire Research in the United States After 1940 (A Historical Account).

Combustion Science and Technology, 39:1–10, 1984.

[6] W.W. Jones. A Review of Compartment Fire Models. NBSIR 83-2684, National Bureau of Standards

(now NIST), Gaithersburg, Maryland, 1983.

[7] J. Quintiere. A Perspective on Compartment Fire Growth. Combustion Science and Technology,

39:11–54, 1984.

[8] G.P. Forney and W.F. Moss. Analyzing and Exploiting Numerical Characteristics of Zone Fire Models.

Fire Science and Technology, 14:49–60, 1994.

[9] S.V. Patankar. Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing, New York, 1980.

[10] 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.

[11] E.S. Oran and J.P. Boris. Numerical Simulation of Reactive Flow. Elsevier Science Publishing Company, New York, 1987.

[12] V. Babrauskas. Ignition Handbook. Fire Science Publishers, Issaquah, Washington USA, 1st edition,

2003. Co-published by the Society of Fire Protection Engineers.

[13] G.P. Forney. User’s Guide for Smokeview Version 5 - A Tool for Visualizing Fire Dynamics Simulation

Data. NIST Special Publication 1017-1, National Institute of Standards and Technology, Gaithersburg,

Maryland, August 2007.

[14] D.A. Anderson, J.C. Tannehill, and R.H. Pletcher. Computational Fluid Mechanics and Heat Transfer.

Hemisphere Publishing Corporation, Philadelphia, Pennsylvania, 1984.

73
[15] R. Peyret and T.D. Taylor. Computational Methods for Fluid Flow. Springer-Verlag, New York, 1983.

[16] J.H. Ferziger and M. Peric. ́ Computational Methods for Fluid Dynamics. Springer-Verlag, Berlin, 2nd

edition, 1999.

[17] R.W. Bilger. Turbulent Diffusion Flames. Annual Review of Fluid Mechanics, 21:101–135, 1989.

[18] J.P. Holman. Heat Transfer. McGraw-Hill, New York, 5th edition, 1989.

[19] F.P. Incropera and D.P. De Witt. Fundamentals of Heat and Mass Transfer. John Wiley and Sons, New

York, 4th edition, 1996.

[20] R. Siegel and J. R. Howell. Thermal Radiation Heat Transfer. Taylor & Francis, New York, 4th

edition, 2002.

[21] P.J. DiNenno, editor. SFPE Handbook of Fire Protection Engineering. National Fire Protection Association, Quincy, Massachusetts, 3rd edition, 2002.

[22] D. Drysdale. An Introduction to Fire Dynamics. John Wiley and Sons, New York, 2nd edition, 2002.

[23] J.G. Quintiere. Principles of Fire Behavior. Delmar Publishers, Albany, New York, 1998.

[24] Verification and Validation of Selected Fire Models for Nuclear Power Plant Applications. NUREG

1824, United States Nuclear Regulatory Commission, Washington, DC, 2007.

[25] 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.

[26] S.B. Pope. Ten questions concerning the large-eddy simulation of turbulent flows. New Journal of

Physics, 6:1–24, 2004.

[27] J. Smagorinsky. General Circulation Experiments with the Primitive Equations. I. The Basic Experiment. Monthly Weather Review, 91(3):99–164, March 1963.

[28] J.W. Deardorff. Numerical Investigation of Neutral and Unstable Planetary Boundary Layers. Journal

of Atmospheric Sciences, 29:91–115, 1972.

[29] 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.

[30] D.K. Lilly. A Proposed Modification of the Germano Subgrid-Scale Closure Method. Physics of Fluids

A, 4(3):633–635, 1992.

[31] 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.

[32] B.E. Poling, J.M. Prausnitz, and J.P. O’Connell. The Properties of Gases and Liquids. McGraw-Hill,

New York, 5th edition, 2000.

[33] 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.

[34] 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.

74
[35] 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.

[36] 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.

[37] G. Heskestad and R.G. Bill. Quantification of Thermal Responsiveness of Automatic Sprinklers Including Conduction Effects. Fire Safety Journal, 14:113–125, 1988.

[38] 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.

[39] 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.

[40] 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.

[41] 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.

[42] G.W. Mulholland and C. Croarkin. Specific Extinction Coefficient of Flame Generated Smoke. Fire

and Materials, 24:227–230, 2000.

[43] R.J. Roby, S.J. Olenick, W. Zhang, D.J. Carpenter, M.S. Klassen, and J.L. Torero. A Smoke Detector

Algorithm for Large Eddy Simulation Modeling. NIST GCR 07-911, National Institute of Standards

and Technology, Gaithersburg, Maryland, 2007.

[44] 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.

[45] 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.

[46] 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.

[47] 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.

[48] T. Ravigururajan and M. Beltran. A Model for Attenuation of Fire Radiation Through Water Droplets.

Fire Safety Journal, 15:171–181, 1989.

[49] A. Tuntomo, C. Tien, and S. Park. Optical Constants of Liquid Hydrocarbon Fuels. Combustion

Science and Technology, 84:133–140, 1992.

75
[50] W.J. Wiscombe. Improved Mie Scattering Algorithms. Applied Optics, 19(9):1505–1509, 1980.

[51] 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.

[52] N. Cheremisinoff. Encyclopedia of Fluid Mechanics, Volume 3: Gas-Liquid Flows. Gulf Publishing

Company, Houston, Texas, 1986.

[53] 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.

[54] K.B. McGrattan, R.G. Rehm, and H.R. Baum. Fire-Driven Flows in Enclosures. Journal of Computational Physics, 110(2):285–291, 1994.

[55] 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.

[56] 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.

[57] 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.

[58] R.A. Sweet. Direct Methods for the Solution of Poisson’s Equation on a Staggered Grid. Journal of

Computational Physics, 12:422–428, 1973.

[59] 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.

[60] 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.

[61] 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.

[62] 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.

76
Appendix A

Nomenclature

As droplet surface area

Aαβ pre-exponential factor for solid phase Arrhenius reaction

B pre-exponential factor for gas phase Arrhenius reaction

C Sprinkler C-Factor

CD drag coefficient

Cs Smagorinsky constant (LES)

cs Solid material specific heat

cp constant pressure specific heat

D diffusion coefficient

dm median volumetric droplet diameter

E activation energy

fb external force vector (excluding gravity)

g acceleration of gravity

g gravity vector, normally (0,0,−g)

H total pressure divided by the density

Hr,αβ heat of reaction for a solid phase reaction

h enthalpy; heat transfer coefficient

hα enthalpy of species α

h

0

α heat of formation of species α

I radiation intensity

Ib radiation blackbody intensity

k thermal conductivity; suppression decay factor

m ̇

00

f

fuel mass flux

m ̇

000

α mass production rate of species α per unit volume

m ̇

00

w water mass flux

m

00

w water mass per unit area

Nu Nusselt number

Pr Prandtl number

p pressure

p0atmospheric pressure profile

pm background pressure of mth pressure zone

p ̃ pressure perturbation

q ̇

00 heat flux vector

q ̇

000 heat release rate per unit volume

77
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

Re Reynolds number

rd droplet radius

rαβ solid phase reaction rate

RTI Response Time Index of sprinkler

s unit vector in direction of radiation intensity

Sc Schmidt number

Sh Sherwood number

Sα solid component production rate

T temperature

t time

U integrated radiant intensity

u = (u,v,w) velocity vector

Wα molecular weight of gas species α

W molecular weight of the gas mixture

We Weber number

x = (x,y,z) position vector

Xα volume fraction of species α

Yα mass fraction of species α

Y

∞

O2 mass fraction of oxygen in the ambient

Y

I

F mass fraction of fuel in the fuel stream

ys soot yield

Z mixture fraction

Zf stoichiometric value of the mixture fraction

γ ratio of specific heats; Rosin-Rammler exponent

∆H heat of combustion

∆HO2

energy released per unit mass oxygen consumed

δ wall thickness

ε dissipation rate

κ absorption coefficient

μ dynamic viscosity

να stoichiometric coefficient, species α

νs yield of solid residue in solid phase reaction

νf yield of gaseous fuel in solid phase reaction

νw yield of gaseous water in solid phase reaction

ρ density

τi j viscous stress tensor

χr radiative loss fraction

σ Stefan-Boltzmann constant; constant in droplet size distribution; surface tension

σd droplet scattering coefficient

σs scattering coefficient

ω = (ωx,ωy,ωz) vorticity vector

78
Appendix B

Derivation of the Velocity Divergence

Randall McDermott, NIST Postdoctoral Fellow

In this appendix we derive the divergence of the velocity field as presented in Eq. (3.17). Note that the

constitutive relationships presented here for the mass diffusion and thermal heat fluxes are valid for direct

numerical simulations (i.e., well-resolved calculations). The minor modifications required of the transport

coefficients for large-eddy simulation are presented in Section 3.1.6. We start the derivation by rearranging

the continuity equation. Next, we differentiate the equation of state to reveal the relationship between

transport equations for mass and energy. We then show how the transport equations may be combined to

yield the velocity divergence constraint. In the last section we present the final result in FDS notation.

Continuity Equation

Let ρ denote the fluid mass density; let u = [u,v,w]

T denote the fluid mass-average velocity; and let  ̇m000

b

denote a bulk source of mass per unit volume (which may come from the evaporation of water droplets, for

example). The continuity equation may be rearranged to yield the following divergence constraint on the

velocity

∇·u =

1

ρ



m ̇

000

b −

Dρ

Dt



(B.1)

where D( )/Dt ≡ ∂( )/∂t +u·∇( ) is the material derivative.

Equation of State

We consider the transport of ns species mass fractions Yα for α = {1,...,ns}, ns−1 of which are independent.

The molecular weight of a given species is denoted Wα and the molecular weight of the mixture, W, is given

by

W =



∑α

Yα

Wα

−1

(B.2)

where as a shorthand notation, which is used throughout this document, we write ∑α for ∑

ns

α=1

. Let pi

(x,t)

denote the hydrostatic pressure in the ith zone of the domain, which in general we take to be a function

of space and time. In practice, however, pi = pi(t) for closed (i.e., sealed or pressurized) domains and

pi = pi(z), where z represents the coordinate aligned with the gravity vector, for large, open domains (e.g.,

forest fires large enough to interact with the stratified atmosphere). The divergence constraint derived below

is based on the ideal gas equation of state (EOS), which, for low-Mach flows, we write as

pi =

ρR T

W

(B.3)

79
where R = 8.3145 kJ/(kmol K) is the gas law constant.

Differentiating the EOS (B.3) we obtain

Dpi

Dt

= ρR T

D

Dt



1

W



+

ρR

W

DT

Dt

+

R T

W

Dρ

Dt

(B.4)

which rearranges to

Dρ

Dt

=

W

R T

Dpi

Dt

−ρW

D

Dt



1

W



−

ρ

T

DT

Dt

(B.5)

Species Transport Equation

The species transport equation plays a role in both the second and third terms on the RHS of (B.5). Including

the bulk mass source, the evolution of species mass fractions is governed by

∂(ρYα)

∂t

+∇·(ρYαu) = −∇· Jα +m ̇

000

α +m ̇

000

b,α

(B.6)

where Jα is the diffusive mass flux vector for species α (relative to the mass-average velocity),  ̇m

000

α

is the

chemical mass production rate of α per unit volume [kg-α produced /(sec m3)], and  ̇m

000

b,α

is the bulk mass

source of α per unit volume [kg-α introduced /(sec m3)]. Note that

∑α

m ̇

000

b,α = m ̇

000

b

(B.7)

and

∑α

m ̇

000

α = 0 (B.8)

Additionally, by construction, the ith component of the species diffusive fluxes sum to zero,

∑α

Jα,i = 0 (B.9)

Thus, as must be the case, summing (B.6) over α yields the continuity equation.

It is convenient to work in terms of the material derivative of the mass fraction. Care must be exercised

in obtaining this expression because the continuity equation is of a non-standard form. Expanding (B.6) we

obtain

ρ

∂Yα

∂t

+Yα

∂ρ

∂t

+ρu·∇Yα +Yα∇·(ρu) = −∇· Jα +m ̇

000

α +m ̇

000

b,α

,

ρ

DYα

Dt

+Yα



∂ρ

∂t

+∇·(ρu)



| {z }

m ̇

000

b

= (B.10)

Thus, the material derivative of the mass fraction can be written as

DYα

Dt

=

1

ρ



m ̇

000

α +m ̇

000

b,α −Yαm ̇

000

b −∇· Jα



=

1

ρ



m ̇

000

α +m ̇

000

b

[Yb,α −Yα]−∇· Jα



(B.11)

where in the second step we use the identity  ̇m

000

b,α = Yαm ̇

000

b with Yb,α being the mass fraction of α in the bulk

prior to its introduction into the fluid mixture.

80
Utilizing (B.2) and (B.11) we obtain

D

Dt



1

W



=

D

Dt



∑α

Yα

Wα



,

= ∑α

1

Wα

DYα

Dt

,

=

1

ρ ∑α

1

Wα



m ̇

000

α +m ̇

000

b

[Yb,α −Yα]−∇· Jα



, (B.12)

which is needed in the second term on the RHS of (B.5).

Enthalpy Transport Equation

The specific sensible enthalpy of species α relative to reference temperature T0 is

hα(T) = Z T

T0

cp,α(T

0

)dT

0

, (B.13)

where the specific heat of α is

cp,α ≡

∂hα

∂T

. (B.14)

The specific sensible enthalpy of the mixture is then given by

hs(Y,T) = ∑αYαhα(T). (B.15)

Neglecting viscous heating and the effect of the fluctuating pressure on dilation work (both assumptions

are valid for low-Mach flows), the transport equation for the sensible enthalpy is

ρ

Dhs

Dt

= −∑α∆h

0

αm ̇

000

α +

Dpi

Dt

−∇·q ̇

00 +m ̇000

b

(hs,b −hs) (B.16)

where ∆h

0

α

is the heat of formation of α at reference temperature T0, hs,b is the specific sensible enthalpy of

the bulk mass source, and q ̇

00 is the heat flux vector which contains contributions from conduction, molecular

diffusion, and radiation,

q ̇

00 = −k∇T +∑α

hs,αJα +q ̇

00

r

(B.17)

Here k is the thermal conductivity of the mixture and q ̇

00

r

is the radiant heat flux.

Relating Enthalpy, Temperature, and Species

Using the chain rule of calculus, we may expand the derivative of the sensible enthalpy hs(Y,T) to obtain

Dhs

Dt

=



∂hs

∂T



DT

Dt

+∑α



∂hs

∂Yα



DYα

Dt

(B.18)

Note that since hs = ∑αYαhα we have

∂hs

∂Yα

=

∂

∂Yα

∑

β

(Yβhβ) = ∑

β

hβ δαβ = hα (B.19)

81
where δαβ is the Kronecker delta. Also,

∂hs

∂T

=

∂

∂T ∑α

Yαhα = ∑αYα



∂hα

∂T



= ∑αYαcp,α ≡ cp (B.20)

defining the specific heat of the mixture. Thus, by rearranging (B.18) and utilizing (B.19) and (B.20) we

obtain

DT

Dt

=

1

cp



Dhs

Dt

−∑αhα

DYα

Dt



(B.21)

Utilizing (B.11) and (B.16) in (B.21) yields

DT

Dt

=

1

ρcp



−∑α∆h

0

αm ̇

000

α +

Dpi

Dt

−∇·q ̇

00 −q ̇000

b +m ̇

000

b

(hs,b −hs)

−∑αhs,α

n

m ̇

000

α +m ̇

000

b

[Yb,α −Yα]−∇· Jα

o

(B.22)

Note that ∑α hs,αm ̇

000

b

[Yb,α −Yα] = m ̇

000

b

(hs,b −hs) and hence these terms cancel, leaving

DT

Dt

=

1

ρcp



−∑α∆h

0

αm ̇

000

α +

Dpi

Dt

−∇·q ̇

00 −q ̇000

b −∑α

hα

n

m ̇

000

α −∇· Jα

o

(B.23)

Assembling Terms

We now have all the pieces we need to construct the divergence constraint which we introduced in Eq. (B.1).

Using (B.5) in (B.1) we obtain

∇·u =

1

ρ



m ̇

000

b −



W

R T

Dpi

Dt

−ρW

D

Dt



1

W



−

ρ

T

DT

Dt



=

1

ρ

m ̇

000

b −

1

pi

Dpi

Dt

+W

D

Dt



1

W



+

1

T

DT

Dt

(B.24)

where in the second step the EOS (B.3) is used to simplify the second term on the RHS. Using (B.12) and

(B.23) in (B.24) yields

∇·u =

1

ρ

m ̇

000

b −

1

pi

Dpi

Dt

+ W



1

ρ ∑α

1

Wα



m ̇

000

α +m ̇

000

b

[Yb,α −Yα]−∇· Jα

	



+

1

T



1

ρcp



−∑α

∆h

0

αm ̇

000

α +

Dpi

Dt

−∇·q ̇

00 −q ̇000

b −∑α

hα



m ̇

000

α −∇· Jα

 (B.25)

Note that W ∑α(Yα/Wα) = 1 and also W ∑α(Yb,α/Wα) = W/Wb, where Wb is the molecular weight of

the bulk mixture prior to its introduction into the fluid mixture. Equation (B.25) thus simplifies to

∇·u =



1

ρcpT

−

1

pi



Dpi

Dt

+

1

ρ



m ̇

000

b

W

Wb

+W ∑α

1

Wα



m ̇

000

α −∇· Jα

	



+

1

ρcpT



−∑α∆h

0

αm ̇

000

α −∇·q ̇

00 −q ̇000

b −∑α

hα



m ̇

000

α −∇· Jα

	



(B.26)

82
FDS Notation

The following relationships are used to rearrange (B.26) into the form shown in the FDS Technical Reference

Guide.

We employ the binary form of Fick’s law using mixture-averaged diffusivities Dα as a constitutive

relation for the diffusive flux,

Jα = −ρDα∇Yα (B.27)

Note that summation is not implied over repeated suffixes. The heat release rate per unit volume is defined

by

q ̇

000 ≡ −∑α

∆h

0

αm ̇

000

α

(B.28)

Taking the z direction to be aligned with the gravity vector we have ∂pi/∂z = −ρig, where g = 9.8 m/s2

and ρiis a specified background density for the ith zone. Thus, the material derivative of the background

pressure may be written as

Dpi

Dt

=

∂pi

∂t

−wρig. (B.29)

Hence, utilizing (B.27), (B.28), and (B.29), and noting 1/(ρT) = R /(W pi) from the EOS, for the ith zone

we may write the divergence (B.26) as

∇·u = D +P

∂pi

∂t

(B.30)

where

P =

1

pi



pi

ρcpT

−1



=

1

pi



R

W cp

−1



(B.31)

and

D =

m ̇

000

b

ρ

W

Wb

+

R

W cp pi



q ̇

000 −q ̇000

b −∇·q ̇

00

−Pwρig

+

W

ρ ∑α

∇·(ρDα∇[Yα/Wα])−

R

W cp pi

∑α

hα∇·(ρDα∇Yα) + 1

ρ ∑α



W

Wα

−

hα

cpT



m ̇

000

α

(B.32)

83
84
Appendix C

A Simple Model of Flame Extinction

Frederick W. Mowrer, Department of Fire Protection Engineering, University of Maryland

A diffusion flame immersed in a vitiated atmosphere will extinguish before consuming all the available

oxygen from the atmosphere. The classic example of this behavior is a candle burning within an inverted jar.

This same concept has been applied within FDS to determine the conditions under which the local ambient

oxygen concentration will no longer support a diffusion flame. In this appendix, the critical adiabatic flame

temperature concept is used to estimate the local ambient oxygen concentration at which extinction will

occur.

Consider a control volume characterized by a bulk temperature, Tm, a mass, m, an average specific heat,

cp, and an oxygen mass fraction, YO2. Complete combustion of the oxygen within the control volume would

release a quantity of energy given by:

Q = mYO2



∆H

rO2



(C.1)

where ∆H/rO2has a relatively constant value of approximately 13100 kJ/kg for most fuels of interest for

fire applications.1 Under adiabatic conditions, the energy released by combustion of the available oxygen

within the control volume would raise the bulk temperature of the gases within the control volume by an

amount equal to:

Q = mcp (Tf −Tm) (C.2)

The average specific heat of the gases within the control volume can be calculated based on the composition

of the combustion products as:

cp =

1

(Tf −Tm) ∑α

Z Tf

Tm

cp,α(T)dT (C.3)

To simplify the analysis, the combustion products are assumed to have an average specific heat of 1.1 kJ/kg/K

over the temperature range of interest, a value similar to that of air, the primary component of the products.

The relationship between the oxygen mass fraction within the control volume and the adiabatic temperature

rise of the control volume is evaluated by equating Eqs. (C.1) and (C.2):

YO2 =

cp(Tf −Tm)

∆H/rO2

(C.4)

If the critical adiabatic flame temperature is assumed to have a constant value of approximately 1700 K for

hydrocarbon diffusion flames, as suggested by Beyler,2

then the relationship between the limiting oxygen

1C. Huggett, “Estimation of the Rate of Heat Release by Means of Oxygen Consumption,” Fire and Materials, Vol. 12, pp. 61-

65, 1980.

2C. Beyler, “Flammability Limits of Premixed and Diffusion Flames,” SFPE Handbook of Fire Protection Engineering (3rd

Ed.), National Fire Protection Association, Quincy, MA, 2003.

85
mass fraction and the bulk temperature of a control volume is given by:

YO2,lim =

cp(Tf,lim −Tm)

∆H/rO2

≈

1.1(1700−Tm)

13100

(C.5)

The relationship represented by Eq. (C.5) is shown in Fig. 3.2. For a control volume at a temperature of

300 K, i.e., near room temperature, the limiting oxygen mass fraction would evaluate to YO2,lim = 0.118.

This value is consistent with the measurements of Morehart, Zukoski and Kubota,3although it is about 10

percent lower than their measured values and the limiting oxygen index values determined by Simmons and

Wolfhard.4 This difference would tend to be conservative in terms of indicating flammability of gas mixtures

at oxygen concentrations lower than might occur in practice. Comparisons with additional experimental data

are needed to further validate the concept. Some improvement might result from a more exact analysis of the

average specific heat of the combustion products or from the use of a slightly higher critical adiabatic flame

temperature in the range of 1700 K to 1800 K, as suggested by the limiting oxygen index data of Simmons

and Wolfhard that is summarized in the SFPE Handbook Chapter of Beyler.

3Morehart, J., Zukoski, E., and Kubota, T., “Characteristics of Large Diffusion Flames Burning in a Vitiated Atmosphere,” Third

International Symposium on Fire Safety Science, Elsevier Science Publishers, pp. 575-583, 1991.

4R.F. Simmons and H.G. Wolfhard, Combustion and Flame, Vol. 1, p. 155, 1957.

86