- Email: [email protected]
Numerical simulation of thermal plumes
Fire Safety Journal 21 (1993) 231-256
~ , ~
Numerical Simulation of Thermal Plumes
Soonil Nam & Robert G. Bill, Jr Factory Mutual Research Corporation, 1151 Boston-ProvidenceTurnpike, Norwood, Massachusetts 02062, USA (Received 22 April 1992; revised version received 28 August 1992; accepted 7 January 1993)
ABSTRACT Thermal plumes generated by pool fires and heptane spray fires were simulated numerically with the aid of the computer code PHOENICS. The flow fields associated with five cases were simulated and compared with experimental results. Modifications were made in the standard k - e model used in the P H O E N I C S code in order to improve the accuracy of that model in the simulation of turbulent thermal plumes. The modification includes an additional term related to the buoyancy generated turbulence kinetic energy, changing constants responsible for the turbulent viscosity and the effective Prandtl number in the energy equation. The results obtained with the modified code were in good agreement with experimental results and significantly improved the results obtained with the standard k-E model.
1 INTRODUCTION This study is the first part of an ongoing research program devoted to numerical modeling of the interaction of fire and sprinkler spray. The long range goal of the program is to develop computer models to predict the effectiveness of sprinkler systems, and reduce the need for costly experiments. One of the first difficulties encountered in simulating fire plumes is selecting an approriate and convenient turbulence model. It is well 231 Fire Safety Journal 0379-7112/93/$06-00 (~) 1993 Elsevier Science Publishers Ltd, England. Printed in Northern Ireland
232
Soonil Nam, Robert G. Bill, Jr
known that the k - e model tends to overpredict the velocities and temperatures along the centerline axis of buoyant plumes and consequently underpredicts the plume width.' 3 The k-E model has been improved for strongly buoyant flows through the addition of an algebraic stress model (ASM). 1~ The addition of an ASM provides a more accurate model for velocity and temperature distributions in a plume. The inclusion of an ASM, however, requires significantly increased computing time. Since the main purpose for our work is the development of a turbulence model which will be adopted for complex, and thus computation intensive, spray-plume interactions, we have chosen instead to modify the standard k-E model by adjusting its constants. This approach seemed reasonable since other studies 5"6 have suggested that the k - e model produces coarse agreement with measurements from fire and smoke flows. Therefore the major focus in this work is to find the optimal combination of some coefficients which are at our disposal in K-e model. The pool fire plume was selected for simulation because it has been extensively studied by many researchers, Morton et al., 7 Heskestad, ~ 1~ and Delichatsios; 12 thus data are available to validate the numerical results. In this study no combustion model was used. The domain of interest for comparison with experimental results, therefore, was above the flame region where the flow behaves as a pure thermal plume. The P H O E N I C S program, ~3 a commercial finite-difference computational fluid dynamics code, was used in this study. It employes the standard k-E model; however the model does not include a buoyancyinduced generation term in its turbulence kinetic energy generation formulation. Thus in addition to optimization of the k-E model constants, another term must be added to the turbulent kinetic energy balance. The modified code has been used to numerically simulate the axi-symmetric plumes generated by large diameter pool fires as reported by Kung and Stavrianidis. 14~5 A schematic of the flow to be simulated is shown in Fig. 1. The results of this study are also compared with plume correlations derived from many experiments. After constructing a turbulence model based on pool fires, the model was further validated with two sets of heptane-spray fire plumes.
2 MATHEMATICAL FORMULATION The plume arising from the pool fire is assumed to be axisymmetric, without swirl. Hence the conservation equations in cylindrical coordinates (with density fluctuation correlations ignored) are:
Numerical simulation of thermal plumes
233
\ \ \ \ . - ---iP-
\ \ \ \ \ \
j \
Vr
=
/ \
\
/
I I
\ \ \
/
Vr
/ /
I
/
\
Air Entrainment
\
\
i,
_71
r
Fig. 1.
The schematic diagram of a plume generated by a pool fire.
(1) Mass continuity:
at
-
(pozr)
r
-4- - - ( p O r t " )
Or
(1)
= 0
where vz and or denote z directional and r directional velocities. (2) z-momentum:
OpVz - -Ot+ r
(PvzrVz)+or((PVr - - ~) -
~z ~ r tZe ~ z J - Or r ~e Or/J
•
(2) _ O (.eO~)+!L S"~=Oz
Oz
f
O~,,
r Or krt~e~z ) -q- g(P=-- p)
(3)
where /z~ represents an effective viscosity which is the sum of a turbulent viscosity/z, and the laminar viscosity /zt, and p~ denotes the ambient density. (3) r-momentum:
apVr 1 1 0 Ot
~---f L ~Z (pllz
rv
O rv r) ~- ~ ( P ~ J r r) .
O {
Ovr~
-~z .~rtZe~z . }
.
O (rtZeOO___rr)] sot or /a
Or \
(4)
Soonil Nam, Robert G. Bill, Jr
234
Sot
OP Or
2/J, el) r -
-
r
+
O( OVz) l O { OVr~ OZ \txe Or / r or krtze--~r ) "
- -
+
(5)
(4) Energy: l{O
Oph + _ (pvzrh) + 0~- r -~z
;
O ( Ix~Oh] O (rlXeOh)] 0" (pv,rh) - Oz r - -= ~rh Oz : Or, oh -O~r
(6)
where h, Crh, q" denote enthalpy, effective Prandtl number for the diffusion of heat, and heat generation rate per unit volume, respectively. The effective Prandtl number in the energy equation, oh, is one of the two parameters which was varied to improve agreement with experimental results. (5) Turbulence kinetic energy k:
Opk~_1[0 Ot
r Oz (pvzrk) +
;
a (rlXtOk) a (rlX,k~l =Sk
(7)
(pv~rk) - -~z \ ~rk ~Z -- ~r \ ~rk Or/1
where ok is the effective Prandtl n u m b e r for the diffusion of turbulence energy and assumed equal to 1.0, and the turbulence production term Sk is composed of Sk = Gk + GB --
pE
(8)
where e, Gk, GB denote turbulence dissipation, shear production, and buoyancy production, respectively, and are defined as
f rl/Oldz\2 (OUrl2 Gk--lx'12[~--~z) + k O r /
(~---~)21 ( ~Uz
+
OVrl2~
+\Or+Oz/
10p
J
(9)
(lO)
GB -= p , , g - -
paz
Note that GB is not included in the standard k-E model and has been included in this study as in Ref. 6. (6) Turbulence dissipation E:
Ot
r Oz(PVzre)+~z(pVrre)-oz
\ro,Oz/-or
\ro,Or/J
where tr, is the effective Prandtl n u m b e r for the diffusion of turbulence dissipation and given as 1.314, and the production of dissipation term
Numerical simulation of thermal plumes
235
S, is defined as
S.=-k [C,(Gk + G.) - CEpE]
(12)
where C1 and C 2 a r e constants given as 1-44 and 1.92, respectively, in the standard model. The k - e model is adopted to obtain the effective turbulent viscosity, which is related to k and E by following formulation
ix, = OCDk2/E
(13)
Here Co is the other constant which was varied in the present study to improve the agreement with experimental results. All the other constants came from the work of Launder and Spalding. 16
3 M O D I F I C A T I O N OF T H E k-E M O D E L The primary task of the simulation is to devise an adequate turbulence model which produces accurate plume velocity and temperature distributions. Among many parameters in the k - e model, two parameters, CD and oh, most directly affect the diffusion of velocity and temperature in the flow. The standard k-E model assigns CD as 0.09, trh as 1.0. A test case with CD ----0.109, o-h----0-614 is reported in Ref. 1. It shows somewhat improved results for a pure jet. The spread rates of buoyant plumes were still significantly underpredicted; the modest improvements observed compared with the standard model was attributed to the inclusion of the buoyancy induced kinetic energy production term, rather than changing of the constants CD and ~rh. Since the standard model always predicts narrower velocity and temperature plume widths than measurements indicate, CD and o'h were varied to broaden the plume widths. It is obvious that a higher diffusion coefficient reduces peak values while increasing the effective region of transport of the property. Therefore we increased the turbulent viscosity by increasing the constant CD until reasonable agreement was obtained between the simulated velocity profile and measurements. The effective Prandtl number, ~rh, strongly affects the temperature distribution, but has little effect on the velocity distribution. A lower trh tends to increase the temperature width, thereby giving lower temperatures. After we obtained a reasonable range of CD based on the velocity distribution with a fixed trh, then both CD and trh were varied until satisfactory results for both temperature and velocity were obtained. The standards of the comparison were the well established plume
Soonil Nam, Robert G. Bill, Jr
236
relations ~°
UZo=Cv%(cpg )l'3qlc/3(z- Zo) 1/3 //
ATo = C To[
T.
~ l/3
"~
! qc-:/3(z _ Zo)-5/3
(14) (15)
\gt~pp~/ C
/ T . \1/2
o
(Z-Zo
(16)
where Vzo is the axial velocity, Vz, at the centerline, Oc is the convective heat relase rate, To is the absolute temperature at the centerline, ATo is the excess temperature above ambient at the centerline, bar is the plume half-width based on the centerline temperature, which is the width to the point where the temperature difference from the ambient becomes one half of the value at the centerline, and where Cv~,, Car,,, and Cb~r are demensionless constants. The constants Co~,., Car,,, C b j for the standards are chosen as 3.4, 9.1 and 0-12 following the recommendation of Heskestad. 1° After numerous trials with different combination of CD and o-h, CD=0"18, O'h=0"85 were chosen as the best combination which provided close agreement between the values of C~o, Car°, C~T and those in Ref. 10. As noted in Section 2, P H O E N I C S does not have a buoyancy produced turbulence kinetic energy term. In addition, P H O E N I C S assigns CD as 0.09 and o-h as 1.0. Therefore proper modifications were made to reflect the changes mentioned above. 4 TEST CONDITIONS The experiments to be simulated were pool fires with large enough diameters to generate turbulent plumes. The general conditions of the experiments reported in Refs 14 and 15 are given in Table 1. TABLE 1 Test Variables in the Experiments 14'17
Test no 1 4 6 8 HS1 HS2
Fuel
Methanol Silicone F1. Hydrocarbon F1. Heptane Heptane Heptane
Pool diam. D (m) 1.219 1.737 1-219 1.219 0.688 0.688
Convective heat flux (lc(Kw)
(l"(kW/m 2)
410 170 792 2 240 163.5 295 -4
351 72 679 1 920 440-1 795.1
Flame height (m) 0.9 0.2 2.2 3-9 ---
Numerical simulation o f thermal plumes
237
Tests 1, 4, 6 and 8 in Refs 14 and 15 were numerically simulated. The numerical simulations reported here are designated by the corresponding test n u m b e r in those references. Since the pure thermal plume behavior should not d e p e n d on any specific character of the fuels but only on the convective heat release rate, simulation of these tests which represent different fuels provides a good test for the validity of a numerical code. In addition to these simulations, two more cases with heptane spray fires ~v were tested for further validation of the model. The plumes were generated by the heptane spray fires at two convective heat release rates of 164kW and 295 kW, which are denoted as HS1 and HS2, respectively. Four heptane spray nozzles were placed at the corners of a 0.61 m square and directed radially toward the center. The conditions of these experiments are also given in Table 1. The flame heights of HS1 and HS2 were not measured.
5 COMPUTATIONAL PROCEDURES 5.1 Numerical domain for the tests
The numerical grids used for the simulations are given in Table 2 and the coordinate system is shown in Fig. 1. Uniform grids were adopted along the z direction, but nonuniform grids were used along the r direction. To place fine grids around r = 0, and coarser grids as r becomes larger, the following formula was adopted: (17)
r ( i ) = r o ( i / n ) ~5
TABLE 2 Numerical Domains for the Tests Test no.
1 4 6 8 HS1 HS2
Pool diam. (m)
1-219 1.737 1.219 1.219 0.688 0-688
Length o f domain (m)
Number o f meshes
r direction
z direction
r direction
z direction
12-52 17.84 12-52 12-52 7.24 7.24
15.24 15-0 30.48 30.84 7.27 7.27
30 60 30 30 30 30
51 51 51 51 51 51
Soonil Nam, Robert G. Bill, Jr
238
Numerical grid distribution
IIIIIIIF Ilfllll I
ililill i ililiil I Ililill I
liltllll Illli i t
.2 []
iliilll iliilll Ililll I Ilillll lillii lilli I Ilili I Illli i IIIli I Illli I Ilrll I iliili Ilill i I ililfl lillll I lillil I I illlii ilillr I ililli I I iliil i liliil f lillil I
H
Radius Fig.
2.
Numerical
grid
distribution.
where r(i) is the radial distance to the mesh i, ro is the maximum length of the radial numerical domain, and n is the total number of meshes in the r direction. The numerical grids used are shown in Fig. 2. In order to avoid a possible numerical instability caused by inappropriate initial guesses for velocities and temperatures with large heat release, an unsteady code was used. To ensure convergence, a small time stip, 0-5 s, was used for the initial stage and then larger time steps, 2.5 s, were employed after 100 time steps. The initial conditions were ambient except for the constant heat flux at the surface of the pool.
5.2 Boundary conditions A convenient way of providing the desired heat release is to assume that the heat release rate is uniform over the area of the pool. This avoids the complex flame shape required for a simulation using a volumetric heat release rate. Since we are only interested in the plume region of the flow, the accuracy of the simulation directly downstream of the area where heat is released is not important to this study.
Numerical simulation of thermal plumes
239
The two sets of heptane spray fires were treated as pool fires of diameter of 0.688 m, which corresponds to the hydraulic diameter of the square on which the nozzles were placed. Boundary conditions for the surfaces r =ro and z = 0 were that pressures are equal to the ambient, and the transport of m o m e n t u m and thermal energy are only allowed in the incoming direction. At each iteration, the direction of moving flow on the surface r =r0 was checked, and if outflow occurred, then the velocity at that point was set to zero until the computation predicted a value with the proper direction at the next iteration. The boundary conditions on the surface z = zoo, the end of the numerical domain in the z direction, were similar to those at r = ro; however, only outgoing instead of incoming flows were allowed. The boundary conditions at r - - 0 were the usual symmetry conditions: O/Or = 0 and v, = O.
6 RESULTS A N D DISCUSSION The variables p, h, Vz, vr, k and • were computed at each node for all time steps. The result for each test was compared with the result computed with approximately twice the number of z- and r-directional grid points reported in Table 2, in order to check the reliability of the numerical solutions. The temperature curve along the radial direction at z = 2zl (explanation is given in the following paragraph) was checked for consistency. All cases but Test 4 showed less than 2% of discrepancy. The heat release rate in Test 4 was at least an order of magnitude smaller than the heat release rates in the other simulated tests. Such a low heat release rate generated very low flows which were not well represented by the selected grid, and there was a practical limitation in number of grid points in the code. For this reason, Test 4 was not included in the analyses of the modified k - • model. Computation of each test took approximately 1 h in SUN Work Station Sparc IPX. The gas enthalpy h was transformed into the corresponding temperature for a comparison with measurements. The half-width based on the centerline velocity is the width to the point where the axial velocity decreases to one half of that of the centerline. The half-width for temperature is based on the excess temperature above ambient. Figures 3 and 4 show the velocity and temperature half-widths. The velocity half-widths are always a little wider than those of temperature. This is consistent with a review of literature data by Heskestad. 1° The plume
Plume -r
half-width
based
I
on velocity
I
I
/
/ ..-'
/ / ..'
/..
,-,
2 I , / . ".
I 1
-
...~" 1..-..~-S '~
0 0 Fig. 3.
5
10 Elevation
16
20
(m)
P l u m e h a l f - w i d t h s b a s e d on the centerline velocities ( - - , Test 1; . . . , - - Test, 8; - . - H S 1 ; - . . . . H S 2 ) . Plume
half-width I
based
on
I
Test 6;
temperature '
'
I
. ..'" /:// ./S 2
/
-
E
..5.'.-¢'" ~.Is
I
t~
1
0 0
5
10 Elevation
Fig. 4.
15
20
(m)
P l u m e half-widths based on the centerline temperatures ( - - , Test 1; • - -, Test 6; , Test 8; - . - , H S I ; . . . . , H S 2 ) .
Numerical simulation o f thermal plumes
241
TABLE 3 P l u m e R e g i o n s O b t a i n e d f r o m Numerical Simulation
Test no. Max. distance in z (m) 1 6 8 HS1 HS2
15-24 30-48 30.48 7.27 7.27
Plume region (m)
zt (m)
Virtual origin: zo (m)
2.2 --
2.2 3.0 4.0 1.2 1.2
-0.2 -0.4 -0"6 -0.2 -0.2
regions of the numerical solution were identified by the elevations where the plume half-widths first exhibited a nearly linear increase with elevation (see eqn 16); their values, z~, are shown in Table 3. In order to compare the numerical results with the plume relations from eqns (14)-(16), the distance to the virtual origin Zo must be found for each test. Equations (14) and (15) are modified for that purpose as follows:
z=C zo;?( g-- ]q<+z0 o o kCpp~TJ | ~ / Z = /C3/5 / baT/---'~-'~_2! L
\gt~pp~l
~.e/5/ATo315 t/c / -I- Zo
j
(18)
(19)
Applying a least-square fitting provides the best fitting constants for Co,o, Zo in eqn (18), and Char, Zo in eqn (19). The values of Zo obtained from the two equations were very close to each other for all tests. The values of Zo given in Table 3 are obtained from eqn (19), based on the centerline temperature, and are used in all the following analyses. After inserting Zo into plume relations from (14) to (16), the coefficients Cozo, Car°, Cb~ are plotted in Figs 5, 6 and 7. All coefficients show quite good agreement with the values quoted previously from Ref. 10: 3.4, 9.1 and 0.12. In order to demonstrate the effect of modified k-E model constants, the results obtained with the standard k - e model for Test 8 are shown in Figs 5-7. The standard k - e model contains the terms which take into account the kinetic energy production due to the buoyancy. The results clearly demonstrate that the modified version predicts the plume behavior far better than the standard k-E model. The standard k - e model overpredicts centerline excess temperature and velocity in Test 8 by 68% and 25%, respectively, which results in a much narrower temperature width for a given plume. Other tests with the standard k - e model also produced similar results, however, they were not included
5
'
1
'
I
I
2
3
Nondimensional
4
(z/zl)
Elevation
Fig. 5. C o m p a r i s o n of the p l u m e coefficient C,~ ( - - , c o n s t a n t in Ref. 10; • • -, Test l; - - , Test 6; - - - , Test 8; . . . . . , HS1; - - , HS2; o . . _ _ Test 8 with the s t a n d a r d k - e model). 20
~
,
12 :'Z.'.:.;.-= -.': . . . . . . .. . . . . . . . . . . .
1
2 Nondimensional
Fig.
--,
Elevation
3 (z/zl)
4
C o m p a r i s o n of the p l u m e coefficient CAr° ( - - , c o n s t a n t in Ref. 10; • •., Test 1; Test 6; - - - , Test 8; . . . . . , HS1; - - , HS2; . . . . , Test 8 with the s t a n d a r d k - e model).
6.
Numerical simulation of thermal plumes
0.2
i
i
I
I
243
0.1
0.0 1
2 3 N o n d i m e n s i o n a l Elevation (z/zz)
4
Fig. 7. Comparison of the plume coefficient Char.
in the figures to avoid being overcrowded. The simulations were compared with the experiments reported in Refs 14 and 15. As stated in Ref. 14, the heat release rates, qc, of the tests were c o m p u t e d based on the measurements at z = 7.62 m. These heat release rate values were used to generate the buoyant plumes in the present numerical simulations. The measured temperatures 14 at z = 7 . 6 2 m and the numerical solutions are shown in Fig. 8 and are in reasonable agreement. The experimental data for Test 8 was not available. Additional comparisons are provided in Figs 9-12 in simulations of the plumes of heptane-spray fires. The overall comparison between numerical results and experimental measurements is good. Figure 9 shows the comparison of temperature and Fig. 10 shows the comparison of velocity with measurements for HS1. Figures 11 and 12 show the same comparisons, respectively, for HS2. Though the earlier comparison based on the plume coefficients is more general and more reliable, the similar comparisons as shown in Figs 9-12 are provided in Figs 13-16 to help to visualize the difference between the modified and the standard model. The numerical results shown in Figs 13-16 are obtained with the standard k-E model. Note the higher plume
P l u m e t e m p e r a t u r e at z = 7 . 6 2 m 1
400
i
J
i
I
i
,
i
i
i
i
i
I
i
i
i
~
]
J
i
i
,
380
O
v
360
t_ c~
340
-
z~
320
".. . ,~
A
O
"'""....
3O0
0.0
0.5
1.0
1.5 Radius
2.0
2.5
(m)
Fig. 8. Comparison of the numerical and experimental results for temperature in pool fires (numerical: - - , Test 1; • • ,, Test 6; experimental: A, Test 1; ~ , Test 6).
T e m p e r a t u r e d i s t r i b u t i o n along radial d i r e c t i o n 1000
I
J
i
i
900
800
700
CA
600
o
500
......... a
o
'...., "'"-.
A
400 -
300
~
0.0
0.2
....
0.4
0.6 Radius
0.8
1.0
(m)
Fig. 9. Comparison of the numerical and experimental results for temperature in HS1 (modified model) (numerical: - - , z = z = l . 4 2 m ; .--, z=2.03m; ---, z=2.64m; experimental: ~ , z = 1.42 m; A, z = 2-03 m; C], z = 2.64 m).
Velocity d i s t r i b u t i o n along c e n t e r l i n e I
10
I
6v . . . . . - ' ............................. @9 0 '-¢ @9
..'"" 4
"'"-...
A
/
2-
0 0
I
l
1
2 Elevation
Fig. 10.
3
(m)
A x i a l variation in centerline velocity in HS1 (modified m o d e l ) (- • -, numerical solution; A , m e a s u r e m e n t ) .
T e m p e r a t u r e d i s t r i b u t i o n along radial d i r e c t i o n I
I
I
I
I
0.2
0.4
0.6
1300
I
1100
@9
900 -
@9
700 @9
300
0.0
- Radius
0.8
_ 1.0
(m)
Fig. 11. Comparison of the numerical and experimental results for temperature in H S 2 (modified m o d e l ) (numerical: - - , z = 1.42 m; • • . , z = 2.03 m; - - - , z = 2.54 m; experimental: ~ , z = 1-42 m; A , z = 2.03 m; V], z = 2.64 m).
Velocity d i s t r i b u t i o n along c e n t e r l i n e I
10
8
I
-
6-
."
A
A ......
/,,,/
.... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
,,-¢ O
/
4-
.,"
A
..........
/'
2-
0
E 1
0
I 2
3 Elevation (m) Axial variation in centerline velocity in HS2 (modified model) ( . . . , numerical solution; A, measurement).
Fig. 12.
1000/
T e m p e r a t u r e d i s t r i b u t i o n along r a d i a l d i r e c t i o n I
I
I
I
J
900 I ~
°°°F\ 00V.< "~
600
I.i'i%,._~kv },\
r
0.0
o
"'~:"2 o
0,2
0.4 0.6 0.B 1.0 Radius (m) Fig. 13. Comparison of the numerical and experimental results for temperature in HS1 (standard model) (numerical: - - , z = z = 1-42 m, • •., z = 2.03 m; - - - , z = 2.64 m; experimental: O, z = 1.42 m; A, Z = 2.03 m; [~, z = 2-64 m).
Velocity
distribution
IO
along
centerline
I
.,
6
..-"
...
..."
.,.. .
.
.
.
,.,.
,-,..
--,.,.
--..
../'"
--...
_
"''"-,.,,....
.. •
v
-,,,. ,...
,..
0
--. . . . . . . . . . .
/
4 /
2-
0 0
I
I
1
2
3
Elevation (m) Fig. 14.
Axial variation in centerline velocity in HS1 (standard model) (. •., numerical solution; O, measurement). Temperature
1300
distribution
i
along
l
radial
direction
i
I
llO0
v
900 L
700 -
",,
,,
5oo!
o
o
',
<>
-?...% 0~
300
I
0.0
0.2
o ~ " ~ - - , " "~" . " ' . - -'. ~ -~ ~_
~ ~
-
-
-
-
-
:
'
,
0.4
.....
0.6
~ - - _ _
I 0.8
1.0
Radius (m) Fig. 15. Comparison of the numerical and experimental results for temperature in HS2 (standard model) (numerical: - - , z = 1.42m; . . . , z = 2-03 m; - - - , z = 2-64m; experimental: ~ , z = 1.42 m, A, z = 2-03 m; O, z = 2-64 m).
248
Soonil Nam, Robert G. Bill, Jr
Velocity distribution 10
along centerline
r
i
..--- ............... "'""'"'"-..,
,, ,,
8" ' " ' -.
, . "'-... " " " " " . . . ,
6-
C9 0
A
. .........
Ix
Z~
4-
2-
0 0
f
I
I
2
3
E l e v a t i o n (m) Fig. 16.
Axial variation in centerline velocity in HS2 (standard model) (. •., numerical solution; A, measurement).
temperature and velocity than the measurements indicate, which prevented the penetration of water particles through the plume region in the earlier fire-sprinkler simulation. ~7 A careful measurement in turbulent axisymmetric buoyant plumes conducted by George et al. ~ predicts the following similar forms for temperature and velocity profiles: T-To - - exp (-657/2) T,,- T~ Pz
- exp ( - 5 5 n 2)
(20) (21)
Vz o
where 7/is a spatial similarity variable defined as r / z . The computational results for Test 8 were compared with the similarity forms. The solid lines in Figs 17 and 18 represent the functional forms in eqns (20) and (21), respectively. Both the calculated temperature and velocity profiles indicate good similarity in terms of the variable 77 which is consistent with the experimental observations. One of the practical interests of plume simulation is the prediction of
T e m p e r a t u r e s i m i l a r i t y for Test 8 1.0
,
,.
'
I
'
'
'
I
'
I
'-
¢D
,...¢
O
"~..,..
0.5
~ ~ :
0.0
,
~
,
,
0.00
I
J
,
,
,
0.05
I
,
,
,
.<.-<...~ .....
,
0,10
I
,
,
0,15
77
Fig. 17. T e m p e r a t u r e similarity in HS2 (experimental similarity, - - ; numerical, . . . , z/D =8; ---,
Velocity
z / D = 12. - . . . . , z / D = 14).
similarity
for
Test
8
1,0
0
0.5 0 0
001 0.00
J
I
l
B
I
l
0.05
l
J
L
I
B
l
0.10
~
I
I
L
L
0.1,5
Fig. 18. Velocity similarity in HS2 (experimental similarity, - - ; numerical, - . . , z/D =8; ---,
z / D = 12; . . . . .
, z / D = 14),
250
Soonil Nam, Robert G. Bill, Jr
Vorticity
contour
for plume
30
(Test #8)
I
15
O
10
~D
5
w
-
I
I
0
I
o
5 Radius (m) Fig. 19.
10
The contour of the vorticity lines to identify r~ for test 8.
air entrainment. The mass flow rate at z, rhz, is given as rhz =
~0rzOVz2~ dr
(22)
where r z denotes the width to the point where vz becomes zero. The axial velocity vz should be zero b e y o n d the plume width, r~, at z. H o w e v e r it is extremely difficult to define the point rz, both experimentally and numerically. The radius rz has been d e t e r m i n e d as the point where the vorticity, Wo, becomes zero. At each z point, Wo was c o m p u t e d along the r direction to identify the plume b o u n d a r y which is d e t e r m i n e d by the point where Wo becomes zero. Figure 19 shows the plume b o u n d a r y obtained through the zero vorticity line. The mass flow, which must be supplied from the e n t r a i n m e n t of air from outside of the plume boundary, has been c o m p u t e d by eqn (22) and c o m p a r e d with the relation developed by Heskestad: 11 mz = 0"071q1'3(z - Zo)5'3[1 + 0"0260c2'3(z - Zo) 5'31
(23)
where rhz is in kg/s, qc is in kW, z and Zo are in m. In Fig. 20, a comparison between the numerical results and eqn (23) indicates good
Numerical simulation of thermal plumes
Air e n t r a i n m e n t
251
in plume
I
I
7 O
O
6
.."
4
~.. ~ . ~ . ~ 7 / "
3
I 2
I 3
4
Elevation (z/zl) Fig. 20.
Air entrainment rates of the plumes in pool fires (numerical: --, Test 1; - •., Test 6, - - - , Test 8; plume formula: O, Test 1; A, Test 6; rl, Test 8).
agreement. All e n t r a i n m e n t rates were nondimensionalized by the mass flow rate at zl for comparison. The values of Zo and zl are provided in Table 3.
7 C E I L I N G JETS The thermal plumes discussed so far were g e n e r a t e d in regions with open boundaries. In o r d e r to simulate the interaction b e t w e e n the plume and sprinkler spray, the main objective of our current research, it is necessary to check the accuracy of plume flows with a ceiling. Numerical solutions of flows g e n e r a t e d with the same heat release rates at the two h e p t a n e spray fires, HS1 and HS2, were obtained u n d e r three different ceiling heights, 2.34, 3-86 and 4.5m. The m a x i m u m radius in the numerical d o m a i n for all cases was 7-24 m. Few experimental data are available to d e t e r m i n e the accuracy of the ceiling jet flow simulation. The scatter a m o n g available data on the m a x i m u m t e m p e r a t u r e and velocity is relatively large. Consequently, it
252
Soonil Nam, Robert G. Bill, Jr
is not easy to determine the accuracy of the numerical results based on the available experimental results. Heskestad and H a m a d a ~9 established that AT/ATo versus r/b gives a better correlation of the m a x i m u m excess t e m p e r a t u r e for strongly b u o y a n t plumes than AT/(Q2/3H 5/3) versus r/H which had been previously used. The two m e t h o d s provide approximately the same level of correlation for weak plumes. H e r e Q is the total heat release rate, AT is the m a x i m u m excess temperature, ATo is the excess t e m p e r a t u r e at the centerline of the undeflected plume at the level of ceiling height H, r is the radius along the ceiling, b is a measure of the width of the undeflected plume at the level of H. T h e y also recomm e n d e d correlation of v~/Vzo versus r/b, where v~,, is the centerline velocity of the undeflected plume at the level of H, instead of v~/(Q]/3H -~/3) versus r/H, for the velocity correlation of the ceiling ,
r
,
'
'
'
I
'
'
~
'
'
'
f [
1.0
<1
x~. "<-~,
[]
<1
[] ~
\
Aa .
A
"\ [] O
.1
I
0.1
I
J
I
I
i
I
L
J
1.0
t
I
I
I
L
I
I
I
'\,A \ O
I
10.0
Fig. 21. C o m p a r i s o n of the numerical solutions and m e a s u r e m e n t on m a x i m u m excess t e m p e r a t u r e (numerical: - - , H S I , H = 2 . 3 4 m ; ..-, HSI, H=3.86m; ---, HSI, H=4.5m; . . . . . , HS2, H = 2 - 3 4 m ; . . . . . , HS2, H = 3 . 8 6 m ; - - , HS2, H = 4 . 5 m ;
experimental: @ l-leskestad and Hamada I~ (propane burners, diameter ranges from 0.15 to 0.61 m, Q ranges from Ii.6 to 382 kW, H ranges from 0.56 to 2.50 m); A, Pickard et al.2° (alcohol pan fires, diameter ranges from 0.15 to 0.9 m, H ranges from 1.2 to 2.4 m, Q ranges from 4-85 to 128kW) E~, Thomas2] (alcohol pan fire, diameter = 0-42 m, Q = 106 kW, H = 1-37 m).
Numerical simulation of thermal plumes
253
1.0 A
A
~o
%
A A
O
A
O
Z~ /X O
\
A x" x
^',
0.1
1.5 Fig. 22.
r/b
lo
Comparison of the numerical solutions and measurement on ceiling jet
velocity: symbols for numerical results are the same as in Fig. 16 (O, King et al.n (double triwall carton with metal liner on wood pallet, growing fire under the quasi-steady assumption, 4¢ ranges from 6.1 to 7.4 MW, H ranges from 1.32 to 5.89 m); A, Pickard et al.2°; ~, Thomas2~). flow. The values for the undeflected plume flow can be obtained from eqns (14)-(16). Figures 21 and 22 show the comparison of m a x i m u m excess temperature and velocity of the ceiling flows b e t w e e n numerical results and m e a s u r e m e n t s . In order to convert the experimental data of Pickard and Thomas, which were plotted b a s e d on r/H, the convective heat release rate qc was assumed 75% of Q, and the location of the virtual origin, Zo, was c o m p u t e d by H e s k e s t a d ' s formula 9 as Zo(m) = - 1 . 0 2 O ( m ) + 0.083qc(kW) 2/5
(24)
w h e r e D is the diameter of the pool. The characteristic plume width, b, was chosen as the half-width of the plume c o m p u t e d from eqn (16) with the coefficient baT = 0.12. The numerical results b e c o m e m o r e scattered as r approaches the limit of the numerical domain due to the coarse grids and the influence of the numerical b o u n d a r y conditions. The overall comparison is quite favorable, except that the data of
Soonil Nam, Robert G. Bill, Jr
254
Velocity distribution I
'
'
~
'
I
'
'
'
'
I
'
along the ceiling
'
'
'
I
'
'
'
'
I
o
•~
:>,
O
O
~'~'-. [] ~"'~.~..
"~'~-~.
~'~ .......
8 8 0
l
0.5
,
,
,
,
I
1.0
,
,
t
,
I
J
,
,
~
1.5
I
2.0
~
,
,
,
I
2.5
,
,
i
i
3.0
r/H Fig. 23. Comparison of the numerical solutions and existing formulas on ceiling jet velocity (--, HS1, H = 2.34 m; • •., HS1, H = 3.86 m; - - - , HS1, H = 4-5 m; . . . . . , HS2, H=2-34m; ..... , HS2, H=3.86m; - - , HS2, H=4.5m; O, Heskestad and Delichatsios's formula; (3, Alpert's formula; both were quoted from Ref. 23 for steady axisymmetric flows). Pickard et al. show substantially higher velocities than the numerical ones in the region close to the centerline. Generally speaking, comparing results of different researchers, the t e m p e r a t u r e data are better correlated than the velocity data in the ceiling jet flows. 23 The velocity curves shown in Fig. 23 which are obtained by the numerical solutions were b e t w e e n H e s k e s t a d ' s correlation and Alpert's correlation as r e p o r t e d by Beyler 23 for steady axisymmetric ceiling jets.
8 CONCLUSION A computational model has been d e v e l o p e d to simulate the thermal plumes arising from large pool fires. B a s e d on plume temperatures, velocities, and spread rates, the standard k - E turbulence model was modified to yield optimal results. Experimental m e a s u r e m e n t s in pool fires and heptane spray fires agree well with the prediction obtained by
Numerical simulation of thermal plumes
255
the model. Ceiling flows were then simulated using the new model and the results were found to be in good agreement with experimentally available data for velocity and temperature. Comparisons of the results between the modified and the standard k - E model show that the present model significantly improves the agreement with the experimental results of buoyant thermal plumes. The new model, which simulates plume and ceiling flow with reasonable accuracy provides a powerful tool for the fire research.
ACKNOWLEDGEMENTS The authors are indebted to Dr Heskestad for many suggestions and insights, and for the experimental data he supplied. They also appreciate stimulating discussions with and experimental data provided by Drs Kung, You, Alpert, de Ris, Tamanini and Delichatsios. This work was supported by the Long Range Sprinkler Technology Program of the FMRC. The support of C. Yao, Manager of the Research Division, is acknowledged.
REFERENCES 1. Hossain, M. S. & Rodi, M., A turbulence model for buoyant flows and its application to vertical buoyant jets. Turbulent Buoyant Jets and Plumes, ed. W. Rodi. Pergamon Press, New York, 1982, pp. 121-78. 2. Davidson, L., Second-order correction of the k - e model to account for non-isotropic effects due to buoyancy. Int. J. Heat Mass Transfer, 33 (1990) 2599-608. 3. Shabbir, A. & Taulbee, D. B., Evaluation of turbulence models for predicting buoyant flows. J. Heat Transfer, ASME, 112 (Nov. 1990) 945-51. 4. Tamanini, F., The effect of buoyancy on the turbulence structure of vertical round jets. J. Heat Transfer, ASME, 100 (Nov. 1978) 659-64. 5. Markatos, N. C., Malin, M. R. & Cox, G., Mathematical modelling of buoyancy-induced smoke flow in enclosures. Int. J. Heat Mass Transfer, 25 (1982) 63-75. 6. Cox, G. & Kumar, S., Field modelling of fire in forced ventilated enclosures. Combust. Sci. Technol., 52 (1987) 7-23. 7. Morton, B. R., Taylor, G. I. & Turner, J. S., Turbulent gravitational convection from maintained and instantaneous sources. Proc. Roy. Soc., 234A (1956) 1-23. 8. Heskestad, G., Similarity relations for the initial convective flow generated by fire. 72-WA/HT-17, The Heat Transfer Division of ASME, Winter Annual Meeting, New York, November, 1972.
256
Soonil Nam, Robert G. Bill, Jr
9. Heskestad, G., Virtual origins of fire plumes. Fire Safety J., 5 (1983) 109-14. 10. Heskestad, G., Engineering relations for fire plumes. Fire Safety J., 7 (1984) 25-32. 11. Heskestad, G., Fire plume air entrainment according to two competing assumptions. Twenty-first Symposium (International) on Combustion, The Combustion Institute, Pittsburgh, PA, 1986, pp. 111-20. 12. Delichatsios, M. A., Air entrainment into buoyant jet flames and pool fires. Combust. Flame, 70 (1987) 33-46. 13. TR/100, The PHOENICS Beginner's Guide. CHAM Ltd, Wimbledon Village, London, 1990. 14. Kung, H. C. & Stavrianidis, P., Plume behavior of large-scale pool fires. FMRC J.I. 0COE6.RA 070(A), Factory Mutual Research Corporation, Norwood, MA, November, 1982. 15. Kung, H. C. & Stavrianidis, P., Buoyant plumes of large-scale pool fires. Nineteenth Symposium (International) on Combustion. The Combustion Institute, Pittsburgh, PA, 1982, pp. 905-12. 16. Launder, B. E. & Spalding, D. B., Mathematical Models of Turbulence. Academic Press, London, 1972. 17. Bill, R. G. & H. C. Kung (eds), FMRC Technical Report, J. I. OQ5NO-8.RA, August, Factory Mutual Research Corporation, Norwood, MA, 1990. 18. George, W. K., Alpert, R. L. & Tamanini, F., Turbulence measurements in an axisymmetric buoyant plume. Int. J. Heat Mass Transfer, 20 (1977) 1145-54.
19. Heskestad, G. & Hamada, T., Ceiling flows of strong fire plumes. FMRC Technical Report, J. I. OKOE1.RU, Factory Mutual Research Corporation, Norwood, MA, 1984. 20. Pickard, R. W., Hird, D. & Nash, P., Note No. 247, Fire Research Station, Boreham Wood, Herts, 1957. 21. Thomas, P. H., Note No. 141, Fire Research Stations, Boreham Wood, Herts, 1955. 22. Kung, H. C., Spaulding, R. D. & You, H. Z., Response of sprinkler links to rack storage fires. FMRC Technical Report, J. I. OG2E7.RA(2), Factory Mutual Research Corporation, Norwood, MA, 1984. 23. Beyler, C. L., Fire plumes and ceiling jets. Fire Safety J., 11 (1986) 53-75.
~ , ~
Numerical Simulation of Thermal Plumes
Soonil Nam & Robert G. Bill, Jr Factory Mutual Research Corporation, 1151 Boston-ProvidenceTurnpike, Norwood, Massachusetts 02062, USA (Received 22 April 1992; revised version received 28 August 1992; accepted 7 January 1993)
ABSTRACT Thermal plumes generated by pool fires and heptane spray fires were simulated numerically with the aid of the computer code PHOENICS. The flow fields associated with five cases were simulated and compared with experimental results. Modifications were made in the standard k - e model used in the P H O E N I C S code in order to improve the accuracy of that model in the simulation of turbulent thermal plumes. The modification includes an additional term related to the buoyancy generated turbulence kinetic energy, changing constants responsible for the turbulent viscosity and the effective Prandtl number in the energy equation. The results obtained with the modified code were in good agreement with experimental results and significantly improved the results obtained with the standard k-E model.
1 INTRODUCTION This study is the first part of an ongoing research program devoted to numerical modeling of the interaction of fire and sprinkler spray. The long range goal of the program is to develop computer models to predict the effectiveness of sprinkler systems, and reduce the need for costly experiments. One of the first difficulties encountered in simulating fire plumes is selecting an approriate and convenient turbulence model. It is well 231 Fire Safety Journal 0379-7112/93/$06-00 (~) 1993 Elsevier Science Publishers Ltd, England. Printed in Northern Ireland
232
Soonil Nam, Robert G. Bill, Jr
known that the k - e model tends to overpredict the velocities and temperatures along the centerline axis of buoyant plumes and consequently underpredicts the plume width.' 3 The k-E model has been improved for strongly buoyant flows through the addition of an algebraic stress model (ASM). 1~ The addition of an ASM provides a more accurate model for velocity and temperature distributions in a plume. The inclusion of an ASM, however, requires significantly increased computing time. Since the main purpose for our work is the development of a turbulence model which will be adopted for complex, and thus computation intensive, spray-plume interactions, we have chosen instead to modify the standard k-E model by adjusting its constants. This approach seemed reasonable since other studies 5"6 have suggested that the k - e model produces coarse agreement with measurements from fire and smoke flows. Therefore the major focus in this work is to find the optimal combination of some coefficients which are at our disposal in K-e model. The pool fire plume was selected for simulation because it has been extensively studied by many researchers, Morton et al., 7 Heskestad, ~ 1~ and Delichatsios; 12 thus data are available to validate the numerical results. In this study no combustion model was used. The domain of interest for comparison with experimental results, therefore, was above the flame region where the flow behaves as a pure thermal plume. The P H O E N I C S program, ~3 a commercial finite-difference computational fluid dynamics code, was used in this study. It employes the standard k-E model; however the model does not include a buoyancyinduced generation term in its turbulence kinetic energy generation formulation. Thus in addition to optimization of the k-E model constants, another term must be added to the turbulent kinetic energy balance. The modified code has been used to numerically simulate the axi-symmetric plumes generated by large diameter pool fires as reported by Kung and Stavrianidis. 14~5 A schematic of the flow to be simulated is shown in Fig. 1. The results of this study are also compared with plume correlations derived from many experiments. After constructing a turbulence model based on pool fires, the model was further validated with two sets of heptane-spray fire plumes.
2 MATHEMATICAL FORMULATION The plume arising from the pool fire is assumed to be axisymmetric, without swirl. Hence the conservation equations in cylindrical coordinates (with density fluctuation correlations ignored) are:
Numerical simulation of thermal plumes
233
\ \ \ \ . - ---iP-
\ \ \ \ \ \
j \
Vr
=
/ \
\
/
I I
\ \ \
/
Vr
/ /
I
/
\
Air Entrainment
\
\
i,
_71
r
Fig. 1.
The schematic diagram of a plume generated by a pool fire.
(1) Mass continuity:
at
-
(pozr)
r
-4- - - ( p O r t " )
Or
(1)
= 0
where vz and or denote z directional and r directional velocities. (2) z-momentum:
OpVz - -Ot+ r
(PvzrVz)+or((PVr - - ~) -
~z ~ r tZe ~ z J - Or r ~e Or/J
•
(2) _ O (.eO~)+!L S"~=Oz
Oz
f
O~,,
r Or krt~e~z ) -q- g(P=-- p)
(3)
where /z~ represents an effective viscosity which is the sum of a turbulent viscosity/z, and the laminar viscosity /zt, and p~ denotes the ambient density. (3) r-momentum:
apVr 1 1 0 Ot
~---f L ~Z (pllz
rv
O rv r) ~- ~ ( P ~ J r r) .
O {
Ovr~
-~z .~rtZe~z . }
.
O (rtZeOO___rr)] sot or /a
Or \
(4)
Soonil Nam, Robert G. Bill, Jr
234
Sot
OP Or
2/J, el) r -
-
r
+
O( OVz) l O { OVr~ OZ \txe Or / r or krtze--~r ) "
- -
+
(5)
(4) Energy: l{O
Oph + _ (pvzrh) + 0~- r -~z
;
O ( Ix~Oh] O (rlXeOh)] 0" (pv,rh) - Oz r - -= ~rh Oz : Or, oh -O~r
(6)
where h, Crh, q" denote enthalpy, effective Prandtl number for the diffusion of heat, and heat generation rate per unit volume, respectively. The effective Prandtl number in the energy equation, oh, is one of the two parameters which was varied to improve agreement with experimental results. (5) Turbulence kinetic energy k:
Opk~_1[0 Ot
r Oz (pvzrk) +
;
a (rlXtOk) a (rlX,k~l =Sk
(7)
(pv~rk) - -~z \ ~rk ~Z -- ~r \ ~rk Or/1
where ok is the effective Prandtl n u m b e r for the diffusion of turbulence energy and assumed equal to 1.0, and the turbulence production term Sk is composed of Sk = Gk + GB --
pE
(8)
where e, Gk, GB denote turbulence dissipation, shear production, and buoyancy production, respectively, and are defined as
f rl/Oldz\2 (OUrl2 Gk--lx'12[~--~z) + k O r /
(~---~)21 ( ~Uz
+
OVrl2~
+\Or+Oz/
10p
J
(9)
(lO)
GB -= p , , g - -
paz
Note that GB is not included in the standard k-E model and has been included in this study as in Ref. 6. (6) Turbulence dissipation E:
Ot
r Oz(PVzre)+~z(pVrre)-oz
\ro,Oz/-or
\ro,Or/J
where tr, is the effective Prandtl n u m b e r for the diffusion of turbulence dissipation and given as 1.314, and the production of dissipation term
Numerical simulation of thermal plumes
235
S, is defined as
S.=-k [C,(Gk + G.) - CEpE]
(12)
where C1 and C 2 a r e constants given as 1-44 and 1.92, respectively, in the standard model. The k - e model is adopted to obtain the effective turbulent viscosity, which is related to k and E by following formulation
ix, = OCDk2/E
(13)
Here Co is the other constant which was varied in the present study to improve the agreement with experimental results. All the other constants came from the work of Launder and Spalding. 16
3 M O D I F I C A T I O N OF T H E k-E M O D E L The primary task of the simulation is to devise an adequate turbulence model which produces accurate plume velocity and temperature distributions. Among many parameters in the k - e model, two parameters, CD and oh, most directly affect the diffusion of velocity and temperature in the flow. The standard k-E model assigns CD as 0.09, trh as 1.0. A test case with CD ----0.109, o-h----0-614 is reported in Ref. 1. It shows somewhat improved results for a pure jet. The spread rates of buoyant plumes were still significantly underpredicted; the modest improvements observed compared with the standard model was attributed to the inclusion of the buoyancy induced kinetic energy production term, rather than changing of the constants CD and ~rh. Since the standard model always predicts narrower velocity and temperature plume widths than measurements indicate, CD and o'h were varied to broaden the plume widths. It is obvious that a higher diffusion coefficient reduces peak values while increasing the effective region of transport of the property. Therefore we increased the turbulent viscosity by increasing the constant CD until reasonable agreement was obtained between the simulated velocity profile and measurements. The effective Prandtl number, ~rh, strongly affects the temperature distribution, but has little effect on the velocity distribution. A lower trh tends to increase the temperature width, thereby giving lower temperatures. After we obtained a reasonable range of CD based on the velocity distribution with a fixed trh, then both CD and trh were varied until satisfactory results for both temperature and velocity were obtained. The standards of the comparison were the well established plume
Soonil Nam, Robert G. Bill, Jr
236
relations ~°
UZo=Cv%(cpg )l'3qlc/3(z- Zo) 1/3 //
ATo = C To[
T.
~ l/3
"~
! qc-:/3(z _ Zo)-5/3
(14) (15)
\gt~pp~/ C
/ T . \1/2
o
(Z-Zo
(16)
where Vzo is the axial velocity, Vz, at the centerline, Oc is the convective heat relase rate, To is the absolute temperature at the centerline, ATo is the excess temperature above ambient at the centerline, bar is the plume half-width based on the centerline temperature, which is the width to the point where the temperature difference from the ambient becomes one half of the value at the centerline, and where Cv~,, Car,,, and Cb~r are demensionless constants. The constants Co~,., Car,,, C b j for the standards are chosen as 3.4, 9.1 and 0-12 following the recommendation of Heskestad. 1° After numerous trials with different combination of CD and o-h, CD=0"18, O'h=0"85 were chosen as the best combination which provided close agreement between the values of C~o, Car°, C~T and those in Ref. 10. As noted in Section 2, P H O E N I C S does not have a buoyancy produced turbulence kinetic energy term. In addition, P H O E N I C S assigns CD as 0.09 and o-h as 1.0. Therefore proper modifications were made to reflect the changes mentioned above. 4 TEST CONDITIONS The experiments to be simulated were pool fires with large enough diameters to generate turbulent plumes. The general conditions of the experiments reported in Refs 14 and 15 are given in Table 1. TABLE 1 Test Variables in the Experiments 14'17
Test no 1 4 6 8 HS1 HS2
Fuel
Methanol Silicone F1. Hydrocarbon F1. Heptane Heptane Heptane
Pool diam. D (m) 1.219 1.737 1-219 1.219 0.688 0.688
Convective heat flux (lc(Kw)
(l"(kW/m 2)
410 170 792 2 240 163.5 295 -4
351 72 679 1 920 440-1 795.1
Flame height (m) 0.9 0.2 2.2 3-9 ---
Numerical simulation o f thermal plumes
237
Tests 1, 4, 6 and 8 in Refs 14 and 15 were numerically simulated. The numerical simulations reported here are designated by the corresponding test n u m b e r in those references. Since the pure thermal plume behavior should not d e p e n d on any specific character of the fuels but only on the convective heat release rate, simulation of these tests which represent different fuels provides a good test for the validity of a numerical code. In addition to these simulations, two more cases with heptane spray fires ~v were tested for further validation of the model. The plumes were generated by the heptane spray fires at two convective heat release rates of 164kW and 295 kW, which are denoted as HS1 and HS2, respectively. Four heptane spray nozzles were placed at the corners of a 0.61 m square and directed radially toward the center. The conditions of these experiments are also given in Table 1. The flame heights of HS1 and HS2 were not measured.
5 COMPUTATIONAL PROCEDURES 5.1 Numerical domain for the tests
The numerical grids used for the simulations are given in Table 2 and the coordinate system is shown in Fig. 1. Uniform grids were adopted along the z direction, but nonuniform grids were used along the r direction. To place fine grids around r = 0, and coarser grids as r becomes larger, the following formula was adopted: (17)
r ( i ) = r o ( i / n ) ~5
TABLE 2 Numerical Domains for the Tests Test no.
1 4 6 8 HS1 HS2
Pool diam. (m)
1-219 1.737 1.219 1.219 0.688 0-688
Length o f domain (m)
Number o f meshes
r direction
z direction
r direction
z direction
12-52 17.84 12-52 12-52 7.24 7.24
15.24 15-0 30.48 30.84 7.27 7.27
30 60 30 30 30 30
51 51 51 51 51 51
Soonil Nam, Robert G. Bill, Jr
238
Numerical grid distribution
IIIIIIIF Ilfllll I
ililill i ililiil I Ililill I
liltllll Illli i t
.2 []
iliilll iliilll Ililll I Ilillll lillii lilli I Ilili I Illli i IIIli I Illli I Ilrll I iliili Ilill i I ililfl lillll I lillil I I illlii ilillr I ililli I I iliil i liliil f lillil I
H
Radius Fig.
2.
Numerical
grid
distribution.
where r(i) is the radial distance to the mesh i, ro is the maximum length of the radial numerical domain, and n is the total number of meshes in the r direction. The numerical grids used are shown in Fig. 2. In order to avoid a possible numerical instability caused by inappropriate initial guesses for velocities and temperatures with large heat release, an unsteady code was used. To ensure convergence, a small time stip, 0-5 s, was used for the initial stage and then larger time steps, 2.5 s, were employed after 100 time steps. The initial conditions were ambient except for the constant heat flux at the surface of the pool.
5.2 Boundary conditions A convenient way of providing the desired heat release is to assume that the heat release rate is uniform over the area of the pool. This avoids the complex flame shape required for a simulation using a volumetric heat release rate. Since we are only interested in the plume region of the flow, the accuracy of the simulation directly downstream of the area where heat is released is not important to this study.
Numerical simulation of thermal plumes
239
The two sets of heptane spray fires were treated as pool fires of diameter of 0.688 m, which corresponds to the hydraulic diameter of the square on which the nozzles were placed. Boundary conditions for the surfaces r =ro and z = 0 were that pressures are equal to the ambient, and the transport of m o m e n t u m and thermal energy are only allowed in the incoming direction. At each iteration, the direction of moving flow on the surface r =r0 was checked, and if outflow occurred, then the velocity at that point was set to zero until the computation predicted a value with the proper direction at the next iteration. The boundary conditions on the surface z = zoo, the end of the numerical domain in the z direction, were similar to those at r = ro; however, only outgoing instead of incoming flows were allowed. The boundary conditions at r - - 0 were the usual symmetry conditions: O/Or = 0 and v, = O.
6 RESULTS A N D DISCUSSION The variables p, h, Vz, vr, k and • were computed at each node for all time steps. The result for each test was compared with the result computed with approximately twice the number of z- and r-directional grid points reported in Table 2, in order to check the reliability of the numerical solutions. The temperature curve along the radial direction at z = 2zl (explanation is given in the following paragraph) was checked for consistency. All cases but Test 4 showed less than 2% of discrepancy. The heat release rate in Test 4 was at least an order of magnitude smaller than the heat release rates in the other simulated tests. Such a low heat release rate generated very low flows which were not well represented by the selected grid, and there was a practical limitation in number of grid points in the code. For this reason, Test 4 was not included in the analyses of the modified k - • model. Computation of each test took approximately 1 h in SUN Work Station Sparc IPX. The gas enthalpy h was transformed into the corresponding temperature for a comparison with measurements. The half-width based on the centerline velocity is the width to the point where the axial velocity decreases to one half of that of the centerline. The half-width for temperature is based on the excess temperature above ambient. Figures 3 and 4 show the velocity and temperature half-widths. The velocity half-widths are always a little wider than those of temperature. This is consistent with a review of literature data by Heskestad. 1° The plume
Plume -r
half-width
based
I
on velocity
I
I
/
/ ..-'
/ / ..'
/..
,-,
2 I , / . ".
I 1
-
...~" 1..-..~-S '~
0 0 Fig. 3.
5
10 Elevation
16
20
(m)
P l u m e h a l f - w i d t h s b a s e d on the centerline velocities ( - - , Test 1; . . . , - - Test, 8; - . - H S 1 ; - . . . . H S 2 ) . Plume
half-width I
based
on
I
Test 6;
temperature '
'
I
. ..'" /:// ./S 2
/
-
E
..5.'.-¢'" ~.Is
I
t~
1
0 0
5
10 Elevation
Fig. 4.
15
20
(m)
P l u m e half-widths based on the centerline temperatures ( - - , Test 1; • - -, Test 6; , Test 8; - . - , H S I ; . . . . , H S 2 ) .
Numerical simulation o f thermal plumes
241
TABLE 3 P l u m e R e g i o n s O b t a i n e d f r o m Numerical Simulation
Test no. Max. distance in z (m) 1 6 8 HS1 HS2
15-24 30-48 30.48 7.27 7.27
Plume region (m)
zt (m)
Virtual origin: zo (m)
2.2 --
2.2 3.0 4.0 1.2 1.2
-0.2 -0.4 -0"6 -0.2 -0.2
regions of the numerical solution were identified by the elevations where the plume half-widths first exhibited a nearly linear increase with elevation (see eqn 16); their values, z~, are shown in Table 3. In order to compare the numerical results with the plume relations from eqns (14)-(16), the distance to the virtual origin Zo must be found for each test. Equations (14) and (15) are modified for that purpose as follows:
z=C zo;?( g-- ]q<+z0 o o kCpp~TJ | ~ / Z = /C3/5 / baT/---'~-'~_2! L
\gt~pp~l
~.e/5/ATo315 t/c / -I- Zo
j
(18)
(19)
Applying a least-square fitting provides the best fitting constants for Co,o, Zo in eqn (18), and Char, Zo in eqn (19). The values of Zo obtained from the two equations were very close to each other for all tests. The values of Zo given in Table 3 are obtained from eqn (19), based on the centerline temperature, and are used in all the following analyses. After inserting Zo into plume relations from (14) to (16), the coefficients Cozo, Car°, Cb~ are plotted in Figs 5, 6 and 7. All coefficients show quite good agreement with the values quoted previously from Ref. 10: 3.4, 9.1 and 0.12. In order to demonstrate the effect of modified k-E model constants, the results obtained with the standard k - e model for Test 8 are shown in Figs 5-7. The standard k - e model contains the terms which take into account the kinetic energy production due to the buoyancy. The results clearly demonstrate that the modified version predicts the plume behavior far better than the standard k-E model. The standard k - e model overpredicts centerline excess temperature and velocity in Test 8 by 68% and 25%, respectively, which results in a much narrower temperature width for a given plume. Other tests with the standard k - e model also produced similar results, however, they were not included
5
'
1
'
I
I
2
3
Nondimensional
4
(z/zl)
Elevation
Fig. 5. C o m p a r i s o n of the p l u m e coefficient C,~ ( - - , c o n s t a n t in Ref. 10; • • -, Test l; - - , Test 6; - - - , Test 8; . . . . . , HS1; - - , HS2; o . . _ _ Test 8 with the s t a n d a r d k - e model). 20
~
,
12 :'Z.'.:.;.-= -.': . . . . . . .. . . . . . . . . . . .
1
2 Nondimensional
Fig.
--,
Elevation
3 (z/zl)
4
C o m p a r i s o n of the p l u m e coefficient CAr° ( - - , c o n s t a n t in Ref. 10; • •., Test 1; Test 6; - - - , Test 8; . . . . . , HS1; - - , HS2; . . . . , Test 8 with the s t a n d a r d k - e model).
6.
Numerical simulation of thermal plumes
0.2
i
i
I
I
243
0.1
0.0 1
2 3 N o n d i m e n s i o n a l Elevation (z/zz)
4
Fig. 7. Comparison of the plume coefficient Char.
in the figures to avoid being overcrowded. The simulations were compared with the experiments reported in Refs 14 and 15. As stated in Ref. 14, the heat release rates, qc, of the tests were c o m p u t e d based on the measurements at z = 7.62 m. These heat release rate values were used to generate the buoyant plumes in the present numerical simulations. The measured temperatures 14 at z = 7 . 6 2 m and the numerical solutions are shown in Fig. 8 and are in reasonable agreement. The experimental data for Test 8 was not available. Additional comparisons are provided in Figs 9-12 in simulations of the plumes of heptane-spray fires. The overall comparison between numerical results and experimental measurements is good. Figure 9 shows the comparison of temperature and Fig. 10 shows the comparison of velocity with measurements for HS1. Figures 11 and 12 show the same comparisons, respectively, for HS2. Though the earlier comparison based on the plume coefficients is more general and more reliable, the similar comparisons as shown in Figs 9-12 are provided in Figs 13-16 to help to visualize the difference between the modified and the standard model. The numerical results shown in Figs 13-16 are obtained with the standard k-E model. Note the higher plume
P l u m e t e m p e r a t u r e at z = 7 . 6 2 m 1
400
i
J
i
I
i
,
i
i
i
i
i
I
i
i
i
~
]
J
i
i
,
380
O
v
360
t_ c~
340
-
z~
320
".. . ,~
A
O
"'""....
3O0
0.0
0.5
1.0
1.5 Radius
2.0
2.5
(m)
Fig. 8. Comparison of the numerical and experimental results for temperature in pool fires (numerical: - - , Test 1; • • ,, Test 6; experimental: A, Test 1; ~ , Test 6).
T e m p e r a t u r e d i s t r i b u t i o n along radial d i r e c t i o n 1000
I
J
i
i
900
800
700
CA
600
o
500
......... a
o
'...., "'"-.
A
400 -
300
~
0.0
0.2
....
0.4
0.6 Radius
0.8
1.0
(m)
Fig. 9. Comparison of the numerical and experimental results for temperature in HS1 (modified model) (numerical: - - , z = z = l . 4 2 m ; .--, z=2.03m; ---, z=2.64m; experimental: ~ , z = 1.42 m; A, z = 2-03 m; C], z = 2.64 m).
Velocity d i s t r i b u t i o n along c e n t e r l i n e I
10
I
6v . . . . . - ' ............................. @9 0 '-¢ @9
..'"" 4
"'"-...
A
/
2-
0 0
I
l
1
2 Elevation
Fig. 10.
3
(m)
A x i a l variation in centerline velocity in HS1 (modified m o d e l ) (- • -, numerical solution; A , m e a s u r e m e n t ) .
T e m p e r a t u r e d i s t r i b u t i o n along radial d i r e c t i o n I
I
I
I
I
0.2
0.4
0.6
1300
I
1100
@9
900 -
@9
700 @9
300
0.0
- Radius
0.8
_ 1.0
(m)
Fig. 11. Comparison of the numerical and experimental results for temperature in H S 2 (modified m o d e l ) (numerical: - - , z = 1.42 m; • • . , z = 2.03 m; - - - , z = 2.54 m; experimental: ~ , z = 1-42 m; A , z = 2.03 m; V], z = 2.64 m).
Velocity d i s t r i b u t i o n along c e n t e r l i n e I
10
8
I
-
6-
."
A
A ......
/,,,/
.... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
,,-¢ O
/
4-
.,"
A
..........
/'
2-
0
E 1
0
I 2
3 Elevation (m) Axial variation in centerline velocity in HS2 (modified model) ( . . . , numerical solution; A, measurement).
Fig. 12.
1000/
T e m p e r a t u r e d i s t r i b u t i o n along r a d i a l d i r e c t i o n I
I
I
I
J
900 I ~
°°°F\ 00V.< "~
600
I.i'i%,._~kv },\
r
0.0
o
"'~:"2 o
0,2
0.4 0.6 0.B 1.0 Radius (m) Fig. 13. Comparison of the numerical and experimental results for temperature in HS1 (standard model) (numerical: - - , z = z = 1-42 m, • •., z = 2.03 m; - - - , z = 2.64 m; experimental: O, z = 1.42 m; A, Z = 2.03 m; [~, z = 2-64 m).
Velocity
distribution
IO
along
centerline
I
.,
6
..-"
...
..."
.,.. .
.
.
.
,.,.
,-,..
--,.,.
--..
../'"
--...
_
"''"-,.,,....
.. •
v
-,,,. ,...
,..
0
--. . . . . . . . . . .
/
4 /
2-
0 0
I
I
1
2
3
Elevation (m) Fig. 14.
Axial variation in centerline velocity in HS1 (standard model) (. •., numerical solution; O, measurement). Temperature
1300
distribution
i
along
l
radial
direction
i
I
llO0
v
900 L
700 -
",,
,,
5oo!
o
o
',
<>
-?...% 0~
300
I
0.0
0.2
o ~ " ~ - - , " "~" . " ' . - -'. ~ -~ ~_
~ ~
-
-
-
-
-
:
'
,
0.4
.....
0.6
~ - - _ _
I 0.8
1.0
Radius (m) Fig. 15. Comparison of the numerical and experimental results for temperature in HS2 (standard model) (numerical: - - , z = 1.42m; . . . , z = 2-03 m; - - - , z = 2-64m; experimental: ~ , z = 1.42 m, A, z = 2-03 m; O, z = 2-64 m).
248
Soonil Nam, Robert G. Bill, Jr
Velocity distribution 10
along centerline
r
i
..--- ............... "'""'"'"-..,
,, ,,
8" ' " ' -.
, . "'-... " " " " " . . . ,
6-
C9 0
A
. .........
Ix
Z~
4-
2-
0 0
f
I
I
2
3
E l e v a t i o n (m) Fig. 16.
Axial variation in centerline velocity in HS2 (standard model) (. •., numerical solution; A, measurement).
temperature and velocity than the measurements indicate, which prevented the penetration of water particles through the plume region in the earlier fire-sprinkler simulation. ~7 A careful measurement in turbulent axisymmetric buoyant plumes conducted by George et al. ~ predicts the following similar forms for temperature and velocity profiles: T-To - - exp (-657/2) T,,- T~ Pz
- exp ( - 5 5 n 2)
(20) (21)
Vz o
where 7/is a spatial similarity variable defined as r / z . The computational results for Test 8 were compared with the similarity forms. The solid lines in Figs 17 and 18 represent the functional forms in eqns (20) and (21), respectively. Both the calculated temperature and velocity profiles indicate good similarity in terms of the variable 77 which is consistent with the experimental observations. One of the practical interests of plume simulation is the prediction of
T e m p e r a t u r e s i m i l a r i t y for Test 8 1.0
,
,.
'
I
'
'
'
I
'
I
'-
¢D
,...¢
O
"~..,..
0.5
~ ~ :
0.0
,
~
,
,
0.00
I
J
,
,
,
0.05
I
,
,
,
.<.-<...~ .....
,
0,10
I
,
,
0,15
77
Fig. 17. T e m p e r a t u r e similarity in HS2 (experimental similarity, - - ; numerical, . . . , z/D =8; ---,
Velocity
z / D = 12. - . . . . , z / D = 14).
similarity
for
Test
8
1,0
0
0.5 0 0
001 0.00
J
I
l
B
I
l
0.05
l
J
L
I
B
l
0.10
~
I
I
L
L
0.1,5
Fig. 18. Velocity similarity in HS2 (experimental similarity, - - ; numerical, - . . , z/D =8; ---,
z / D = 12; . . . . .
, z / D = 14),
250
Soonil Nam, Robert G. Bill, Jr
Vorticity
contour
for plume
30
(Test #8)
I
15
O
10
~D
5
w
-
I
I
0
I
o
5 Radius (m) Fig. 19.
10
The contour of the vorticity lines to identify r~ for test 8.
air entrainment. The mass flow rate at z, rhz, is given as rhz =
~0rzOVz2~ dr
(22)
where r z denotes the width to the point where vz becomes zero. The axial velocity vz should be zero b e y o n d the plume width, r~, at z. H o w e v e r it is extremely difficult to define the point rz, both experimentally and numerically. The radius rz has been d e t e r m i n e d as the point where the vorticity, Wo, becomes zero. At each z point, Wo was c o m p u t e d along the r direction to identify the plume b o u n d a r y which is d e t e r m i n e d by the point where Wo becomes zero. Figure 19 shows the plume b o u n d a r y obtained through the zero vorticity line. The mass flow, which must be supplied from the e n t r a i n m e n t of air from outside of the plume boundary, has been c o m p u t e d by eqn (22) and c o m p a r e d with the relation developed by Heskestad: 11 mz = 0"071q1'3(z - Zo)5'3[1 + 0"0260c2'3(z - Zo) 5'31
(23)
where rhz is in kg/s, qc is in kW, z and Zo are in m. In Fig. 20, a comparison between the numerical results and eqn (23) indicates good
Numerical simulation of thermal plumes
Air e n t r a i n m e n t
251
in plume
I
I
7 O
O
6
.."
4
~.. ~ . ~ . ~ 7 / "
3
I 2
I 3
4
Elevation (z/zl) Fig. 20.
Air entrainment rates of the plumes in pool fires (numerical: --, Test 1; - •., Test 6, - - - , Test 8; plume formula: O, Test 1; A, Test 6; rl, Test 8).
agreement. All e n t r a i n m e n t rates were nondimensionalized by the mass flow rate at zl for comparison. The values of Zo and zl are provided in Table 3.
7 C E I L I N G JETS The thermal plumes discussed so far were g e n e r a t e d in regions with open boundaries. In o r d e r to simulate the interaction b e t w e e n the plume and sprinkler spray, the main objective of our current research, it is necessary to check the accuracy of plume flows with a ceiling. Numerical solutions of flows g e n e r a t e d with the same heat release rates at the two h e p t a n e spray fires, HS1 and HS2, were obtained u n d e r three different ceiling heights, 2.34, 3-86 and 4.5m. The m a x i m u m radius in the numerical d o m a i n for all cases was 7-24 m. Few experimental data are available to d e t e r m i n e the accuracy of the ceiling jet flow simulation. The scatter a m o n g available data on the m a x i m u m t e m p e r a t u r e and velocity is relatively large. Consequently, it
252
Soonil Nam, Robert G. Bill, Jr
is not easy to determine the accuracy of the numerical results based on the available experimental results. Heskestad and H a m a d a ~9 established that AT/ATo versus r/b gives a better correlation of the m a x i m u m excess t e m p e r a t u r e for strongly b u o y a n t plumes than AT/(Q2/3H 5/3) versus r/H which had been previously used. The two m e t h o d s provide approximately the same level of correlation for weak plumes. H e r e Q is the total heat release rate, AT is the m a x i m u m excess temperature, ATo is the excess t e m p e r a t u r e at the centerline of the undeflected plume at the level of ceiling height H, r is the radius along the ceiling, b is a measure of the width of the undeflected plume at the level of H. T h e y also recomm e n d e d correlation of v~/Vzo versus r/b, where v~,, is the centerline velocity of the undeflected plume at the level of H, instead of v~/(Q]/3H -~/3) versus r/H, for the velocity correlation of the ceiling ,
r
,
'
'
'
I
'
'
~
'
'
'
f [
1.0
<1
x~. "<-~,
[]
<1
[] ~
\
Aa .
A
"\ [] O
.1
I
0.1
I
J
I
I
i
I
L
J
1.0
t
I
I
I
L
I
I
I
'\,A \ O
I
10.0
Fig. 21. C o m p a r i s o n of the numerical solutions and m e a s u r e m e n t on m a x i m u m excess t e m p e r a t u r e (numerical: - - , H S I , H = 2 . 3 4 m ; ..-, HSI, H=3.86m; ---, HSI, H=4.5m; . . . . . , HS2, H = 2 - 3 4 m ; . . . . . , HS2, H = 3 . 8 6 m ; - - , HS2, H = 4 . 5 m ;
experimental: @ l-leskestad and Hamada I~ (propane burners, diameter ranges from 0.15 to 0.61 m, Q ranges from Ii.6 to 382 kW, H ranges from 0.56 to 2.50 m); A, Pickard et al.2° (alcohol pan fires, diameter ranges from 0.15 to 0.9 m, H ranges from 1.2 to 2.4 m, Q ranges from 4-85 to 128kW) E~, Thomas2] (alcohol pan fire, diameter = 0-42 m, Q = 106 kW, H = 1-37 m).
Numerical simulation of thermal plumes
253
1.0 A
A
~o
%
A A
O
A
O
Z~ /X O
\
A x" x
^',
0.1
1.5 Fig. 22.
r/b
lo
Comparison of the numerical solutions and measurement on ceiling jet
velocity: symbols for numerical results are the same as in Fig. 16 (O, King et al.n (double triwall carton with metal liner on wood pallet, growing fire under the quasi-steady assumption, 4¢ ranges from 6.1 to 7.4 MW, H ranges from 1.32 to 5.89 m); A, Pickard et al.2°; ~, Thomas2~). flow. The values for the undeflected plume flow can be obtained from eqns (14)-(16). Figures 21 and 22 show the comparison of m a x i m u m excess temperature and velocity of the ceiling flows b e t w e e n numerical results and m e a s u r e m e n t s . In order to convert the experimental data of Pickard and Thomas, which were plotted b a s e d on r/H, the convective heat release rate qc was assumed 75% of Q, and the location of the virtual origin, Zo, was c o m p u t e d by H e s k e s t a d ' s formula 9 as Zo(m) = - 1 . 0 2 O ( m ) + 0.083qc(kW) 2/5
(24)
w h e r e D is the diameter of the pool. The characteristic plume width, b, was chosen as the half-width of the plume c o m p u t e d from eqn (16) with the coefficient baT = 0.12. The numerical results b e c o m e m o r e scattered as r approaches the limit of the numerical domain due to the coarse grids and the influence of the numerical b o u n d a r y conditions. The overall comparison is quite favorable, except that the data of
Soonil Nam, Robert G. Bill, Jr
254
Velocity distribution I
'
'
~
'
I
'
'
'
'
I
'
along the ceiling
'
'
'
I
'
'
'
'
I
o
•~
:>,
O
O
~'~'-. [] ~"'~.~..
"~'~-~.
~'~ .......
8 8 0
l
0.5
,
,
,
,
I
1.0
,
,
t
,
I
J
,
,
~
1.5
I
2.0
~
,
,
,
I
2.5
,
,
i
i
3.0
r/H Fig. 23. Comparison of the numerical solutions and existing formulas on ceiling jet velocity (--, HS1, H = 2.34 m; • •., HS1, H = 3.86 m; - - - , HS1, H = 4-5 m; . . . . . , HS2, H=2-34m; ..... , HS2, H=3.86m; - - , HS2, H=4.5m; O, Heskestad and Delichatsios's formula; (3, Alpert's formula; both were quoted from Ref. 23 for steady axisymmetric flows). Pickard et al. show substantially higher velocities than the numerical ones in the region close to the centerline. Generally speaking, comparing results of different researchers, the t e m p e r a t u r e data are better correlated than the velocity data in the ceiling jet flows. 23 The velocity curves shown in Fig. 23 which are obtained by the numerical solutions were b e t w e e n H e s k e s t a d ' s correlation and Alpert's correlation as r e p o r t e d by Beyler 23 for steady axisymmetric ceiling jets.
8 CONCLUSION A computational model has been d e v e l o p e d to simulate the thermal plumes arising from large pool fires. B a s e d on plume temperatures, velocities, and spread rates, the standard k - E turbulence model was modified to yield optimal results. Experimental m e a s u r e m e n t s in pool fires and heptane spray fires agree well with the prediction obtained by
Numerical simulation of thermal plumes
255
the model. Ceiling flows were then simulated using the new model and the results were found to be in good agreement with experimentally available data for velocity and temperature. Comparisons of the results between the modified and the standard k - E model show that the present model significantly improves the agreement with the experimental results of buoyant thermal plumes. The new model, which simulates plume and ceiling flow with reasonable accuracy provides a powerful tool for the fire research.
ACKNOWLEDGEMENTS The authors are indebted to Dr Heskestad for many suggestions and insights, and for the experimental data he supplied. They also appreciate stimulating discussions with and experimental data provided by Drs Kung, You, Alpert, de Ris, Tamanini and Delichatsios. This work was supported by the Long Range Sprinkler Technology Program of the FMRC. The support of C. Yao, Manager of the Research Division, is acknowledged.
REFERENCES 1. Hossain, M. S. & Rodi, M., A turbulence model for buoyant flows and its application to vertical buoyant jets. Turbulent Buoyant Jets and Plumes, ed. W. Rodi. Pergamon Press, New York, 1982, pp. 121-78. 2. Davidson, L., Second-order correction of the k - e model to account for non-isotropic effects due to buoyancy. Int. J. Heat Mass Transfer, 33 (1990) 2599-608. 3. Shabbir, A. & Taulbee, D. B., Evaluation of turbulence models for predicting buoyant flows. J. Heat Transfer, ASME, 112 (Nov. 1990) 945-51. 4. Tamanini, F., The effect of buoyancy on the turbulence structure of vertical round jets. J. Heat Transfer, ASME, 100 (Nov. 1978) 659-64. 5. Markatos, N. C., Malin, M. R. & Cox, G., Mathematical modelling of buoyancy-induced smoke flow in enclosures. Int. J. Heat Mass Transfer, 25 (1982) 63-75. 6. Cox, G. & Kumar, S., Field modelling of fire in forced ventilated enclosures. Combust. Sci. Technol., 52 (1987) 7-23. 7. Morton, B. R., Taylor, G. I. & Turner, J. S., Turbulent gravitational convection from maintained and instantaneous sources. Proc. Roy. Soc., 234A (1956) 1-23. 8. Heskestad, G., Similarity relations for the initial convective flow generated by fire. 72-WA/HT-17, The Heat Transfer Division of ASME, Winter Annual Meeting, New York, November, 1972.
256
Soonil Nam, Robert G. Bill, Jr
9. Heskestad, G., Virtual origins of fire plumes. Fire Safety J., 5 (1983) 109-14. 10. Heskestad, G., Engineering relations for fire plumes. Fire Safety J., 7 (1984) 25-32. 11. Heskestad, G., Fire plume air entrainment according to two competing assumptions. Twenty-first Symposium (International) on Combustion, The Combustion Institute, Pittsburgh, PA, 1986, pp. 111-20. 12. Delichatsios, M. A., Air entrainment into buoyant jet flames and pool fires. Combust. Flame, 70 (1987) 33-46. 13. TR/100, The PHOENICS Beginner's Guide. CHAM Ltd, Wimbledon Village, London, 1990. 14. Kung, H. C. & Stavrianidis, P., Plume behavior of large-scale pool fires. FMRC J.I. 0COE6.RA 070(A), Factory Mutual Research Corporation, Norwood, MA, November, 1982. 15. Kung, H. C. & Stavrianidis, P., Buoyant plumes of large-scale pool fires. Nineteenth Symposium (International) on Combustion. The Combustion Institute, Pittsburgh, PA, 1982, pp. 905-12. 16. Launder, B. E. & Spalding, D. B., Mathematical Models of Turbulence. Academic Press, London, 1972. 17. Bill, R. G. & H. C. Kung (eds), FMRC Technical Report, J. I. OQ5NO-8.RA, August, Factory Mutual Research Corporation, Norwood, MA, 1990. 18. George, W. K., Alpert, R. L. & Tamanini, F., Turbulence measurements in an axisymmetric buoyant plume. Int. J. Heat Mass Transfer, 20 (1977) 1145-54.
19. Heskestad, G. & Hamada, T., Ceiling flows of strong fire plumes. FMRC Technical Report, J. I. OKOE1.RU, Factory Mutual Research Corporation, Norwood, MA, 1984. 20. Pickard, R. W., Hird, D. & Nash, P., Note No. 247, Fire Research Station, Boreham Wood, Herts, 1957. 21. Thomas, P. H., Note No. 141, Fire Research Stations, Boreham Wood, Herts, 1955. 22. Kung, H. C., Spaulding, R. D. & You, H. Z., Response of sprinkler links to rack storage fires. FMRC Technical Report, J. I. OG2E7.RA(2), Factory Mutual Research Corporation, Norwood, MA, 1984. 23. Beyler, C. L., Fire plumes and ceiling jets. Fire Safety J., 11 (1986) 53-75.