ARTICLE IN PRESS
Fire Safety Journal 42 (2007) 127–138 www.elsevier.com/locate/firesaf
Validation of FDS for the prediction of medium-scale pool fires J.X. Wen, K. Kang, T. Donchev, J.M. Karwatzki Faculty of Engineering, Kingston University, Friars Avenue, Roehampton Vale, London SW15 3DW, UK Received 10 January 2006; received in revised form 26 June 2006; accepted 10 August 2006
Abstract Fire dynamics simulator (FDS) has been applied to simulate a medium-scale methanol pool fire. The simulation used predominantly the existing features in FDS except that an additional sub-grid-scale combustion model based on the laminar flamelet approach of Cook AW and Riley JJ [Combust and Flame 1998;112:593–606] was used alongside the default mixture fraction combustion model for comparison. The predictions of the two different combustion models for temperature and axial velocity distributions were found to be in reasonably good agreement with each other and the experimental data. The pulsating nature of air entrainment was demonstrated by the air entrainment velocity fluctuations and the instantaneous velocity vectors, which revealed formation and shedding of vortices and the well-known ‘‘neck-in’’ at a distance of approximately one diameter from the pool surface. The predicted variations of air entrainment at different heights agreed well with some published data and correlation. Although the limitation of the code in predicting the puffing frequency was noticed as the spectra of temperature fluctuations failed to demonstrate any dominant frequency, the present study has demonstrated the capability of FDS to deliver reliable predictions on most important parameters of pool fires. r 2006 Elsevier Ltd. All rights reserved. Keywords: Pool fires; Large eddy simulation (LES); SGS turbulence modelling; SGS combustion modelling
1. Introduction The fire dynamics simulator (FDS), developed by the National Institute of Standards and Technology, is receiving increasingly wide applications within the fire community. Following a recent paper on small pool fire studies also conducted with FDS [1], this paper reports on some validation studies of the code for the predictions of medium-scale pool fires. The pool fire scenario is chosen as it represents a simplified version of many fires that occur in practice, such as furniture fires. It is one of the most basic forms of fuel combustion and is often present in accidental fires. Most early investigations on pool fires were conducted in the fully turbulent far field. Gebhart and co-workers [2] pioneered the study of buoyant plume/jet instabilities and flow transitions in wall-bounded and unbounded buoyant flows. More recently, these phenomena have been experiCorresponding author. Tel.: +44 020 8547 7836; fax: +44 020 8547 7992. E-mail address:
[email protected] (J.X. Wen).
0379-7112/$ - see front matter r 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.firesaf.2006.08.007
mentally investigated by Cetegen and co-workers for buoyant plumes [3], and by Weckman and Sobiesiak [4], Cegen and Ahmed [5], Malasekera et al. [6], Hamins et al. [7] and Mandin and Most [8] for pool fires. These investigations have led to the development of correlation that relates oscillation frequency to the convection time scale for the flame. More detailed reviews on pool fire study can be found in Pagni [9] and Joulain [10]. Numerous investigations have also been conducted using CFD-based techniques. Blunsdon et al. [11] used the RANS-based CFD approach to simulate the transient behaviour of buoyant turbulent diffusion flames. They considered a flame above a honeycomb mesh, and assigned ‘‘inlet’’ values (to those quantities) which are thought to be rather high [12]. The results did not capture the large, toroidal vortices that are closely related to puffing. The recent CFD simulation of Sinai [12] did not assume any extraneous supply of turbulence energy at the pool and captured the torodial vortices at a frequency which agreed with experiments but the predicted temperature was found to decay too slowly with height. Liu and Wen [13] also used the RANS approach to model the behaviour of pool fires
ARTICLE IN PRESS J.X. Wen et al. / Fire Safety Journal 42 (2007) 127–138
128
Nomenclature A Cs D Dl F Frf h I _ ent m _f m P Pe q_ 000 R r Re S Sc s T t U
cross-sectional area Smagorinsky constant diameter diffusivity external force Froude number enthalpy radial intensity entrained mass flow rate mass burning rate pressure Peclect number heat source ideal gas constant stoichiometric coefficient Reynolds number strain stress Schmidt number stoichiometric number temperature time velocity
but the focus of their work was on the effect of turbulence models. Ghoniem et al. [14] conducted vortex-based simulation of an axis-symmetric fire plume. Their analysis suggested that puffing results from an intrinsic flow instability which was a Kelvin–Helmholtz-type mechanism associated with the formation of an axis-symmetric vortex sheet along the boundary between the inner hot mixture and the outer colder gases. Mell et al. [15] and Jiang and Luo [16] carried out DNS on axis-symmetric and planar plumes, but similar to the vortex-based study of Ghoniem et al. [14], the applications of the simulations is limited by the 2D characteristics of the calculations. Furthermore, DNS is still limited to small Reynolds numbers and the RANS or FANS-based approaches are unable to simulate the interaction between the many length and time scales which exist in turbulent reacting flows. Large eddy simulation (LES) is a promising alternative. It is computationally less expensive than DNS, but still capable of tackling the scale-dependent dynamic behaviour. Baum et al. [17] simulated fire plumes using an earlier version of FDS where the fire itself was prescribed as an ensemble of thermal elements ejecting from the burning surface at an initial velocity and radiation was prescribed as a fixed percentage heat loss. Recently, the present authors have reported on some small pool fire studies using FDS where the focus was on some unique characteristic of small pool fires such as flame anchoring and double flame [1]. In the current paper, we examine the predictions for a mediumscale pool fire.
Uent air entrainment velocity u, v, w three components of velocity vectors w_ reaction rate Y mass fraction Z height Others r rN x w s z k
density density of air mixture fraction scalar dissipation rate SGS item in momentum equation SGS item in mixture fraction equation SGS item in energy equation
Subscripts i, j, k f O p sf w
three directions in co-ordinates fuel oxygen combustion products stoichiometric surface wall
2. Numerical modelling 2.1. Governing equations The filtered equations can be written as qr¯ qr¯ u~ j þ ¼ 0, qt qxj
(1)
qtij qu~ i qu~ i u~ j q¯p þ r¯ þ rg , ¼ ¯ þfi þ qxi qt qxj qxj
(2)
qðr¯ Y~ l Þ qðr¯ u~ j Y~ l Þ q þ ¼ ðruj Y l r¯ u~ j Y~ l Þ qt qxj qxj q qY l rDl þ þ w¯_ l , qxj qxj ~ ~ qðr¯ u~ j hÞ qðr¯ hÞ D¯p qqrj q qT~ þ ¼ þ k qxj qt Dt qxj qxj qxj X q qY~ l ~ rD þ , ¯ l hl qxj qxj l
ð3Þ
ð4Þ
where r¯ is the filtered density, k is the thermal conductivity, ~ Y~ l Cp is the specific heat, Dl is the material diffusivity, u, and h~ are the Favre-filtered velocity, mass species concentration and enthalpy, qr is the radiative heat flux, ~_ l is the Favre-filtered Dp/Dt is the material derivative and w species source term. A Favre-filtered quantity denoted by a ~ ¼ rf=r. tilde is defined by f ¯
ARTICLE IN PRESS J.X. Wen et al. / Fire Safety Journal 42 (2007) 127–138
2.2. SGS turbulence modelling The Smagorinsky’s [18] eddy viscosity model is used as the SGS turbulence model with a coefficient (Cs) of 0.17. Smagorinsky [18] originally proposed a value of 0.23. Previous researchers have used values ranging from 0.1 to 0.23 for different applications. The value of 0.17 was chosen after referencing previously reported FDS smoke movement studies which used a value of 0.14 [19] and the value which is generally used for indoor air quality study which is 0.16. In FDS version 4 Technical Reference Guide [20], a default value of 0.2 is suggested following some numerical tests of different Smagorinsky parameters carried out by Zhang et al. [21]. A value of 0.5 was specified for the turbulent Schmidt and Prandtl numbers. It was found in the tests of Zhang et al. [21] that these two numbers have much less influence on the results than the Smagorinsky constant. 2.3. SGS combustion modelling Two combustion models are used in the current simulation. One is developed by Mell et al. [15] for the FDS code (referred to as ‘‘mixture fraction model’’) and the other is the laminar flamelet approach of Cook and Riley [22] which has been slightly modified by the present authors [1] (hereby referred as MLFM). Both models use mixture fraction to describe the reaction process. Chemical reaction is assumed to take place only at the stoichiometric surface where the mixture fraction reached the stoichiometric value xsf ¼
Y O2 , rY f1 þ Y O2
(5)
where Yf1and Y O2 are the mass fraction of fuel and oxidant, respectively. r is the stoichiometric coefficient. The thin flame located at the stoichiometric surface is treated as one-dimensional. The whole flame is treated as an assemble of those thin flames (usually referred to as ‘‘flamelet’’). The evolution of mixture fraction was controlled by the governing equation generated from the combination of its definition and the governing equation of species concentration. 2.3.1. The mixture fraction model The mixture fraction that represents the chemical reaction is solved by the governing equation [15,19]. Dx ¼ rrDrx. (6) Dt Based on the assumption of ‘‘infinitely fast reaction’’, the oxygen and fuel cannot co-exist at any local point. This leads to the ‘‘state relation’’ between the oxygen mass fraction YO and the mixture fraction: x 1 Y O ðxÞ ¼ Y O 1 (7) when xoxsf , xsf r
Y O ðxÞ ¼ 0
when xoxsf .
(8)
129
A look-up table is pre-constructed containing the species concentrations of oxygen, fuel and other products as functions of mixture fraction. Once the mixture fraction at each location is calculated in LES, the corresponding species concentration is obtainable by referring back to the look-up table. The heat release rate is calculated from the mass consumption rate and the heat of combustion of the oxygen _ O. q_ 000 ¼ DH O m
(9)
Combining with the radiant heat loss, the heat release rate is then put back to the conservation equation of energy (Eq. (3)) to calculate the temperature and other flow variables. 2.3.2. The modified laminar flamelet model Following Cook and Riley [22], the localized value of mixture fraction is obtained by solving the following governing equation in LES: qrx qðrxuÞ 1 q qx þ ¼ m (10) þ Bx . qt qxi ReSc qxi qxi The scalar dissipation rate can be calculated either by its governing equation or modelled by mixture fraction. Some simplification was introduced to the governing equation of the mixture fraction to further increase computational efficiency. x, x2 , w are chosen as inputs. x is obtained through its transport equation (Eq. (10)), in which zx represents the sub-grid-scale quantities zx ¼
q ðrxu rxuÞ. qxi
(11)
The species concentration of fuel is then calculated by solving the following second-order differential equation: wðx; tÞ
d2 Y f ¼ w_ f , dx2
(12)
where w(x, t) is the scalar dissipation rate, w_ f is the reaction rate. In analogue to SGS turbulence modelling, we re-write Eq. (10) as follows. Here the modelled turbulent viscosity mt is used to replace m and hence eliminates the need to include the SGS term zx qrx qðrxuÞ 1 q qx þ ¼ mt . (13) qt qxi ReSc qxi qxi x2 can be obtained either by solving a transport equation or by modelling. In the present study, the following model of Kops et al. [23] was used to get x2 : 2
x2 ¼ x þ C x rD2 rx rx,
(14)
where Cx is a model constant and D is the filtered length in LES.
ARTICLE IN PRESS J.X. Wen et al. / Fire Safety Journal 42 (2007) 127–138
130
w was modelled using the following expression of Gao et al. [24] : 2 m m qx þ t w¼ . (15) Pe Sct qxi
simple approach through some band-mean absorption coefficients.
As both x2 and w can be obtained from x, the original 3D look-up table in the flamelet calculation could be mapped to a 1D look-up table. The species concentration can be located by the mixture fraction only. This simplification brings significant advantages in computational efficiency. Further details about the implementation of this model into FDS can be found in Kang and Wen [1]. Although both SGS combustion models have used mixture fraction to describe the combustion process in fires, there are some differences between them.
The methanol pool fire tested by Weckman and Strong [26] was considered. Methanol was supplied at a constant rate to give a total heat release rate of 24.6 kW on a circular burner with a diameter of 30.5 cm. The fire was set in the open air at room temperature.
(1) The first difference is in the governing equation of mixture fraction. In the mixture fraction model, the mixture fraction is governed by Eq. (8), which has been filtered to ensure that all the items left in the equation are large-scale quantities. The quantities with a scale smaller than the filter width are simply filtered out. In the modified laminar flamelet model, the compensation of energy from the small scales has been taken into consideration in the governing equation of mixture fraction (Eq. (10)). (2) The second difference is the way of obtaining species concentration from the mixture fraction. In FDS, the mass fraction of oxygen was calculated from the socalled ‘‘ideal state relation’’, which was generated from the definition of mixture fraction. The corresponding mass fraction of fuel and products were obtained from the oxygen mass fraction. While in the MLFM, the mass fraction of fuel was calculated from Eq. (12). (3) In the mixture fraction model, the heat release rate is calculated from the mass consumption rate and the heat of combustion of the oxygen. In MLFM, the calculation of heat release rate is based on both the reactants and the products, (4) X ðrDÞ ¯ iþð1=2Þ;j;k ðxiþ1;j;k xi;j;k =dzÞ dY n 000 , DH n q_ ijk ¼ dx dy dz dx xox n f
(16) where n ¼ f, o and p, respectively.
2.4. Radiation and soot The default radiation model in FDS is used. It solves the finite-volume-based radiation transport equations [25] and considers soot as the most important combustion product controlling the thermal radiation from the fire and hot smoke. The contribution of soot is accounted for using a
3. Experiment considered
4. Computational domain and grid sensitivity study The fire is basically modelled as the ejection of pyrolyzed fuel from a liquid surface that burns when mixed with oxygen. The stoichiometry of the reaction is set to identify the fuel type. All species associated with the combustion process are accounted for by way of the mixture fraction variable. The heat of vaporization is specified so the burning rate of the fuel is dependent on the heat feedback from the fire. A rectangular computational domain was set with dimensions of 1.6 m(W) 1.6 m(D) 3.2 m(H) in a Cartesian co-ordinate. The ‘‘stair stepping’’ method is used to construct the grid around the circular burner. To lessen the impact of stair stepping on the flow field near the wall. FDS has built-in facility to prevent vorticity from being generated at sharp corners and smooth out the jagged steps that make up the obstruction [20]. Except the ground, all other boundaries were set open. The surrounding air was still at the beginning of the simulation. Such a computational domain was chosen from a series of sensitivity studies which found that narrow width/depth (less than 4 times of burner diameter) limited the development of eddies and gave incomplete description of air entrainment. The use of large domain (more than 8 times burner diameter in width/depth), however, would necessitate the use of coarser grids due to the limit of the computational resource. During the analysis of the results, no visible eddies were observed at the edges of the computational domain implying that the use of pressure boundary there is reasonable. Around two million cells (128 128 144) were employed in the simulation giving an average mesh resolution of 1.3 cm. Cells in the centre of the fire were further stretched to give a finer resolution of 1 cm. Fig. 1 shows the comparison of the predicted mean temperature distributions using three different grid resolutions, i.e., the coarse grid of 64 64 96, the medium gird of 108 108 128 and the find grid of 128 128 144. The predictions with the coarse grid were found to have the largest discrepancies with the experimental data. The predictions with the fine grid, as employed in the present study, were closest to the data while the predictions with the medium grid were found to be reasonably close to that of the fine grid and the data. Hence the fine grid was chosen for the present study.
ARTICLE IN PRESS Mean temperature at z = 6 cm 1.60E+03 1.40E+03 1.20E+03 1.00E+03 8.00E+02 6.00E+02 4.00E+02 2.00E+02 0.00E+00 0.00E+00 5.00E-02 1.00E-01 1.50E-01 2.00E-01 Radius (m)
Modified laminar
108×108×128
1600
64×64×96
1400
Weckman [29]
1200 1000 800 600 400 0 0.00
1600
128×128×144
1400
108×108×128 Weckman [29]
1000 800 600 400 200 0 0.00E+00
5.00E-02
1.00E-01 1.50E-01 Radius (m)
2.00E-01
0.03
0.05
0.08 0.11 Radius (m)
0.13
0.16
Mean temperature at z = 12 cm
64×64×96
Temperature (K)
Temperature (K)
flamelet model model Weckman [98]
200
Mean temperature at z = 12 cm
1200
131
Mean temperature at z = 6 cm
128×128×144
Temperature (k)
Temperature (K)
J.X. Wen et al. / Fire Safety Journal 42 (2007) 127–138
1400
Modified laminar
1200
flamelet model model Weckman [98]
1000 800 600 400 200 0 0.00
Mean temperature at z = 20 cm
1200
0.03
128×128×144
0.05
0.08 0.11 Radius (m)
0.13
0.16
108×108×128
Mean temperature at z = 20 cm
64×64×96 Weckman [29]
800
Modified laminar
1200
flamelet model model Weckman [98]
1000 600
Temperature (k)
Temperature (K)
1000
400 200 0 0
0.05
0.1 Radius (m)
0.15
800 600 400 200
0.2 0 0
0.05
Fig. 1. Grid sensitivity study.
0.1 Radius (m)
0.15
0.2
Mean Temperature at z = 30 cm Modified laminar
1200 Temperature (K)
Using the correlation of Malalasekera [6], the pulsation period of the fire was found to be around 0.32 s. The computational time was set to 20 s to give enough time for the fire to reach the fully developed stage. During the simulation, the time step was dynamically adjusted by the instantaneous velocities (DTomin(Dx/u, Dy/v, Dz/w)). The instantaneous results were monitored regularly. The fire was assumed to have reached quasi-steady state when the instantaneous values appeared periodically. All the time-averaging processes were carried out within the converged period over 10 puffing cycles.
flamelet model model Weckman [98]
1000 800 600 400 200 0 0.00
0.03
0.05 0.08 0.11 Radius (cm)
0.13
0.16
Fig. 2. Temperature distributions at different heights.
5. Results and discussion 5.1. Mean quantities 5.1.1. Temperature and velocity distributions The mean temperature at different heights (6,12, 20 and 30 cm) and the centreline temperature distributions are
plotted in Figs. 2 and 3. The predictions with the two different combustion models are in reasonably good agreement with each other and with the experimental data. The distribution of axial velocities is presented in Fig. 4. Both set of predictions have captured the correct trends of axial velocity distributions. The magnitude of axial velocity
ARTICLE IN PRESS J.X. Wen et al. / Fire Safety Journal 42 (2007) 127–138
132
Central line temperature 1400
Modified laminar flamelet model mixture fraction model Weckman [98]
1.2 1
1200
Velocity (m/s)
Temperature (k)
Axial velocity at z = 6 cm
Modified laminar flamelet model model Weckman [98]
1600
1000 800 600 400
0.8 0.6 0.4 0.2
200
0
0 0.01
0.02
0.03
0.04 Z /Q
0.05
0.06
0.07
0.08
0
0.05
(2 /5)
0.15
0.2
Axial velocity at z = 12 cm
Velocity (m/s)
Fig. 3. Centerline temperature.
2 1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0
Modified laminar flamelet model mixture fraction model Weckman [98]
0
0.05
0.1 0.15 Radius (m)
0.2
Axial velocity at z = 20 cm Modified laminar flamelet model mixture fraction model Weckman [98]
3.50
velocity (m /s)
3.00 2.50 2.00 1.50 1.00 0.50 0.00 0.00
0.05
0.10 Radius (m)
0.15
0.20
Axial velocity at z = 30 cm 3.5
Modified laminar flamelet model mixture fraction model Weckman [98]
3 velocity (9m /s)
was determined by vortex movement due to buoyancy produced by heat release. At lower heights, the flow was dominated by the chemical reaction. Eddies were generated at the burner rim then entrained towards the fire centre. Driven by the combined effect of molecular mixing and buoyancy, the vortices moved upwards while evolving, resulting in a relatively high axial velocity in the centre. The axial velocity decreased along the radius as the buoyancy became weaker at larger radius. In the region above the reaction zone, where the majority of the flow was in the fully developed turbulent regime, the intensity of turbulent mixing and the buoyancy effect caused a rapid rise of the axial velocity. At lower heights, where the flow was believed to be in the transition region, the mixture fraction model under-predicted the axial velocity while in the reaction zone and in fully developed turbulent region, the mixture fraction model over-predicted the velocity in the centre of the fire and underestimated at the outside region. Overall for all the heights listed in Fig. 4, the two sets of predictions are in reasonably good agreement. The predictions of the MLFM are only marginally closer to the data. The lined contours of the radial velocity are shown in Fig. 5. The distribution patterns of radial velocities from both set of predictions were similar to the experimental data. The highest radial velocity occurred near the burner rim at a radial position of 13 cm, where the air entrainment rate was expected to reach the maximum. The radial velocities in the centre of the fire were also high due to chemical reaction. The magnitude of radial velocity decreased with the height and the peak value moved towards the centre at the same time. This was thought to be the result of the completion of combustion and the development of vortices. Although the predictions of the MLFM were found to be slightly closer to the experimental measurement, the magnitudes of both sets of predicted radial velocities were much larger than the measured values. While the exact cause is not clear, it is thought that this may be linked to some simplifications in the momentum equation where the hydrostatic pressure gradient was subtracted and the baroclinic torque due to the non-alignment of density and pressure gradients was
0.1 Radius (m)
2.5 2 1.5 1 0.5 0 0
0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 Radius (m)
Fig. 4. Distributions of axial velocity at different heights.
neglected. Those simplifications were believed to have more impact on the predictions of radial velocities than the axial velocities. The simplified pressure treatment using the Possion equation may also have some influence.
ARTICLE IN PRESS J.X. Wen et al. / Fire Safety Journal 42 (2007) 127–138 Radial velocity contour
133
Radial velocity contour
Modified laminar flamalet model mixture fraction model
0.30
0.30
-0.2
0.28
0.28
0.26
-0.3
0.0
0.22 -0.1
0.20
-0.3
-0.1
Height (m)
Height (m)
-0.1
-0.5 -0.4
0.16 0.14 0.1
0.12
0.1
0.18
0.0
0.16 -0.4
0.14 -0.5
0.12
0.10
0.2
0.2
0.10 0.3 0.2
-0.1
-0.6
0.08
-0.1
-0.4
0.06
0.06
-0.5
0.02
-0.1
0.24
-0.1
0.20
0.04
0.3
0.2
0.22
0.08
-0.2
0.26
0.24
0.18
-0.1
0.0
0.04
-0.6
0.4 0.5 -0.3 0.0
0.3 0.50.4 0.70.6
-0.6
0.02
-0.2 -0.1
-0.4
-0.5 -0.3
-0.6 -0.7 -0.2
-0.1 0.00 -0.04 -0.02 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16
0.00 -0.04 -0.02 0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 Radius (m)
Radius (m) Experimental measurement (Weckman [98])
30 -0
02
-0.
2
2
.0
-0.0
28 26 24 22
-0.02 -0.05
-0
0
18
.1
.0
5
-0
0.05
16 15
-0.
-0.02
height, z [cm]
20
14 12 0 .1
-0
10
-0.20
0.05
-0
8
5 10 -0.
5 .2
-0.0
6 4 2
-4
-2
0
2
4 C
L
6
8
radius, r [cm]
10
12
14
16
nm
Fig. 5. Contours of the predicted and measured radial velocities.
5.1.2. Air entrainment Following Zhou and Gore [27] and also the above observation that the highest radial velocity occurred near the burner rim, we took the radial velocity at the radial position of 15 cm (radius of the burner) as the air entrainment velocity and plotted this in Fig. 6. Both SGS combustion models gave similar distribution of air entrainment velocities. However, the mixture fraction
model under-predicts the air entrainment near the burner base and over-predicts it at the top. In Fig. 7, the predicted mean air entrainment rates were compared with the correlation of Delichatsios [28], the data fit of Delichatsios and Orloff [29] and the experimental data of Beyler [30]. The related radial velocity was also taken at a radial position of 15 cm. According to Delichatsios [28], air entrainment in pool fires can be
ARTICLE IN PRESS J.X. Wen et al. / Fire Safety Journal 42 (2007) 127–138
134
Air entrainment velocity
Modified laminar flamelet model mixture fraction model
0.2
Velocity (m /s)
0.1 0 -0.1 -0.2 -0.3 -0.4 0. 00 0. 07 0. 13 0. 20 0. 27 0. 33 0. 40 0. 47 0. 53 0. 60 0. 67 0. 73 0. 80 0. 87
-0.5 Height (m) Fig. 6. Air entrainment velocity against height.
3 Temperature
2
3
200
20 1
100
10
0 0
20
(a) 1 Isothermal layer growth development rate (mm/min) 2 Thickness (cm)
300 Temperature (°C)
30
0
(b)
400
1 Isothermal layer growth development rate 2 Thickness
40 Time (min)
60
80
Run #2 300
40 3
30 200 2
1 Isothermal layer growth development rate 2 Thickness 3 Temperature
20
100
10
Temperature (°C)
1 Isothermal layer growth development rate (mm/min) 2 Thickness (cm)
Run #1 40
1
0
0 0
20
40 Time (min)
60
80
Fig. 7. Comparison of the predicted air entrainment rates with published correlation.
divided into three different regions. The parameters that could affect the air entrainment behave differently in each region. Delichatsios suggested that the dimensionless air entrainment rate meant (normalized by the fuel supply rate mf) up to a vertical distance Z above the fuel source was determined by the stoichiometric number s, the Froude
number Frf and the normalized height Z/D. He proposed that ððment =ðs þ 1Þmf ÞÞFrf ðZ=DÞ. Since the air entrainment rate was proportional to r1 AU ent , Delichatsios [28] further argued that the correlation between air entrainment rate and the height varied at different regions. Near the burner base, where AD2, the air entrainment rate was proportional to (Z/D)1/2. At the ‘‘neck-in’’ area, where most reaction took place, the air entrainment rate was proportional to (Z/D)3/2. In this region, air entrainment rate was strongly dependent on height. The air entrainment rate increased quickly with the height and most of the fuel was burned here, resulting in the so-called ‘‘neck-in’’ area. At the upper part of the fire, where the air entrainment was dominated by the turbulent mixing and had little to do with burner diameter, the air entrainment rate was proportional to (Z/D)5/2. In this region, less fuel was left and height-related buoyancy became the dominant factor for air entrainment. Both sets of predictions behaved similarly while in lower region the predictions with the MLFM were slightly closer to the correlation. In general, the FDS predictions agreed well with Delichatsios’s analysis for all three regions. 5.1.3. Heat release rate per unit volume (HRRPUV) The predicted values of HRRPUV are shown in Fig. 8. In the prediction with the MLFM as shown in Fig. 8(a), the highest heat release rate per unit volume was found in the outer region of the flame near the pool surface, where the mixture contained abundant fuel gas/vaporized fuel. As revealed in the current prediction as well as that of Zhou et al. [1], the oxygen concentration was also higher there because of the strong inflow of air towards the centre near the base induced by the sharp pressure gradient between the flame and the surroundings. Near the base, the heat release rates in the centre were relatively low due to lack of oxygen. Further up, with more fresh air being entrained by the rising flame and brought into the centre by the evolving vortices, the heat release rate reached the highest in the fire centre. Consistent with previous findings, the combustion intensities were strongest in the ‘‘neck-in’’ region around half a diameter above the pool surface. Following that, combustion continued but the high heat release regions were becoming gradually smaller before the combustion processes almost fully completed at a height of approximately 1.8 m, i.e., 6 diameters above the pool surface. Along the radial direction apart from immediately above the pool, the intensity of combustion decreased gradually from centre outwards, leading to the reduction of volumetric heat release rates as well as mean temperature and axial velocity (Figs. 2 and 4). 5.2. The dynamic behaviour Four sequential plots of instantaneous temperature contours are shown in Fig. 9. They demonstrate the dynamic behaviour of the fire after it reached the quasistable state. The time interval between two adjacent plots is
ARTICLE IN PRESS J.X. Wen et al. / Fire Safety Journal 42 (2007) 127–138
135
2.5
2.5
2.0
2.0 4500
4500 4500
1.5
1.5 Height (m)
Height (m)
4500 4500 4500 4500
5000 5000
1.0
1.0
5000 5500
4500
4500
6000
4500
6500
0.5 7500
4500 5000
7000
7500
4000 1003000 0 2000 5000 3000 4000 5000 6000
0.5
5000 6500 5500
2000 100 0
5500
0.0 -0.32 -0.24
6500 6000 800850 900 1000 1050 0000000 950 7000 4500
6500
6000
9000 10000 8000 60007000 4500
7000 7500900 000 9500 950 5508000 0 850 6000 5000
-0.16 -0.08 0.00 0.08 Radius (m)
0 0004000 1200 1000 0 7000 8000 9000 03000 600
0.0 0.16
0.24
0.32
-0.32
-0.24
-0.16
5000 5000
4000
3000 5000 6000 4000 10000 7000 8000 9000 2000 100 0 3000
-0.08 0.00 0.08 Radius (m)
0.16
0.24
0.32
Fig. 8. Contours of heat release rate per unit volume.
0.08 s, which was determined by the puffing frequency. As calculated from Malalasekera’s correlation [6], the puffing frequency of the methanol fire considered was 2.56 Hz, leading to a puffing period of 0.32 s. To give four sets of contours within one pulsation period, 0.08 s was chosen as an interval. As shown in the temperature contours in Fig. 9, hot eddies formed at the edge near the base of the burner. Small eddies rolled into the centre of the fire with the air entrainment and rose upwards, continuing to evolve during the process. At the height of approximately one pool diameter, they started to separate from the plume and gradually vanished in the turbulent fire plume. The wellknown phenomena of pool fire ‘‘neck-in’’, where the air entrainment velocities reached peak values and the visible flame there shrank as a neck, was identifiable and the fluctuation of the ‘‘neck-in’’ height from about half to one pool diameter above the base can also be seen. Most of the burning completed before the separation point at a height between one to two pool diameters and just after the fire had ‘‘necked-in’’. Some isolated pockets of hot eddies can still be seen above that. This is believed to be linked to the existence of fuel-rich pockets. In the mean time, new eddies started to form near the base. The periodical motion of hot eddy generation and shedding is thought to be the dominant cause of the fire pulsation. The above findings are consistent with previously reported experimental
observations [5,31] and numerical predictions of Ghoniem et al. [31] using a vorticity-based approach. To further demonstrate the ‘‘neck-in’’ area, four line plots of the instantaneous air entrainment velocities as a function of height are presented in Fig. 10, representing four different moments within a pulsation period. The velocities were taken at the edge of fire (radial position ¼ 15 cm). The negative values of velocity represent air being entrained into the fire centre. The maximum air entrainment velocity occurred near the base of the burner, with the second peak entrainment value occurring between half and one diameter above the base. At the same height, the entrainment velocity changes periodically. The results suggest that the ‘‘neck-in’’ is linked with the shedding of large flame puffs and small vortices near the pool surface. Previous studies suggested that most combustion reaction took place near and within the ‘‘neck-in’’ area [21,27]. The spectra of temperature fluctuations are plotted in Fig. 11. They failed to show clearly any dominant frequency at different heights, implying that the code cannot predict the puffing frequency. It is thought that this limitation might result from the simplified pressure treatment using a Poisson equation. Some separate numerical tests with an in-house LES code without such simplification have demonstrated that it is possible to capture the puffing frequency with the LES approach. However, this code runs much more slowly than FDS and
ARTICLE IN PRESS 136
J.X. Wen et al. / Fire Safety Journal 42 (2007) 127–138
Fig. 9. Instantaneous temperature contours.
hence it is not suitable for routine fire simulations in engineering applications. In Fig. 12 where the frequency of temperature crossing a range of reference values is plotted, the experimental measurements of Cox and Chitty [32] for a methane fire on a 30 cm square burner are also shown. This data set was chosen partly because of its availability and partly because the adiabatic flame temperatures of methane and methanol are very close. It is seen that the predictions are consistent with the measurements. The peak frequency of crossing 1000 1C occurs at Z0 ¼ 0.06, where is believed to be in the
reaction zone. The peak frequency of crossing 500 and 750 1C occurs within the intermittent region (0.08oZ0 o0.2) while the peak frequency for crossing 250 1C occurs in the fully developed plume region (Z0 40.2). 6. Concluding remarks This study examines in detail the predictions using FDS for a medium-scale pool fire. The simulations were carried out predominantly using the existing features in FDS
ARTICLE IN PRESS J.X. Wen et al. / Fire Safety Journal 42 (2007) 127–138
Air entrainment
0
0 Velocity (m /s)
-0.05 -0.1 -0.15 -0.2 -0.25
0 0.03 0.05 0.08 0.1 0.13 0.15 0.18 0.2 0.23 0.25 0.28 0.3
-0.05 -0.1 -0.15 -0.2 -0.25
-0.3
-0.3 Height (m)
Height (m)
Air entrainment
Air entrainment
0
3 0.
0. 2 22 5 0. 25 0. 27 5 0.
1 0.
12 5 0. 15 0. 17 5
0
0.
02 5 0. 05 0. 07 5
0. 3
0. 2 0. 22 5 0. 25 0. 27 5
5 07
0.
5
05 0.
02
0. 1 0. 12 5 0. 15 0. 17 5
-0.1 Velocity (m/s)
Veloicty (m /s)
-0.05
0.
0
0 0.
Velocity (m /s)
Air entrainment
0.05
0 0.03 0.05 0.08 0.1 0.13 0.15 0.18 0.2 0.23 0.25 0.28 0.3
137
-0.1 -0.15 -0.2 -0.25 -0.3
-0.2 -0.3 -0.4 -0.5
Height (m)
Height (m)
Fig. 10. Instantaneous air entrainment velocity.
Temperature Spectra on the Plume Axis
20 18 16 Frequency (Hz)
2 Power spectral density (C /Hz)
1e+5
1e+4
1e+3
14
250°C
500°C
750°C
1000°C
250°C [143]
500°C [143]
750°C [143]
1000°C [143]
12 10 8 6 4
1e+2
2
1e+1
-2
0 0.01 0.023 0.057 0.08 0.1 0.13 0.17 Normalised height z'
1e+0 No 0.25 rm 0.20 0.15 (m alise / K d h 0.10 0.05 w2 e / 5 ight 0.00 ) z'
0.2
0.3
Fig. 12. Frequency of temperature crossing a range of reference values. 100 10 1
y (Hz)
c Frequen
Fig. 11. Temperature spectra on the plume axis.
except that an additional sub-grid-scale combustion model based on the laminar flamelet approach of Cook and Riley [22] was used for comparison. The results have shown that FDS with its existing features can deliver accurate predictions for most important parameters of pool fires that are of significance in the fire safety context. These parameters include mean values of temperature and axial velocity distributions, as well as the air entrainment ratios. Although the results have shown the limitation of the code in predicting the puffing frequency, this is thought to be a
consequence of some approximation in FDS which needs to be there in order to obtain the high efficiency of the FFT-based fast direct solver for the Poisson equation for total pressure. Nevertheless, the predictions have captured the pulsating nature of air entrainment and the well-known ‘‘neck-in’’ at a height of approximately half to one diameter above the pool surface. The overall agreement between the two sets of predictions using different combustion models should help to demonstrate the robustness of the existing mixture fraction combustion model for medium-scale pool fire predictions. Further work would be desirable to explore the issue of radial velocity predictions and soot modelling in FDS. In the later aspect, it is encouraging to note the recent work of Lautenberger et al. [33].
ARTICLE IN PRESS J.X. Wen et al. / Fire Safety Journal 42 (2007) 127–138
138
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11]
[12] [13] [14] [15] [16] [17] [18]
Kang Y, Wen JX. Combust Sci Technol 2004;176(12):2193–223. Gebhart B. Annu Rev Fluid Mech 1973;5(213). Cetegen BM. Phys Fluids 1997;9(12):3742–52. Weckman EJ, Sobiesiak A. Proc Combust Inst 1988;22:1299–310. Cetegen BM, Ahmed TA. Combust Flame 1993;93:157–84. Malalasekera WMG, Versteeg HK, Gilchrist K. Fire Mater 1996;20: 261–71. Hamins A, Yang JC, Kashiwagi T. Proc Combust Inst 1992;24: 1695–702. Mandin O, Most JM. In: Proceedings of the sixth international symposium on fire safety science, 1999. p. 1137–1148. Pagni PJ. Appl Mech Rev 1990;43:153–70. Joulain P. Proc Combust Inst 1998:272691–706. Blunsdon CA, Beeri Z, Malalasekera WMG, Dent JC. HTD, vol. 296, Fire, combustion and hazardous waste processing, New York: ASME; 1994. Sinai YL. Fire Safety J 2000;35(1):51–62. Liu F, Wen JX. Fire Safety J 2002;37(2):125–51. Ghoniem AF, Lakkis I, Steriou M. Proc Combust Inst 1996;26: 1531–9. Mell WE, Mcgrattan KB, Baum HR. Proc Combust Inst 1996;26: 1523–30. Jiang X, Luo KH. Proc Combust Inst 2000;28:1989–96. Baum HA, Mcgrattan KB, Rehm RG. In: Proceedings of the fifth international symposium on fire safety science, 1997. p. 511–522. Smagorinsky J. Monthly Weather Rev 1963;91:99–165.
[19] McGrattan KB, Forney GP, Floyd JE. Fire dynamics simulator— technical guide (version 2). National Institute of Standards and Technology, USA, June 2001. [20] McGrattan KB. Fire dynamics simulator—technical guide (version 4). National Institute of Standards and Technology, USA, July 2004. [21] Zhang W, Hamer A, Klassen M, Carpenter D, Roby R. Turbulence statistics in a fire room model by large eddy simulation. Fire Safety J 2002;37(8):721–52. [22] Cook AW, Riley JJ. Combust Flame 1998;112:593–606. [23] Kops DB, Riley SM, Kosaly G, Cook AW. Flow, Turbulence Combust 1998;60:105–22. [24] Kang MY, Wen JX, McGrattan KB, Baum HR. Interflam 2001:743–54. [25] Huggett C. Fire Mater 1980;4(2):61–5. [26] Gao F, O’Brien EE. A large eddy simulation scheme for turbulent reacting flows. Phys Fluids 1993;5(6):1282–4. [27] Cetegen BM. Phys. Fluids 1997;9(12):3842–52. [28] Delichatsios MA. Combust Flame 1987;70:33–46. [29] Delichatsios MA, Orloff L. Proc Combust Inst 1984;20:367–75. [30] Beyler CL. PhD thesis, Harvard University, 1983. [31] Cetegen BM, Zukoski EE, Kubota T. Proc Combust Inst 1984; 39:305. [32] Cox G, Chitty R. Fire Mater 1982;6:127–34. [33] Lautenberger CW, de Ris JL, Dembsey NA, Barnett JR, Baum HR. A simplified model for soot formation and oxidation in CFD simulation of non-premixed hydrocarbon flames. Fire Safety J 2005; 40(22):141–76.