Fire Safety Journal 40 (2005) 213–244

The critical ventilation velocity in tunnel fires—a

computer simulation

C.C. Hwang, J.C. Edwards

National Institute for Occupational Safety and Health, Pittsburgh Research Laboratory, Pittsburgh,

PA 15236-0070, USA

Received 18 July 2003; received in revised form 20 September 2004; accepted 25 November 2004

Available online 18 January 2005

Abstract

In ventilated tunnel fires, smoke and hot combustion products may form a layer near the

ceiling and flow in the direction opposite to the ventilation stream. The existence of this reverse

stratified flow has an important bearing on fire fighting and evacuation of underground mine

roadways, tunnels and building corridors. In the present study, conducted by the National

Institute for Occupational Safety and Health, a CFD program (fire dynamics simulator) based

on large eddy simulations (LES) is used to model floor-level fires in a ventilated tunnel.

Specifically, the critical ventilation velocity that is just sufficient to prevent the formation of a

reverse stratified layer is simulated for two tunnels of different size. The computer code is

verified by checking the computed velocity profile against experimental measurements. The

CFD results show the leveling-off of the critical ventilation velocity as the heat release rate

surpasses a certain value. At this critical ventilation, the ceiling temperature above the fire

reaches a maximum for both tunnels. The velocity leveling-off can be explained from this

observation. An extended correlation of Newman (Combust. Flame 57 (1984) 33) is applied to

the temperature profiles obtained by CFD. At the critical ventilation, temperature

stratification exists downstream from the fire. The computed critical ventilation velocity

shows fair agreement with available experimental data taken from both horizontal and

inclined fire tunnels. The CFD simulations indicate that the Froude modeling is an

approximation for tunnel fires. The Froude-scaling law does not apply to two geometrically

similar fire tunnels. The CFD results are compared with two simple theories of

critical ventilation by Kennedy et al. (ASHRAE Trans. Res. 102(2) (1996) 40) and Kunsch

ARTICLE IN PRESS

www.elsevier.com/locate/firesaf

0379-7112/$ - see front matter Published by Elsevier Ltd.

doi:10.1016/j.firesaf.2004.11.001

Corresponding author

E-mail address: ckh9@cdc.gov (C.C. Hwang).
ARTICLE IN PRESS

Nomenclature

A tunnel cross-sectional area, m2

cp specific heat, kJ kg1 K1

Fr Froude number, dimensionless

g acceleration of gravity, ms2

H tunnel height, m

H ̄ hydraulic tunnel height, defined as the ratio of 4 times the tunnel crosssectional area to the tunnel perimeter, m

k thermal conductivity, Wm1 K1

m_ mass flow rate, kgs1

P pressure, Pa

q_00 heat flux, Wm2

Q chemical heat release rate, kW

Qc convective heat release rate, kW

Q dimensionless heat release rate

t time, s

T temperature, K

u ventilation velocity, ms1

u dimensionless ventilation velocity

V velocity, ms1

W tunnel width, m

x,y,z Cartesian coordinates, m

Greck letters

r air density, kg m3

f tunnel inclination angle, degree

O a parameter defined in Eq. (9)

DT temperature difference associated with stratification, K

DT

0 constant with a value equal to 6.13, dimensionless

Subscript

avg average

b blackbody

c convective

cf difference between ceiling and floor

cr critical

f floor

h temperature rise above ambient near the ceiling

in tunnel inlet

n nth band of radiation intensity

O ambient

214 C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244
(Fire safety J. 37 (2002) 67).

Published by Elsevier Ltd.

Keywords: Tunnel Fire; Critical ventilation velocity; Computer simulation; Stratified flow; Froude

modeling; Large eddy simulation

1. Introduction

1.1. Critical ventilation velocity

When a fire is started on the floor of a straight tunnel without a ventilation (cross)

flow, a hot plume rises above the fire and entrains the surrounding cold air into the

plume. The plume, upon reaching the ceiling, forms two gas streams flowing in

opposite directions along the ceiling. When a cross ventilation current exists, the

symmetry in the rising plume and in the ceiling gas streams is broken. The ventilation

current bends the plume and the length of the ceiling layer flowing against the

ventilation current is reduced. This situation is depicted in Fig. 1, where the ceiling

layer is represented by the gas temperature contours. In Fig. 1, the velocity vectors

are also shown to indicate the general flow pattern in the tunnel.

In the event of a tunnel fire or smoke emergency, a main concern is to maintain

an evacuation path that is free of smoke and hot gases. The existence of reverse

stratified layer (also called back-layer) of hot combustion products has an

important bearing on fire fighting and evacuation of underground mine roadways,

tunnels and building corridors. Consider a scenario involving a stopped vehicle

on fire in a tunnel, disrupting traffic and requiring worker and passenger

evacuation. Another scenario involves a conveyer-belt fire in an underground

mine entry, producing smoke and toxic combustion products. The consideration

ARTICLE IN PRESS

ref reference

s material

Fig. 1. Stratified ceiling layer of combustion products against the ventilation current from left, with the

fire source at the floor.

C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244 215
for emergency planning must focus on determining the ventilation required to

maintain a single evacuation path from the fire source clear of smoke and hot

gases.

Experimental data show that the length of a back-layer upstream of the fire source

is a function of uin for fixed heat release rate from the fire [1–3]. As the value of uin

increases, the length of the back-layer decreases. The ‘‘critical ventilation velocity’’

ucr is defined as the value of the ventilation velocity uin that is just able to prevent the

formation of a back-layer. In the present CFD study, the ventilation velocity for a

given Q is ucr when the front end of the back-layer is at x 1⁄4 0 (front end of the fire

source).

The risk from accidental fires as well as the subsequent smoke movement depends

largely on this applied ventilation current uin: It is of practical importance to

understand the physical parameters and flow conditions under which the reverse

stratified flow can be prevented.

1.2. Previous works

Previous authors employed simple empirical models to determine the critical

ventilation velocity needed to prevent upstream movement of smoke from a fire in a

tunnel [1,4–7]. In general, these models considered the buoyancy head and the

dynamic head in the system, and deduced appropriate quantities for correlation.

Hwang et al. [8] and Charters et al. [9] employed phenomenological models

that provided more detail than the simple empirical approaches. A recent review of

tunnel fires by Grant et al. [10] pointed out that existing experimental data still

show an inadequate fundamental understanding of the interaction between

buoyancy-driven combustion products and forced ventilation, the validity of

extrapolation of small-scale results to large scales, the influence of slope on smoke

movement, and the effect of tunnel geometry. Wu and Bakar [11] carried out

experimental tests of tunnel fires using tunnels of different cross sections. They

used the hydraulic diameter as the characteristic length in the dimensionless

groups for correlation. The correlation was able to encompass their own data and

large tunnel data of other workers. The correlation shows that at large values of Q;

ucr levels off. Kunsch [12] derived an expression that shows the decreasing effect of Q

on uin as Q increases. His equation shows that the aspect ratio of the tunnel crosssection is also a parameter.

1.3. Objectives of the present study

The present study addresses the problem of critical ventilation velocity in

longitudinally ventilated tunnels. Specifically, the leveling-off of ucr for large

values of Q is analyzed. Computer simulations are made to see whether the levelingoff of ucr can be observed. If leveling-off of ucr can be simulated, a possible cause

is searched and studied. The simulation results are compared with available

experimental data and simple theories.

ARTICLE IN PRESS

216 C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244
2. Theoretical considerations

2.1. Experimental data

Fig. 2 shows a plot of existing experimental data on the critical ventilation velocity

as a function of the heat release rate in fire tunnels. The tunnel size is expressed by

the hydraulic tunnel height, H ̄ ; defined as 4  (cross-sectional area)/(perimeter). The

values of H ̄ range from 0.18 to 7.72 m, a factor of 43. Note that the scaling factor H ̄

is not included in the correlation of Fig. 2. The factor H ̄ will be included in the

correlation in the next section. If we ignore the local variations in the plot, the global

trend can be represented by an equation

ucr 1⁄4 aQ1=5

with ucr in m/s, Q in kW, and a; a constant. An explanation for this 1/5th power

correlation will be attempted in Sections after the flow field in the fire plume zone

have been studied. If one studies the details of an individual set of experimental data

in Fig. 2, especially those of Wu and Bakar [11] and Memorial tunnel [2], one notices

the leveling-off of ucr values as Q increases. Other authors also reported this trend

ARTICLE IN PRESS

X X XXXX XX

Q, kW

ucr, m/s

100 101 102 103 104 105

0.5

1

1.5

2

2.5

3

3.5

Memorial Tunnel

Buxton Lab.

Wu & Bakar, D

Wu & Bakar, A

Hwang & Wargo

1/5th power line

X

Fig. 2. Experimental data of critical ventilation velocity in tunnel fires. Memorial tunnel: H ̄ 1⁄4 7:72 m;

Buxton lab: H ̄ 1⁄4 2:38 m; Wu and Bakar, tunnel A: H ̄ 1⁄4 0:176 m; Wu and Bakar, tunnel D: H ̄ 1⁄4 0:40 m;

Hwang and Wargo: H ̄ 1⁄4 0:343 m:

C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244 217
[7,10]. A more satisfactory way to treat existing experimental data will be presented

in the next Section.

2.2. Scaling parameters

McCaffrey and Quintiere [13] studied buoyancy driven countercurrent flows

generated by a fire source, and found that the general nature of the flow is not

dependent on scale. This can be achieved by assuming geometric similarity and

assuming dynamic similarity holds if the Froude number is maintained constant. The

Froude number criterion leads to equal corresponding temperatures between the

scaled model (MODEL) and full-scale (FS) systems. The relationship

V

ffiffiffiffi

Hp



FS

1⁄4 V

ffiffiffiffi

Hp



MODEL

holds for velocity V where H is height, provided

Q

H5=2



FS

1⁄4 Q

H5=2



MODEL

;

where Q is the heat generated at the fire source. In the correlations of experimental

data for tunnel fires, the following dimensionless groups have been used [7,11]

u 1⁄4 uin

ffiffiffiffiffiffiffi

gH p and Q 1⁄4 Q

r0T0cp

ffiffiffiffiffiffiffiffiffi

gH5 p : (1)

To accommodate for the tunnel cross-sectional geometry, Wu and Bakar [11] used

the hydraulic tunnel height H ̄ in place of H: As shown in Fig. 3, plots of u

cr versus

Q encompass the tunnel sizes from 0.18 to 7.7 m. In Fig. 3, a line representing

u / ðQ

Þ

1=3 is drawn. This relationshipappears to hold for small values of Q

: The

dimensionless critical ventilation velocity u

cr is defined as the dimensionless value of

the ventilation velocity that is just sufficient to prevent the formation of a reverse

stratified layer upstream of the fire source. For low values of Q; u

cr increases with

Q: For high values of Q; u

cr is approximately constant and becomes independent of

the value of Q: This behavior of the critical-velocity leveling-off has been explained

as fire blockage to longitudinal ventilation [7]. Wu and Bakar [11] observed that the

temperature distribution in the fire plume changed as the value of Q was increased.

When the fire heat release rate is increased to certain level, the ‘intermittent’ flames

reach the tunnel ceiling and occupy the upper part of the tunnel. The intermittent

flames have the feature of constant flow speed; therefore the buoyancy force in the

back-layering is not sensitive to the heat release rate.

2.3. Temperature field in tunnel fires

The previous authors attribute the leveling-off of ucr to fire blocking or

intermittent flames. These explanations are related to the flow and temperature

fields above and downstream from the fire source. Therefore, a detailed study of the

temperature field downstream from the fire may reveal the cause for ucr leveling-off.

ARTICLE IN PRESS

218 C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244
Newman [14] studied experimentally tunnel fires and expressed his results as shown

in Figs. 4 and 5. It is noted that there is a similarity between the line drawn in Fig. 3

and a line connecting points in Fig. 5. In both diagrams, the dependent quantity

behaves differently as the independent quantity increases beyond a certain value.

This suggests that a quantity associated with the temperature field above the fire may

play a role in ucr leveling-off. Our first guess is the temperature stratification under

critical ventilation. Newman [14] described the temperature stratification in tunnel

fires by

DT

DTref

1⁄4 f ðFrÞ 1⁄4 f Vref

ðDT=TrefÞgH ̄  1=2

( ); (2)

where DT is a temperature difference associated with the stratification; Tref is the

reference temperature; Vref is the reference flow velocity; g is the acceleration due to

gravity; H ̄ is the characteristic dimension of the duct. The main features of

Newman’s results (Figs. 4 and 3a in [14]) are shown in Figs. 4 and 5. In Fig. 4, a

Froude number is used based on

DT 1⁄4 DTcf;Tref 1⁄4 Tavg; Vref 1⁄4 Vavg

ARTICLE IN PRESS

X

X

XX X X X X

Q*

u*

10-3 10-2 10-1 100

0.1

0.2

0.3

0.4

0.5

Memorial Tunnel

Buxton Lab.

Wu & Bakar, D

Wu & Bakar, A

Hwang & Wargo

1/3rd power line

X

cr

Fig. 3. Experimental data of critical ventilation velocity in tunnel fires using non-dimensional parameters.

u

cr appears to depend on ðQÞ

1=3 for small values of Q:

C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244 219
where DTcf is the difference between the gas temperature near the ceiling and the gas

temperature near the floor, and DTavg is the average temperature rise relative to

ambient. The correlation can be represented by the relationship

DTcf

DTavg

1⁄4 1:5 Vavg

DTcf=Tavg 	 
gH ̄  1=2

( )1

: (3)

Three regions can be defined in Fig. 4:

Region I: For Frp0:9;DTcf=DT avgX1:7; i.e., buoyancy dominating temperature

stratification.

Region II: For 0:9pFrp10; 1:7XDTcf=DTavgX0:12; i.e., significant interaction of

the ventilation velocity with the fire-induced buoyancy occurs.

Region III: For FrX10; DTcf=DTavgp0:12; i.e., stratification is insignificant.

Newman [14] also showed (Fig. 5) that for region I, DTcf=DTh 1⁄4 1:0 independent

of DTcf=DTavg; where DTh is the gas temperature rise above ambient near the ceiling.

This implies that the gas near the floor is essentially at ambient conditions and the

maximum temperature stratification is established. In a ventilated tunnel fire, an

increase in the heat generation Q increases the temperature stratification that can be

ARTICLE IN PRESS

Fr = Vavg/ [(∆Tcf/Tavg)gH]0.5

∆

Tcf/ ∆

Tavg

10-1 100 101

10-1

100

101

Region III

Region I Region II

Line passing

through data points

Fig. 4. Gas temperature stratification versus Froude number in tunnel fires. This diagram is based on Fig.

4 of Ref. [14].

220 C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244
expressed by an increase in the value of DTcf=DTh: Eventually, a condition is

reached at which DTcf=DTh 1⁄4 1:0; and additional increase in Q cannot change the

value of DTcf=DTh:

2.4. u2Q relation

When a fire exists in a tunnel of uniform cross-section, the change in the average

gas velocity is caused by the change in the gas density. Consider a control volume

enclosed by two tunnel cross-sections and the connecting internal walls. The

conservation of mass for the control volume gives

Vavg 1⁄4 rin

ravg

uin 
 Tavg

Tin

uin; (4)

where the subscripts ‘‘in’’ designate the condition at the tunnel entrance. Solving for

uin from Eqs. (3) and (4) yields

uin 1⁄4 1:5

Tin

Tavg  DTavgDTcf  DTcfTavg gH ̄

 1=2

(5)

ARTICLE IN PRESS

∆Tcf /∆Tavg

∆

Tcf /∆

Th

10-1 100 101

0.5

1

1.5

2

Line passing

through data points

Region II Region I

Fig. 5. Gas temperature stratification heptane data in tunnel fires. This diagram is based on Fig. 3a of Ref.

[14].

C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244 221
Neglecting changes in the kinetic energy and the potential energy, an energy balance

for the control volume gives:

Qc 1⁄4 cpm_ inDTavg 1⁄4 cprinuinADTavg 1⁄4 cprinuin

p

4 H ̄ 2

DTcf

DTavg

DTcf

; (6)

where Qc is the convective heat release rate. It is assumed that the mass addition

from the fire source is negligible compared with m_ in: Solving for DTcf yields

DTcf 1⁄4 DTcf

DTavg  4Qc

princpuinH ̄ 2

 !: (7)

Substituting Eq. (7) into Eq. (5) and solving for uin gives

uin 1⁄4 1:5

4

p

1=2 Tin

Tavg 3=2 DTavgDTcf 1=2 gQc

rincpTinuinH ̄

1=2

: (8)

Eq. (8) can be cast into a non-dimensional form

uin

ffiffiffiffiffiffiffi

gH ̄ p 1⁄4 1:42 Tin

Tavg  DTavgDTcf 1=3 Qc

rincpTin

ffiffiffiffiffiffiffiffiffi

gH ̄ 5

q

0

B@

1

CA

1=3

or in terms of the dimensionless groupdefined by Eq. (1)

u 1⁄4 1:42OðQ

c Þ

1=3 (9)

where

O 1⁄4 Tin

Tavg  DTavgDTcf 1=3

Eq. (9) which relates u to Q is an extended correlation of Eq. (3). In theory,

Eq. (9) applies to a cross-section downstream from a tunnel fire, and should apply

to the situation when u is the critical ventilation velocity u

cr: This must be verified in

the present simulation. Since the values of Tin=Tavg and DTavg=DTcf are unknown at

this point, a full discussion of Eq. (9) will be delayed until this information is

available.

3. CFD simulation

3.1. Turbulence model

In our previous study of the reverse stratified flow generated by a floor-level fire,

the standard k   turbulence model was employed for simulation [15]. The

advantages of this model are its simplicity and cost effectiveness. For fire

applications, one of fundamental limitations of this model is the averaging

procedure at the root of the model equations. Since the k   model was developed

ARTICLE IN PRESS

222 C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244
as a time-averaged approximation to the conservation equations of fluid dynamics,

the results of fire simulations appear smooth. The evolution of eddy structures

characteristic of most fire plumes is lost with such an approach [16]. Other known

drawbacks of k   model are:

 The model assumes high Reynolds number and has to be adapted for low

Reynolds number flows.

 The model assumes an isotropic eddy viscosity and has to be modified to account

for the effect of the mean streamline curvature.

 Wall functions (algebraic laws) are generally required to describe flow fields in the

vicinity of walls.

Workers dealing with buoyancy-driven flows have used modified versions of the

k   model [11,17–19]. The modification in general was made by inclusion of a

buoyancy production term in the governing equations. The present study employed a

fire dynamics simulator (FDS Version 2) code based on large eddy simulations (LES)

techniques [20]. The application of LES techniques to fire is aimed at extracting

greater temporal and spatial fidelity from simulation of fire. This model explicitly

calculates the turbulent large scales and models the effects of smaller ones using sub

grid closure rules. Compared to direct numerical simulations (DNS), the description

of the unresolved small scales is lost. Compared to Reynolds ensemble averaging,

LES provides the instantaneous resolved field. LES of reacting flows can resolve the

instantaneous position of a large-scale flame, so that LES captures the low-frequency

variations of flow parameters. The approach based on LES has a particular

advantage over the Reynolds-averaging procedures in that only the effects of smallscale turbulence motion have to be modeled [21]. The FDS code adopts the ‘‘low

Mach number’’ combustion equations that describe the low speed motion of a gas

driven by chemical heat release and buoyancy forces. The balance equations for LES

are obtained by filtering the instantaneous balance equations. In the FDS code, the

Smagorinsky subgrid-scale model is employed [22]. Although the Smagorinsky

model is known as being too dissipative [23], the LES model is believed to be a better

model than a k   model for fire simulation. The governing equations for this model

can be found in Ref. [16].

3.2. Combustion model

The FDS model assumes that the combustion is mixing-controlled and that all

species of interest can be described in terms of a mixture fraction [16]. The mixture

fraction is a conserved quantity representing the fraction of material at a given

point that originated as fuel. The local heat release rate is computed from the local

oxygen consumption rate at the flame surface. The flame surface can be located by

requiring that fuel and oxidizer simultaneously vanish. The thermal radiation is

computed by solving the radiative transport equation (RTE) for a non-scattering gas

[16]. In practice, the radiation intensity integrated over a band is solved. Over the

wavelength of 1–200 mm, the number of bands is from 6 to 10.

ARTICLE IN PRESS

C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244 223
3.3. Description of fire tunnels

In the present study, simulations were performed for two tunnels for which

experimental data were available. Fig. 6 shows the schematics of the tunnels with

their fire-source geometries. The first (small) tunnel was 4.9-m long with a uniform

cross-section 0.4 m wide and 0.3 m high. A fire tunnel of this size was employed to

measure distributions of the temperature and velocity in reverse stratified layers

generated by floor-level fires [3]. Assigning a volumetric flow rate at the tunnel

inlet and atmospheric pressure at the outlet simulated the ventilation flow. In the

experiment, natural gas was the fire source burning stoichiometrically at a floor area

of 0.15 m axial length and 0.4 m width. In the simulation, methane and propane were

used as the fuel. The heat generation per unit area at the fire source (HRRPUA in

FDS2) provided the fire intensity for simulation.

ARTICLE IN PRESS

Fig. 6. Schematics of fire tunnels and fire-source configurations. (a) small tunnel; (b) large tunnel. The

direction of g (acceleration of gravity) as shown is for ascensional ventilation (uphill slope).

224 C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244
The second (large) tunnel used for simulation was 853 m long, 8.6 m wide and

7.0 m high. The tunnel size was selected to simulate experiments conducted at

Memorial tunnel [2]. The tunnel used in the experiments has an arched ceiling. The

height of the tunnel for simulation was selected so that the cross-sectional area

(60.2 m2

) is equal to that of the experimental tunnel. In the experiments, the

fire source was No. 2 fuel oil in various-size pans, with the oil surface at 0.9 m

from the floor. In simulation, the fire source was liquid heptane burning

stoichiometrically from a plane at a height 1.0 m from the floor. The width

of the pan was fixed at 2.6 m wide, so that the length of the pan varied with

the total heat generation. For heptane fuel, the heat of vaporization is

specified so that the burning rate of fuel is dependent on heat feed back. The data

are given in the database of FDS code. Propane was also used in the simulation. For

example, the length of the pan for a 50 MW fire is approximately 12 m. In the

experiments, jet fans near the tunnel exit provided longitudinal ventilation. In the

simulation, a volumetric flow rate was assigned at the exit and atmospheric pressure

was assigned at the inlet.

It is seen that the fire-source geometry for these two tunnels are entirely

different. The small tunnel has the fire source spread over the entire width

of the floor, and the fire intensity is varied by varying the heat generation

per unit area. The fire source for the large tunnel, on the other hand, does

not occupy the entire width of the floor. It has a fixed heat generation per unit area,

and the total fire intensity is linearly proportional to the area of the fire source.

This implies that as the fire intensity increases, the longitudinal extent of the fire

source increases.

3.4. Boundary conditions

The surface material of the tunnel is assumed thermally thick. A one-dimensional

heat conduction equation for the material temperature, Tsðn;tÞ; is applied in

the direction ~n pointing into air/solid interface (n 1⁄4 0; where n is the normal

coordinate)

rscs

qTs

qt 1⁄4 ks

q2Ts

qn2 ; ks

qTs

qn

ð0;tÞ 1⁄4 q_

00

c þ q_

00

r

where rs; cs and ks are the (constant) density, the specific heat and the thermal

conductivity of the material; and q_00

c and q_

00

r are the convective and radiative heat

fluxes at the surface [15].

For both tunnels, a uniform velocity is assigned at the inlet and atmospheric

pressure is assigned at the exit. In some computations for the large tunnel,

atmospheric pressure is assigned at the inlet and volumetric flow rate is assigned at

the exit. There was no difference in the computed value of the critical ventilation

velocity when the latter boundary conditions were used. At the tunnel wall, the gas

velocity has slipas given by the code. This implies that no wall function is used near

a wall to accommodate the no-slipboundary condition.

ARTICLE IN PRESS

C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244 225
3.5. Verification of the computer code

The computer code employed was FDS2 that has been verified in many

applications [24–28]. To see how well the code predicts the reverse-stratified flow

in a tunnel fire, CFD simulations were made to predict the experimental data of

Hwang and Wargo (small tunnel) [3] and the data obtained at the Memorial tunnel

(large tunnel) [2]. The descriptions of these tunnels are given in the previous section.

The coordinate system (see Fig. 6 for symbols) for the small tunnel was x from –3.05

to 1.85 m, y from –0.2 to 0.2 m, and z from 0 to 0.3 m; the coordinate system for the

large tunnel was x from –615 to 238 m, y from 4.3 to 4.3 m, and z from 0 to 7 m. To

reduce the computation time, in the computations to determine the critical

ventilation velocity, the range of x was 100 m upstream from the fire and

50–90 m downstream from the fire. Compared with the former configuration, it was

assumed that the inlet flow development length and an additional length downstream

of 50 m had small effect on the results. These assumptions were checked with actual

computations. Simulation times up to 300 s were run for both size tunnels. A steady

ARTICLE IN PRESS

u, m/s

z, m

-0.5 0 0.5

0

0.1

0.2

0.3

u, m/s

-0.5 0 0.5

0

0.1

0.2

0.3

u, m/s

-0.5 0 0.5

0

0.1

0.2

0.3

u, m/s

-0.5 0 0.5

0

0.1

0.2

0.3

u, m/s

-0.5 0 0.5

0

0.1

0.2

0.3

T, °C

z, m

20 40 60 80 100

0

0.1

0.2

0.3

x=0m

T, °C

20 40 60 80 100

0

0.1

0.2

0.3

x=-0.31m

T, °C

20 40 60 80 100

0

0.1

0.2

0.3

x=-0.61m

T, °C

20 40 60 80 100

0

0.1

0.2

0.3

x=-0.92m

T, °C

20 40 60 80

0

0.1

0.2

0.3

x=-1.23m

Fig. 7. Comparison of the experimental data (’) of the small tunnel [3] and its prediction (—) by CFD at

various cross sections. Upper diagram: velocity profiles. Lower diagram: temperature profiles. Heat

generation 1⁄4 3.3 kW.

226 C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244
state was established by visualizing the layer growth using the ‘Smokeview’ package

in FDS2. Since the measurements were made under steady conditions and an output

from LES simulations had flow oscillations, the results from the computations were

the mean values from the final 10- to 20-second of the computer outputs.

Fig. 7 shows the velocity and temperature profiles at various stations along the

small fire tunnel. Experimental results of Hwang and Wargo [3] are also plotted for

comparison. The computations predict the velocity profiles quite well close to the fire

zone. In the region upstream of the fire zone, the computations under-predict the

ceiling layer velocity and consequently under-predict the layer thickness. Overall,

the predicted temperatures agree well with the experiment. Fig. 8 shows the velocity

profiles and the temperature contours for the large fire tunnel. The computations

over-predict the ceiling layer velocity and the layer thickness. The discrepancy may

be attributed to the difference in the configurations of the ceiling region used in the

simulation and the actual tunnel. The shape of the temperature contours agrees well

with the experimental measurements [2].

Varying the number of grids for both tunnels checked the influence of number of

grid cells. For the tunnel having dimensions of 4.9 m  0.4 m  0.3 m, two sizes of

ARTICLE IN PRESS

u, m/s

-5 0 5

0

1

2

3

4

5

6

7

x=-189m

u, m/s

-5 0 5

0

1

2

3

4

5

6

7

x=-107m

u, m/s

-5 0 5

0

1

2

3

4

5

6

7

x=-62m

u, m/s

z, m

z, m

z, m

z, m

-5 0 5

0

1

2

3

4

5

6

7

x=-292m

30

3 60

6

60

90

90

90

9

120

120

120

120

15

150

150

150

150

21

210

210

210

33

330

330

390

390

39

51

510 510

x, m

z, m

-500 -400 -300 -200 -100 0 100 200

0

1

2

3

4

5

6

7

plane at y=0

30

60

60

30

0 210 390

(a)

(b)

Fig. 8. Results of simulations for the large fire tunnel. (a) Velocity profiles at various cross-section, y 1⁄4 0

plane. Experimental data: ’; CFD prediction: —. (b) Temperature-contour plot in y 1⁄4 0 plane. Heat

generation 1⁄4 50 MW.

C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244 227
grids were used. In the first case, the number of grids was 108  20  16, and in the

second case, the number of grids was 162  36  32. Table 1 shows a comparison of

these two cases. The conditions for both cases are close to the critical ventilation

conditions as indicated by a short length of the back layering. Both cases have

almost identical plots of temperature contours and the velocity vector fields. For the

large tunnel, efforts were made to decrease the size of the cells. As a guide,

approximately the cell size used by McGrattan et al. [27] in their simulation of a 2 m,

5-MW heptane fire was employed. The cell size for Case 2 (Table 2) was of the order

of 0.3 m. The cell size for Case 1 was 25% larger. Because of the scheme with which

the heat generation is computed in the code, the value of Q for Case 2 is larger by

approximately 4%, and the reverse ceiling-layer length is longer (11.1 versus 6.38 m).

Because of slightly different mixing patterns in the fire zone, Case 1 actually has a

higher temperature at the ceiling than in Case 2. It is judged that the cell sizes in

Tables 1 and 2 are sufficient in the simulation to predict the critical ventilation

velocities in tunnel fires. Both cell sizes (Cases 1 and 2) given in Tables 1 and 2 were

used in the computations reported in the Results Section.

4. Results

4.1. Critical ventilation velocity, u

cr

For a fixed heat generation rate Q; the program was run with a selected ventilation

velocity (volumetric flow rate) and the formation of the back-flow was checked.

ARTICLE IN PRESS

Table 1

Effect of cell size for the small tunnel

Case 1 Case 2

No. of grids 108  20  16 162  36  32

Cell size 0.05  0.02  0.0188 m 0.033 0.011  0.0094 m

uin (m/s) 0.620 0.620

Q 10.82 kW 10.73 kW

Layer length (m) 0:0oxo0:200 0:067oxo0:233

Highest temperature at the ceiling (1C) 139 at (0.350, 0, 0.30) 141 at (0.367, 0, 0.30)

Table 2

Effect of cell size for the large tunnel

Case 1 Case 2

No. of grids 240  24  24 300  30  30

Cell size 0.3750.3580.292 m 0.300  0.287 0.233 m

uin (m/s) 3.17 3.15

Q 16680 kW 17320 kW

Layer length (m) 6:38oxo9:00 11:1oxo9:90

Highest temperature at the ceiling (1C) 266 at (9.75, 0, 7 m) 256 at (8.70, 0, 7 m)

228 C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244
Runs were repeated until a range of ventilation velocities encompassed the

conditions exhibiting back-flow and no back-flow along the ceiling relative to the

upwind edge of the fire zone. Figs. 9 and 10 show the plots of ventilation velocity uin

versus the heat generation rate Q for the two tunnels used in the simulation. As

shown in Figs. 9 and 10, the demarcation line of the two regions separating back flow

and no back flow defines the critical ventilation values. The critical ventilation

velocity ucr can then be obtained from these plots as a function of Q: The uncertainty

in the values of ucr is 73% for the large tunnel and 71.5% for the small tunnel. The

uncertainty appears to arise from continuous oscillations in the flow conditions. As

mentioned previously, the ‘Smokeview’ package was employed to study the

approach of the steady flow conditions. The relevant flow variables in the last

15–20 s were then averaged for plotting.

The maximum critical velocity is less than approximately 0.9 m/s for the small

tunnel and less than approximately 3.5 m/s for the large tunnel. For both tunnels, ucr

varies approximately as Q1=3 for low values of Q; but ucr levels off for high values of

Q: In the small tunnel, methane and propane were used as the fuel. In the large

tunnel, propane and heptane were used. Figs. 9 and 10 show that fuel types have

negligible effects on the value of ucr: This may be attributed to the computer code

which employs the mixture fraction model for combustion, and to the fact that the

heat of combustion is approximately same for the hydrocarbon fuels used in the

ARTICLE IN PRESS

Q, kW

uin, m/s

100 101 102

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

methane, no back flow

methane, back flow

propane, no back flow

propane, back flow

critical values

Fig. 9. A diagram used for determination of critical velocities for the small fire tunnel. CFD results. The

power of the line at the lower end is 0.299.

C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244 229
simulation. This suggests that the main parameters affecting ucr are Q and a tunnel

dimension H (or H ̄ ) as given by Eq (1). Fig. 11 shows non-dimensional plots u

cr

versus Q using the values of ucr versus Q obtained from Figs. 9 and 10. As seen in

Fig. 11, the scaling parameters introduced to non-dimensionalize the plots collapsed

the data from two different size tunnels into curves sitting side by side. If the scaling

were perfect, the two curves would merge into one. It is noted that in Froude

modeling, the Reynolds number (viscous force) is ignored and the radiation and

conduction groups may not be preserved. At high heat generation rates, i.e., QX0:3;

the value of u

cr is approximately independent of Q: Other workers have observed

this trend [7,11]. Fig. 11 also shows experimental data from small-scale tunnels [11]

and large-scale tunnels (given in [11]) that are included in Fig. 2.

4.2. Froude modeling

As observed in Fig. 11, various experimental data and CFD results cannot be

collapsed into a single curve by applying the scaling parameter H ̄ : This indicates that

the Froude modeling is only an approximation for these data sets. This is expected

because there is no geometric and dynamic similarity among these data sets.

To assess the effects of tunnel size and fire-source geometry on the Froude

modeling, a series of CFD runs are performed. Tunnels I(a) and I(b) in Table 3 are

geometrically similar with a scale factor of 21.5. The width of tunnel I(a) is equal to

ARTICLE IN PRESS

Q, kW

uin, m/s

103 104 105

1

2

3

4

5

heptane, no back flow

heptane, back flow

propane, no back flow

propane, back flow

critical values

Fig. 10. A diagram used for determination of critical ventilation velocities for the large fire tunnel. CFD

results. The power of the line at the lower end is 0.321.

230 C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244
that of the small tunnel and the width of tunnel I(b) is equal to that of the large

tunnel. The tunnel length of I(b) is smaller than that of the large tunnel used in the

computation. The values of Q are adjusted to be approximately equal for Cases I(a)

ARTICLE IN PRESS

Table 3

The effects of tunnel size and fire-source geometry on the Froude modeling

Case H ̄ ðmÞ Q(kW) ucrðm=sÞ Q u

cr Fire source (m) Fuel Remarks

I (a) 0.343 11.1 0.605 0.145 0.330 0.15  0.4 CH4 a

(b) 7.37 24,000 3.41 0.147 0.401 3.23  8.6 CH4 a

(c) 0.343 11.6 0.66 0.152 0.360 0.5  0.12 CH4

II (a) 0.343 4.7 0.463 0.062 0.252 0.15  0.4 CH4 a

(b) 0.686 26.6 0.71 0.062 0.274 0.30  0.8 CH4 a

(c) 0.343 4.8 0.51 0.063 0.278 0.5  0.12 CH4

III (a) 7.72 15,100 3.04 0.082 0.349 9.92  2.6 C3H8

(b) 7.72 15,400 2.79 0.084 0.321 3.0  8.6 C3H8 a

IV (a) 7.72 44,800 3.60 0.244 0.414 10.8  2.6 C7H16

(b) 7.72 46,300 3.96 0.252 0.455 3.3  8.6 C7H16 a

V (a) 0.343 119 0.875 1.557 0.477 0.15  0.4 CH4 a

(b) 0.343 121 0.99 1.584 0.540 0.5  0.12 CH4

a

Fire source spans across the tunnel floor.

Q*

u*

10-2 10-1 100 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

large tunnel, CFD

small tunnel, CFD

Wu & Bakar, Tunnel A

Wu & Bakar, Tunnel D

Hwang & Wargo

Memorial Tunnel

cr

Fig. 11. Non-dimensional plots of the critical ventilation velocity versus the fire heat generation, showing

a comparison of CFD and experimental results.

C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244 231
and I(b). The value of u

cr for I(b) is 21% larger than that for I(a). This observation

alone may account for the separation of the lines for the CFD results in Fig. 11.

Cases I(a) and I(c) compare the effect of the fire-source geometry on u

cr: Case I(a)

has a fire source spanning across the cross-section, while Case I(c) has a rectangular

fire source along the axis of the tunnel. Both fire-source areas are equal. The

value of Q for Case I(c) is 5% greater than that of Case I(a), the value of u

cr is 9%

greater. Case II is similar to Case I with different value of Q and different ratio

of the tunnel size. The trends are similar. Case III and IV show the effects of

the fire-source geometry on the value of u

cr: The effects are opposite

depending on the value of Q: The fuel type is shown to have a negligible

effect on u

cr; the effect appears to be the fire-source geometry or other effect such as

radiation which is not considered in the Froude Modeling. In Case V, very high

values of Q are considered. The effect of fire-source geometry is similar to that

observed in Cases I and II.

In a geometrically similar tunnel, altering the fire-source geometry alters the

kinematics of flow. Dynamical similarity can no longer be preserved. Fig. 12 are

vector plots showing the velocity fields along a cross-sectional plane for Cases I(a)

and I(c) in Table 3. Fig. 12(a) has the fire source spanning across the floor,

showing upward flows along the entire floor. The upward flow is bent over at

about 1/3rd of the tunnel height by the ventilation flow. The ventilation-flow vector

is approximately perpendicular to the cross-sectional plane, and cannot be

represented in Fig. 12. Fig. 12(b) has the fire source running the center region of

the floor, producing a rising fire plume in the center region of the cross-section.

The vector plots in a cross-sectional plane at x 1⁄4 2:15 m for Case I(b)

(not shown) also show upward flow along the entire floor, but bends over at about

1/6th of the tunnel height. At the upper corners, strong streams toward the

center region are visible. There is no such stream in Fig. 12(a). Case I(a)

and I(b) are geometrically similar tunnels, yet the flow patterns are different at

the geometrically similar locations. These CFD simulations confirm that the

Froude-scaling law expressed by the parameters in Eq. (1) is an approximation for

tunnel fires.

4.3. Effects of tunnel inclination on u

cr

For the same fire-source conditions, the critical ventilation velocity is dependent

on channel inclination angles. Fig. 6 shows that when the tunnel inclination angle is

positive (f40), the ventilation flow is ascensional. Experiments show that a tunnel

with rising air-flow (f40; uphill slope) requires a smaller critical ventilation velocity

than that with horizontal air-flow [3]. Fig. 13 shows the results of the present CFD

and the experimental data of Hwang and Wargo [3] both with Q 1⁄4 4:5 kW: The

CFD results are in good agreement with the experimental results. Also shown are

CFD results for Q 1⁄4 35kW: It may be noted that the ucr  f relationshipis nonlinear.

Atkinson and Wu obtained experimentally expressions for the critical velocity in

tunnels that slope downhill with an angle between 01 and 101 [29]. The expressions

ARTICLE IN PRESS

232 C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244
show that the critical velocity ucr varies linearly with the inclination angle f:

ucr 1⁄4 ðgH ̄ Þ

1=2

u

max

Q

0:12 1=3

ð1  0:014fÞ for Qo0:12; (10)

ucr 1⁄4 ðgH ̄ Þ

1=2

u

maxð1  0:014fÞ for Q40:12 (11)

It is noted that f is negative for a downhill slope. (Note that Eq. (5) of their paper

has an error). Their results are also shown in Fig. 13. For Q

o0:12 (corresponding

to Q 1⁄4 4:5 kW), Eq. (10) checks well with the results of Hwang and Wargo [3]. For

ARTICLE IN PRESS

y, m

z, m

-0.2 -0.1 0 0.1 0.2

0

0.1

0.2

0.3

(b)

y, m

z, m

-0.2 -0.1 0 0.1 0.2

0

0.1

0.2

0.3

(a)

Fig. 12. Vector plots of the flow fields downstream of the fire source. (a) Case I(a) in Table 3 with

Q 1⁄4 11:1 kW and ucr 1⁄4 0:605 m=s: Tunnel cross-section at x 1⁄4 0:1 m (inside of fire source). (b) Case I(c)

with Q 1⁄4 11:6MW and ucr 1⁄4 0:66 m=s: Tunnel cross-section at x 1⁄4 0:1 m (inside of fire source).

C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244 233
Q

40:12 (corresponding to Q 1⁄4 35 kW), Eq. (11) under-predicts the results of the

present computation.

Kennedy et al. [6] have an expression for the effect of tunnel inclination that is a

part of the equations for computing ucr:

ucr 1⁄4 K1Kg

gHQ

rcpATf

 1=3

(12)

with

K1 1⁄4 Fr1=3 cr ;

Kg 1⁄4 1 þ 0:0374ðgradeÞ

0:8

Tf 1⁄4 Q

rcpAucr

þ Tin

where A is the tunnel cross-sectional area, Frcr is the critical Froude number, taken

as 4.5, ‘‘grade’’ is the absolute value of the tunnel downgrade in %, and Tin is

absolute temperature of approaching air. Eq. (12) is employed to predict ucr for

the small tunnel with downgrades from 01 to 101. The results are shown in Fig. 13.

ARTICLE IN PRESS

Tunnel Inclination Angle, φ, degree

ucr, m/s

-20 -10 0 10 20

0

0.2

0.4

0.6

0.8

1

1.2

1.4

Q=4.5 kW, Hwang & Wargo

Q=4.5 kW, CFD

Q=35 kW, CFD

Q=4.5 kW, Atkinson & Wu

Q=35 kW, Atkinson & Wu

Q=4.5kW, Kennedy et al.

Q=35kW, Kennedy et al.

downhill slope uphill slope

4.5 kW

35 kW

Fig. 13. Effects of channel inclination on the critical ventilation velocity. Positive inclination angle

corresponds to ascensional ventilation air-flow.

234 C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244
Eq. (12) slightly under-predicts the experimental results. The factor representing the

effect of tunnel inclination Kg appears to over-predict the effect.

4.4. Effects of ambient temperature on u

cr

A series of runs were made to see any effect of the ambient temperature on the

critical ventilation velocity. The ambient temperature was varied from 5 to 35 1C for

given values of Q and uin: Table 4 shows the results of this study using the small

tunnel. As the ambient temperature was varied, the ceiling-layer position and its

layer length stayed approximately unchanged. This implies that within the variation

of the ambient temperature studied, the ambient temperature has negligible effect on

the critical ventilation velocity.

4.5. Temperature maximum at critical ventilation

The trend that the critical ventilation velocity u

cr is approximately independent of

the fire heat generation Q beyond approximately Q 
 0:3 will be discussed. The

present CFD simulation found that as Q was increased maintaining the criticalventilation conditions, the temperature field at the ceiling region above the fire

source approached a maximum. Fig. 14 shows the temperature contours at y 1⁄4 0

plane for the small tunnel under the critical ventilation conditions. Referring to Fig.

9, the three values of Q’s selected represent the values at 1/3rd-power region of the

line, at the transitional region, and at the constant-ucr region. In Fig. 14(a), with a

low value of Q the temperature at the ceiling region is low at approximately 40 1C. In

Fig. 14(b), Q is increased to 53 kW, and the ceiling temperature is now 100–350 1C.

As Q is increased to 89 kW in Fig. 14(c), the ceiling temperature stays at

approximately 200–350 1C. In Fig. 15, the temperature contours at y 1⁄4 0 plane are

for the large tunnel. In reference to Fig. 10, the three values of Q selected represent

those at 1/3rd-power region of the line, at the transitional region, and at the

constant-ucr: The ceiling temperatures in Figs. 15(a)–(c) are approximately 100,

200–250, and 200–250 1C, respectively. It may be noted that, despite a difference in

the tunnel size and the heat release rate, the temperatures above the fire source as

shown in Fig. 14(b), (c) and 15(b), (c) do not exceed about 350 1C. At the critical

ventilation, the back layering terminates at x 1⁄4 0 and the fire plume tilts forward. In

ARTICLE IN PRESS

Table 4

Effect of ambient temperature on the critical ventilation velocity

Temperature (1C) Q(kW) uinðm=sÞ Layer length (m)

Min. x Max. x

5 4.53 0.50 0.08 0.32

15 4.54 0.50 0.04 0.28

25 4.52 0.50 0.08 0.28

35 4.52 0.50 0.04 0.32

C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244 235
a ventilated tunnel, the forward layer carries more heat energy than the back layer.

At low Q’s, to maintain the back layer at x 1⁄4 0; ucr increases as Q increases. As Q

increases close to a value for which ucr levels off, the fire plume bends to the extent

that the back layer receives the same amount of heat energy. Any additional heat

energy goes to the forward layer. At this point, ucr is unchanged and sustains the

back layer at x 1⁄4 0; even though Q keeps increasing. For Figs. 14(b), (c) and 15(b),

ARTICLE IN PRESS

9

1

x, m

z, m

-0.1 0 0.1 0.2 0.3 0.4 0.5 0.6

0

0.1

0.2

0.3

(a)

0

0 15

45

35

25

055

x, m

z, m

-0.1 0.1 0.2 0.4 0.5 0.6

0.1

0.2

0.3

0.3

(b)

2 5

5

50 5

x, m

z, m

-0.1 0.1 0.2 0.3 0.4 0.5 0.6

0

0

0.1

0.2

0.3

(c)

Fig. 14. Temperature contour at y 1⁄4 0 plane for the small tunnel, showing the conditions at the critical

ventilation. (a) 3 kW; (b) 53 kW; (c) 89 kW.

236 C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244
15(c) the maximum temperature is reached for critical ventilation. The attainment of

a temperature maximum at the critical ventilation provides a demarcation at which

ucr no longer changes as Q increases. It may be noted from Figs. 14 and 15 that the

forward-layer temperature may reach much higher level than the maximum

temperature.

ARTICLE IN PRESS

50

50

50

50

50

100

100

100

100

100

150

150

200

x, m

z, m

0 10 20 30 40

0

2

4

6

(a)

50

50

50

100

100

00

100

100

150

150

15

0

150

20

200

200

200

200

250

250

250

250

250

300

300

300

300

300

300

350

350

350

350

350

400

400

450

450

0 10 20 30 40

0

2

4

6

z, m

(b) x, m

50

50

50

50

100

100

100

150

150

150

150

200

200

200

200

15

250

250

300

300

300

300

300

3

350

350

530

4

400

400

400

450

450

450

450

500 550

0 10 20 30 40

0

2

4

6

z, m

(c) x, m

Fig. 15. Temperature contours at y 1⁄4 0 plane for the large tunnel, showing the conditions at the critical

ventilation. (a) Q 1⁄4 13MW; (b) Q 1⁄4 67MW; (c) Q 1⁄4 144MW:

C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244 237
4.6. Temperature stratification

An analysis of Eq. (9) requires temperature profiles downstream from the fire.

Because of a large variation in a given cross-section in the region downstream from

the fire source, the temperature profile at a cross-section is obtained by averaging in

the y-direction for a given value of z: The value of Tavg for a cross-section A is

computed by the equation

Tavg 1⁄4

Z ruT dA

ru dA

Fig. 16 shows the temperature profiles for the cases shown in Figs. 14 and 15. The

cases (a)B, (a)C (small tunnel) and (b)B, (b)C (large tunnel) correspond to the

leveling-off of u

cr in Fig. 11. As shown in Table 5, the estimated ucr agrees well with

the CFD value except for Case (a)A in which Eq (9) over-estimates the CFD value.

Except in Case (a)A, there is a range of x-values at which the temperature profiles

yield u

cr approximately equals to the CFD vlalue.

The parameter used by Newman [14] for temperature stratification DTcf=DTavg is

not consistent for Case (b)A. A more appropriate parameter appears to be the ratio

DTcf=Tavg: By comparing the values in Table 5 with the temperature profiles in

Fig. 16, the degree of temperature stratification is consistent with the value of

DTcf=Tavg: It appears that when this value reaches approximately 0.75 in critical

ventilation, the temperature stratification is complete, and leveling-off of u

cr takes

place. It is noted that the value of 0.75 applies to both the small and the large tunnels.

Based on the forgoing analysis, Eq. (9) is shown to predict the relationship of

u

cr  Q in Fig. 11. In summary, the following statement may be made for Eq. (9):

As Qðor Q

c Þ ! 0; Tin

Tavg

! 1 and DTavg

DTcf

! 0

0

; indeterminate:

But the right-hand side of Eq. (9) becomes u

cr as Q ! 0; so that DTavg=DTcf is a

large constant. Therefore,

O ! const: as Q ! 0;

and

Q ! large values; O / ðQÞ

1=3

:

The last statement is difficult to argue; however, the CFD values for u

cr in Table 5

appear to verify the statement.

4.7. Explanation using a free plume

Oka and Atkinson [7] tried to explain the observed trend using results for a free

plume [30]. This argument is presented in the following. For a very small fire where

ARTICLE IN PRESS

238 C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244
the flame length is a small fraction of the tunnel height,

Dr /

Q2=3

z5=3 ;

ARTICLE IN PRESS

C

T, °C

z, m

0 200 400 600 800

0

0.1

0.2

0.3

A B

(a)

T, °C

z, m

0 200 400 600

0

1

2

3

4

5

6

7

A B C

(b)

Fig. 16. Temperature profiles downstream of the fire source. (a) Small tunnel: Profile A) Fig. 14(a),

x 1⁄4 1:85 m; Profile B) Fig. 14(b), x 1⁄4 0:61 m; Profile C) Fig. 14(c), x 1⁄4 1:05 m: (b) Large

tunnel: Profile A )Fig. 15(a), x 1⁄4 10 m; Profile B) Fig. 15(b), x 1⁄4 37 m; Profile C) Fig. 15(c),

x 1⁄4 51 m:

C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244 239
where Dr is the density difference above ambient. The buoyancy head Dp can then be

estimated as

Dpbuoyancy 1⁄4

Z H

0

Drg dz /

Q2=3

H2=3 (13)

For a very large fire with continuous flaming to roof level, the temperature is roughly

constant above the fire and the density difference becomes a constant value Drf for

all z and Q: Figs. 14 and 15 show that this statement may not be true. The buoyancy

head becomes

Dpbuoyancy 1⁄4

Z H

0

Drfg dz / H (14)

Oka and Atkinson [7] argued that, at the critical ventilation velocity, the ratio of the

dynamic head of the inflow to the buoyancy head above the fire must exceed a

critical value, or

Dpdynamic

Dpbuoyancy " #

cr

4Ccr (15)

For small fires this becomes (using Eq. (13) and Dpdynamic 1⁄4 rinðu2

cr=2)

ucr4C1=2 cr

Q

H

 1=3

(16)

This implied that ucr is proportional to Q1=3 for a small Q: For large fires we get,

using Eq. (14),

ucr4C1=2 cr H1=2 (17)

Thus, ucr is independent of Q for a large Q:

4.8. Comparison with simple theories

Kennedy et al. [6] developed a theory for the critical ventilation velocity in a fire

tunnel. The theory is based on an overall heat balance in one-dimensional flow of a

fire tunnel, and is given by Eq. (12). It is noted that this set of equations in Q and ucr

ARTICLE IN PRESS

Table 5

Estimation of u

cr using Eq. (9)

Q (kW) T-profile Q u

cr Eq. (9) ucr CFD DTcf=DTavg DTcf=Tavg H ̄ ðmÞ

3 (a)A 0.039 0.392 0.268 0.928 0.078 0.343

53 (a)B 0.693 0.503 0.480 1.708 0.836 0.343

89 (a)C 1.170 0.506 0.482 1.146 0.776 0.343

13500 (b)A 0.073 0.310 0.325 3.47 0.400 7.07

67200 (b)B 0.366 0.444 0.406 1.906 0.745 7.07

143800 (b)C 0.783 0.475 0.406 1.257 0.771 7.07

240 C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244
are non-linear. A FORTRAN program was written to obtain a set of Q and ucr for

the large tunnel.

By using a two-dimensional approach, Kunsch developed a model for the critical

ventilation velocity in a ventilated tunnel [12]. The critical ventilation velocity is

given by

u

cr 1⁄4 C3

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

C1DT

0

p

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

1 þ ð1  C2=C1ÞDT

0Q2=3

q

1 þ DT

0Q2=3 Q1=3 (18)

with

C1 1⁄4 1  0:1ðH=WÞ

1 þ 0:1ðH=WÞ 1 þ 0:10 HW

  0:015 H

W

2 " # ffi 1  0:10 H

W

;

C2 1⁄4 1  0:10 H=W 	 


1 þ 0:10 H=W 	 
 0:574 1  0:20 H

W

;

C3 1⁄4 0:613;

where DT

0 is a dimensionless constant with a value equal to 6.13, and W is the

tunnel width in meter.

ARTICLE IN PRESS

Q, kW

ucr, m/s

103 104 105 106

1

1.5

2

2.5

3

3.5

4

4.5

5

CFD, this study

Kennedy et al. (6)

Kunsch (12)

Memorial tunnel data

1/3 - power line

Fig. 17. Comparison of two simple theories with the results of the present CFD simulations for the large

fire tunnel.

C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244 241
Fig. 17 shows the plots of Eqs. (12) and (18) with data that applied to the large tunnel.

The CFD results are also shown for comparison. The lines representing these simple

theories and the CFD results cross at Q 
 105 kW: For Q below 105 kW, the CFD

simulation predict the values of ucr approximately the average of those predicted by

these two theories. Although the simple theories predict leveling-off of u

cr as Q

increases, CFD results give better agreement with the experimental trend. It is noted that

the simple theories are based on analyses of one-dimensional or two-dimensional flows

while the CFD results are based on three-dimensional tunnel flow. A simple theory does

not include the complex flow patterns associated with the rising hot plume. Neither does

it include the fire-source geometry or gas radiation. The present CFD results provide

details of 3-D flow patterns with reasonable computer time.

5. Discussions and conclusions

In the present investigation, a CFD code FDS2 was employed to predict the

critical ventilation velocity in fire tunnels. Tunnels of different sizes and fire-source

geometries were selected for simulations. The following observations are made from

the present study.

(1) When the critical ventilation velocity uin is plotted against the fire heat generation

rate Q as shown in Fig. 2, ucr is roughly proportional to the 1/5 power of Q: This

plot encompasses all available data of tunnel sizes. It is noted that there are

subgroups of points for tunnels of different size. Each subgroup has its maximum

ucr; the critical velocity at leveling-off. For each tunnel size, if we take the

points in Fig. 2 where ucr begins to level off and connect these points with a

straight line, the line has approximately 1/5th power. For example, for the

Memorial tunnel the point is (3  104 MW; 2.9 m/s), and for the tunnel A, Wu

and Bakar the point is (5 kW, 0.48 m/s). The straight line connecting these two

points in Fig. 2 has approximately 1/5th power. According to Quintiere [31], the

reference velocity V for a fire plume is

V 1⁄4 ffiffiffig p Q

r0

ffiffiffi

g p cpT0

1=5

It can be argued that the fire-plume velocity is the base velocity of the back-layer

velocity when the plume reaches the ceiling region. The back-layer velocity

eventually becomes zero (or turns around) at x 1⁄4 0 for critical ventilation. This

process involves the ventilation velocity in contact with the back-layer. Therefore,

ðucrÞmax / V / Q1=5

where ðucrÞmax is the value of ucr at leveling-off. In Fig. 3, the subgroupslides along

the 1/5th power line when the scaling factor H ̄ is applied through the Froude

Modeling. Within the subgroup, ucr / Q1=3 for small Q’s.

(2) The CFD results clearly show the leveling-off of the critical ventilation velocity

as the fire heat generation increases. Previous workers have established similar

ARTICLE IN PRESS

242 C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244
trend experimentally. The leveling-off of the critical ventilation velocity can be

explained by the temperature maximum above the fire source.

(3) The concept of temperature stratification in fire tunnel [14] is applied to the

temperature profiles obtained by the present CFD simulations for both the small

and the large tunnels. A new parameter for temperature stratification appears to be

more consistent for the present application than the one defined by Newman [14].

(4) The CFD results show that the fuel type and the ambient temperature have

negligible effects on the value of ucr:

(5) The tunnel configuration (size and fire-source configuration) appears to have

effects on the u

cr  Q relationship. It may be noted that the Froude modeling

used in the correlation requires geometric similarity. When such similarity is not

maintained as in our small channel and large tunnels, a perfect Froude modeling

cannot be expected. The results of the present simulations show that maintaining

the geometrical similarity for tunnels cannot maintain the kinematical similarity.

This indicates that the Froude-scaling law is an approximation for tunnel fires.

(6) In the simple theories using one- or two-dimensional approach, the leveling-off

of the critical ventilation is also predicted, though not as sharply as the results

shown in experiments and CFD. The parameter, layer width W; appears in the

two-dimensional theory, and does not appear in the Froude modeling [12]. The

introduction of the hydraulic tunnel height H ̄ effectively absorbs the parameter

W [11]. As long as the channel cross-section does not deviate too far from a

square, approximate geometrical similarity is maintained. Since the aspect ratio

of the tunnels used in the practical application are not too large, an extrapolation

of small-scale results is approximately valid.

(7) The FDS code using LES turbulence model [20] is capable of predicting the

critical ventilation velocity in channels of various size and configuration. The

effect of tunnel inclination on the critical ventilation velocity is also predicted.

The results compare well with existing experimental data. Comparisons with the

results of previous authors are also made. The agreement with these results and

the present CFD results is fair. As expected, in an inclined tunnel with

ascensional ventilation, a smaller ventilation velocity is required to prevent back

layering of the hot gas than in a horizontal tunnel. The change in the ventilation

velocity is a non-linear function of the inclined angle of the tunnel.

References

[1] Thomas PH. Movement of smoke in horizontal corridors against an air flow. Inst Fire Engrs Q

1970;30:45–53.

[2] Massachusetts Highway Department. Boston, MA: Memorial tunnel fire ventilation test program,

Test Report, 1995.

[3] Hwang CC, Wargo JD. Experimental study of thermally generated reverse stratified layers in a fire

tunnel. Combust Flame 1986;66:171–80.

[4] Heselden AJM. Studies of fire and smoke behavior relevant to tunnels. In: Proceedings of the second

international symposium on the aerodynamics and ventilation of vehicle tunnels. Cambridge, UK:

BHRA Fluid Engineering; 1976, pp. 23–5.

ARTICLE IN PRESS

C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244 243
[5] Guelzim A, Souil JM, Vantelon JP, Sou DK, Gabay D, Dallest D. Modelling of a reverse layer of fireinduced smoke in a tunnel. Fire safety science proceeding of the fourth international symposium;

1994, pp. 277–88.

[6] Kennedy WD, Gonzales JA, Sanchez JG. Derivation and application of the SES critical velocity

equations. ASHRAE Trans: Research 1996;102(2):40–4.

[7] Oka Y, Atkinson GT. Control of smoke flow in tunnel fires. Fire Safety J 1996;25:305–22.

[8] Hwang CC, Chaiken RF, Singer JM, Chi DNH. Reverse stratified flow in duct fires: a twodimensional approach. In: 16th symposium international on combustion. Pittsburgh: The

Combustion Institute; 1977, pp. 1385–95.

[9] Charters DA, Gray WA, McIntosh AC. A computer model to assess fire hazards in tunnels (FASIT).

Fire Technol 1994;30:134–54.

[10] Grant GB, Jagger SF, Lea GJ. Fires in tunnels. Philos Trans R Soc London A 1998;356:2873–906.

[11] Wu Y, Bakar MZA. Control of smoke flow in tunnel fires using longitudinal ventilation systems—a

study of the critical velocity. Fire Safety J 2000;35:363–90.

[12] Kunsch JP. Simple model for control of fire gases in a ventilated tunnel. Fire Safety J 2002;37:67–81.

[13] McCaffrey BJ, Quintiere JG. Buoyancy driven countercurrent flows generated by a fire source. In:

Spalding DB, Afgan N, editors. Heat Transfer and Turbulent Buoyant Convection. Washington,

USA: Hemisphere Publishing co.; 1977. p. 457–72.

[14] Newman JS. Experimental evaluation of fire-induced stratification. Combust Flame 1984;57:33–9.

[15] Hwang CC, Edwards JC. CFD Modeling of smoke reversal. In: Proceedings of the international

conference on engineering fire protection design. Bethesda, Maryland, USA: Society of Fire

Protection Engineers; 2001, pp. 376–87.

[16] McGrattan KB, Baum HR, Rehm RG, Hamins A, Forney GP, Floyd JE, Hostikka S. Fire dynamics

simulator (version 2)—technical reference guide. NIST report no NISTIR 6783. Gaithersburg,

Maryland: National Institute of Standards and Technology; 2001.

[17] Xue H, Hihara E, Saito T. Turbulence model of fire-induced air flow in a ventilated tunnel. Int.

J. Heat Mass Transfer 1993;36:1739–48.

[18] Sinai YL, Owens MP. Validation of CFD modeling of unconfined pool fire with cross-wind: flame

geometry. Fire Safety J 1995;24:1–34.

[19] Woodburn PJ, Britter RE. CFD simulations of a tunnel fire—Part I. Fire Safety J 1996;26:35–62.

[20] McGratten KB, Forney GP, Floyd JE, Hostikka S. Fire dynamics simulator –—user’s guide (version

2). NIST report no NISTIR 6784. Gaithersburg, Maryland: National Institute of Standards and

Technology, 2001.

[21] Givi P. Model-free simulations of turbulent reactive flows. Prog. Energy Combut. Sci. 1989;15:1–107.

[22] Smagorinsky J. General circulation experiments with the primitive equations. I. the basic experiment.

Monthly Weather Rev 1963;91:99–164.

[23] Poinsot T, Veynante D. Theoretical and numerical combustion. Philadelphia, PA, USA: Edwards

Inc.; 2001.

[24] Baum HR, McGrattan KB, Rehm RG. Simulation of smoke plumes from large pool fires. In:

Proceedings of the 25th symposium (international) on combustion. Pittsburgh: The Combustion

Institute; 1994, pp. 1463–69.

[25] Baum HR, Ezekoye OA, McGrattan KB, Rehm RG. Mathematical modeling and computer

simulation of fire phenomena. Theor Comput Fluid Dynam 1994;6:125–39.

[26] Baum HR, McGrattan KB, Rehm RG. Three dimensional simulations of fire plume dynamics. J Heat

Transfer Soc of Japan 1997;35:452–5.

[27] McGrattan KB, Baum HR, Rehm RG. Large eddy simulations of smoke movement. Fire Safety

J 1998;30:161–78.

[28] Emmerich S, McGrattan KB. Application of large eddy simulation model to study room air flow.

ASHRAE Trans 1998;104(1B):1128–40.

[29] Atkinson GT, Wu Y. Smoke control in sloping tunnels. Fire Safety J 1996;27:335–41.

[30] MaCaffery BJ. Purely buoyant diffusion flames: some experimental results. NBSIR 79-1910.

Gaithersburg, Maryland, USA: National Bureau of Standards; 1979.

[31] Quintiere JG. Scaling applications in fire research. Fire Safety J 1989;15:3–29.

ARTICLE IN PRESS

244 C.C. Hwang, J.C. Edwards / Fire Safety Journal 40 (2005) 213–244