Numerical Studies on Heat Release Rate in Room Fire on Liquid Fuel under Different Ventilation Factors

Heat release rate (HRR) of the design fire is the most important parameter in assessing building fire hazards. However, HRR in room fire was only studied by computational fluid dynamics (CFD) in most of the projects determining fire safety provisions by performance-based design. In contrast to ten years ago, officers in the Far East are now having better knowledge of CFD. Two common questions are raised on CFD-predicted results on describing free boundaries; and on computing grid size. In this work, predicting HRR by the CFD model was justified with experimental room pool fire data reported earlier. The software fire dynamics simulator (FDS) version 5 was selected as the CFD simulation tool. Prescribed input heating rate based on the experimental results was used with the liquid fuel model in FDS. Five different free boundary conditions were investigated to predict HRR. Grid sensitivity study was carried out using one stretched mesh and multiple uniform meshes with different grid sizes. As it is difficult to have the entire set of CFD predicted results agreed with experiments, macroscopic flow parameters on the mass flow rate through door opening predicted by CFD were also justified by another four conditions with different ventilation factors.


Introduction
There are many big construction projects [1] with difficulties to comply with the fire codes while developing in the Far East. Performance-based design (PBD) was then applied to determine fire safety provisions. However, budget allocated on fire safety used to be small, resulting in a lack of full-scale burning tests to measure heat release rate. There are even no resources allocated in some places like Hong Kong to compile a database of heat release rate for local products including train compartments of the railway system [2] as in Japan, China, and Korea. Consequently, wrong concept was adopted in estimating the heat release rate, taking "average value" as "peak value" in many projects submitted to be evaluated by the principal author [1]. At most, computational fluid dynamics (CFD) was applied [3,4] in hazard assessment. Fire behavior involves complex dynamics driven by critical events, such as the ignition of secondary items, flashover, window breakage, and falling down of glass systems. All these phenomena have not yet been modeled realistically without using empirical parameters. Authorities having jurisdictions (AHJ) are now more knowledgeable in fire science and engineering. Taking Hong Kong as an example, a huge percentage of senior officers approving fire projects are well trained and have a master degree. CFD-predicted results are, therefore, evaluated more in depth. Two questions [5,6] on free boundaries and grid size are commonly raised in using CFD. The hazard assessment is normally rejected, and this is very different from fifteen years ago in accepting highcost projects but the fire engineering reports had no indepth experimental justification released to the public [7][8][9]. Therefore, very tight inspection scheme was implemented or proposed to implement in the Far East on those existing projects accepted 20 years ago. Most of the new project applications based on CFD are only accepted by AHJ for designing smoke control design in big and tall buildings [10]. Even so, field tests with hot smoke are required to justify the CFD predictions.
In fire safety design, the most important parameter is the heat release rate (HRR), which is the single most important 2 International Journal of Chemical Engineering variable in characterizing the "flammability" of products and their consequent fire hazard. It gives information on fire size, fire growth rate, available egress time, and suppression system impact [11]. The potential for ignition of nearby items, flashover potential in a room, and the rate of water needed to extinguish the fire can be estimated [12]. The evolution of HRR with time becomes the most important input variable [13] which must be estimated properly to fire simulations. However, HRR used to be predicted by CFD without experimental justification.
The CFD software fire dynamics simulator (FDS) was selected to study fire driven fluid flow. It was developed [4] by National Institute of Standards and Technology (NIST), used in solving practical fire problems. The Navier-Stokes equations were derived to study low-speed, thermally driven flow with an emphasis on smoke and heat transport from fires. Large eddy simulation (LES) was used in simulating turbulence. The FDS version 5 was applied in this study.
These two questions on predicting HRR by CFD will be justified in this paper. HRR is predicted by FDS under different free boundary conditions and grid sensitivities. Results measured from full-scale burning test in gasoline fuel in a room calorimeter [14,15] were used. Five boundary conditions were used to evaluate free open boundary conditions in applying CFD to predict HRR. Effect of grid size was investigated using one stretched mesh and multiple uniform meshes with different grid sizes.
As raised by Chen [16], it is difficult to have the whole set of CFD-predicted results agreed with experiments. However, macroscopic flow parameters predicted by CFD are very useful [17]. Four other fire scenarios with different ventilation conditions were investigated in this paper with empirical equations on HRR with ventilation factor. Four ventilation factors were used to predict HRR.

Earlier Experimental Results
Experimental data with heat release rates measured in a room calorimeter at Lanxi, Harbin, Heilongjiang, China [14,15] were taken used in justifying CFD predictions. The room was constructed of brick with a cement finish of length 3.6 m, width 2.4 m, and height 2.4 m as shown in Figure 1. The thickness of the brick wall was 0.25 m, with the cement finish of thickness 20 mm. The ceiling was constructed with 0.2 m reinforced concrete. There was a door of height 2.0 m and width 0.8 m.
Three thermocouple trees, each with six thermocouples were placed at the center of the room, near the wall, and the center of the door as shown in Figure 1. The uppermost thermocouple was placed below the ceiling and the other five thermocouples were placed at 0.4 m intervals. A heat flux meter was placed 0.6 m from the door to measure the radiactive heat flux received at the floor level. The ambient temperature was 15 • C.
A gasoline pool fire of 0.46 m diameter was placed at the center of the room. The HRR was measured by oxygen consumption calorimeter. The measured HRR shown in Figure 2 is used to carry out grid sensitivity study in this paper.

Pyrolysis Model for Liquid Fuel
Pyrolysis of solids and liquids can be predicted in FDS by different models [4] depending on the availability of material properties. Evaporation rate of liquid fuel is determined by the Clausius-Clapeyron equation [18] through the liquid temperature and the concentration of fuel vapor above the pool surface. The volume fraction of the fuel vapor above the surface (X f ) is taken as a function of the liquid boiling temperature (T b ), the heat of vaporization (h v ) and the molecular weight (W f ) by: At the start of simulation, the fuel vapor mass flux m i was generated first by the initial vapor volume flux V i specified by the user: The evaporation mass flux will be updated in the simulation by inspecting the difference between the predicted close-tothe surface volume fraction of fuel vapor with the equilibrium value estimated by (1). For simplicity, the liquid fuel itself is treated as a thermally thick solid in simulating thermal conduction. Convection within the liquid pool is not considered. The fuel mass flux cannot be expected as an explicit function of temperature, the value can only be estimated by iteration. Temperature and flow conditions would affect the result. Furthermore, estimation of the evaporation rate depends strongly on grid size and distribution. The material properties cannot be traced with measurement.
In applying FDS to study smoke spread and heat transfer, the heat release rate can be taken as an input parameter by the user. The desired value of HRR is transformed [18] into a mass flux for fuel (m f ) at a given solid surface given in terms of the heat release rate per unit area (q) and a specified time ramp function f (t) by: The mass loss rate can then be computed according to the input q in FDS simulation. Typical values concerned are shown in the example calculation of the manual as listed in the following sections.

Free Open Boundary
There is an open boundary on the exterior boundary of the computational domain to let bi-directional flow with hot gas flowing out and cool air coming into the room through the openings. In FDS [18], the geometry as shown in Figure   and density ρ 1 to a point P ∞ far away from the outflow line at pressure P, velocity 0, and density ρ ∞ as: The pressure P is set to ambient pressure P ext , ρ ∞ is the ambient density.
x y X Ghost cell The Bernoulli theorem was applied to the inflow condition at an OPEN vent by taking the flow as inviscid, steady, and incompressible. Taking the point A 2 to A 3 of the inflow streamline, the fluid element on the boundary accelerated from point A 2 at pressure P 2 , density ρ 2 , and velocity u 2 to point A 3 at state pressure P 3 , density ρ 3 , and velocity u 3 :  The kinetic energy at A 2 is (1/2)ρ 2 u 2 2 with ambient pressure P = P ext . The fluid accelerates to u 3 at point A 3 which is on an inflow boundary.
The hydrodynamic pressure (head) H 3 at point 3 can be obtained by: Substituting (6) to (5) would give: Rearranging gives The density ρ 3 is taken as the average density between the gas-phase and ghost cells. Note that such ghost cells are assigned with the same size as their neighbor cells to establish gradients of various quantities at the boundary as shown in the figure.
Taking I as the index of the last gas phase cell in the x direction, u I, jk , v I, jk , and w I, jk are the x, y, and z components of velocity at the boundary, respectively. At open boundaries (say i = I), H depends on whether the flow is incoming or outgoing. Defining u at the Ith cell as: The density at the boundary wall is denoted by ρ w . The background density used in defining H is ρ ∞ . u is the velocity component at the boundary. The boundary condition on H is: International Journal of Chemical Engineering  The external dynamic pressure H ext is: The initial dynamic pressure H 0 is:  Both external and internal dynamic pressures are specified by the user (P ext is set with DYNAMIC PRESSURE for a specific VENT in FDS). The condition that (8) follows from prescribing the fluctuating hydrodynamic pressure P ext at an outflow vent (zero by default) and assuming the Bernoulli equation applies to an inflow vent where the fluid accelerates from the state {P ext ; ρ 0 ; u 0 ; v 0 ; w 0 } along a streamline. Note that the background density ρ ∞ is usually equal to the initial ambient density ρ 0 at the start of a calculation, but may change in time if the baroclinic correction term is included.
There are different views on specifying free open boundary conditions in applying CFD and these views should be watched carefully [19]. There have been some attempts to International Journal of Chemical Engineering solve the problem with free boundary. Markatos and coworkers [20,21] extended the flow domain to the "free boundary" region outside the doorway when studying the smoke flow in enclosures and obtained results that agreed reasonably with experimental data. Schaelin and co-workers [22] pointed out that extending the computing domains outside was a better approach when simulating plume flow. Galea and Markatos [23] pointed out in their case study on simulating fire development in an aircraft that it is desirable to extend the solution domain outside the fire compartment in order to find physically realistic behaviour in the vicinity of the open doors. Some pioneering work on fire modelling [19][20][21]23] demonstrated that the flow pattern in the vicinity of doorway was entirely different if the free boundary had not been extended sufficiently. However, in applying FDS, Hadjisophocleous and Ko [24] suggested that the impact of the open boundary at the exterior of the computational domain was minor when the boundary had been extended up to 2 m outside a geometry of width 10 m. Therefore, it may not be necessary to extend the computational domain to some distance beyond the opening to obtain good results while using FDS version 4.07. They also pointed out that

Numerical Simulations
The FDS version 5 was used to predict the HRR in a room fire. A three-dimensional Cartesian coordinate system was assumed with length along the x-direction, width along the y-direction, and height along the z-direction. Free boundary conditions were imposed on the outside part between the computational domain and the external environment. Fluid can enter or leave the computational domain freely. Pressure at the boundary was taken as the same as the ambient pressure.
The physical properties of the brick wall and cement were taken as thermal conductivity 1.0 Wm −1 K −1 , density 1700 kgm −3 , specific heat capacity 1000 Jkg −1 K −1 . The physical properties of the ceiling were taken as constant values [26] with a thermal conductivity of 1.6 Wm −1 K −1 , density of 2400 kgm −3 , and specific heat capacity of 920 Jkg −1 K −1 .
The molecular formula of gasoline is C 8 H 18 , complete combustion with oxygen gives water vapor (H 2 O) and carbon dioxide (CO 2 ) as: Input parameters in FDS are the boiling temperature, heat of reaction and specific heat, values are found to be 155 • C, 338.9 kJkg −1 and 2.22 kJkg −1 K −1 , respectively, in the literature [5].
The pool fire was simulated by a polygon with equally calculated area to the circular pan in the experiment. The volume of gasoline was 8 liters. The initial air temperature was 15 • C.
The Smagorinsky subgrid-scale turbulence model was used in large eddy Simulation (LES) of FDS. In this study, the value of the Smagorinsky coefficient Cs is 0.2, the value of the turbulent Prandtl number Prt is 0.5, and the value of the turbulent Schmidt number Sct is 0.5 as recommended [4].
A scenario SC1 with the same experimental condition with a door of width 0.8 m and height 2.0 m as shown in Figure 4 Table 1. Simulation results were compared with experimental data. Appropriate boundary conditions were then used for sensitivity analysis on grid. The total time for each simulation is taken to be 1500 s. The time step is determined by the Courant-Friedrichs-Lewy (CFL) condition to satisfy the stability criteria [4,18].
In order to quantify this comparison precisely, functional analysis proposed by Peacock et al. [27] on zone modeling was applied to evaluate the CFD results [6,28]. Transient predicted and measured data are expressed as vectors P and M. The Euclidean norm and secant inner product cosine between P and M are calculated: International Journal of Chemical Engineering Values of norm and cosine are used to compare CFDpredicted results with measured data. Values of norm should be 0, and cosine should be close to 1 for good agreement.

Open Boundary Conditions
Predicted HRR results for all five simulations OB1, OB2, OB3, OB4, and OB5 are compared with each other in Figure 9. The functional analysis results on norm and cosine are shown in Table 1. Small norm values of 0.03 indicate that the predicted HRR curves are almost identical in magnitude with the experiment curve. Cosine values lying between 0.95 to 0.97 indicate that the curve shapes are very similar. Velocity vectors, temperature contour, and the pressure contour at the central vertical plane of each condition at 600 s are shown from Figures 6, 7, and 8. The flow pattern of OB1 is very different from those of other four free boundary conditions. This suggests that extending the free Taking both computing accuracy and the computing time required into consideration, the boundary condition OB2 was selected for further simulation in this paper. Computing domain was extended to 7.2 m length, 4.8 m width and 4.8 m height in the Cartesian co-ordinate system to simulate air movement and distribution of combustible gases. The size of computing domain was twice of the chamber along the x, y and z directions respectively as shown in Figure 10.

Grid Size Variation
In LES, transient larger eddies are solved and smaller unresolvable eddies are modeled with a time averaged component and a fluctuating perturbation about that average [3]. The fineness of the numerical grid would determine the size of eddies that can be solved. The nominal size of the mesh cell δx and the characteristic fire diameter D * given by fire poweṙ Q, air density ρ ∞ , temperature T ∞ , specific heat of air C p , and gravitational acceleration g are important in simulating buoyant plumes [18]: The ratio D * /δx can be taken as the number of computational cells spanning the characteristic diameter of the fire. A refined grid system can improve the accuracy of results of LES. One stretched mesh (stretched in xand y-directions and uniform in z-direction) and multiple uniform meshes (finer mesh in the vicinity of the fire and door) as shown in Figure 11 were applied to study grid sensitivity in this paper. Different nominal grid size systems varying from relatively coarse mesh to fine mesh were tested to study the effect of different grid systems on predicting HRR. Systems T1, T2, T3 are results predicted under stretched mesh; and M1, M2, M3 are results predicted under stretched mesh with different grid sizes as shown in Table 2. Under each mesh system, the grid size was reduced gradually. The predicted HRR was compared with experimental data [14,15].
At first, the curve of heat release rates from the experiment was taken as the input function. The experimental HRR profile was compared with calculated HRR profiles for different grid sizes of stretched mesh as shown in Figure 12. As shown in Figure 12(a), there were significant changes when the grid system had changed from T1 to T2. Changing the grid system from T2 to T3 gives little change in the predicted HRR profile. However, a longer computation time of over 180 hours is required. The same conclusion can be applied to multiple mesh system file M2 as shown in Figure 12(b). Functional analysis results of the point-topoint comparison for grid size file T1, T2, T3, M1, M2, and M3 are presented in Table 2.
Grid system T1 give lower values than experiment as shown in Figure 12(a). Grid systems T2 and T3 have identical shape compared with experimental curve(s). Functional analysis results give a cosine of T1 (of 0.74). Values of cosine from T2 and T3 are 0.96 and 0.99, respectively. These values are very close to unity. Fining the system from T2 to T3 gives little changes in both curve shape and magnitude. Based on the functional analysis results, fining grid size than a certain value might not result in a better prediction.
Grid system M2 in the core area is the same as that of grid size file T2. Though the comparison, as shown in Figure 12, suggests that agreement would not be better for both the multiple meshes system and the stretched mesh system, the functional analysis, as shown in Table 2, confirms that grid sytem M2 agrees better with experiment. The results from finer grid file M3 for multiple meshes system agree well with experiment, as shown in Figure 12(b). Functional analysis in M3 confirms this. However, M3 required a 430hour computing time, which is much longer than M2. And the changes in results triggered by M3 are not significant.
Within the grid size file T1 and T2, HRR was simulated by the liquid fuel model in FDS. The fuel was allowed to burn by itself. Predicted HRR by both systems did not agree with experiment as shown in Figure 13. The fuel was burnt up in a much shorter period than measured experimentally. Functional analysis indicates that better results can be predicted by finer grid system T2. However, predicted curves of both systems did not agree with experiment. The burning rate in a room fire is strongly affected by radiative heat feedback which might not be well simulated in the FDS International Journal of Chemical Engineering   [34], and Thomas et al. [35]. Therefore, it is suggested that liquid fuel model should be applied carefully to engineering application. Further investigation to work on the accuracy of CFD prediction should be conducted. Therefore, the heat release rate of a fire is better taken as an input parameter by ramped heating rate based on the experimental results. This will be shown in the following simulations with different ventilation factors. As pointed out by the program developers [4], inherited limitations in the FDS liquid fuel model would generate problems.

Ventilation Factors
Different ventilation scenarios were investigated by adjusting the door height from sill to upper boundary. In addition to SC1 with door height 2.0 m from the floor, four different ventilation conditions SC2 to SC5 with a fixed door of width 0.8 m but different heights, as shown in Figure 4 were used in this study: In SC1, the four simulations SC2 to SC5 on different ventilation factors were conducted for 1500 s, using the grid size profile M2 and boundary condition OB2.
Predicted curve of the HRR in different scenarios by the FDS are compared with experiment in Figure 14. It is observed that, except SC4, all predicted HRR curves agree with the input HRR. SC4 has the lowest ventilation factor. Velocity vectors, temperature contour, and pressure contour of the central vertical plane at 600 s are shown in Figures 15,16 and 17. Prediction in oxygen supply with the decreased opening area would decrease the heat release rate.
The macroscopic flow parameter transient net air flow of mass throw the door predicted by FDS is used to justify CFD predictions. Four different ventilation factors V f defined in terms of the openings A 0 and the opening height H o are changed: An empirical equation was derived by Babrauskas and Williamson [36] to express the air flow rateṁ a through V f : Based on the predicted results by FDS, integrated mass flux through the door, and averaged (noun) over the computing period can give the intake airflow rateṁ a . Figure 18 shows the intake airflow rateṁ a for all ventilation scenarios. With the fire load located in the center of the room, increasing the size of the ventilation openings would increase the averaged intake airflow rateṁ a . Under the same ventilation condition, the averaged intake airflow rateṁ a would be increased by higher HRR value. This relationship can be described by linear correlation: This equation is very close to the empirical equation given by (17). It has to be pointed out that, for scenario SC5, though the ventilation area was similar to SC3, the amount of outtake air caused by convection was not as much as in scenario SC3 due to the different locations of opening. Buoyancy of the hot air might be a factor. Velocity vectors in Figures 15(c) and 15(e) show that ventilation rate increases if the opening is closer to the top of the room.

Conclusion
Heat release rate of a room pool fire was studied by FDS version 5 with different free open boundary conditions. A sensitivity study was carried out to compare predicted results with the experiment results under different grid systems. Predictions with the coarse grid have deviated more away from the experimental data. The predictions with a medium grid were found to agree with the experiment well, and selected for this study. In general, the prediction of HRR by FDS is good. Although the computational domain outside did not give much difference to the predicted HRR, predicted velocity patterns are very different.
Free opening boundary condition should be evaluated before being applied in CFD simulations, especially when the combustion process is included while simulating fires in tall or supertall buildings. Extending the computational domain to a sufficient distance beyond the opening is recommended.
The effects of ventilation factor on HRR and mass flow rate through the door opening were analyzed by FDS. On the five fire scenarios with different ventilation conditions investigated in this paper, air intake rate might be higher to give more oxygen for combustion under higher ventilation factor. Linear correlations can be fitted for the relationship between the intake airflow rate and the ventilation factor, with results very close to those equations derived from simple hydraulic theory. Averaged predicted HRR by the liquid fuel model in FDS did not agree with the experimental results. The burning rate in a room fire is strongly affected by radiative heat feedback which might not be well simulated in the FDS model. This result was incorporated with other studies [4,[33][34][35]. It is suggested that liquid fuel model should be applied carefully to engineering application. Further related work on validation and verification of liquid fuel model in FDS should be conducted.