This document has been superseded 

by a later version. For the latest 

version go to the web site: 

http://fire.nist.gov/fds
NISTIR 6783, 2002 Ed.

Fire Dynamics Simulator (Version 3) – Technical Reference Guide

Kevin B. McGrattan

Howard R. Baum

Ronald G. Rehm

Anthony Hamins

Glenn P. Forney

Jason E. Floyd

Simo Hostikka

Kuldeep Prasad
NISTIR 6783, 2002 Ed.

Fire Dynamics Simulator (Version 3) – Technical Reference Guide

Kevin B. McGrattan

Howard R. Baum

Ronald G. Rehm

Anthony Hamins

Glenn P. Forney

Kuldeep Prasad

Fire Research Division

Building and Fire Research Laboratory

Jason E. Floyd

Hughes Associates, Inc.

Baltimore, Maryland

Simo Hostikka

VTT Building and Transport

Finland

November 2002

UNITED STATES OF AMERICA

DEPARTMENT OF COMMERCE

U.S. Department of Commerce

Donald L. Evans, Secretary

National Institute of Standards and Technology

Arden L. Bement, Director
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 analyses

performed using these tools.

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

thermodynamics, combustion, and heat transfer, 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.

i
Contents

1 Introduction 1

2 Hydrodynamic Model 3

2.1 Conservation Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2.2 State, Mass and Energy Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2.3 The Momentum Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.4 Diffusive Terms (LES) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.5 Diffusive Terms (DNS) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

3 Combustion 8

3.1 Mixture Fraction Combustion Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8

3.2 Enhancements to the Mixture Fraction Model . . . . . . . . . . . . . . . . . . . . . . . . . 10

3.3 Finite-rate Reaction (DNS) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

4 Thermal Radiation 13

5 Thermal Boundary Conditions 16

5.1 Convective Heat Transfer to Walls . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

5.2 Pyrolysis Model, Thermally-Thick Solid . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

5.3 Pyrolysis Model, Thermally-Thin Solid . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

5.4 Pyrolysis Model, Liquid Fuels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

6 Sprinklers 19

6.1 Sprinkler Activation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

6.2 Sprinkler Droplet Size Distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

6.3 Sprinkler Droplet Trajectory in Air . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

6.4 Sprinkler Droplet Transport on a Surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

6.5 Mass and Energy Transfer from Droplets . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

6.6 Interaction of Droplets and Radiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

6.7 Fire Suppression by Water . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

7 Numerical Method 25

7.1 Simplified Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

7.2 Temporal Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

7.3 Spatial Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

7.4 Large Eddy vs. Direct Numerical Simulation . . . . . . . . . . . . . . . . . . . . . . . . . 27

7.5 The Mass Transport Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

7.5.1 Convective and Diffusive Transport . . . . . . . . . . . . . . . . . . . . . . . . . . 28

7.6 Discretizing the Divergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

7.6.1 Heat Release Rate (Mixture Fraction) . . . . . . . . . . . . . . . . . . . . . . . . . 29

7.6.2 Heat Release Rate (Finite-Rate Reaction) . . . . . . . . . . . . . . . . . . . . . . . 30

7.7 Thermal Radiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

7.8 Interaction of Droplets and Radiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

7.9 Thermal and Material Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . 34

7.10 The Momentum Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

7.10.1 Force Terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

7.10.2 Time Step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

ii
7.11 The Pressure Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

7.12 Particle Tracking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

8 Conclusion 41

9 Nomenclature 42

References 44

iii
1 Introduction

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

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

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

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

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

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 necessary 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 cool 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. An excellent description of the physical

and mathematical assumptions behind the zone modeling concept is given by Quintiere [2], who chronicles

developments through 1983. Model development since then has progressed to the point where documented

and supported software implementing these models are widely available [3].

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

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

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 governing equations, in particular

the k−ε turbulence model pioneered by Patankar and Spalding [4]. 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. The k − ε model was 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

1
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 the k − ε model 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

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

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

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

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

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

is ultimately limited by the computing power available. Present day 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 roughly 100 ∼ 200.

Unfortunately, the range of length scales that need to be accounted for if all relevant fire processes are to be

simulated is roughly 104 ∼ 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 meters or 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.

2
2 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 [5]. This gives the equations an elliptic character, consistent with

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

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

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

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

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

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

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

processes to some extent.

2.1 Conservation Equations

First, consider the conservation equations of mass, momentum and energy for a thermally-expandable,

multi-component mixture of ideal gases [5]:

Conservation of Mass

∂ρ

∂t

+∇· ρu = 0 (1)

Conservation of Species

∂

∂t

(ρYl) +∇· ρYlu = ∇· ρDl∇Yl +m ̇

000

l

(2)

Conservation of Momentum

ρ



∂u

∂t

+ (u·∇)u



+∇p = ρg+f+∇· τ (3)

Conservation of Energy

∂

∂t

(ρh) +∇· ρhu =

Dp

Dt

−∇·qr +∇· k∇T +∑

l

∇· hlρDl∇Yl (4)

Note that the external force on the fluid, represented by the term f in Eq. (3), consists of the drag exerted by

water droplets emanating from sprinklers plus other external forces. The term Dp/Dt = ∂p/∂t +u ·∇p is a

material derivative. All other symbols are listed in the Nomenclature (Section 9).

2.2 State, Mass and Energy Equations

The conservation equations are supplemented by an equation of state relating the thermodynamic quantities.

An approximation to the ideal gas law is made by decomposing the pressure into a “background” component,

a hydrostatic component, and a flow-induced perturbation

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

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

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

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

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

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

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

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

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

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

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

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

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

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

the expression for the flow divergence, an important quantity in the analysis to follow.

The divergence of the flow is obtained by taking the material derivative of the equation of state, and

then substituting terms from the mass and energy conservation equations. First, define the constant-pressure

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

l. Next, define the enthalpy h = ∑l hl Yl where

hl(T) = h

0

l +

Z T

T

0

cp,l(T

0

) dT0(7)

and h

0

l

is the heat of formation of species l. Now the divergence can be written

∇·u =

1

ρcp T

 

∇· k∇T +∇·∑

l

R

cp,ldT ρDl∇Yl −∇·qr

!

+

M

ρ ∑

l

∇· ρDl∇(Yl/Ml)−

1

ρcp T ∑

l

R

cp,ldT ∇· ρDl∇Yl +

1

ρ ∑

l



M

Ml

−

hl

cp T



m ̇

000

l +



1

ρcp T

−

1

p0



d p0

dt

(8)

This expression can be simplified by making some approximations. First, assume that Rcp,ldT ≈ cp,l T.

Further, assume the specific heat can be expressed in terms of the number of internal degrees of freedom νl

active in the molecule.

cp,l =



2+νl

2



R

Ml

=



γl

γl −1



R

Ml

(9)

If the ratio of specific heats γl for each species is assumed to be constant, the second line of Eq. (8) disappears, and the only term left from the production term in line 3 is

1

ρ ∑

l



M

Ml

−

hl

cp T



m ̇

000

l = −

1

ρcp T ∑

l

h

0

l m ̇

000

l

(10)

which can be regarded as the energy due to the reaction. From here on, the reaction energy release rate per

unit volume will be written  ̇q

000 = −∑l h0

l m ̇

000

l

.

The approximate form of the divergence used in the calculation is

∇·u =

1

ρcp T

 

∇· k∇T +∇·∑

l

R

cp,ldT ρDl∇Yl −∇·qr +q ̇

000!

+



1

ρcp T

−

1

p0



d p0

dt

(11)

4
Notice that the assumption of a temperature-independent specific heat was made only to eliminate minor

terms in the divergence expression, and thereby reduce the cost of the calculation. In general, it is not

assumed that the specific heat is independent of temperature. The pressure rise term on the right hand side

of the divergence expression is non-zero only if it assumed that the enclosure is tightly sealed, in which case

the background pressure p0 can no longer be assumed constant due to the increase (or decrease) in mass

and thermal energy within the enclosure. The evolution equation for the pressure is found by integrating

Eq. (11) over the entire domain Ω

d p0

dt

=

Z

Ω

1

ρcp T



∇· k∇T +...dV −

Z

∂Ω

u· dS

,Z

Ω



1

p0

−

1

ρcp T



dV (12)

2.3 The Momentum Equation

The momentum equation is simplified by subtracting off the hydrostatic pressure gradient from the momentum equation (3), and then dividing by the density to obtain1

∂u

∂t

−u×ω+

1

2

∇|u|

2 +

1

ρ

∇p ̃ =

1

ρ

((ρ−ρ∞)g+f+∇· τ) (13)

To simplify this equation further, a substitution is made

∇H ≈

1

2

∇|u|

2 +

1

ρ

∇p ̃ (14)

The basis for this approximation is seen in the evolution equation for the circulation, obtained by integrating

Eq. (13) over a closed loop moving with the fluid (in the absence of any external force)

dΓ

dt

=

I

1

ρ

(−∇p ̃ + (ρ−ρ∞)g+∇· τ)· dx (15)

There are three sources of vorticity: the baroclinic torque due to the non-alignment of the density and

pressure gradients, buoyancy due to horizontal density gradients, and viscosity. Buoyancy is the dominant

source of vorticity, and the approximation above is equivalent to neglecting the baroclinic torque.2

Neglecting the baroclinic torque simplifies the elliptic partial differential equation obtained by taking

the divergence of the momentum equation

∇

2H = −

∂(∇·u)

∂t

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

1

ρ

((ρ−ρ∞)g+f+∇· τ) (16)

The linear algebraic system arising from the discretization of Eq. (16) has constant coefficients and can be

solved to machine accuracy by a fast, direct (i.e. non-iterative) method that utilizes fast Fourier transforms.

No-flux or forced-flow boundary conditions are specified by asserting that

∂H

∂n

= −Fn −

∂un

∂t

(17)

1Note the use of the vector identity (u·∇)u =

1

2∇|u|

2 −u×ω.

2 An option exists in the code to restore the baroclinic torque by decomposing the pressure term

∇p ̃

ρ

=

∇p ̃

ρ ̄

+



1

ρ

−

1

ρ ̄



∇p ̃

and evaluating the second term on the right hand side at the previous time step. The expression ρ ̄ is an average density, equal to

2ρminρmax/(ρmin +ρmax). For most large-scale applications, the baroclinic torque is relatively small compared to buoyancy.

5
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 (18)

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.

2.4 Diffusive Terms (LES)

The viscous stress tensor in the momentum equation is given by

τ = μ



2 def u−

2

3

(∇·u)I



(19)

where I is the identity matrix and the deformation tensor is defined

def u ≡

1

2



∇u+ (∇u)

t



=





∂u

∂x

1

2



∂u

∂y +

∂v

∂x



1

2



∂u

∂z +

∂w

∂x



1

2



∂v

∂x +

∂u

∂y



∂v

∂y

1

2



∂v

∂z +

∂w

∂y



1

2



∂w

∂x +

∂u

∂z



1

2



∂w

∂y +

∂v

∂z



∂w

∂z





(20)

In the numerical model, there are two options for treating the dynamic viscosity μ. For a Large Eddy

Simulation (LES) where the grid resolution is not fine enough to capture the mixing processes at all relevant

scales, a sub-grid scale model for the viscosity is applied. Following the analysis of Smagorinsky [7], the

viscosity can be modeled as

μLES = ρ(Cs ∆)

2



2(def u)·(def u)−

2

3

(∇·u)

2

1

2

(21)

where Csis an empirical constant, ∆ is a length on the order of the size of a grid cell, and the deformation

term is related to the Dissipation Function

Φ ≡ τ ·∇u ≡ μ



2(def u)·(def u)−

2

3

(∇·u)

2



(22)

= μ

"

2



∂u

∂x

2

+2



∂v

∂y

2

+2



∂w

∂z

2

+ (23)



∂v

∂x

+

∂u

∂y

2

+



∂w

∂y

+

∂v

∂z

2

+



∂u

∂z

+

∂w

∂x

2

−

2

3



∂u

∂x

+

∂v

∂y

+

∂w

∂z

2

#

(24)

The dissipation function is the rate at which kinetic energy is transferred to thermal energy. It is a source

term in the energy conservation equation that is usually neglected because it is small – an approximation

consistent with the low Mach number equations.

In an LES calculation, the thermal conductivity and material diffusivity are related to the turbulent

viscosity by

kLES =

μLES cp

Pr

; (ρD)l,LES =

μLES

Sc

(25)

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

There have been numerous refinements of the original Smagorinsky model [8, 9, 10], but it is difficult to

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

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

gives results almost identical to those obtained using the Smagorinsky model [11]. 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.

2.5 Diffusive Terms (DNS)

For a Direct Numerical Simulation (DNS), the viscosity, thermal conductivity and material diffusivity are

approximated from kinetic theory. The viscosity of the lth species is given by

μl =

26.69×10−7(Ml T)

1

2

σ

2

l Ωv

kg

m s

(26)

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

kl =

μl cp,l

Pr

W

m K

(27)

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

μDNS = ∑

l

Yl μl; kDNS = ∑

l

Yl kl (28)

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

Dlm =

2.66×10−7 T

3/2

M

1

2

lm σ

2

lm ΩD

m2

s

(29)

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

−1

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

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

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

the given species diffusing into nitrogen

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

where species 0 is nitrogen.

7
3 Combustion

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

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

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

is used.

3.1 Mixture Fraction Combustion Model

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

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

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

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

(often limited) understanding of the phenomena involved.

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

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

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

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

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

conserved quantity representing the fraction of material at a given point that originated as fuel. The relations

between the mass fraction of each species and the mixture fraction are known as “state relations”. The state

relation for the oxygen mass fraction provides the information needed to calculate the local oxygen mass

consumption rate. The form of the state relation that emerges from classical laminar diffusion flame theory

is a piecewise linear function. This leads to a “flame sheet” model, where the flame is a two dimensional

surface embedded in a three dimensional space. The local heat release rate is computed from the local

oxygen consumption rate at the flame surface, assuming that the heat release rate is directly proportional

to the oxygen consumption rate, independent of the fuel involved. This relation, originally proposed by

Huggett [13], is the basis of oxygen calorimetry.

Start with the most general form of the combustion reaction

νF Fuel+νO O2 → ∑

i

νP,i Products (31)

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

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

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

m ̇

000

F

νFMF

=

m ̇

000

O

νOMO

(32)

The mixture fraction Z is defined as:

Z =

sYF −(YO −Y

∞

O

)

sYI

F +Y

∞

O

; s =

νOMO

νFMF

(33)

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

on its undepleted ambient value, Y

∞

O

. Note that Y

I

F

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

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

conservation law

ρ

DZ

Dt

= ∇· ρD∇Z (34)

8
FIGURE 1: State relations for propane.

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

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

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

surface as:

Z(x,t) = Zf; Zf =

Y

∞

O

sYI

F +Y

∞

O

(35)

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

fraction YO and Z

YO(Z) = 

Y

∞

O

(1−Z/Zf) Z < Zf

0 Z > Zf

(36)

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

a hydrocarbon fuel:

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

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

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

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

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

be obtained from the right hand side of Eq. (37). State relations for the ideal reaction of propane and air is

shown in Fig. 1.

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

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

9
of the oxygen consumption

q ̇

000 = ∆HO m ̇000

O (38)

Here, ∆HO is the heat release rate per unit mass of oxygen consumed. The oxygen mass conservation

equation

ρ

DYO

Dt

= ∇· ρD∇YO +m ̇

000

O (39)

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

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

−m ̇

000

O = ∇·



ρD

dYO

dZ

∇Z



−

dYO

dZ

∇· ρD∇Z = ρD

d

2YO

dZ2

|∇Z|

2

(40)

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

. However, an expression for the

oxygen consumption rate per unit area of flame sheet can be derived from Eq. (40)

−m ̇

00

O =

dYO

dZ









Z<Zf

ρD ∇Z ·n (41)

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

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

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

cell, consistent with all other gas phase quantities.

3.2 Enhancements to the Mixture Fraction Model

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

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

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

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

∗/δx, where D∗

is a characteristic fire diameter

D

∗ =



Q ̇

ρ∞ cp T∞

√g

2

5

(42)

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

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

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

∗

is small relative

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

Z = Zf will underestimate the observed flame height [14]. It has been found empirically that a good estimate

of flame height can be found for crude grids if a different value of Z is used to define the combustion region

Zf,eff

Zf

= min 1 , C

D

∗

δx



(43)

HereC is an empirical constant to be used for all fire scenarios. As the resolution of the calculation increases,

the Zf,eff approaches the ideal value, Zf. The benefit of the expression is that it provides a quantifiable

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

, where D is the physical

diameter of the fire.

10
measure of the grid resolution that takes into account not only the size of the grid cells, but also the size of

the fire.

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

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

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

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

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

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

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

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

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

H/D = 3.7 Q

∗2/5 −1.02 ; A = πr

p

r

2 +h2

(44)

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

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

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

the upper bound is redistributed over the entire flame volume.

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

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

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

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

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

these phenomena are complex, and even simplified models still rely on a reasonably accurate estimation of

the temperature and local strain rate in the neighborhood of the flame sheet. Sub-grid scale modeling of

gas phase suppression and extinction is still an area of active research in the combustion community. Until

reliable models can be developed for building-scale fire simulations, simple empirical rules can be devised

that prevent burning from taking place when the atmosphere immediately surrounding the fire cannot sustain

the combustion. Based on the work of Quintiere, Mowrer and others, a model for flame extinction has been

implemented in FDS. The mixture fraction continues to be used to track the progress of the fuel mixing with

the surrounding air, but now the surrounding volume is assessed to determine if it is more or less likely to

support combustion. Figure 2 shows values of temperature and oxygen concentration for which burning can

and cannot take place. Note that once the combustion region falls in the “No Burn” zone, the state relations

(Fig. 1) are no longer valid for values of Z below stoichiometric, since now some fuel may be mixed with

the other combustion products. To account for the deviation from the ideal state relations, at least one other

scalar quantity in addition to the mixture fraction would have to be tracked in the calculation.

3.3 Finite-rate Reaction (DNS)

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

fuel

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

The reaction rate is given by the expression

d[CxHy]

dt

= −B[CxHy]

a

[O2]

b

e

−E/RT (46)

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

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

11
FIGURE 2: Oxygen-temperature phase space showing where

combustion is allowed and not allowed to take place.

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

provides a valid description of the combustion process.

12
4 Thermal Radiation

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

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

4π

Z

4π

Φ(s,s

0

) Iλ(x,s

0

) dΩ

0

(47)

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

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

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

non-scattering gas the RTE becomes

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

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

in the gas phase. The interaction of radiation with droplets is explained in Section 6.6.

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

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

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

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

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

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

/π (50)

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

From a series of numerical experiments it has been found that six bands are usually enough (N = 6). If the

absorption of the fuel is known to be important, separate bands can be reserved for fuel, and the total number

of bands is increased to ten (N = 10). For simplicity, the fuel is assumed to be CH4. The limits of the bands

are shown in Table 1.

TABLE 1: Limits of the spectral bands.

9 Band Model 1 2 3 4 5 6 7 8 9

Soot CO2 CH4 Soot CO2 H2O H2O Soot Soot Major Species

H2O, Soot Soot Soot Soot CH4, Soot

ν (1/cm) 10000 3800 3400 2800 2400 2174 1429 1160 1000 50

λ (μm) 1.00 2.63 2.94 3.57 4.17 4.70 7.00 8.62 10.0 200

6 Band Model 1 2 3 4 5 6

Soot CO2 CH4 CO2 H2O, CH4, Soot Soot Major Species

H2O, Soot Soot Soot

Even with a reasonably small number of bands, the solution of N RTE’s is very time consuming. Fortunately, in most large-scale fire scenarios soot is the most important combustion product controlling the

13
thermal radiation from the fire and hot smoke. As the radiation spectrum of soot is continuous, it is possible

to assume that the gas behaves as a gray medium. The spectral dependence is lumped into one absorption

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

Ib(x) = σT(x)

4

/π (52)

In optically thin flames, where the amount of soot is small compared to the amount of CO2 and water, the

gray gas assumption may produce significant overpredictions of the emitted radiation.

For the calculation of the gray or band-mean absorption coefficients, κn, a narrow-band model, RadCal [19], has been implemented in FDS. At the start of a simulation the absorption coefficient(s) are tabulated

as a function of mixture fraction and temperature. During the simulation the local absorption coefficient is

found by table-lookup.

In calculations of limited spatial resolution, the source term, Ib, in the RTE requires special treatment

in the neighborhood of the flame sheet because the temperatures are smeared out over a grid cell and are

thus considerably lower than one would expect in a diffusion flame. Because of its dependence on the

temperature raised to the fourth power, the source term must be modeled in those grid cells cut by the flame

sheet. Elsewhere, there is greater confidence in the computed temperature, and the source term can assume

its ideal value there

κ Ib =



κσT

4/π Outside flame zone

χr q ̇

000/4π Inside flame zone (53)

Here,  ̇q

000 is the chemical heat release rate per unit volume and χr

is the local fraction of that energy emitted

as thermal radiation. Note the difference between the prescription of a local χr and the resulting global

equivalent. For a small fire (D < 1 m), the local χris approximately equal to its global counterpart. However,

as the fire increases in size, the global value will typically decrease due to a net re-absorption of the thermal

radiation by the increasing smoke mantle.

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

Iw(s) = εIbw +

1−ε

π

Z

s

0

·nw<0

Iw(s

0

) |s

0

·nw| dΩ (54)

where Iw(s) is the intensity at the wall, ε is the emissivity, and Ibw is the black body intensity at the wall.

The radiative transport equation (49) is solved using techniques similar to those for convective transport

in finite volume methods for fluid flow [20], thus the name given to it is the Finite Volume Method (FVM).

To obtain the discretized form of the RTE, the unit sphere is divided into a finite number of solid angles. In

each grid cell a discretized equation is derived by integrating equation (48) over the cell i jk and the control

angle δΩl, to obtain

Z

Ωl

Z

Vi jk

s·∇In(x,s)dV dΩ =

Z

Ωl

Z

Vi jk

κn(x) [Ib,n(x)−In(x,s)]dV dΩ (55)

The volume integral on the left hand side is replaced by a surface integral over the cell faces using the

divergence theorem. Assuming that the radiation intensity I(x,s) is constant on each of the cell faces, the

surface integral can be approximated by a sum over the cell faces. More detail on the discretization and

solution of the RTE can be found in Section 7.7.

The radiant heat flux vector qris defined

qr(x) = ZsI(x,s) dΩ (56)

The radiative loss term in the energy equation is

−∇·qr(x) = κ(x) [U(x)−4πIb(x)] ; U(x) = Z

4π

I(x,s)dΩ (57)

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

that which is emitted.

15
5 Thermal Boundary Conditions

The type of thermal boundary condition applied at any given surface depends on whether that surface is to

heat up and burn, whether the burning rate will simply be prescribed, or whether there is to be any burning

at all.

5.1 Convective Heat Transfer to Walls

The heat fluxes to a solid surface consist of gains and losses from convection and radiation. The radiative

flux at the surface is obtained from the boundary condition for the radiation equation, Eq. (54).

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

(58)

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

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.43 for a horizontal surface and 0.95 for

a vertical surface) [21], 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.

5.2 Pyrolysis Model, Thermally-Thick Solid

If the surface material is assumed to be thermally-thick, a one-dimensional heat conduction equation for the

material temperature, Ts(n,t), is applied in the direction n pointing into the air/solid interface (n = 0)

ρs cs

∂Ts

∂t

= ks

∂

2Ts

∂n

2

; −ks

∂Ts

∂n

(0,t) = q ̇

00

c +q ̇

00

r −m ̇

00 ∆Hv (60)

where ρs, cs and ks are the (constant) density, specific heat and conductivity of the material;  ̇q

00

c

is the

convective and  ̇q

00

r

is the (net) radiative heat flux at the surface,  ̇m

00 is the mass loss rate of fuel and ∆Hv is

the heat of vaporization. It is assumed that fuel pyrolysis takes place at the surface, thus the heat required

to vaporize fuel is extracted from the incoming energy flux. The pyrolysis rate is given by an Arrhenius

expression

m ̇

00 = A e−E/RT (61)

The value of the pre-exponential factor, A, is fixed, R is the universal gas constant, and E is adjusted so that

the material burns in the neighborhood of a prescribed temperature. The actual burning rate is governed

by the overall energy balance in the solid. These parameters are often difficult to obtain for real fuels; the

intent of using the given expression for the mass loss rate is to mimic the behavior of burning objects when

details of their pyrolysis mechanisms are unknown. Figure 3 displays a few plots of the mass flux versus

temperature for two different temperatures.

16
FIGURE 3: Fuel pyrolysis rates for two given “ignition” temperatures. See Eq. 61.

5.3 Pyrolysis Model, Thermally-Thin Solid

If the surface material is assumed to be thermally-thin, that is, its temperature is assumed uniform across its

width, Ts(t) is affected by gains and losses due to convection, radiation and pyrolysis. The thermal lag of

the material is a function of the product of its density, specific heat and thickness δ

dTs

dt

=

q ̇

00

c +q ̇

00

r −m ̇

00 ∆Hv

ρscsδ

(62)

The convective and radiative fluxes are summed over the front and back surfaces of the thin fuel. Unless

otherwise specified, the back surface is assumed to face an ambient temperature void. Note that the individual values of the parameters ρs

, cs and δ are not as important as their product, thus often in the literature

and in the computer program, the three values are lumped together as a product. The pyrolysis rate for a

thermally-thin fuel is the same as for a thermally-thick; see Eq. 61.

5.4 Pyrolysis Model, Liquid Fuels

The rate at which liquid fuel evaporates when burning is a function of the liquid temperature and the concentration of fuel vapor above the pool surface. Equilibrium is reached when the partial pressure of the fuel

vapor above the surface equals the Clausius-Clapeyron pressure

pcc = p0 exp−

hvMf

R



1

Ts

−

1

Tb

 (63)

17
where hv is the heat of vaporization, Mf

is the molecular weight, Tsis the surface temperature, and Tb is the

boiling temperature of the fuel [22].

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.

18
6 Sprinklers

Simulating the effects of a sprinkler spray involves a number of pieces: predicting activation, computing the

droplet trajectories and tracking the water as it drips onto the burning commodity.

6.1 Sprinkler Activation

The temperature of the sensing element of a given sprinkler is estimated from the differential equation put

forth by Heskestad and Bill [23], with the addition of several terms to account for radiative heating and

cooling 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| (64)

Here Tlis the link temperature, Tg is the gas temperature in the neighborhood of the link, Tm is the temperature of the sprinkler mount (assumed ambient), and β is the volume fraction of (liquid) water in the gas

stream. The sensitivity of the detector is characterized by the value of RTI. The amount of heat conducted

away from the link by the mount is indicated by the “C-Factor”, C. The constant C2 has been empirically

determined by DiMarzo [24] to be 6×106 K/(m/s) 1

2 , and its value is relatively constant for different types

of sprinklers.

6.2 Sprinkler Droplet Size Distribution

Once activation is predicted, a sampled set of spherical water droplets is tracked from the sprinkler to either

the floor or the burning commodity. In order to compute the droplet trajectories, the initial size and velocity

of each droplet must be prescribed. This is done in terms of random distributions. The initial droplet

size distribution of the sprinkler spray is expressed in terms of its Cumulative Volume Fraction (CVF),

a function that relates the fraction of the water volume (mass) transported by droplets less than a given

diameter. Researchers at Factory Mutual have suggested that the CVF for an industrial sprinkler may be

represented by a combination of log-normal and Rosin-Rammler distributions [25]

F(d) =







√

1

2π

Z d

0

1

σd

0 e

−

[ln(d

0/dm)]2

2σ2 dd0

(d ≤ dm)

1−e

−0.693(

d

dm

)

γ

(dm < d)

(65)

where dm is the median droplet diameter (i.e. half the mass is carried by droplets with diameters of dm

or less), and γ and σ are empirical constants equal to about 2.4 and 0.6, respectively 4. The median drop

diameter is a function of the sprinkler orifice diameter, operating pressure, and geometry. Research at

Factory Mutual has yielded a correlation for the median droplet diameter [26]

dm

D

∝ We−

1

3 (66)

where D is the orifice diameter of the sprinkler. The Weber number, the ratio of inertial forces to surface

tension forces, is given by

We =

ρwU

2D

σw

(67)

where ρw is the density of water, U is the water discharge velocity, and σw is the water surface tension

(72.8 × 10−3 N/m at 20 ◦C). The discharge velocity can be computed from the mass flow rate, which is a

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

√

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

19
FIGURE 4: Cumulative Volume Fraction and Cumulative Number Fraction functions of the droplet

size distribution from a typical industrial-scale sprinkler. The median diameter dm is 1 mm, σ = 0.6

and γ = 2.4.

function of the sprinkler’s operating pressure and K-Factor. FM reports that the constant of proportionality in

Eq. (66) 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 [26].

In the numerical algorithm, the size of the sprinkler droplets are chosen to mimic the Rosin-Rammler/lognormal distribution. A Probability Density Function (PDF) for the droplet diameter is defined

f(d) = F

0

(d)

d

3

Z ∞

0

F

0

(d

0

)

d

03

dd0(68)

Droplet diameters are randomly selected by equating the Cumulative Number Fraction of the droplet distribution with a uniformly distributed random variable U

U(d) = Z d

0

f(d

0

)dd0(69)

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

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

Typically, 1,000 droplets per sprinkler per second are tracked (50 droplets every 0.05 s, depending on user

preference). The procedure for selecting droplet sizes is as follows: Suppose water is leaving the sprinkler at

a mass flow rate of  ̇m. Suppose also that the time interval for droplet insertion into the numerical simulation

is δt, and the number of droplets inserted each time interval is N. Choose N uniformly distributed random

numbers between 0 and 1, call them Ui, obtain N droplet diameters di based on the given droplet size

distribution (Eq. (69), and then compute a weighting constant C from the mass balance

m ̇ δt = C

N

∑

i=1

4

3

πρw



di

2

3

(70)

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

20
6.3 Sprinkler Droplet Trajectory in Air

For a sprinkler spray, the force term f in Eq. (3) represents the momentum transferred from the water droplets

to the gas. It is obtained by summing the force transferred from each droplet in a grid cell and dividing by

the cell volume

f =

1

2

∑ρCDπr

2

d

(ud −u)|ud −u|

δxδyδz

(71)

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

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

CD =







24/Re Re < 1

241+0.15Re0.687

/Re 1 < Re < 1000

0.44 1000 < Re

(73)

Re =

ρ|ud −u|2rd

μ

(74)

where μ is the dynamic viscosity of air.

6.4 Sprinkler Droplet Transport on a Surface

When a water droplet hits a solid horizontal surface, it is assigned a random horizontal direction and moves

at a fixed velocity until it reaches the edge, at which point it drops straight down at the same fixed velocity.

This “dripping” velocity has been measured to be on the order of 0.5 m/s [27]. Penetration of water into

porous materials is handled very crudely by assigning a fraction of the water droplets that strike a solid

horizontal surface to go straight through the solid at a slow velocity. Neither the fraction nor the velocity

has been validated.

6.5 Mass and Energy Transfer from Droplets

The evaporation of water droplets is handled semi-empirically. A water droplet suspended in air will evaporate as a function of the droplet equilibrium vapor mass fraction, the local gas phase vapor mass fraction,

the heat transfer to the droplet, and the droplet’s motion relative to the gas. A correlation for the mass loss

rate of a droplet that involves these parameters is given here [28]

dmd

dt

= −2π rd Sh ρD (Yd −Yg) (75)

The subscripts d and g refer to the droplet and gas, respectively, md is the droplet mass, D is the diffusion

coefficient for water vapor into air, Y is the water vapor mass fraction, and Sh is the droplet Sherwood

number, given by a correlation involving the Reynolds and Schmidt numbers

Sh = 2+0.6 Re

1

2 Sc

1

3 (76)

The vapor mass fraction of the gas, Yg, is obtained from the overall set of mass conservation equations and

the vapor mass fraction of the droplet is obtained from the Clausius-Clapeyron equation

Xd = exp

hv Mw

R



1

Tb

−

1

Td

 ; Yd =

Xd

Xd(1−Ma/Mw) +Ma/Mw

(77)

21
where Xd is the droplet water vapor volume fraction, hv is the heat of vaporization, Mw is the molecular

weight of water, Ma is the molecular weight of air, R is the gas constant, Tb is the boiling temperature of

water and Td is the droplet temperature.

In addition to calculating the mass transfer due to evaporation, the transfer of energy must also be

calculated. The droplet heats up due to the convective heat transfer across the surface of the droplet minus

the energy required to evaporate water

md cp,w

dTd

dt

= Ad hd (Tg −Td)−

dmd

dt

hv (78)

Here cp,w is the specific heat of water, Ad = 4πr

2

d

is the surface area of the droplet, hd is the heat transfer

coefficient, given by

hd =

Nu k

2rd

; Nu = 2+0.6 Re

1

2 Pr

1

3 (79)

Nu is the Nusselt number, k is the thermal conductivity of air, and the Prandtl number, Pr, is about 0.7 for

air. The Sherwood number, Sh, is analagous to the Nusselt number, with the Schmidt number about 0.6

compared to the 0.7 for the Prandtl number.

Finally, the exchange of mass and energy between the droplets and the gas results in an additional term

that must be added to the expression for the divergence, Eq. (8)

∇·u = ... +

R

γ p0



ρ ∑(Yi/Mi)

∂T

∂t

+

T

Mw

m ̇

000

w



(80)

where  ̇m

000

w

is the water evaporation rate per unit volume. The liquid water droplets are assumed to occupy

no volume, simplifying the analysis.

6.6 Interaction of Droplets and Radiation

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

mist systems [29]. Water droplets attenuate thermal radiation through a combination of scattering and absorption [30]. 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. (47) are temporarily neglected for simplicity, the radiative transport equation becomes

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

4π

Z

4π

Φ(s,s

0

) Iλ(x,s

0

) dΩ

0

(81)

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

where r is the droplet radius and Ca and Cs are absorption and scattering cross sections, respectively, given

by Mie theory. The droplet number density function f(r,dm) is assumed to have the same form as the droplet

size distribution, but a different mean.

An accurate computation of the in-scattering intergal on the right hand side of Eq (81) would be extremely time consuming. It is here approximated by dividing the total 4π solid angle to a “forward angle”

22
δΩl

and “ambient angle” δΩ∗ = 4π−δΩl. For compatibility with the FVM solver, δΩlis set equal to the

control angle given by the angular discretization. However, it is assumed to be symmetric around the center

of the control angle. Within δΩlthe intensity is Iλ(x,s) and elsewhere it is approximated as

U

∗

(x,λ) = U(x,λ)−δΩl

Iλ(x,s)

δΩ∗

(83)

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)] (84)

where χf = χf(r,λ) is a fraction of the total intensity originally within the solid angle δΩlthat is scattered

into the same angle δΩl. Defining now an effective scattering coefficient section

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

4π−δΩl

Z ∞

0

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

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

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

the same as what is used for the gas phase RTE. After the band integrations, the spray RTE for band n

becomes

s·∇In(x,s) = −[κd,n(x) +σ ̄ d,n(x)]In(x,s) +κd,n(x) Ib,d,n(x) + σ ̄ d(x,λ)

4π

Un(x) (87)

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

scattering cross sections and the scattering phase function are calculated using the MieV code developed

by Wiscombe [31]. Both κd and σ ̄ d are averaged over the possible droplet radii and wavelength before

the actual simulation. A single constant temperature is used in the wavelength averaging. This “radiation

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

used by default. The averaged quantities, being now functions of the droplet mean diameter only, are saved

in one-dimensional arrays. During the simulation, the local properties are calculated as a table lookup using

the local mean droplet diameter [32]. Details of the computation are given in Section 7.8.

6.7 Fire Suppression by Water

The above two sections describe heat transfer from a droplet of water to a hot gas, a hot solid, or both.

Although there is some uncertainty in the values of the respective heat transfer coefficients, the fundamental

physics are fairly well understood. However, when the water droplets encounter burning surfaces, simple

heat transfer correlations become more difficult to apply. The reason for this is that the water is not only

cooling the surface and the surrounding gas, but it is also changing the pyrolysis rate of the fuel. If the

surface of the fuel is planar, it is possible to characterize the decrease in the pyrolysis rate as a function of

the decrease in the total heat feedback to the surface. Unfortunately, most fuels of interest in fire applications

are multi-component solids with complex geometry at scales unresolvable by the computational grid.

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

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

(88)

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

(89)

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

(90)

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

m ̇

00

f

(t) = m ̇

00

f,0

(t) e

−

R

k(t)dt (91)

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

(92)

Note that a is an empirical constant.

24
7 Numerical Method

This section presents the details of the numerical algorithm. First the equations that are being solved are

presented. Each of the conservation equations emphasize the importance of the velocity divergence and

vorticity fields, as well as the close relationship between the thermally expandable fluid equations [5] and the

Boussinesq equations for which the authors have developed highly efficient solution procedures [34, 35]. All

spatial derivatives are approximated by second order central differences and the flow variables are updated

in time using an explicit second order predictor-corrector scheme.

7.1 Simplified Equations

Regardless of whether one is performing an LES or a DNS calculation, the overall solution algorithm is the

same. The equations derived in Section 2 that are to be solved numerically are listed again here.

Conservation of Mass

∂ρ

∂t

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

Conservation of Species

∂ρYl

∂t

+u·∇ρYl = −ρYl∇·u+∇· ρD∇Yl +m ̇

000

l

(94)

Conservation of Momentum

∂u

∂t

+u×ω+∇H =

1

ρ

((ρ−ρ∞)g+f+∇· τ) (95)

Divergence Constraint

∇·u =

1

ρcp T

 

∇· k∇T +∇·∑

l

R

cp,ldT ρDl∇Yl −∇·qr +q ̇

000!

+



1

ρcp T

−

1

p0



d p0

dt

(96)

Equation of State

p0(t) = ρTR ∑

l

Yl/Ml (97)

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.

7.2 Temporal Discretization

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

n

,

Y

n

i

, u

n

, H n, and p

n

0

are known. All other quantities can be derived from them. Note that the superscript

(n+1)e refers to an estimate of the value of the quantities at the (n+1)st time step.

1. The thermodynamic quantities ρ, Yi, and p0 are estimated at the next time step with an explicit Euler

step. For example, the density is estimated

ρ

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

·∇ρn +ρ

n∇·un

) (98)

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.

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

∇

2Hn = −

(∇·u)

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

δt

−∇·F

n

(99)

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

These will be described in detail below. Then the velocity is estimated at the next time step

u

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

) (100)

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

δx

u

,

δy

v

,

δz

w



(101)

If the time step is too large, it is reduced so that it satisfies the CFL condition and the procedure starts

from the beginning of the time step. If the time step satisfies the stability condition, the procedure

continues.

3. The thermodynamic quantities ρ, Yi, and p0 are corrected at the next time step. For example, the

density is corrected

ρ

n+1 =

1

2



ρ

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

·∇ρ(n+1)e +ρ

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

)



(102)

The divergence (∇·u)

(n+1)

is derived from the corrected thermodynamic quantities.

4. The pressure is recomputed using estimated quantities

∇

2H(n+1)e = −

2(∇·u)

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

δt

−∇·F

(n+1)e

(103)

The velocity is then corrected

u

n+1 =

1

2

h

u

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



F

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

i (104)

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

divergence.

7.3 Spatial Discretization

Spatial derivatives in the governing equations are written as second order accurate finite differences on a

rectilinear grid. The overall domain is a rectangular box that is divided into rectangular grid cells. Each cell

is assigned indices i, j and k representing the position of the cell in the x, y and z directions, respectively.

Scalar quantities are assigned in the center of each grid cell, thus ρ

n

i jk is the density at the nth time step

in the center of the cell whose indices are i, j and k. Vector quantities like velocity are assigned at cell

faces, thus the x component of velocity u is defined at the faces whose normals are parallel to the x-axis,

the y component v is defined at the faces whose normals are parallel to the y-axis, and the z component w is

defined at the faces whose normals are parallel to the z-axis. The quantity u

n

i jk is the x component of velocity

at the forward pointing face of the i jkth cell; u

n

i−1, jk is at the backward pointing face of the i jkth cell.

26
7.4 Large Eddy vs. Direct Numerical Simulation

The major difference between an LES and a DNS calculation is the form of the viscosity, and the thermal

and material diffusivities. For a Large Eddy Simulation, the dynamic viscosity is defined at cell centers

μi jk = ρi jk (Cs ∆)

2

|S| (105)

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

(106)

The quantity |S| consists of second order spatial differences averaged at cell centers. The thermal conductivity and material diffusivity of the fluid are related to the viscosity by

ki jk =

cp,0 μi jk

Pr

; (ρD)i jk =

μi jk

Sc

(107)

where Pr is the Prandtl number and Sc is the Schmidt number, both assumed constant. Note that the specific

heat cp,0 is that of the dominant species of the mixture. Based on simulations of smoke plumes, Csis 0.14,

Pr and Sc are 0.2. There is no rigorous justification for these choices.

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

ki jk = ∑

l

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

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

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

27
7.5 The Mass Transport Equations

Due to the low Mach number approximation being used in the model, the mass and energy equations are

combined by way of the divergence. The divergence of the flow field contains much of the fire-specific

source terms described above.

7.5.1 Convective and Diffusive Transport

The density at the center of the i jkth cell is updated in time with the following predictor-corrector scheme.

In the predictor step, the density at the (n+1)st time level is estimated based on information at the nth level

ρ

(n+1)e

i jk −ρ

n

i jk

δt

+ (u·∇ρ)

n

i jk = −ρ

n

i jk(∇·u)

n

i jk (111)

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

The species conservation equations are differenced the same way

(ρYl)

(n+1)e

i jk −(ρYl)

n

i jk

δt

+ (u·∇ρYl)

n

i jk = −(ρYl)

n

i jk(∇·u)

n

i jk + (∇· ρD∇Yl)

n

i jk +m ̇

000

i jk (113)

at the predictor step, and

(ρYl)

(n+1)

i jk −

1

2



(ρYl)

n

i jk + (ρYl)

(n+1)e

i jk 

1

2

δt

+(u·∇ρYl)

(n+1)e

i jk = −(ρYl)

(n+1)e

i jk (∇·u)

(n+1)e

i jk +(∇·ρD∇Yl)

(n+1)e

i jk +m ̇

000

i jk

(114)

at the corrector step.

The convective terms are written as upwind-biased differences in the predictor step and downwindbiased differences in the corrector step. In the expressions to follow, the symbol ± means + in the predictor

step and − in the corrector step. The opposite is true for ∓.

(u·∇ρ)i jk =

1∓εu

2

ui jk

ρi+1, jk −ρi jk

δx

+

1±εu

2

ui−1, jk

ρi jk −ρi−1, jk

δx

+

1∓εv

2

vi jk

ρi, j+1,k −ρi jk

δy

+

1±εv

2

vi, j−1,k

ρi jk −ρi, j−1,k

δy

+

1∓εw

2

wi jk

ρi j,k+1 −ρi jk

δz

+

1±εw

2

wi j,k−1

ρi jk −ρi j,k−1

δz

(115)

(u·∇ρYl)i jk =

1∓εu

2

ui jk

(ρYl)i+1, jk −(ρYl)i jk

δx

+

1±εu

2

ui−1, jk

(ρYl)i jk −(ρYl)i−1, jk

δx

+

1∓εv

2

vi jk

(ρYl)i, j+1,k −(ρYl)i jk

δy

+

1±εv

2

vi, j−1,k

(ρYl)i jk −(ρYl)i, j−1,k

δy

+

1∓εw

2

wi jk

(ρYl)i j,k+1 −(ρYl)i jk

δz

+

1±εw

2

wi j,k−1

(ρYl)i jk −(ρYl)i j,k−1

δz

(116)

Note that without the inclusion of the ε’s, these are simple central difference approximations. The ε’s are

local CFL numbers, εu = uδt/δx, εv = vδt/δy, and εw = wδt/δz, where the velocity components are those

that immediately follow. Their role is to bias the differencing upwind. Where the local CFL number is near

unity, the difference becomes nearly fully upwinded. Where the local CFL number is much less than unity,

the differencing is more centralized [36].

28
7.6 Discretizing the Divergence

The divergence in both the predictor and corrector step is discretized

(∇·u)i jk =

1

ρT cp

 

q ̇

000

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

l

(∇· hlρD∇Yl)i jk!+



1

ρT cp

−

1

p0



d p0

dt

(117)

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



(118)

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

Ti jk =

p0

ρi jkR ∑

N

l=0

(Yl,i jk/Ml)

(119)

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

N

∑

l=0

Yl,i jk

Ml

=

1

M0

+

N

∑

l=1



1

Ml

−

1

M0



Yl (120)

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

ρi jk =

p0

T∞R ∑

N

l=0

Yl,i jk/Ml

(121)

Again, because only species 1 through N are explicitly computed, this expression can be written

ρi jk =

M0 p0

T∞R

+

N

∑

l=1



1−

M0

Ml



(ρYl)i jk (122)

7.6.1 Heat Release Rate (Mixture Fraction)

Energy from combustion is released into those grid cells through which the flame sheet (Z = Zf) passes.

The analytical expression for the heat release rate per unit area of flame sheet is

q ̇

00 = ∆HO

dYO

dZ









Z<Zf

(ρD)∇Z ·n (123)

where n is the outward facing unit normal. Note that both dYO/dZ and ∇Z · n are negative. To convert the

analytical expression for the HRR per unit area into a discretized expression for the HRR per unit volume,

all cells through which the flame sheet passes must be identified. Then, the normal derivative of Z must be

computed component by component. For example, suppose the flame sheet passes between cell i jk and cell

29
i+1, jk, and also suppose that Zi jk > Zf > Zi+1, jk. The contribution to the heat release rate per unit volume

in the i jkth and i+1, jkth cells from the small flame sheet is

q ̇

000 = −∆HO

dYO

dZ









Z<Zf

h

(ρD)

i+

1

2

, jk (Zi jk −Zi+1, jk)/δx+ui jk ρi+

1

2

, jk (Z ̄ −Zf)

i

δy δz

δx δy δz

(124)

The distribution of energy is based on a linear interpolation of the value of the mixture fraction. A similar

expression can be derived for all other possible cuts through the cell by the flame sheet. Note that the second

term in the expression above corrects for error associated with locating the flame sheet. The terms above are

a consequence of writing the convective and diffusive flux of fuel (mixture fraction) across the cell interface

in conservative form:

ui jk

(ρZ)i+1, jk + (ρZ)i jk

2

+ (ρD)

i+

1

2

, jk

Zi jk −Zi+1, jk

δx

=

ui jk Zf ρi+

1

2

, jk +ui jk (Z ̄ −Zf) ρi+

1

2

, jk + (ρD)i+

1

2

, jk

Zi jk −Zi+1, jk

δx

(125)

where

Z ̄ =

(ρZ)i+1, jk + (ρZ)i jk

ρi+1, jk +ρi jk

; ρi+

1

2

, jk =

ρi+1, jk +ρi jk

2

(126)

The second two terms on the right hand side constitute the “diffusive” flux, from which the heat release rate

is derived.

7.6.2 Heat Release Rate (Finite-Rate Reaction)

In a DNS calculation (usually), a one-step, finite-rate reaction of a hydrocarbon fuel is assumed

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

For each grid cell, at the start of a time step where t = t

n

and Y

n

CxHy,i jk ≡ YF(t

n

) and Y

n

O2,i jk ≡ YO(t

n

), the

following ODE is solved numerically with a 2nd order Runge-Kutta scheme

dYF

dt

= −

Bρ

a+b−1

i jk

Mb

O Ma−1

F

YF(t)

aYO(t)b

e

−E/RTi jk (128)

dYO

dt

= −

νO MO

νF MF

dYF

dt

(129)

The temperature Ti jk and density ρi jk are fixed at their values at time t

n

and the ODE is iterated from t

n

to

t

n+1

in about 100 time steps. The pre-exponential factor B, the activation energy E, and the exponents a and

b are input parameters. The average heat release rate over the entire time step is given by

q ̇

000n

i jk = ∆H ρ

n

i jk

YF(t

n

)−YF(t

n+1

)

δt

(130)

where δt = t

n+1 −tn

. The species mass fractions are adjusted at this point in the calculation (before the

convection and diffusion update)

Y

n

l,i jk = Yl(t

n

)−

νl Ml

νF MF

(YF(t

n

)−YF(t

n+1

)) (131)

30
7.7 Thermal Radiation

The discretized RTE is derived by integrating equation (48) over the grid cell i jk and the control angle δΩl

Z

Ωl

Z

Vi jk

s·∇I(x,s)dV dΩ =

Z

Ωl

Z

Vi jk

κ(x) [Ib(x)−I(x,s)]dV dΩ (132)

The volume integral on the left hand side is replaced by a surface integral over the cell faces using the

divergence theorem. Assuming that the radiation intensity I(x,s) is constant on each of the cell faces, the

surface integral can be approximated by a sum over the cell faces. Assuming further that I(x,s) is constant

within the volume Vi jk and over the angle δΩl we obtain

6

∑

m=1

Am I

l

m

Z

Ωl

(s·nm)dΩ = κi jk hIb,i jk −I

l

i jki

Vi jk δΩl(133)

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

It must be noticed, that while the intensity is assumed constant within the angle δΩl, its direction covers the

angle δΩlexactly.

In Cartesian coordinates5, the normal vectors nm are the base vectors of the coordinate system and the

integrals over the solid angle do not depend on the physical coordinate, but the direction only. The intensities

on the cell boundaries, I

l

m, are calculated using a first order upwind scheme. If the physical space is swept in

the direction s

l

, the intensity I

l

i jk can be directly solved from an algebraic equation. This makes the numerical

solution of the FVM very fast. Iterations are needed only to account for the reflective boundaries. However,

this is seldom necessary in practice, because of the small time step set by the flow solver.

The spatial discretization for the RTE solver is the same as for the fluid solver. The coordinate system

used to discretize the solid angle is shown in Figure 5. 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φ(θ) 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) (134)

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

Ω

(135)

rounded to the nearest even integer. The number of φ-angles on each band is

Nφ(θ) = max4,0.5NΩ



cos(θ

−)−cos(θ+)

	 (136)

5

In the axisymmetric case equation (133) 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.

31
φ

s

x

z

y

θ

FIGURE 5: Coordinate system of the angular discretization.

rounded to the nearest integer that is divisible by 4. θ

− and θ+ are the lower and upper bounds of 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 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 [37], 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 (133) 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 (133) 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 (137)

where

a

l

i jk = Ax|D

l

x

|+Ay|D

l

y

|+Az|D

l

z

|+κi jk Vi jkδΩl(138)

a

l

x = Ax|D

l

x

| (139)

a

l

y = Ay|D

l

y

| (140)

a

l

z = Az

|D

l

z

| (141)

b

l

i jk = κi jk Ib,i jk Vi jk δΩl

(142)

δΩl =

Z

Ωl

dΩ =

Z

δφ

Z

δθ

sinθ dθ dφ (143)

32
D

l

x =

Z

Ωl

(s

l

·i)dΩ (144)

=

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Ω (145)

=

Z

δφ

Z

δθ

sinφsinθsinθ dθdφ

=

1

2



cosφ

− −cosφ+

 ∆θ−

cosθ

+ sinθ+ −cosθ− sinθ−



D

l

z =

Z

Ωl

(s

l

·k)dΩ (146)

=

Z

δφ

Z

δθ

cosθsinθ dθdφ

=

1

2

∆φhsinθ

+

2 −

sinθ

−

2

i

Here i, j and k are the base vectors of the Cartesian coordinate system. θ

+, θ−, φ+ and φ+ are the upper and

lower boundaries of the control angle in the polar and azimuthal directions, respectively, and ∆θ = θ

+ −θ−

and ∆φ = φ

+ − φ−. The solution method of (137) is based on an explicit marching sequence [38]. The

marching direction depends on the the propagation direction of the radiation intensity. As the marching is

done in the “downwind” direction, the “upwind” intensities in all three spatial directions are known, and the

intensity I

l

i jk can be solved directly. Iterations may be needed only with the reflective walls and optically

thick situations. Currently, no iterations are made.

The boundary condition on a solid wall is given as

I

l

w = ε

σT

4

w

π

+

1−ε

π ∑

Dl

0

w<0

I

l

0

w

|D

l

0

w

| (147)

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 radiative heat flux on the wall is

qw =

NΩ

∑

l=1

I

l

w

Z

δΩl

(s·nw)dΩ =

NΩ

∑

l=1

I

l

wD

l

n

(148)

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

beforehand.

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

(149)

Computationally intensive integration over all the incoming directions is avoided by keeping the solid angle

discretization symmetric x, y and z planes. The connection matrix associates one incoming direction l

0

to

each mirrored direction on each wall cell.

33
The local incident radiation intensity is

Ui jk =

NΩ

∑

l=1

I

l

i jkδΩl

(150)

7.8 Interaction of Droplets and Radiation

The computation of χf for a similar but simpler situation has been derived in Ref. [39]. 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(151)

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

(152)

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Ω

(153)

The limits of the innermost integral are

μd,0 = μμ0 +

p

1−μ

2

p

1−μ

02

; μd,π = μμ0 −

p

1−μ

2

p

1−μ

02

(154)

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

7.9 Thermal and Material Boundary Conditions

Four types of thermal boundary conditions are applied at solid surfaces. The first, and simplest, is an

adiabatic boundary condition that states that there is no temperature gradient normal to the surface. It is

implemented by assigning to the grid cell that is embedded in the solid (the ghost cell) the same temperature

as the first cell in the gas (the gas cell).

The second type of boundary condition is where the solid surface has a prescribed temperature (usually

this prescribed temperature is a function of time).

The third type of boundary condition assumes the solid to be thermally-thin. The surface temperature is

updated in time according to

T

n+1

w = T

n

w +δts

q ̇

00

c +q ̇

00

r

ρscsδ

(155)

where Tw is the wall temperature, δtsis the time step used when updating the thermal boundary conditions

(usually greater than the hydrodynamic time step δt), and ρs, cs, δ are the input density, specific heat and

thickness of the wall. In a DNS calculation where the boundary layer is resolved, the convective flux to the

wall is given by

q ̇

00

c = −k

Tgas −Tw

δn/2

(156)

34
where δn is the size of a grid cell in the normal direction to the wall. In an LES calculation where the

boundary layer is not resolved,

q ̇

00

c = C|Tgas −Tw|

1

3 (Tgas −Tw) W/m2

(157)

whereC is an empirical coefficient (0.95 for vertical surface; 1.43 for horizontal), and Tgas is the temperature

of the gas in the cell bordering the wall.

The fourth type of thermal boundary condition is for a thermally-thick solid. In this case, a one dimensional heat transfer calculation is performed at each boundary cell designated as thermally-thick. The width

of the solid δ is partitioned into N cells, clustered near the front face. The cell boundaries are located at

points xi

xi = f(ξi) = δ

e

sξi/δ −1

e

s −1

(158)

where 0 ≤ i ≤ N, ξi = iδξ, δξ = δ/N, and 0 < s ≤ 1 is a measure of the degree of clustering of the cells

at the front face. The width of each cell is δxi = f

0

(ξi−

1

2

)δξ, 1 ≤ i ≤ N where ξi−

1

2

= (i −

1

2

)δξ. The

temperature at the center of the ith cell is denoted Ts,i. These temperatures are updated in time using an

implicit Crank-Nicholson scheme

T

n+1

s,i −T

n

s,i

δt

=

α

2δxi

 

T

n

s,i+1 −T

n

s,i

δxs,i+

1

2

−

T

n

s,i −T

n

s,i−1

δxs,i−

1

2

+

T

n+1

s,i+1 −T

n+1

s,i

δxi+

1

2

−

T

n+1

s,i −T

n+1

s,i−1

δxi−

1

2

+

!

(159)

for 1 ≤ i ≤ N. The boundary condition is discretized

−ks

T

n+1

s,1 −T

n+1

s,0

δx 1

2

= q ̇

00

c +q ̇

00

r −4 ε σ T

n

3

s,

1

2



T

n+1

s,

1

2

−T

n

s,

1

2



(160)

where Ts,

1

2

= (Ts,1 +Ts,0)/2 is the temperature at the front face. Notice that the radiative emission term has

been linearized

T

(n+1)

4

s,

1

2

−T

n

4

s,

1

2

≈ 4 T

n

3

s,

1

2



T

n+1

s,

1

2

−T

n

s,

1

2



(161)

The wall temperature is defined Tw ≡ Ts,

1

2

= (Ts,0 +Ts,1)/2.

Regardless of how the wall temperature is determined, there are two ways of coupling the wall temperature with the fluid calculation. Gas phase temperatures are defined at cell centers; the wall is defined at

the boundary of the bordering gas phase cell and a “ghost” cell inside the wall. As far as the gas phase calculation is concerned, the normal temperature gradient at the wall is expressed in terms of the temperature

difference between the “gas” cell and the “ghost” cell. The wall temperature affects the gas phase calculation

through the prescription of the ghost cell temperature. This ghost cell temperature has no physical meaning

on its own. Only the difference between ghost and gas cell temperatures matters, for this defines the heat

transfer to the wall. In a DNS calculation, the wall temperature is assumed to be an average of the ghost cell

temperature and the temperature of the first cell in the gas, thus the ghost cell temperature is defined

Tghost = 2Tw −Tgas (162)

For an LES calculation, the heat lost to the boundary is equated with an empirical expression

k

Tgas −Tghost

δn

= C|Tgas −Tw|

1

3 (Tgas −Tw) (163)

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

35
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 (164)

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

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

(166)

where  ̇m

00

l

is the mass flux of species l per unit area, un is the normal component of velocity at the wall

pointing into the flow domain, and δn is the distance between the center of the ghost cell and the center of

the gas cell. Notice that the last term on the right hand side is subtracted at the predictor step and added at

the corrector step, consistent with the biased 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/Ml)

(167)

36
7.10 The Momentum Equation

The three components of the momentum equation are

∂u

∂t

+Fx +

∂H

∂x

= 0 ; Fx = wωy −vωz −

1

ρ



fx +

∂τxx

∂x

+

∂τxy

∂y

+

∂τxz

∂z



(168)

∂v

∂t

+Fy +

∂H

∂y

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

1

ρ



fy +

∂τyx

∂x

+

∂τyy

∂y

+

∂τyz

∂z



(169)

∂w

∂t

+Fz +

∂H

∂z

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

1

ρ



fz +

∂τzx

∂x

+

∂τzy

∂y

+

∂τzz

∂z



(170)

The spatial discretization of the momentum equations take the form

∂u

∂t

+Fx,i jk +

Hi+1, jk −Hi jk

δx

= 0 (171)

∂v

∂t

+Fy,i jk +

Hi, j+1,k −Hi jk

δy

= 0 (172)

∂w

∂t

+Fz,i jk +

Hi j,k+1 −Hi jk

δz

= 0 (173)

where Hi jk is taken at center of cell i jk, ui jk and Fx,i jk are taken at the side of the cell facing in the forward x

direction, vi jk and Fy,i jk at the side facing in the forward y direction, and wi jk and Fz,i jk at the side facing in

the forward z (vertical) direction. In the definitions to follow, the components of the vorticity (ωx,ωy,ωz) are

located at cell edges pointing in the x, y and z directions, respectively. The same is true for the off-diagonal

terms of the viscous stress tensor: τzy = τyz, τxz = τzx, and τxy = τyx. The diagonal components of the stress

tensor τxx, τxx, and τxx; the external force components (fx, fy, fz); and the upwinding bias terms εu, εv, and

εw are located at the respective cell faces.

Fx,i jk =



1∓εw

2

wi+

1

2

, jk ωy,i jk +

1±εw

2

wi+

1

2

, j,k−1 ωy,i j,k−1



−



1∓εv

2

vi+

1

2

, jk ωz,i jk +

1±εv

2

vi+

1

2

, j−1,k ωz,i, j−1,k



−

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



(174)

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



(175)

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



(176)

ωx,i jk =

wi, j+1,k −wi jk

δy

−

vi j,k+1 −vi jk

δz

(177)

37
ωy,i jk =

ui j,k+1 −ui jk

δz

−

wi+1, jk −wi jk

δx

(178)

ωz,i jk =

vi+1, jk −vi jk

δx

−

ui, j+1,k −ui jk

δy

(179)

τxx,i jk = μi jk 2

ui jk −ui−1, jk

δx

−

2

3

(∇·u)i jk≡ μi jk 

4

3

(∇·u)i jk −2

vi jk −vi, j−1,k

δy

−2

wi jk −wi j,k−1

δz



(180)

τyy,i jk = μi jk 2

vi jk −vi, j−1,k

δy

−

2

3

(∇·u)i jk≡ μi jk 

4

3

(∇·u)i jk −2

ui jk −ui−1, jk

δx

−2

wi jk −wi j,k−1

δz



(181)

τzz,i jk = μi jk 2

wi jk −wi j,k−1

δz

−

2

3

(∇·u)i jk≡ μi jk 

4

3

(∇·u)i jk −2

ui jk −ui−1, jk

δx

−2

vi jk −vi j−1,k

δy



(182)

τ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



(183)

τ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



(184)

τ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



(185)

εu =

uδt

δx

(186)

εv =

vδt

δy

(187)

εw =

wδt

δz

(188)

The variables εu, εv and εw are local CFL numbers evaluated at the same locations as the velocity component immediately following them, and serve to bias the differencing of the convective terms in the upwind

direction. The subscript i+

1

2

indicates that a variable is an average of its values at the ith and the (i+1)th

cell. The divergence defined in Eq. (117) 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

(189)

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.

7.10.1 Force Terms

The external force term components, in addition to including the effects of buoyancy, may also include the

drag force from sprinkler droplets.

fx,i jk =

1

2

∑ρCDπr

2

d

(ud −ui jk)|ud −u|

δxδyδz

−(ρi+

1

2

, jk −ρ∞)gx (190)

fy,i jk =

1

2

∑ρCDπr

2

d

(vd −vi jk)|ud −u|

δxδyδz

−(ρi, j+

1

2

,k −ρ∞)gy (191)

fz,i jk =

1

2

∑ρCDπr

2

d

(wd −wi jk)|ud −u|

δxδyδz

−(ρi j,k+

1

2

−ρ∞)gz (192)

where g = (gx,gy,gz) is the gravity vector, rd is the radius of a droplet, u = (ud,vd,wd) the velocity of a

droplet, CD the drag coefficient, and δxδyδz the volume of the i jkth cell. The summations represent all

droplets within a grid cell centered about the x, y and z faces of a grid cell respectively.

38
7.10.2 Time Step

The time step is determined by the CFL condition, and in cases of high viscosity, a parabolic stability

criterion typical of explicit second order accurate schemes

δt < min

δx

ui jk

,

δy

vi jk

,

δz

wi jk

,

ρi jk δx

2

8μi jk

,

ρi jk δy

2

8μi jk

,

ρi jk δz

2

8μi jk 

(193)

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 above

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 parabolic stability criterion is only invoked

for a DNS calculation.

7.11 The Pressure Equation

The divergence of the momentum equation yields a Poisson equation for the pressure

Hi+1, jk −2Hi jk +Hi−1, jk

δx

2

+

Hi, j+1,k −2Hi jk +Hi, j−1,k

δy

2

+

Hi j,k+1 −2Hi jk +Hi j,k−1

δz

2

= −

Fx,i jk −Fx,i−1, jk

δx

−

Fy,i jk −Fy,i, j−1,k

δy

−

Fz,i jk −Fz,i j,k−1

δz

−

∂

∂t

(∇·u)i jk (194)

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 that is part of a library

of routines for solving elliptic PDEs called CRAYFISHPAK [40]. To ensure that the divergence of the fluid

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

∂

∂t

(∇·u)i jk =

(∇·u)

(n+1)e

i jk −(∇·u)

n

i jk

δt

(195)

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

(196)

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

Direct Poisson solvers are most efficient if the domain is a rectangular region, although other geometries

such as cylinders and spheres can be handled almost as easily. For these solvers, the no-flux condition (197)

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

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 (197) cannot be directly

prescribed at the boundaries of these blocked cells. Fortunately, it is possible to exploit the relatively small

changes in the pressure from one time step to the next to enforce the no-flux condition. At the start of a time

step, the components of the convection/diffusion term F are computed at all cell faces that do not correspond

to walls. At those cell faces that do correspond to solid walls, prescribe

Fn = −

∂H

∂n

∗

+βun (198)

39
where Fn is the normal component of F at the wall, and β is a relaxation factor empirically determined to be

about 0.8 divided by the time step δt. The asterisk indicates the most recent value of the pressure. Obviously,

the pressure at this particular time step is not known until the Poisson equation is solved. Equation (198)

asserts that following the solution of the Poisson equation for the pressure, the normal component of velocity

un will be driven closer to zero according to

∂un

∂t

≈ −βun (199)

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

(200)

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

7.12 Particle Tracking

Thermal elements are introduced into the flow field as a means of introducing heat and as a way to visualize

the flow. The position xp of each thermal element is governed by the equations

dxp

dt

= u (202)

The thermal element positions are updated according to the same predictor-corrector scheme that is applied

to the other flow quantities. Briefly, the position xp of a given thermal element is updated according to the

two step scheme

x

(n+1)e

p = x

n

p +δt u

n

(203)

x

n+1

p =

1

2



x

n

p +x

(n+1)e

p +δt u

(n+1)e



(204)

where the bar over the velocity vector indicates that the velocity of the fluid is interpolated at the element’s

position.

40
8 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 analyze industrial-scale fires, it can be used reliably when the

fire size is specified and the building is relatively large in relation to the fire. In these cases, the model predicts flow velocities and temperatures to an accuracy of 10 to 20% compared to experimental measurements.

Currently, research is focussed on improving both the gas phase and solid phase descriptions of combustion

in the model so that simulations involving fire growth and suppression, especially in residential sized rooms,

can be improved.

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.

41
9 Nomenclature

As water droplet surface area

B pre-exponential factor for Arrhenius reaction

C Sprinkler C-Factor

CD drag coefficient

Cs Smagorinsky constant (LES)

cp constant pressure specific heat

D diffusion coefficient

D

∗

characteristic fire diameter

dm median volumetric droplet diameter

E activation energy

f external force vector (excluding gravity)

g acceleration of gravity

H total pressure divided by the density

h enthalpy; heat transfer coefficient

hi enthalpy of ith species

h

0

i

heat of formation of ith species

I radiation intensity

Ib radiation blackbody intensity

k thermal conductivity; suppression decay factor

M molecular weight of the gas mixture

Mi molecular weight of ith gas species

m ̇

00

f

fuel mass flux

m ̇

00

w water mass flux

m

00

w water mass per unit area

m ̇

00

O

oxygen consumption rate per unit area

Nu Nusselt number

Pr Prandtl number

p pressure

p0 background pressure

p ̃ pressure perturbation

qr radiative heat flux vector

q ̇

000 heat release rate per unit volume

q ̇

00

r

radiative flux to a solid surface

q ̇

00

c

convective flux to a solid surface

Q ̇ total heat release rate

Q

∗

characteristic fire size

R universal gas constant

Re Reynolds number

rd water droplet radius

RTI Response Time Index of sprinkler

s unit vector in direction of radiation intensity

Sc Schmidt number

Sh Sherwood number

T temperature

t time

tb thermal element burn-out time (LES)

42
U integrated radiant intensity

u = (u,v,w) velocity vector

m ̇

000

i

production rate of ith species per unit volume

We Weber number

x = (x,y,z) position vector

Xi volume fraction of ith species

Yi mass fraction of ith species

Y

∞

O mass fraction of oxygen in the ambient

Y

I

F mass fraction of fuel in the fuel stream

ys soot yield

Z mixture fraction

γ ratio of specific heats; Rosin-Rammler exponent

∆H heat of combustion

∆HO energy released per unit mass oxygen consumed

δ wall thickness

κ absorption coefficient

μ dynamic viscosity

νi stoichiometric coefficient, species i

Φ dissipation function

ρ density

τ viscous stress tensor

χr radiative loss fraction

σ Stefan-Boltzmann constant

ω = (ωx,ωy,ωz) vorticity vector

43
References

[1] H.C. Hottel. Stimulation of Fire Research in the United States After 1940. Combustion Science and

Technology, 39:1–10, 1984.

[2] J. Quintiere. A Perspective on Compartment Fire Growth. Combustion Science and Technology,

39:11–54, 1984.

[3] G.P. Forney and W.F. Moss. Analyzing and Exploiting Numerical Characteristics of Zone Fire Models.

Fire Science and Technology, 14:49–60, 1994.

[4] S.V. Patankar. Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing, New York, 1980.

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

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

1999.

[7] J. Smagorinsky. General Circulation Experiments with the Primitive Equations. I. The Basic Experiment. Monthly Weather Review, 91:99–164, 1963.

[8] J.W. Deardorff. Numerical Investigation of Neutral and Unstable Planetary Boundary Layers. Journal

of Atmospheric Sciences, 29:91–115, 1972.

[9] M. Germano, U. Piomelli, P. Moin, and W.H. Cabot. A Dynamic Subgrid-Scale Eddy Viscosity Model.

Physics of Fluids A, 3:1760–1765, 1991.

[10] D.K. Lilly. A Proposed Modification of the Germano Subgrid-Scale Closure Method. Physics of Fluids

A, 4:633–635, 1992.

[11] 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:45–52, 1997.

[12] R.C. Reid, J.M. Prausnitz, and B.E. Poling. Properties of Gases and Liquids. McGraw-Hill, Inc., New

York, 4th edition, 1987.

[13] C. Huggett. Estimation of the Rate of Heat Release by Means of Oxygen Consumption Measurements.

Fire and Materials, 4:61–65, 1980.

[14] T. Ma. Numerical Simulation of an Axi-symmetric Fire Plume: Accuracy and Limitations. Master’s

thesis, University of Maryland, 2001.

[15] G. Heskestad. SFPE Handbook, chapter Fire Plumes. National Fire Protection Association, Quincy,

Massachusetts, 2nd edition, 1995.

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

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

44
[18] R. Siegel and J. R. Howell. Thermal Radiation Heat Transfer. Hemisphere Publishing Corp., Philadelphia, 3 edition, 1992.

[19] 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 20899, 1993.

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

[21] J.P. Holman. Heat Transfer. McGraw-Hill, New York, 5th edition, 1989.

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

[23] G. Heskestad and R.G. Bill. Quantification of Thermal Responsiveness of Automatic Sprinklers Including Conduction Effects. Fire Safety Journal, 14:113–125, 1988.

[24] M. DiMarzo. The Effect of Minute Water Droplets on a Simulated Sprinkler Link Thermal Response.

Technical report, University of Maryland, 1999.

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

[26] 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, International

Association For Fire Safety Science, pages 1165–1176, 1986.

[27] A. Hamins and K.B. McGrattan. Reduced-Scale Experiments to Characterize the Suppression of Rack

Storage Commodity Fires. Technical Report NIST Internal Report (NISTIR 6439), National Institute

of Standards and Technology, Gaithersburg, Maryland 20899, 1999.

[28] N. Cheremisinoff. Encyclopedia of Fluid Mechanics, Volume 3: Gas-Liquid Flows. Gulf Publishing

Co., Houston, Texas, 1986.

[29] T. Ravigururajan and M. Beltran. A Model for Attenuation of Fire Radiation Through Water Droplets.

Fire Safety Journal, 15:171–181, 1989.

[30] A. Tuntomo, C. Tien, and S. Park. Optical Constants of Liquid Hydrocarbon Fuels. Combustion

Science and Technology, 84:133–140, 1992.

[31] W. Wiscombe. Mie Scattering Calculations: Advances in Technique and Fast, Vector-Speed Computer

Codes. Technical report, National Center for Atmospheric Research (NCAR), 1979.

[32] G.M. Hale and M.R. Quarry. Optical constants of water in the 200 nm to 200 mm wavelength region.

Applied Optics, 12:555–563, 1973.

[33] 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, International Association For Fire Safety

Science, pages 901–912, 1994.

[34] K.B. McGrattan, R.G. Rehm, and H.R. Baum. Fire-Driven Flows in Enclosures. Journal of Computational Physics, 110(2):285–292, 1994.

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

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

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

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

[39] W. Yang, T. Parker, H.D. Ladouceur, and R.J. Kee. The Interaction of Thermal Radiation and Water Mists in Fire Suppression. In Proceedings of the Second NRIFD Symposium. National Research

Institute of Fire and Disaster (Japan), 2002.

[40] Green Mountain Software, Boulder, Colorado. CRAYFISHPAK User’s Guide, Cray Version 1.1, 1990.

46