Fire Dynamics Simulator –

Technical Reference Guide

Kevin B. McGrattan

Howard R. Baum

Ronald G. Rehm

Anthony Hamins

Glenn P. Forney

NISTIR 6467
NISTIR 6467

Fire Dynamics Simulator -

Technical Reference Guide

Kevin B. McGrattan

Howard R. Baum

Ronald G. Rehm

Anthony Hamins

Glenn P. Forney

Building and Fire Research Laboratory

January 2000

U.S. Department of Commerce

William M. Daley, Secretary

Technology Administration

Dr. Cheryl L. Shavers, Under Secretary of Commerce for Technology

National Institute of Standards and Technology

Raymond G. Kammer, 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.

ii
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

3 Combustion 7

3.1 Thermal Elements . ....................................... 7

3.2 One-step, Finite-rate Combustion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

4 Thermal Boundary Conditions 9

5 Sprinklers 11

5.1 Sprinkler Activation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

5.2 Sprinkler Droplet Size Distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

5.3 Sprinkler Droplet Trajectory in Air . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

5.4 Sprinkler Droplet Transport on a Surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

5.5 Fire Suppression by Water . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

6 Numerical Method 20

6.1 Simplified Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

6.2 Temporal Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

6.3 Spatial Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

6.4 Large Eddy vs. Direct Numerical Simulation . . . . . . . . . . . . . . . . . . . . . . . . . 22

6.5 The Mass Transport Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

6.5.1 Convective and Diffusive Transport . . . . . . . . . . . . . . . . . . . . . . . . . . 23

6.5.2 Heat Release Rate (LES) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

6.5.3 Heat Release Rate (DNS) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

6.5.4 Thermal Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

6.5.5 Species Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

6.5.6 Density Boundary Condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

6.6 The Momentum Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

6.6.1 Force Terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

6.6.2 Time Step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

6.7 The Pressure Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

6.8 Particle Tracking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

7 Conclusion 33

8 Nomenclature 34

References 36

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 = ∇  (ρD)l∇Yl +W ̇ 000

l (2)

Conservation of Momentum

ρ

∂u

∂t + (u  ∇)u+ ∇p = ρg + f + ∇  τ (3)

Conservation of Energy

∂

∂t (ρh) + ∇  ρhu 

Dp

Dt = q ̇000 + ∇  k∇T + ∇ 

∑

l

hl (ρD)l∇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. The energy driving the system is represented by the heat release

rate  ̇q000 in Eq. (4). The term Dp=Dt = ∂p=∂t + u  ∇p is a material derivative. All other symbols are listed

in the Nomenclature (Section 8).

2.2 State, Mass and Energy Equations

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

density, pressure and enthalpy; ρ, p and h. The pressure is decomposed into three components

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

The first term on the right hand side is the “background” pressure, the second is the hydrostatic contribution,

and the third is the flow-induced perturbation pressure. For most applications, p0 is constant. However, if

3
the enclosure is tightly sealed, p0 is allowed to increase (or decrease) with time as the pressure within the

enclosure rises due to thermal expansion or falls due to forced ventilation. Also, if the height of the domain

is on the order of a kilometer, 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.

A further assumption about the thermodynamic variables is that the constant-pressure specific heat of

the ith species cp;i is assumed to be independent of temperature. Under this assumption, the enthalpy can be

written as:

h = ∑

l

hlYl = T ∑

l

cp;lYl (7)

The specific heat for each species can be expressed in terms of the number of internal degrees of freedom νl

active in that molecule.

cp;l =



2 + νl

2 

R

Ml =



γl

γl  1

R

Ml

(8)

If the ratio of specific heats γl for each species is assumed to be that of a diatomic molecule (ν = 5, γ = 7=5),

the equation of state can be rewritten in the form1

p0(t ) = γ  1

γ ρh (9)

Taking the material derivative of Eq. (9) and using the mass and energy conservation equations, the divergence of the velocity field, ∇  u, can be written in terms of the thermodynamic quantities

∇  u =

γ  1

γ p0 

q ̇000 + ∇  k∇T + ∇ 

∑

l

cp;lTρD∇Yl 

1

γ  1

d p0

dt !

(10)

If the enclosure is tightly sealed, 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. (10) over the entire domain Ω

d p0

dt =

γ  1

V 

Z

Ω

q ̇000 dV +

Z

∂Ω

k∇T  dS + ∑

l

Z

∂Ω

cp;lTρD∇Yl  dS!



γp0

V Z

∂Ω

u  dS (11)

where V is the volume of the enclosure.

1The basis of this approximation is that nitrogen will be the dominant species in most fire scenarios.

4
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 obtain2

∂u

∂t

+

1

2

∇juj

2  u ω +

1

ρ

∇p ̃ =

1

ρ ((ρ  ρ∞)g + f + ∇  τ) (12)

To simplify this equation further, a substitution is made

∇H 

1

2

∇juj

2 +

1

ρ

∇p ̃ (13)

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

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

dΓ

dt =

I

1

ρ (∇p ̃ + (ρ  ρ∞)g + ∇  τ)  dx (14)

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.

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 + ∇  τ) (15)

The linear algebraic system arising from the discretization of Eq. (15) 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 (16)

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 = juj

2=2 outgoing

H = 0 incoming (17)

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.

The components of the viscous stress tensor are given by

τi j = μ

∂ui

∂xj

+

∂uj

∂xi  δi j

2

3

∂uk

∂xk 

(18)

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

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

2∇juj

2  u ω.

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

viscosity can be modeled as

μLES = max 

μDNS; ρ (Cs ∆)

2 jSj

(19)

where Cs is an empirical constant, ∆ is a length on the order of the size of a grid cell, and jSj is the magnitude

of the deformation tensor

jSj

2 = 2

∂u

∂x 2

+ 2

∂v

∂y2

+ 2

∂w

∂z 2

+



∂u

∂y

+

∂v

∂x2

+



∂u

∂z

+

∂w

∂x 2

+



∂v

∂z

+

∂w

∂y 2

(20)

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

kLES =

μLES cp

Pr ; (ρD)l;LES =

μLES

Sc (21)

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.

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 107(Ml T )

1

2

σ2

l Ωv

kg

m s

(22)

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

kl =

μl cp;l

Pr

W

m K (23)

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

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

Dlm =

2:66 107 T3=2

M 1

2

lm σ2

lm ΩD

m2

s

(25)

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

where species 0 is nitrogen.

6
3 Combustion

In the model, the combustion processes are modeled in two ways. In a DNS calculation where the diffusion of fuel and oxygen can be modeled directly, a global one-step, finite-rate chemical reaction is usually

used. In cases where the computational grid is not fine enough to resolve the diffusion of fuel and oxygen,

Lagrangian particles or “thermal elements” are used to introduce the thermal energy of the fire into the

calculation, and also to visualize the movement of the smoke and hot gases. The thermal elements carry

the heat released by the fire, providing a self-consistent description of the smoke transport at all resolvable

length and time scales.

The terms LES and DNS refer to solution methodologies of the equations of fluid dynamics. The two

relatively simple combustion models described below are not inherent to LES or DNS. In fact, the terms

“Direct Numerical Simulation” and “Large Eddy Simulation” are more difficult to define when used in the

context of a reacting flow calculation.

3.1 Thermal Elements

In an LES calculation, combustion is occuring at length scales well below the resolution limits of the underlying numerical grid, thus the mixing of fuel gases and air cannot be calculated directly. Instead, the

fire is represented by discrete Lagrangian particles (or thermal elements) that originate at solid surfaces and

release heat at a specified rate. Surfaces heat up due to both convective and radiative heat transfer from some

external source, like an ignitor. When a surface heats up to its prescribed ignition temperature, thermal elements are ejected from the surface, and burned at a prescribed rate. The thermal elements are introduced at

a rate of  ̇n00 particles per unit time per unit area with a small initial velocity. This initial velocity is a function

of the mass loss rate per unit area of the fuel bed, which can be obtained from the heat release rate per unit

area  ̇q00

f , the density of the fuel gases ρf , and the heat of combustion ∆H

un =

q ̇00

f

ρf ∆H (27)

The heat release rate of a single thermal element is given by

q ̇p =

q ̇00

f

n ̇00

1

tb

(t t0 < tb) (28)

where tb is the burn-out time of the thermal element, and t

0 is the time the element is ejected from the burning

surface. A fraction of the energy assigned to each thermal element is emitted as radiation and potentially reabsorbed by surrounding smoke and hot gases. The volumetric heat release rate term in the energy equation

is thus given by

q ̇000 =

∑i q ̇p;i(1  χr)

δxδyδz

+ κ(x)

Np

∑

i=1

χr q ̇p;i

4πjxp;i  xj

2 e

R

κ(l) dl (29)

The first summation is over all active thermal elements in the grid cell whose volume is δxδyδz. The

second summation is over all other active thermal elements outside of the grid cell whose radiative energy

is potentially re-absorbed by the soot contained within the grid cell. The fraction of the elements’ energy

converted into thermal radiation is constant and given by χr. The position of the ith element is xp;i, the heat

release rate is  ̇qp;i, and the segment of the line connecting the points xp;i and x is given by dl. Since there are

hundreds of thousands of thermal elements in a typical calculation, the summation is made over a sampling

of the elements that are still burning. The absorption coefficient κ is given by

κ = 1186 ρYs T

ρsoot

m1 (30)

7
where ρYs is the mass of soot per unit volume and ρsoot is the density of soot itself [13]. More simply,

ρYs=ρsoot  fv the soot volume fraction.

This simple ray-tracing technique to compute the transport of radiant energy from the fire to its surroundings works well in cases where the fire itself is the source of the radiation. However, once walls

and the smoke attain temperatures above 400

C, the fire is no longer the sole source of radiation and the

use of ray-tracing become more expensive and cumbersome. A better methodology would be to solve the

governing equation for thermal radiation directly. An effort is underway to implement this solver.

The heat release rate per unit area  ̇q00

f is prescribed by the user and is based on experimental measurements for a configuration resembling as much as possible that being modeled. Typically, a calorimetry

experiment will produce a time history of the total heat release rate of the fuel array. Dividing the heat

release rate by the estimated area of the burning surface yields a time history of the heat release rate per unit

area  ̇q00

f (t ). Often this function can be approximated as a constant, but a time-dependent function can also

be used in the calculation.

The burn-out time tb is obtained from the plume correlation of Baum and McCaffrey [14]. It is assumed

that a thermal element burns out somewhere in the intermittent region of the plume, 1:32D < z < 3:30D

,

where z is the height above the fuel surface, D = (Q ̇ =cpρ∞T∞p

g)

2=5 is the characteristic diameter of the

fire, and Q ̇ is the total heat release rate of the fire. An estimate of the burn-out time can be made

Z

1:32D

0

dz

w(z) < tb <

Z

3:30D

0

dz

w(z)

(31)

where

w(z) = 

2:18pgz (z < 1:32D ) flame region

2:45p

gD (1:32D < z < 3:30D ) intermittent region (32)

The burn-out time falls somewhere between 1:05pD=g < tb < 1:86pD=g and is usually a few tenths of

a second.

The burn-out time of any thermal element will vary based on the concentration of oxygen in the gas

surrounding it3. Oxygen is consumed in any given grid cell based on the amount of heat generated in that

grid cell. The source term in the oxygen transport equation becomes

W ̇ 000

O2 = 

q ̇000

∆HO2

(33)

where ∆HO2 is the amount of heat liberated per unit mass of oxygen consumed (by default 13,100 kJ/kg).

When the oxygen mass fraction YO2 falls to a certain prescribed lower limit, combustion is assumed to stop,

and the unburned fuel associated with the thermal elements remains unburned until more oxygen can be

found. This combustion model is preliminary. In cases where a fire in a room is spreading to the point

of flashover, it can no longer be assumed that there is a constant flux of fuel emanating from the burning

surfaces, nor can it be assumed that the fuel even burns. In simulations approaching flashover, this simple

model breaks down, and an improved combustion model is being developed.

Another useful input parameter associated with a given solid fuel is the available potential energy per

unit volume of the unburned fuel. If prescribed, the computational cell representing a piece of the fuel will

disappear once the potential energy has been liberated. This parameter is denoted by q000 or Heat Released

Per Unit Volume, and is expressed in units of kJ/m3. If prescribed, this parameter permits fire spread through

and consumption of the solid fuel, as in the burning of a boxed commodity or a wood crib.

3By default, oxygen transport is not included in an LES calculation, and the burn-out time of any thermal element is based on

standard flame height correlations.

8
3.2 One-step, Finite-rate Combustion

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

νCxHy CxHy + νO2 O2 ! νCO2 CO2 + νH2O H2O (34)

The reaction rate is given by the expression

d [CxHy]

dt = B [CxHy]

a [O2]b eE=RT (35)

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

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

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

provides a valid description of the combustion process.

4 Thermal Boundary Conditions

There are four types of thermal boundary conditions at solid surfaces: adiabatic, fixed temperature, thermallythin or thermally-thick. The choice of boundary condition depends on the surface material and whether or

not it plays a role in the fire. The simplest boundary condition is adiabatic, that is, no heat transfer at all. In

this case, the solid is assumed to be at the same temperature as the surrounding gas. The next type of boundary condition is where a temperature is prescribed at the solid surface. If the surface material is assumed to

be thermally-thick, the one-dimensional heat conduction equation is applied in the direction n normal to the

solid surface

ρs cs

∂Ts

∂t = ks

∂2Ts

∂n2 ; ks

∂Ts

∂n









n=0 = q ̇00

c + q ̇00r  q ̇00rr (36)

where ρs, cs and ks are the density, specific heat and conductivity of the material; and  ̇q00

c ,  ̇q00r ,  ̇q00rr are the

convective, radiative and re-radiative fluxes at the surface. If the surface material is assumed to be thermallythin, its temperature is governed by its density, specific heat and thickness δ

dTs

dt =

q ̇00

c + q ̇00r  q ̇00rr

ρscsδ (37)

In this case, the individual values of the parameters ρs, cs and δ are not as important as their product.

The heat fluxes to the surface for either a thermally-thick or thermally-thin material consist of gains and

losses from convection and radiation. The convective heat flux to the surface  ̇q00

c is either obtained directly

from the gradient of the gas temperature at the boundary in a DNS calculation

q ̇00

c = k

∂T

∂n

(38)

or it is obtained from a correlation of the form

q ̇00

c = h ∆T W/m2 ; h = C j∆T j

1

3 W/m2/K (39)

in an LES calculation. Here ∆T is the difference between the wall and gas temperature, andC is an empirical

constant with a default value of 1.43 for a horizontal surface and 0.95 for a vertical surface [17].

9
The radiative heat flux to the surface,  ̇q00

r , is calculated based on the assumption that a prescribed fraction

of the heat released from the burning thermal elements is radiated away, and this energy is absorbed by the

surrounding surfaces, as well as by the smoke-laden gas. At a given point xs on a solid surface, the radiative

flux is given as

q ̇00

r =

Np

∑

i=1

cos(φi)

χr q ̇p;i

4πjxp;i  xsj

2 e

R

κ(l) dl (40)

xp;i is the position and  ̇qp;i is the heat release rate of the ith thermal element, dl is an element of the line

segment connecting the points xp;i and xs, and φi is the angle formed by the normal to the surface and the

vector xp;i  xs.

Energy from the heated wall is lost as thermal radiation according to

q ̇00

rr = σε (T4s  T4∞ ) (41)

where σ is the Stefan-Boltzmann constant σ = 5:67 108 W/m2/K4, and ε is the emissivity of the surface.

10
5 Sprinklers

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

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

5.1 Sprinkler Activation

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

forth by Heskestad and Bill [18], 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

juj

RTI (Tg  Tl ) 

C

RTI (Tl  Tm) 

C2

RTIβjuj (42)

Here Tl is the link temperature, Tg is the gas temperature in the neighborhood of the link, Tm is the temperature of the sprinkler mount, 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 [19]

to be 6 106 K/(m/s) 1

2 , and its value is relatively constant for different types of sprinklers.

5.2 Sprinkler Droplet Size Distribution

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

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

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

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

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

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

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

F (d) =

8

<

:

p

1

2π Z

d

0

1

σd 0

e [ln(d 0=dm)]2

2σ2 dd 0 (d  dm)

1  e0:693(

d

dm )

γ (dm < d)

(43)

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

m or

less), and γ and σ are empirical constants equal to about 2.4 and 0.6, respectively. 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 [21]

dm

D ∝ We 1

3 (44)

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 =

ρwU2D

σw

(45)

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

(72:8 103 N/m at 20

C). The discharge velocity can be computed from the mass flow rate, which is a

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

Eq. (44) appears to be independent of flow rate and operating pressure. Three different sprinklers were tested

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

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)

d3 Z

∞

0

F 0(d 0)

d 03 dd 0

(46)

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

U (d) = Z

d

0

f (d 0) dd 0(47)

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

Figure 1: 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:43.

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

Typically, 1,000 particles per sprinkler per second are tracked. The total number of droplets represented

by each computed droplet is given by  ̇mw=(n ̇ md ), where  ̇mw is the mass flow rate of water from a single

sprinkler,  ̇n is the number of droplets tracked per sprinkler per second, and md is the average mass of a

12
droplet. The average mass of a droplet can be expressed in terms of the diameter PDF

md =

Z

∞

0

f (d 0 )

4

3

πρw 

d 0

2 3

dd 0 =

1

6

πρw Z

∞

0

F 0(d 0) dd 0 Z

∞

0

F 0(d 0)

d 03 dd 0

(48)

where ρw is the density of water. The number of droplets tracked per sprinkler per second,  ̇n, is controlled

by the user.

5.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πr2

d (ud  u)jud  uj

δxδyδz

(49)

where Cd is a drag coefficient, 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 πr2

d (ud  u)jud  uj (50)

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

CD = 8

<

:

24=Re Re < 1

24 

1 + 0:15Re0:687 =Re 1 < Re < 1000

0:44 1000 < Re

(51)

Re =

ρ jud  uj 2rd

μ

(52)

where μ is the dynamic viscosity of air. The sprinkler spray droplet temperature Td and mass md are governed

by the following equations

dTd

dt =

Ashd (Tg  Td )

cw md

Td < Tv (53)

dmd

dt = 

Ashd (Tg  Td )

hv

Td = Tv (54)

where cw is the specific heat of water, hv is the energy of vaporization of water, As = 4πr2

d is the surface area

of the droplet, Tg is the gas temperature, Tv is the vaporization temperature of water, hd is the heat transfer

coefficient, given by

hd =

Nuk

2rd

(55)

Nu is the Nusselt number given by

Nu = 2 + 0:6 Re 1

2 Pr 13 (56)

A good approximation for the Prandtl number is about 0.7, and k is the thermal conductivity of air.

13
5.4 Sprinkler Droplet Transport on a Surface

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

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

This “dripping” velocity has been measured to be on the order of 0.5 m/s [22]. Both the temperature of nonburning surfaces and the heat release rate of burning surfaces are affected by the droplets. At non-burning

surfaces, water absorbs heat from the hot surface and radiant energy from the fire. The equations are slightly

different than those governing heat transfer to the airborne drops.

dTd

dt =

0:5As(hs(Ts  Td ) + q ̇00

r )

cw md

Td < Tv (57)

dmd

dt = 

0:5As(hs(Ts  Td ) + q ̇00

r )

hv

Td = Tv (58)

where Ts is the temperature of the surface,  ̇q00

r is the incoming radiant energy flux, and hs is the heat transfer

coefficient between the solid surface and the water droplet, assumed to be constant at 3000 W/m2/K [22].

5.5 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. [23]. The authors consider dozens of rack storage commodity fires of different

geometries and water application rates, and characterize the suppression rates in terms of a few global

parameters. Their analysis yields an expression for the total heat release rate from a rack storage fire after

sprinkler activation

Q ̇ = Q ̇ 0 ek(tt0 )

(59)

where Q ̇ 0 is the total heat release rate at the time of application t

0, and k is a fuel-dependent constant. For

the FMRC Standard Plastic commodity k is given as

k = 0:716  ̇m00

w  0:0131 s1 (60)

where  ̇m00

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  ̇m00

w  0:0040 s1 (61)

Unfortunately, this analysis is based on global water flow and burning rates. Equation (59) 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 heat release rate of the fuel

can be expressed in the form [22]

q ̇00

f (t ) = q ̇00f ;0(t ) 

e

R

k1 dt + k2(t t0)

(62)

14
Here  ̇q00

f ;0(t ) is the heat release rate per unit area of the fuel when no water is applied and k

1 and k2 are

functions of the local water mass per unit area, m00

w, expressed in units of kg/m2.

k1 = a1 m00

w s1 (63)

k2 = a2 m00

w + b2 s1 (64)

The linear term in Eq. (62) is based on the observation that for a boxed commodity, it is possible for the

local heat release rate to increase as the fire burns into the box and is protected from the water droplets

by material overhead, thus often a gradual increase in the heat release rate is observed following the initial

decrease after water is applied.

To develop the suppression model for the FMRC Standard Plastic commodity, 19 experiments were

conducted at UL under a 2 MW calorimeter [22]. These experiments were designed as small-scale RDD

(Required Delivered Density) tests. The fuel/sprinkler arrangement consisted of four boxes of the FMRC

Plastic Commodity. The boxes were stacked two high. The two stacks were positioned 15 cm (6 in) apart,

the same separation that is commonly used in full-scale tests. A water applicator was positioned above the

boxes to deliver a uniform water flux onto the tops of the boxes. The applicator consisted of four nozzles

that were 60 cm (2 ft) apart and 30 cm (1 ft) above the plane of the box tops. Several nozzle sizes were used,

depending on the desired water flow. Table 4 lists the average water application rate per unit area and the

time of water application. The time of water application was varied from 30 s to 200 s. The water flux at the

box top was varied from 0.03 kg/m2/s (0.05 gpm/ft2) to 0.66 kg/m2/s (0.97 gpm/ft2). The ignition source

Table 1: Time and Rate of Water Application

Test Application Total Water Flow Average Water Flux

No. Time (s) (L/s) (gpm) (L/m2/s) (gpm/ft2)

1 380 0.98 15.5 0.66 0.97

2 470 0.57 9.0 0.38 0.56

3 65 0.41 6.5 0.28 0.41

4 106 0.41 6.5 0.28 0.41

5 115 0.11 1.8 0.074 0.11

6 122 0.11 1.8 0.074 0.11

7 150 0.079 1.3 0.053 0.08

8 93 0.11 1.8 0.074 0.11

9 93 0.21 3.3 0.14 0.20

10 110 0.21 3.3 0.14 0.20

11 205 0.21 3.3 0.14 0.20

12 116 0.16 2.5 0.11 0.16

13 63 0.16 2.5 0.11 0.16

14 64 0.28 4.5 0.19 0.28

15 71 0.079 1.3 0.053 0.08

16 62 0.047 0.9 0.032 0.05

17 104 0.047 0.9 0.032 0.05

18 58 0.079 1.3 0.053 0.08

19 30 0.079 1.3 0.053 0.08

was a propane igniter that consisted of two parallel 12.5 mm diameter copper tubes each 30 cm long.

The heat release rate histories for the experiments and the simulations are given in Figs. 2–4. The decay,

and in some cases re-growth, of the fire is captured reasonably well by the simulations. A weakness of the

15
suppression algorithm, however, is its reliance on 5 empirical coefficients that are not easily measured. It

is hoped that further work in this area will provide more insight into fire suppression, and the numerical

algorithm will reflect this improved understanding.

16
Figure 2: Simulated (solid lines) and experimental (dashed lines) heat release rates for Tests 1, 3–7.

17
Figure 3: Simulated (solid lines) and experimental (dashed lines) heat release rates for Tests 8–13.

18
Figure 4: Simulated (solid lines) and experimental (dashed lines) heat release rates for Tests 14–19.

19
6 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 [24, 25]. 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.

6.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 (65)

Conservation of Species

∂ρYl

∂t

+ u  ∇ρYl = ρYl∇  u + ∇  ρD∇Yl +W ̇ 000

l (66)

Conservation of Momentum

∂u

∂t

+ u ω + ∇H =

1

ρ ((ρ  ρ∞)g + f + ∇  τ) (67)

Divergence Constraint

∇  u =

γ  1

γ p0 

q ̇000 + ∇  k∇T + ∇ 

∑

l

cp;lTρD∇Yl 

1

γ  1

d p0

dt !

(68)

Equation of State

p0(t ) = ρTR ∑

l

Yl=Ml (69)

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.

6.2 Temporal Discretization

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

n,

Yn

i , un, H n, and pn

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

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.

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

∇2

H n = 

(∇  u)

(n+1)e  (∇  u)

n

δt  ∇  Fn (71)

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 + ∇H n) (72)

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 

(73)

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

The divergence (∇  u)

(n+1)

is derived from the corrected thermodynamic quantities.

4. The pressure is recomputed using estimated quantities

∇2

H (n+1)e = 

2(∇  u)

n+1  (∇  u)

(n+1)e  (∇  u)

n

δt  ∇  F(n+1)e (75)

The velocity is then corrected

un+1 =

1

2 h

un + u(n+1)e  δt F(n+1)e + ∇H (n+1)ei (76)

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

divergence.

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

ijk 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

ijk is the x component of velocity

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

n

i1; jk is at the backward pointing face of the i jkth cell.

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

μijk = ρijk (Cs ∆)

2 jSj (77)

where Cs is an empirical constant, ∆ = (δxδyδz)

1

3 , and

jSj

2 = 2

∂u

∂x 2

+ 2

∂v

∂y2

+ 2

∂w

∂z 2

+



∂u

∂y

+

∂v

∂x2

+



∂u

∂z

+

∂w

∂x 2

+



∂v

∂z

+

∂w

∂y 2

(78)

The quantity jSj 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

kijk =

cp;0 μijk

Pr ; (ρD)ijk =

μijk

Sc (79)

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, Cs is 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

μijk = ∑

l

Yl ;ijk μl (Tijk) (80)

kijk = ∑

l

Yl ;ijk kl (Tijk) (81)

Dl;ijk = Dl0(Tijk) (82)

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.

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

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

ijk  ρn

ijk

δt + (u  ∇ρ)

n

ijk = ρnijk(∇  u)

n

ijk (83)

Following the prediction of the velocity and background pressure at the (n + 1)st time level, the density is

corrected

ρ(n+1)

ijk  1

2 

ρn

ijk + ρ(n+1)e

ijk 

1

2 δt + (u  ∇ρ)

(n+1)e

ijk = ρ(n+1)eijk (∇  u)

(n+1)e

ijk (84)

The species conservation equations are differenced the same way

(ρYl )

(n+1)e

ijk  (ρYl )

n

ijk

δt + (u  ∇ρYl )

n

ijk = (ρYl )

n

ijk(∇  u)

n

ijk + (∇  ρD∇Yl )

n

ijk +W ̇ 000

ijk (85)

at the predictor step, and

(ρYl )

(n+1)

ijk  1

2 

(ρYl )

n

ijk + (ρYl )

(n+1)e

ijk 

1

2 δt + (u  ∇ρYl)

(n+1)e

ijk = (ρYl )

(n+1)e

ijk (∇  u)

(n+1)e

ijk + (∇  ρD∇Yl)

(n+1)e

ijk +W ̇ 000

ijk

(86)

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  ∇ρ)ijk =

1  εu

2 uijk

ρi+1; jk  ρijk

δx

+

1  εu

2 ui1; jk

ρijk  ρi1; jk

δx

+

1  εv

2 vijk

ρi; j+1;k  ρijk

δy

+

1  εv

2 vi; j1;k

ρijk  ρi; j1;k

δy

+

1  εw

2 wijk

ρi j;k+1  ρijk

δz

+

1  εw

2 wi j;k1

ρijk  ρi j;k1

δz

(87)

(u  ∇ρYl )ijk =

1  εu

2 uijk

(ρYl )i+1; jk  (ρYl )ijk

δx

+

1  εu

2 ui1; jk

(ρYl )ijk  (ρYl )i1; jk

δx

+

1  εv

2 vijk

(ρYl )i; j+1;k  (ρYl )ijk

δy

+

1  εv

2 vi; j1;k

(ρYl )ijk  (ρYl )i; j1;k

δy

+

1  εw

2 wijk

(ρYl )i j;k+1  (ρYl )ijk

δz

+

1  εw

2 wi j;k1

(ρYl )ijk  (ρYl )i j;k1

δz

(88)

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

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

The divergence in both the predictor and corrector step is discretized

(∇  u)ijk =

γ  1

γp0 

q ̇000

ijk + (∇  k∇T )ijk + ∑

l

(∇  Tcp;lρD∇Yl )ijk 

1

γ  1 

d p0

dt !

(89)

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

(∇  k∇T )ijk =

1

δx 

ki+ 1

2 ; jk

Ti+1; jk  Tijk

δx  ki 1

2 ; jk

Tijk  Ti1; jk

δx  +

1

δy 

ki; j+ 1

2 ;k

Ti; j+1;k  Tijk

δy  ki; j 1

2 ;k

Tijk  Ti; j1;k

δy  +

1

δz 

ki j;k+1

2

Ti j;k+1  Tijk

δz  ki j;k1

2

Tijk  Ti j;k1

δz 

(90)

(∇  cp;lTρD∇Yl )ijk =

cp;l

δx 

Ti+ 1

2 ; jk ρDl ;i+ 12 ; jk

Yl;i+1; jk Yl ;ijk

δx  Ti 1

2 ; jk ρDl;i 12 ; jk

Yl ;ijk Yl;i1; jk

δx  +

cp;l

δy 

Ti; j+ 1

2 ;k ρDl ;i; j+12 ;k

Yl ;i; j+1;k Yl;ijk

δy  Ti; j 1

2 ;k ρDl;i; j12 ;k

Yl;ijk Yl ;i; j1;k

δy  +

cp;l

δz 

Ti j;k+1

2

ρDl ;i j;k+1

2

Yl;i j;k+1 Yl ;ijk

δz  Ti j;k1

2

ρDl;i j;k1

2

Yl ;ijk Yl;i j;k1

δz 

(91)

(∇  ρD∇Yl )ijk =

1

δx 

ρDl;i+ 1

2 ; jk

Yl;i+1; jk Yl;ijk

δx  ρDl ;i 1

2 ; jk

Yl ;ijk Yl ;i1; jk

δx  +

1

δy 

ρDl;i; j+1

2 ;k

Yl;i; j+1;k Yl ;ijk

δy  ρDl;i; j1

2 ;k

Yl ;ijk Yl ;i; j1;k

δy  +

1

δz 

ρDl;i j;k+1

2

Yl;i j;k+1 Yl;ijk

δz  ρDl ;i j;k1

2

Yl ;ijk Yl;i j;k1

δz 

(92)

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

Tijk =

p0

ρijkR ∑N

l=0(Yl;ijk=Ml )

(93)

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

N

∑

l=0

Yl;ijk

Ml =

1

M0

+

N

∑

l=1

1

Ml 

1

M0 

Yl (94)

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

ρijk =

p0

T∞R ∑N

l=0Yl;ijk=Ml

(95)

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

ρijk =

M0 p0

T∞R +

N

∑

l=1

1 

M0

Ml 

(ρYl )ijk (96)

24
6.5.2 Heat Release Rate (LES)

For an LES calculation, heat is added to the flow domain through the use of heat releasing Lagrangian

particles, or thermal elements. The thermal elements are ejected from burning surfaces with a normal

velocity that is either specified by the user or automatically determined based on the specified heat release

rate per unit area, the heat of combustion, and the density of the fuel gases

un =

q ̇00

r

∆H ρf

(97)

The thermal elements release energy at a constant rate. A specified fraction of the energy is emitted as

thermal radiation, and this energy can be re-absorbed by the smoke-laden gases or by the walls if desired.

The non-radiated energy from the thermal elements is interpolated on the computational grid. As input, the

user specifies the fraction of that energy lost as thermal radiation, χr, and the burn-out time of the elements

tb. The heat release rate per unit volume of the i jkth grid cell is given by the non-radiated energy

q ̇000

ijk =

∑m(1  χr) q ̇p;m

δxδyδz

;  ̇qp;m =

q ̇00

f

tb n ̇00 (98)

where  ̇q00

f is the heat release rate per unit area assigned to the surface from which the mth element originated,

n ̇00 is the number of thermal elements introduced per unit time per unit area at this same surface, and the

summation is over all thermal elements within the grid cell whose indices are i jk.

If desired, the radiated fraction of the energy from the thermal elements can be re-absorbed by the

smoke-laden gases, in which case an additional contribution to the heat release rate per unit volume in a

given grid cell is given by

q ̇000

ijk = κijk

Np

∑

m=1

χr q ̇p;m

4πjxp;m  xijk j

2 e

R

κ(l) dl (99)

Note that here the summation is carried out over all (or a sampling) of the thermal elements, not just those

within the i jkth cell. The absorption coefficient κ is computed at the center of the grid cell. It is based on

the mass of particulate matter and the temperature within that cell

κijk = 1186Tijk fv (100)

Here fv is the soot volume fraction, given by

fv =

∑mp;m

ρsoot δxδyδz

(101)

where ρsoot is the density of the soot, and mp;m is the particulate mass carried by the mth thermal element

mp;m =

χs q ̇p;m tb

∆H

max

t t0;m

tb ; 1

(102)

Here χs is the soot yield of the given fuel, and t

0;m is the time when the mth thermal element was introduced

at the burning surface.

If sprinklers are flowing, the water droplets can cool the hot gases. The heat release rate per unit volume

is decreased

q ̇000

ijk = 

∑Adhd (Tijk  Td )

δxδyδz

(103)

where the summation is carried out over all sprinkler droplets in a grid cell of volume δxδyδz. Here T is the

gas temperature, Td is the droplet temperature, Ad = 4πr2

d is the surface area of a droplet and hd = Nuk=2rd

is a heat transfer coefficient.

25
6.5.3 Heat Release Rate (DNS)

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

νCxHy CxHy + νO2 O2 ! νCO2 CO2 + νH2O H2O (104)

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

n and Yn

CxHy ;ijk  YF (t

n) and Yn

O2 ;ijk  YO(t

n), the

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

dYF

dt = 

Bρa+b1

ijk

Mb

O Ma1

F

YF (t )

aYO(t )b eE=RTijk (105)

dYO

dt = 

νO MO

νF MF

dYF

dt (106)

The temperature Tijk and density ρijk are fixed at their values at time t

n and the ODE is iterated from tn 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 ̇

000 n

ijk = ∆H ρnijk

YF (t

n) YF (tn+1)

δt (107)

where δt = t

n+1  tn. The species mass fractions are adjusted at this point in the calculation (before the

convection and diffusion update)

Yn

l;ijk = Yl (t

n) 

νl Ml

νF MF

(YF (t

n) YF (tn+1)) (108)

6.5.4 Thermal 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

Tn+1 w = Tn

w + δts

q ̇00

c + q ̇00r  q ̇00rr

ρscsδ (109)

where Tw is the wall temperature, δts is 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 (110)

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 = CjTgas  Twj

1

3 (Tgas  Tw) W/m2 (111)

26
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 radiative flux to the wall is given by

q ̇00

r =

Np

∑

m=1

χr q ̇p;m cosφm

4π jxp;m  xj

e

R

κ(l) dl (112)

where  ̇qp;m is the heat release rate of the mth thermal element and φm is the angle formed by the normal to

the surface and the line connecting the thermal element and the point on 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) = sξi +

1  s

δ2 ξ3

i (113)

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 12 = (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

Tn+1 s;i  Tn

s;i

δt =

α

2δxi 

Tn

s;i+1  Tn

s;i

δxs;i+ 1

2 

Tn

s;i  Tn

s;i1

δxs;i 1

2 +

Tn+1 s;i+1  Tn+1 s;i

δxi+ 1

2 

Tn+1 s;i  Tn+1 s;i1

δxi 1

2 +!

(114)

for 1  i  N. The boundary condition is discretized ks

Tn+1 s;1  Tn+1 s;0

δx 1

2 = q ̇00

c + q ̇00r  εσ h

Tn4

s;

1

2  T4

∞ + 4Tn3

s;

1

2 

Tn+1

s;

1

2  Tn

s;

1

2 i (115)

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

εσ h

T (n+1)

4

s;

1

2  T4

∞ i  εσ h

Tn4

s;

1

2  T4

∞ + 4Tn3

s;

1

2 

Tn+1

s;

1

2  Tn

s;

1

2 i (116)

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

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

k

Tgas  Tghost

δn = CjTgas  Twj

1

3 (Tgas  Tw) (118)

where δn is the distance between the center of the ghost cell and the center of the gas cell. This equation is

solved for Tghost, so that when the conservation equations are updated, the amount of heat lost to the wall is

equivalent to the empirical expression on the right hand side. Note that Tghost is purely a numerical construct.

It does not represent the temperature within the wall, but rather establishes a temperature gradient at the wall

consistent with the empirical correlation.

27
6.5.5 Species Boundary Conditions

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

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  ̇m00

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

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

(121)

where  ̇m00

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.

6.5.6 Density Boundary Condition

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 )

(122)

28
6.6 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 

(123)

∂v

∂t

+ Fy +

∂H

∂y = 0 ; Fy = uωz  wωx 

1

ρ 

fy +

∂τyx

∂x

+

∂τyy

∂y

+

∂τyz

∂z 

(124)

∂w

∂t

+ Fz +

∂H

∂z = 0 ; Fz = vωx  uωy 

1

ρ 

fz +

∂τzx

∂x

+

∂τzy

∂y

+

∂τzz

∂z 

(125)

The spatial discretization of the momentum equations take the form

∂u

∂t

+ Fx;ijk +

Hi+1; jk  Hijk

δx = 0 (126)

∂v

∂t

+ Fy;ijk +

Hi; j+1;k  Hijk

δy = 0 (127)

∂w

∂t

+ Fz;ijk +

Hi j;k+1  Hijk

δz = 0 (128)

where Hijk is taken at center of cell i jk, uijk and Fx;ijk are taken at the side of the cell facing in the forward x

direction, vijk and Fy;ijk at the side facing in the forward y direction, and wijk and Fz;ijk 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 ( f

x ; fy ; fz); and the upwinding bias terms εu, εv, and

εw are located at the respective cell faces.

Fx;ijk =



1  εw

2 wi+ 1

2 ; jk ωy;ijk +

1  εw

2 wi+ 1

2 ; j;k1 ωy;i j;k1





1  εv

2 vi+ 1

2 ; jk ωz;ijk +

1  εv

2 vi+ 1

2 ; j1;k ωz;i; j1;k



1

ρi+ 1

2 ; jk 

fx;ijk +

τxx;i+1; jk  τxx;ijk

δx

+

τxy;ijk  τxy;i; j1;k

δy

+

τxz;ijk  τxz;i; j;k1

δz 

(129)

Fy;ijk =



1  εu

2 ui; j+ 1

2 ;k ωz;ijk +

1  εu

2 ui1; j+1

2 ;k ωz;i1; jk





1  εw

2 wi; j+ 1

2 ;k ωx;ijk +

1  εw

2 wi; j+ 1

2 ;k1 ωx;i j;k1



1

ρi; j+ 1

2 ;k 

fy;ijk +

τyx;ijk  τyx;i1; jk

δx

+

τyy;i; j+1;k  τyy;ijk

δy

+

τyz;ijk  τyz;i; j;k1

δz 

(130)

Fz;ijk =



1  εv

2 vi j;k+1

2

ωx;ijk +

1  εv

2 vi; j1;k+1

2

ωx;i; j1;k





1  εu

2 ui j;k+1

2

ωy;ijk +

1  εu

2 ui1; j;k+1

2

ωy;i1; jk



1

ρi j;k+1

2 

fz;ijk +

τzx;ijk  τzx;i1; jk

δx

+

τzy;ijk  τzy;i; j1;k

δy

+

τzz;i j;k+1  τzz;ijk

δz 

(131)

ωx;ijk =

wi; j+1;k  wijk

δy 

vi j;k+1  vijk

δz

(132)

29
ωy;ijk =

ui j;k+1  uijk

δz 

wi+1; jk  wijk

δx

(133)

ωz;ijk =

vi+1; jk  vijk

δx 

ui; j+1;k  uijk

δy

(134)

τxx;ijk = μijk 

2

uijk  ui1; jk

δx 

2

3 (∇  u)ijk  μijk 

4

3 (∇  u)ijk  2

vijk  vi; j1;k

δy  2

wijk  wi j;k1

δz 

(135)

τyy;ijk = μijk 

2

vijk  vi; j1;k

δy 

2

3 (∇  u)ijk  μijk 

4

3 (∇  u)ijk  2

uijk  ui1; jk

δx  2

wijk  wi j;k1

δz 

(136)

τzz;ijk = μijk 

2

wijk  wi j;k1

δz 

2

3 (∇  u)ijk  μijk 

4

3 (∇  u)ijk  2

uijk  ui1; jk

δx  2

vijk  vi j1;k

δy 

(137)

τxy;ijk = τyx;ijk = μi+ 1

2 ; j+ 12 ;k 

ui; j+1;k  uijk

δy

+

vi+1; jk  vijk

δx 

(138)

τxz;ijk = τzx;ijk = μi+ 1

2 ; j;k+ 12 

ui j;k+1  uijk

δz

+

wi+1; jk  wijk

δx 

(139)

τyz;ijk = τzy;ijk = μi; j+ 1

2 ;k+ 12 

vi j;k+1  vijk

δz

+

wi; j+1;k  wijk

δy 

(140)

εu =

uδt

δx

(141)

εv =

vδt

δy

(142)

εw =

wδt

δz

(143)

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. (89) is identically equal to the divergence defined by

(∇  u)ijk =

uijk  ui1; jk

δx

+

vijk  vi; j1;k

δy

+

wijk  wi j;k1

δz

(144)

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.

6.6.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;ijk =

1

2

∑ρCdπr2

d (ud  uijk)jud  uj

δxδyδz  (ρi+ 1

2 ; jk  ρ∞)gx (145)

fy;ijk =

1

2

∑ρCdπr2

d (vd  vijk)jud  uj

δxδyδz  (ρi; j+ 1

2 ;k  ρ∞)gy (146)

fz;ijk =

1

2

∑ρCdπr2

d (wd  wijk)jud  uj

δxδyδz  (ρi j;k+1

2  ρ∞)gz (147)

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.

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

uijk ;

δy

vijk ;

δz

wijk ;

ρijk δx2

8μijk ;

ρijk δy2

8μijk ;

ρijk δz2

8μijk 

(148)

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.

6.7 The Pressure Equation

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

Hi+1; jk  2Hijk + Hi1; jk

δx2 +

Hi; j+1;k  2Hijk + Hi; j1;k

δy2 +

Hi j;k+1  2Hijk + Hi j;k1

δz2 = 

Fx;ijk  Fx;i1; jk

δx 

Fy;ijk  Fy;i; j1;k

δy 

Fz;ijk  Fz;i j;k1

δz 

∂

∂t (∇  u)ijk (149)

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 [27]. To ensure that the divergence of the fluid

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

∂

∂t (∇  u)ijk =

(∇  u)

(n+1)e

ijk  (∇  u)

n

ijk

δt (150)

at the predictor step, and then

∂

∂t (∇  u)ijk =

2(∇  u)

n+1

ijk  (∇  u)

(n+1)e

ijk  (∇  u)

n

ijk

δt (151)

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

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

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

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 (152) 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 (153)

31
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 (153)

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

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 = (u2

I ; jk + v2

I; j 1

2 ;k + w2

I; j 1

2 ;k )=2 uI; jk > 0

HI+ 1

2 ; jk = 0 uI; jk < 0 (155)

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

6.8 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 (157)

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 = xn

p + δt un (158)

xn+1

p =

1

2 

xn

p + x(n+1)e

p + δt u(n+1)e

(159)

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

position.

32
7 Conclusion

The equations and numerical algorithm described in this document form the core of an evolving fire model.

As research into specific fire-related phenomena continues, the relevant parts of the model can be improved.

Because the model was originally designed to 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 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.

33
8 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

k thermal conductivity; suppression decay factor

M molecular weight of the gas mixture

Mi molecular weight of ith gas species

m ̇ 00

w water flux per unit area

m00

w water mass per unit area

Nu Nusselt number

Pr Prandtl number

p pressure

p0 background pressure

p ̃ pressure perturbation

q ̇000 heat release rate per unit volume

q ̇00

f fire heat release rate per unit area

q ̇00

r radiative flux to a solid surface

q ̇00

c convective flux to a solid surface

q ̇00

rr radiative loss from a solid surface

q ̇00

w water cooling per unit area

R universal gas constant

Re Reynolds number

rd water droplet radius

RTI Response Time Index of sprinkler

S deformation tensor

Sc Schmidt number

T temperature

t time

tb thermal element burn-out time (LES)

u = (u; v; w) velocity vector

W ̇ 000

i production rate of ith species per unit volume

We Weber number

x = (x; y; z) position vector

Yi mass fraction of ith species

γ ratio of specific heats

34
∆H heat of combustion

δ wall thickness

ε emissivity

κ absorption coefficient

μ dynamic viscosity

νi stoichiometric coefficient

ρ density

τ viscous stress tensor

χr radiative loss fraction

χs smoke or soot yield

σ Stefan-Boltzmann constant

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

35
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] A. Atreya and S. Agrawal. Effect of Heat Loss on Diffusion Flames. Combustion and Flame, 115:372–

382, 1998.

[14] H.R. Baum and B.J. McCaffrey. Fire Induced Flow Field – Theory and Experiment. In Fire Safety Science – Proceedings of the Second International Symposium, pages 129–148. International Association

for Fire Safety Science, 1989.

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

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

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

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

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

Technical report, University of Maryland, 1999.

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

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

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

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

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

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

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

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

37