Abstract

In fire hazard calculations, knowledge of the heat-release rate (HRR) of a burning item is imperative. Typically, room-scale calorimetry is conducted to determine the HRRs of common combustible items. However, this process can be prohibitively expensive. In this work, a method is proposed to invert for the HRR of a single item burning in a room using transient heat flux measurements at the walls and ceiling near the item. The primary device used to measure heat flux is the directional flame thermometer (DFT). The utility of the inverse method is explored on both synthetically generated and experimental data using two so-called forward models in the inversion algorithm: fire dynamics simulator (FDS) and the consolidated model of fire and smoke transport (CFAST). The fires in this work have peak HRRs ranging from 200 kW to 400 kW. It was found that FDS outperformed CFAST as a forward model at the expense of increased computational cost and that the error in the inverse reconstruction of a 400 kW steady fire was on par with room-scale oxygen consumption calorimetry.

1 Introduction

The heat-release rate (HRR) of a burning item is key to understanding the thermal effects of a fire on its surroundings [1]. HRR is the rate of exothermic energy release from the fire and drives the hazardous consequences of the fire. The HRR can be used in conjunction with a numerical model to predict quantities of interest such as the damage caused by the fire and the tenability conditions in a building.

The HRR is usually measured with one of many calorimetry methods. Bench scale calorimeters that take milligram size samples such as the bomb calorimeter, differential scanning calorimeter, and microscale combustion calorimeter are used to determine the heat of combustion and/or heat of gasification of a material which can then be used to calculate the HRR given a mass burning rate [24]. At the bench scale for items weighing of the order of a several grams, the oxygen consumption (cone) calorimeter is commonly used to determine the HRR [5]. The cone calorimeter relies on oxygen consumption measurement wherein all of the products of combustion are collected and the HRR is calculated based on the concentration of oxygen in the exhaust gases [6,7]. At the room scale, oxygen consumption calorimetry can be cost-prohibitive for large volumetric flow rates as the gas collection process requires a large hood, duct work, blower, and gas analysis equipment. Measuring HRR under the well-ventilated conditions of an oxygen consumption calorimeter does not necessarily reflect the HRR for burning scenarios in compartments for which radiative feedback and/or ventilation limits the burning rates. A typical gram scale sample of polyurethane burns with a peak HRR of about 200 W/g, and the oxygen consumption calorimeter fume hood requires an air flow of 0.035 m3/s [8]. For a furniture item like a sofa, the peak HRR could be as large as 3 MW with a required airflow of 3.5 m3/s [9,10].

Several studies have sought to estimate the HRR of a fire in a compartment without gas measurements by using inverse modeling techniques. Each of the following studies estimated the HRR that best matched the observed data using a fire model to reduce data. For these problems, a so-called “zone” model was used which simplifies the description of a compartment in to a hot upper zone and cooler lower zone with uniform temperature and properties. Richards et al. sought to determine the HRR and location of a fire using temperature data from ceiling mounted sensors to develop a system for automatic detection of growing fires [11,12]. Neviackas assessed the size of a steady fire in a compartment to aid firefighting operations using upper gas layer (UGL) temperatures as the observed data source [13]. Price et al. continued the work of Neviackas by using gas temperature and smoke layer height to predict the HRR of steady fires [14]. This work succeeded in providing reasonable predictions with simple models and sparse data but tended to overshoot the peak HRR. Using gas temperature measurements to estimate the hot gas layer temperature, Overholt and Ezekoye employed a fire size correlation and zone model in a predictor/corrector scheme to estimate the transient HRR of a fire in a compartment [15]. This method captured the magnitude of the HRR but struggled with response time and overshoot, even for steady cases.

In this work, we propose an inversion algorithm that utilizes the transient heat flux (as opposed to temperature) from an array of sensors and two distinct forward models to predict the transient HRR of a fire in a burn chamber that is approximately the size of a common residential compartment. We seek to build upon the HRR inversion techniques found in the literature by providing the framework for reconstructing HRR profiles that are comparable to those from oxygen consumption calorimetry, while using heat flux sensors and a higher fidelity forward model. We compare the performance of a zone model and a coarsely resolved computational-fluid-dynamics (CFD) model, when used as forward models in the inversion scheme, both in accuracy of the reconstructed HRR and computational cost. First, synthetic observed heat flux data are generated from a CFD model, fire dynamics simulator (FDS), for two HRR profiles with varying peak HRRs to test the performance of the inversion framework. Experimental data are collected from three sets of experiments each containing three replicates with a different HRR profile for each set. In the reminder of this paper, we will present the experimental system, measurements and procedures, fire evolution models, and HRR inversion model which will be used on synthetic and experimental data.

2 Experimental System, Conditions and Measurements

2.1 Experimental System.

The experimental structure is a room measuring 5.7 × 4.6 m with a ceiling height of 2.17 m. In this work, the two windows in the structure were always closed while the door was always open. The door is located on the west wall in the southwest corner and measures 0.91 m wide by 2.0 m tall. All interior walls and the ceiling were lined with 1.6 cm thick gypsum board and the floor is constructed from 2.54 cm thick concrete pavers over a sand and wood support structure.

The fires in the experiments and models were produced by two propane burners with horizontal dimensions of 0.32 m by 0.32 m and a height of 0.4 m. Three HRR profiles were used in the experiments, each with a duration of 5 min: a 200 kW steady HRR, a 400 kW steady HRR, and a “triangular” HRR profile. The final profile began at 40 kW, increased linearly to 400 kW by 90 s, and decreased linearly to 40 kW by 5 min before being shut off. The experimental HRR was estimated as the mass flowrate measured by the flow controller multiplied by the heat of combustion of propane assuming complete combustion. The fire location was held constant at 1.5 m from the north wall and 3.0 m from the west wall to the center of the first burner. When the second burner was added for the 400 kW tests, the first burner was not moved.

Heat flux and temperature measurements were taken during each experiment with typical profiles for sensors located approximately 2 m from a 200 kW steady fire shown in Fig. 1. In this figure, thermocouple 01 is located near the ceiling and the remaining thermocouples are spaced evenly in the vertical direction showing the transition from the hot upper gas layer to cooler lower layer. For thermocouples closer to the fire and for larger 400 kW fires, peak temperatures were of the order of 400–500 °C. The heat fluxes in Fig. 1 are taken from the front side of two directional flame thermometers (DFTs), where the front side refers to the side facing down toward the fire. Typical peak heat fluxes were less than 10 kW/m2, with the exception of DFTs that saw flame impingement where heat fluxes could reach 25 kW/m2. Previous work on this DFT configuration estimated experimental uncertainties between 9% and 13% for the range of heat fluxes observed in this experiment [16].

Fig. 1
Typical temperatures (a) and heat fluxes (b) for a 200 kW steady fire at sensors approximately 2 m from the fire
Fig. 1
Typical temperatures (a) and heat fluxes (b) for a 200 kW steady fire at sensors approximately 2 m from the fire
Close modal

2.2 Heat Flux and Temperature Measurements.

In the experimental structure, heat flux data were collected with DFTs constructed following the methods outlined by Kokel et al. [17]. DFTs are composed of ceramic fiber insulation sandwiched between two square steel plates measuring 11.43 cm wide by 0.1905 cm thick. This insulation was initially 2.54 cm thick and it has been compressed to 2.22 cm by four bolts and spacers. A K-type thermocouple is fixed in the center of the plate between the plate and insulation with a thin square of steel shim stock. The net heat flux at each DFT is calculated by solving the 1D transient heat conduction equation in the insulation and steel plates approximating the plate temperatures as the thermocouple temperatures.

There are a total of 28 DFTs in the structure and they are divided between wall and ceiling stands as shown in Fig. 2. There are eight ceiling stands with a total of 16 horizontally mounted DFTs forming an approximately 1 m grid on the ceiling. Each stand contains four thermocouples spaced evenly along the supporting vertical bar with the top thermocouple located at the ceiling. The two wall stands each contain six vertically mounted DFTs arranged in two rows of three DFTs at heights 0.96 and 1.9 m with 0.45 m between the DFTs on each row. One thermocouple is located above the center DFT on each row.

Fig. 2
Layout of experimental sensor locations in the structure (a). Dots represent heat flux sensor locations where an increase in darkness indicates a decrease in elevation. Picture of the experimental structure taken from the door (b).
Fig. 2
Layout of experimental sensor locations in the structure (a). Dots represent heat flux sensor locations where an increase in darkness indicates a decrease in elevation. Picture of the experimental structure taken from the door (b).
Close modal

3 Fire Models and Data Reduction Modeling

3.1 Forward Models.

This work employs two forward models to issue predictions for the transient heat flux data: consolidated model of fire and smoke transport version 7 (CFAST) and fire dynamics simulator version 6. CFAST is a two-zone fire model primarily used for simulating the behavior of compartment fires where the zones are a simplified representation of the hot upper gas layer and cooler lower gas layer in a compartment. It relies on a combination of physics and correlations to predict the bulk mass and energy transport processes that occur in a compartment or series of connected compartments (Fig. 3).

Fig. 3
A zone model approximates the gases in a compartment as a hot upper layer at Tu and cold lower layer at Tl. The fire (Q˙) is modeled as a bulk energy flow in to the room. Mass flows (m˙) are modeled with simplified physics and correlations.
Fig. 3
A zone model approximates the gases in a compartment as a hot upper layer at Tu and cold lower layer at Tl. The fire (Q˙) is modeled as a bulk energy flow in to the room. Mass flows (m˙) are modeled with simplified physics and correlations.
Close modal

Fire dynamics simulator is a computational fluid dynamics code that uses large eddy simulation models to numerically solve the Navier–Stokes equations [18]. We employ FDS as a forward model for the inversion scheme with a cell size of 20 cm. FDS, even in these coarse simulations predicts the flow field in the compartment whereas CFAST only models the flow at the compartment openings and fire plume.

To determine the maximum acceptable grid size, a scaling study would need to be conducted to find the smallest cell size beyond which the solution does not change within an acceptable tolerance. We choose a relatively coarse mesh to gain speed at the expense of accuracy as the forward model is called of the order of hundreds of times to reconstruct the HRR curve. In FDS at this mesh resolution, 15 s of simulation time takes approximately 60 s of wall clock time for a 1000 kW fire, whereas CFAST requires less than 2 s of wall clock time for 300 s of simulation time for the same fire.

A finely resolved FDS model was used to generate the synthetic data for which the inversion methodology was tested and refined prior to being used on experimental data. To generate synthetic observations of the heat flux at each sensor, FDS was run with a mesh size of 5 cm. In the synthetic data sets, virtual net heat flux sensors were spaced at a 1.0 m grid (excluding the windows) along the west, north, and east walls and ceiling of the main compartment resulting in a total of 36 sensors. The virtual net heat flux sensors were placed on the surfaces of the gypsum boards in the simulated compartment.

3.2 Directional Flame Thermometer Modeling.

To use the DFTs in an inversion scheme for compartment fire cases, they must be represented in the two forward models: CFAST and FDS. The DFT itself also requires a model to invert for the net heat flux from the temperature data collected at each plate. The net heat flux at each DFT is calculated by solving the 1D transient heat conduction equation in the insulation by approximating the plate temperatures as the measured thermocouple temperatures and inverting for the net heat flux at the outside face of each plate. The DFTs constructed for these experiments followed the work of Kokel et al. with a few revisions to the plate width and thickness as well as the insulation thickness [17]. A full description of the 1D transient conduction code can be found in Kokel's thesis [19]

3.2.1 Consolidated Model of Fire and Smoke Transport Directional Flame Thermometer Model.

Consolidated model of fire and smoke transport does not allow for sensors made of multiple materials to be included in the simulation or have temperature dependent thermal properties, so the 1D transient conduction code was coupled with the boundary conditions from surrogate sensors in CFAST. The energy balance for the thermally lumped steel plates is
(1)

where lp, ρ, cp, α, ϵ, and T are the thickness, density, specific heat, absorptivity, emissivity, and temperature of the steel plate, respectively. The incident radiant heat flux per unit area from the surroundings and conduction to the insulation are qinc and qcond, respectively, and Tinf is the surrounding gas temperature. The final parameters, h and σ, are the heat transfer coefficient at the plate and Stefan-Boltzmann constant. From the surrogate sensor in CFAST we take the incident radiant heat flux, surrounding gas temperature, and heat transfer coefficient to define the boundary condition at the front side of each DFT. CFAST does not provide back side conditions, so these are approximated to be purely convective with the same heat transfer coefficient and gas temperature.

These boundary conditions were applied in two ways in this work: entirely from CFAST and as hybrid of CFAST, experiments, and the Alpert correlation [20]. The hybrid method takes the incident radiant heat flux from CFAST, the surrounding gas temperature from the nearest thermocouple in the experiments, and the heat transfer coefficient from the Alpert correlation for a given HRR.

3.2.2 Fire Dynamics Simulator Directional Flame Thermometer Model.

Fire dynamics simulator contains a 1D transient conduction solver with the ability to specify multiple layers of materials with temperature dependent properties. A virtual DFT was constructed in FDS by specifying material properties for the steel plates and cerablanket with the thermal property versus temperature tables from Kokel [19]. There are two notable differences between the FDS DFT model and 1D transient conduction model from Kokel: the steel plate is lumped in the Kokel model while in FDS the plate is discretized within the FDS conduction model, and the thermal properties in the Kokel model are calculated from a curve fit to the table values while FDS interpolates properties from the tables.

Preliminary comparisons between FDS with known HRRs and experiments showed that FDS was consistently under-predicting the net heat flux at the ceiling sensors by a factor of one half (Fig. 4). The heat transfer coefficients at the ceiling DFTs ranged from 7 to 15 W/m2/K which is indicative of natural convection [21]. However, heat transfer at the ceiling is driven by the fire plume impinging on the ceiling and spreading out into what is called a ceiling jet. The well known Alpert ceiling jet correlation predicts heat transfer coefficients ranging from 17 to 111 W/m2/K for these same sensors for HRRs ranging from 200 to 400 kW [20]. This correlation calculates the heat transfer to a ceiling as a function of the fire size in compartment geometry. We also note that this correlation was derived for unconfined ceilings where a hot gas layer does not develop. When the heat transfer coefficient predicted by the Alpert correlation as a function of HRR is applied directly in FDS, the heat flux predictions improve as seen in Fig. 4. However, the magnitude of the heat flux leaving the DFTs after the fire has extinguish is over-predicted as the correlation was not designed for a no-fire scenario.

Fig. 4
Measured (Exp.) versus predicted (FDS) heat flux at the front side of DFTs 01 and 02 for the 200 kW steady fire using the default FDS heat transfer coefficient (a) and Alpert correlation heat transfer coefficient (b)
Fig. 4
Measured (Exp.) versus predicted (FDS) heat flux at the front side of DFTs 01 and 02 for the 200 kW steady fire using the default FDS heat transfer coefficient (a) and Alpert correlation heat transfer coefficient (b)
Close modal
The Alpert correlation will not work with the inversion algorithm as FDS does not allow a user specified heat transfer coefficient to vary as a function of heat release rate. Fortunately, FDS allows the user to supply their own correlation for the Nusselt number (Nu)
(2)
where L is the convective length scale defined by the user (default 1 m) and k is the thermal conductivity of the fluid. The Nusselt number correlation in FDS follows the form:
(3)

were C1, C2, n, and m are correlation constants, Re is the Reynolds number, and Pr is the Prandtl number. The default parameters are from the correlation for turbulent flow over an isothermal flat plate from Incropera and De-Witt [21]. The values for this correlation are C1=0,C2=0.037, n  = 4/5, m  = 1/3, and Pr = 0.7.

For this work, a correlation was derived with data from experiments and simulations of the 200 kW and 400 kW steady HRR fires to fit this form. The Prandtl number and exponent, m, are assumed to be constant at 0.7 and 0.33, respectively. The constant C1 remains set to zero. The characteristic length scale of the Reynolds number and Nusselt number was set to be the distance from the fire in the horizontal plane. We fit log(C2) and n from Eq. (3) using least squares regression by using a mixture of experimental and FDS data to calculate the Reynolds and Nusselt numbers.

To calculate the Reynolds and Nusselt number at each DFT, the following properties must be evaluated: gas velocity, heat transfer coefficient, film temperature, and the thermal conductivity, specific heat, density, and viscosity of the gas at the film temperature. As the gas velocity was not measured, the velocity from the 20 cm grid FDS simulations at the front side of the virtual DFTs was used. This fit was designed around the 20 cm case to facilitate its use as a forward model in future inversion algorithms. All thermal properties are evaluated at the average film temperature that was estimated from the experiments by using the DFT surface temperature and the closest thermocouple to each DFT during the second half of the test from 2.5 min after ignition to 30 s before extinction (2 min total). Additionally, DFT 3 was left out of the fit due to its proximity to the fire in relation to the location of the closest thermocouple.

To calculate the experimental Nusselt number, we require an estimate of the heat transfer coefficient at each DFT for both HRRs. We approximate the convective heat flux at each DFT (qconv) by taking the net heat flux from the experiments and subtracting the net radiative flux from FDS
(4)
The heat transfer coefficient can then be estimated at each DFT as follows:
(5)

The Nusselt number was calculated using the average heat transfer coefficients from the 400 kW and 200 kW experiments. The values for C and n from the least squares regression are 1.017 and 0.665, respectively. The resulting Nusselt number fit is 6–7 times larger than default correlation for the range of Reynolds numbers in this scenario. With the new fit, the same DFTs from Fig. 4 are shown in Fig. 5. This fit shows good prediction of the heat fluxes during the fire as well as after the fire has extinguished. This adaptability will be beneficial when the model is used to invert for fires with time-varying HRRs.

Fig. 5
Measured (Exp.) versus predicted (FDS) heat flux at the front side of DFTs 01 and 02 for the 200 kW steady fire using the new Nusselt number correlation
Fig. 5
Measured (Exp.) versus predicted (FDS) heat flux at the front side of DFTs 01 and 02 for the 200 kW steady fire using the new Nusselt number correlation
Close modal

The estimated experimental Nusselt number fit was applied to the FDS input files and run for the 5 cm and 20 cm grid resolutions and the 200 kW and 400 kW steady fires. The l2 error between the simulation and experiment from 10 s after ignition to 10 s before extinction for only the ceiling DFTs is summarized in Table 1. In each case where both the default FDS and Nusselt number fit simulations were run, the error in the heat flux predictions decreased by approximately a factor of two. These results indicate that it will be feasible to use the 20 cm simulation as a forward model to invert for HRRs between 200 kW and 400 kW.

Table 1

l2 error in the heat flux at all ceiling DFTs for the default heat transfer coefficient and new fit

HRRResolution (cm)Default FDS errorC, n fit error
200 kW20417169
5368193
400 kW20640284
5N/A329
HRRResolution (cm)Default FDS errorC, n fit error
200 kW20417169
5368193
400 kW20640284
5N/A329

While comparisons to experimental data showed good agreement between the simulated FDS DFTs and the experiments at the ceiling DFTs, FDS with a grid resolution of 20 cm did not adequately predict the experimental heat flux at the DFTs mounted on the wall stand. It was found that the radiant heat flux predicted by FDS at the midlevel wall DFTs closely followed the net heat flux measured in the experiments indicating that the convective heat flux was being predicted incorrectly. In this work, the midlevel wall DFTs were brought into the cost function in the form of experimental net heat flux and FDS radiant heat flux. The implications of this choice will be explored in Sec. 4.

3.3 Data Processing Model.

Two methods of postprocessing the observed data were used in this work, a floating average smoother and a median filter. The floating average smoother computes the floating average within a specified window around the current point. The median filter sorts the points within a specified window and assigns the median to the current point which has the advantage of preserving the effects of relatively quick transient events such as the ignition of a fire while still smoothing turbulent fluctuations like the floating average smoother.

Both the synthetic sensors in FDS and the experimental DFTs were set to record data at a sampling rate of 1 Hz. The filter windows were manually set to 9 s and 61 s for the floating average smoother and median filter, respectively, in an attempt to preserve the balance between detecting ignition events and smoothing out turbulent fluctuations. The two data processing methods are shown in Fig. 6 on data measured by a ceiling mounted DFT located approximately 1 m from a fire with a steady HRR of 200 kW burning for 300 s. The window smoother decreases the slope of the heat flux measured at the ignition event and tracks the measured heat flux during larger fluctuations while the median filter preserves the slope of the heat flux after ignition and tracks the slowly decreasing trend of the measured heat flux. If the window of the floating average was increased to be the same as the median filter, the small fluctuations during the experiment would be smoothed out but the sharp increase in heat flux after ignition would not be preserved.

Fig. 6
Example DFT heat flux profile from 200 kW steady fire experiment before and after smoothing with the window smoother (a) and median filter (b)
Fig. 6
Example DFT heat flux profile from 200 kW steady fire experiment before and after smoothing with the window smoother (a) and median filter (b)
Close modal

3.4 Inversion Algorithm.

The inversion algorithm seeks to find the HRR of a burning item at n equally spaced points in the time-line of the experiment by sequentially minimizing a cost function (Cn) at each point. Brent's method (as implemented in Scipy) is used to minimize the cost function composed of the squared l2 error in the transient heat flux profile (d) at specified points on or adjacent to the solid surfaces [22]. This method is an improvement on the golden section search wherein the first three iterations of Brent's method are identical to a golden section search. After the first three iterations, inverse parabolic interpolation is used to propose the next iteration and the cost function is evaluated at the proposed point if it is within the search bounds and less than one half the distance of the second to last step. Otherwise, the golden section search is used for the next point. The second criteria is needed to prevent the parabolic steps from oscillating around the minimum. In this particular inverse framework, Brent's method shaves 3–4 iterations off of the optimization step when compared to the golden section search alone which takes 15–16 iterations for each point in the reconstructed transient HRR.

The cost function is specified as
(6)

where Nw and Nc are the number of wall and ceiling sensors, respectively. The heat flux data (d) are divided between sensors located on or adjacent to the walls and ceiling, subscripts w and c, respectively, and between observed and predicted values denoted by subscripts obs and pre where the observed data are either the synthetic or experimental data and the predicted data are from the forward models. The hyper-parameter β is used to weight the influence of the ceiling and wall sensors where β can be set between zero and one such that a value approaching one weights the ceiling sensors higher in the cost function (Eq. (6)).

The inversion algorithm begins at time n  = 0 where the HRR (Q˙) is set to zero kW. Brent's method is used to optimize for the HRR at time n  + 1 by iteratively calculating a new point (Q˙n+1), running the forward model, and evaluating the cost function. The total simulation/experimental time is discretized in to nend points where the HRR will be reconstructed.

The following convergence criterion as implemented in the optimization toolbox in Scipy version 0.19.0 is checked after each cost function evaluation [23]:
(7)
where Q˙ is the most recent proposal for Q˙n+1, a and b are the moving lower and upper bounds, respectively, of the search bracket, and τ is given by
(8)

In this case, the stopping tolerance, Q˙tol, is set to 1 kW and ϵ is set to 2.2 × 10−16 [23]. Once the convergence criterion is met, the algorithm saves the best HRR and proceeds to the next point in time keeping all previous HRR values fixed. This process is repeated until the full HRR has been reconstructed up to the time-step n=nend.

4 Results and Discussion

The inverse reconstructions are judged on three error categories when compared to the true HRR: the l2 error, the mean absolute error, and the mean l1 error as a percentage of the peak HRR. This section details the results of the inverse HRR reconstruction algorithm on synthetic and experimental transient heat flux data using both CFAST and FDS as the forward model in the algorithm. We limit the exploration of the inversion algorithm to two HRR profiles: a steady fire and a triangular fire. The steady HRR profile begins at 0 kW and increases to the specified HRR within 1 s of ignition and burns with a duration of 300 s. The triangular HRR fire is a simplified approximation of a real fuel such as a chair or bag of combustible material parameterized by a peak HRR, time to peak HRR, and total burning duration. The triangular fire begins and ends at a fixed value, increases linearly to the peak HRR at the specified time to peak, and then decreases linearly to the ending HRR at the total burning time.

The HRR solution was reconstructed with CFAST as the forward model with a reconstruction interval of 15 s for slowly changing HRRs (the triangle HRR) and 5 s for HRRs with relatively fast changes (approximately greater than 100 kW per second) to minimize overshoot in the solution. However, it was found that using FDS as the forward model produced an oscillatory solution at reconstruction intervals of 15 s and 5 s. This issue was mostly alleviated by fixing the reconstruction interval at 5 s and running the simulation 4 s beyond the end of the reconstruction interval. The HRR was projected during the extra 4 s based on the slope of the HRR curve from the previous reconstruction time-step following Eq. (9). Additional data were generated for Δtre1 seconds so as not to generate a new “restart” file in FDS which occurs every Δtre seconds
(9)

The projected HRR (Q˙) is a function of the proposed HRR at time n  + 1 (Q˙), the HRR saved at the last reconstruction time (Q˙n), and the reconstruction interval (Δtre). While including the projected net heat flux data in the cost function increased the computational time required to complete the inversion algorithm, the final reconstructed solution is greatly improved by the projection method.

4.1 Synthetic Data.

Synthetic net heat flux data were generated by FDS with a 5 cm grid for two triangular HRR profiles with peak HRRs of 200 kW and 1000 kW. Each triangular fire began and ended at 0 kW with a time to peak HRR of 90 s and total burning time of 300 s. In these synthetic examples, the data from all ceiling and wall mounted sensors were used in computing the cost function with a weighting of β following Eq. (6), and the optimal β was found via a grid search with a resolution of 0.05 by minimizing the l2 error between the reconstructed HRR and “true” synthetic HRR. A summary of these errors for the best values of β can be found in Table 2.

Table 2

Errors in the HRR reconstructions of synthetic triangular fires for the best values of β

Peak HRRForward modelβl2 error (kW)Mean absolute error (kW)Percent of mean HRR (%)
200 kWCFAST0.7538.46.616.61
FDS0.21089.549.54
1000 kWCFAST0.595715931.8
FDS0.1533135.47.08
Peak HRRForward modelβl2 error (kW)Mean absolute error (kW)Percent of mean HRR (%)
200 kWCFAST0.7538.46.616.61
FDS0.21089.549.54
1000 kWCFAST0.595715931.8
FDS0.1533135.47.08

The window smoother was used for processing the synthetic data before inverting with both CFAST and FDS as forward models. Figure 7 shows the inversion results for a 200 kW triangular fire using CFAST and FDS with 15 and 5 s reconstruction intervals, respectively. The CFAST solution is stable and tracks the synthetic observed HRR up until about 180 s after which it begins to over-predict the HRR. A β value of 0.75 was found to yield the CFAST solution with the lowest l2 error of 38.4 kW. The inverse FDS reconstruction over-predicts the HRR in the growth period but tracks the synthetic HRR better after the peak with an l2 error of 108.4 kW.

Fig. 7
Inversion for a 200 kW triangular fire on synthetic data using CFAST (a) and FDS (b)
Fig. 7
Inversion for a 200 kW triangular fire on synthetic data using CFAST (a) and FDS (b)
Close modal

In the case of the 1000 kW triangular fire, CFAST significantly over-predicts the HRR above approximately 500 kW (Fig. 8). This is a result of the plume temperature correlation used by CFAST having an upper limit of ΔT=900 K, where this is the temperature difference between the temperature at the point where the plume impinges on the ceiling and the ambient temperature. At about 60 s of simulation time, the plume temperature reaches this threshold which in turn propagates to the ceiling jet temperature through the Alpert correlation which drives the convective heat transfer to the ceiling sensors. As the ceiling jet temperature is effectively constant at HRRs above approximately 500 kW, the only hope for the inversion algorithm to increase heat transfer to the ceiling is via radiation from the upper gas layer and surfaces of the compartment. Therefore the HRR that minimizes the cost function becomes increasingly higher than the true synthetic HRR that generated the data.

Fig. 8
Inversion for a 1000 kW triangular fire on synthetic data using CFAST (a) and FDS (b)
Fig. 8
Inversion for a 1000 kW triangular fire on synthetic data using CFAST (a) and FDS (b)
Close modal

Fire dynamics simulator performs better than CFAST on the 1000 kW triangular HRR fire with an optimal β value of 0.15 (Fig. 8(b)). Similar to the 200 kW triangular HRR case, the inversion algorithm with FDS over-predicts during most of the growth period, however it under-predicts the HRR during the decay period. FDS was expected to perform better than CFAST as a more finely resolved FDS simulation was used to generate the synthetic data, but the performance of CFAST in the 200 kW triangular HRR case was surprising.

The results of using the inverse HRR reconstruction algorithm on synthetic data represent a lower limit on the potential error of the current algorithm when used on experimental data. The performance of the inversion algorithm can be compared across different fire scenarios by taking the mean absolute error as a percentage of the mean HRR. Based on Table 2, we aspire to reach a mean absolute error on experimental data of 6.61% and 9.54% of the peak HRR using CFAST and FDS, respectively, for HRRs of the order of 200 kW. We expect CFAST will not work for peak HRRs of 500 kW and larger, while using FDS as the forward model can potentially produce results with errors as low as 7.08% of the mean HRR.

4.2 Experimental Directional Flame Thermometer Data.

The inverse HRR reconstruction algorithm was tested on net heat flux data gathered by DFTs from three sets of experiments: a 200 kW steady HRR, a 400 kW steady HRR, and a triangular fire with a starting and ending HRR of 40 kW and a peak HRR of 400 kW at 90 s (each test burned for a duration of 300 s). Initially, only data from the ceiling DFTs were used (β = 1.0) with Consolidated model of fire and smoke transport and FDS, as significant effort was put in to understanding the heat transfer at the ceiling. The midlevel wall stand DFTs were only brought in to the cost function with FDS as the forward model as CFAST did not adequately predict the heat flux at these DFTs.

The median filter was used on the experimental data as its use results in smoother HRR reconstructions with CFAST as the forward model. Using the median filter versus the window smoother has little effect on FDS HRR reconstructions as the projection to future time steps effectively stabilizes the resulting HRR profile. The CFAST inversion process takes around 30 min of wall clock time on a single processor core for most of the fires in this study. Inverse HRR reconstruction with FDS generally requires 40× the computational time of CFAST. However, FDS is parallelizable and each simulation has been setup to use four cores resulting in approximately 300 min of wall clock time per inverse FDS reconstruction.

4.2.1 Consolidated Model of Fire and Smoke Transport Inversion.

Using only the ceiling DFTs, CFAST over-predicts consistently in both the 200 kW and 400 kW steady fire cases (with the 400 kW case shown in Fig. 9(a)). The driving mechanism behind this trend is the over-prediction of the gas temperature near the DFT while the heat transfer coefficient is greatly under-predicted resulting in CFAST under-predicting the net heat flux at most of the ceiling sensors.

Fig. 9
Inversion for a 400 kW steady fire on experimental data using only ceiling DFTs with CFAST providing all DFT boundary inputs (a) and the hybrid approach where CFAST provides radiation only (b)
Fig. 9
Inversion for a 400 kW steady fire on experimental data using only ceiling DFTs with CFAST providing all DFT boundary inputs (a) and the hybrid approach where CFAST provides radiation only (b)
Close modal

Modifying the DFT boundary conditions to include the experimentally measured gas temperatures and the heat transfer coefficient from the Alpert correlation as detailed in Sec. 3.2.1 results in a reconstruction that under-predicts the experimental HRR (Fig. 9(b)). Again, the same trend was observed in both the 200 kW and 400 kW steady fire cases. Based on the results of Sec. 3.2.1, the inversion algorithm was not able to recover from CFAST over and under-predicting the net heat flux at different sensors during the first 2 min of the fire. The errors for CFAST on steady HRR experiments are reported in Table 3, and the mean absolute error as a percentage of the mean HRR is around 20% which is significantly worse than the synthetic data cases.

Table 3

Errors in the HRR reconstructions of experimental fires with CFAST and β= 1.0

Peak HRRBoundary conditionl2 error (kW)Mean absolute error (kW)Percent of mean HRR (%)
200 kW steadyCFAST53267.633.8
Hybrid45253.927.0
400 kW steadyCFAST56269.017.3
Hybrid90610125.3
400 kW triangularCFAST44551.923.6
Hybrid96011250.9
Peak HRRBoundary conditionl2 error (kW)Mean absolute error (kW)Percent of mean HRR (%)
200 kW steadyCFAST53267.633.8
Hybrid45253.927.0
400 kW steadyCFAST56269.017.3
Hybrid90610125.3
400 kW triangularCFAST44551.923.6
Hybrid96011250.9

4.2.2 Fire Dynamics Simulator Inversion.

Investigation of the performance of FDS as the forward model began with only including the ceiling DFTs in the cost function. Figure 10 shows the inverse reconstruction on the 200 kW and 400 kW steady HRR fires with only the ceiling DFTs (i.e., β= 1.0). The reconstructed HRR is within 25% of the experimental values for most of the duration of the fire, however the reconstructed HRR initially under-predicts the experimental HRR and increases over time by around 100 kW in each case. The result is an over-prediction of the experimental HRR in the 200 kW steady fire case for the final 150 s of the experiment while the 400 kW steady solution appears to fluctuate around the experimental HRR in the final 150 s.

Fig. 10
Inversion for a 200 kW (a) and 400 kW (b) steady fire on experimental data using FDS with only the ceiling DFTs
Fig. 10
Inversion for a 200 kW (a) and 400 kW (b) steady fire on experimental data using FDS with only the ceiling DFTs
Close modal

Based on investigation of the heat flux at the wall DFTs, the net heat fluxes at these DFTs are not predicted correctly by FDS at a 20 cm grid resolution. It was determined that the incorrect predictions were the result of incorrect convective heat transfer at the walls, and when the components of the net heat flux are broken down as in Fig. 11 it can be seen that the radiative component of the net heat flux from FDS provides a good approximation to the net experimental heat flux. We also note that height of the hot UGL is under-predicted in the 20 cm FDS simulations such that the midlevel DFTs are predicted to be in the UGL when they are below the UGL in the experiment, resulting in over-prediction of the convective heat flux.

Fig. 11
Radiative (Rad) and convective (Conv) components of the net heat flux from FDS versus the net experimental heat flux measured by the midlevel DFTs on stand I for the 400 kW steady fire
Fig. 11
Radiative (Rad) and convective (Conv) components of the net heat flux from FDS versus the net experimental heat flux measured by the midlevel DFTs on stand I for the 400 kW steady fire
Close modal

The six midlevel wall DFTs were included in the cost function and a value of β that minimized the l2 error was found with a partial grid search that had a resolution of 0.05 for the 200 kW and 400 kW steady HRR fires as well as the 400 kW triangular HRR fire. In each case, there is a significant valley in the l2 error near the best β values which are either 0.05 or 0.1. We note that as there are 16 ceiling DFTs and 6 wall DFTs included in the cost function, and equal weighting corresponds to a β value of 0.25.

The l2 and mean absolute errors are reported in Table 4 at the best β value for each experimental HRR profile as these values will be referenced in the proceeding discussion. The mean absolute error as a percentage of the mean HRR will be used for comparing the errors of fires with different profiles and maximum HRRs. When using the inversion code for calorimetry of an item with an unknown HRR, this method of finding β will not be applicable. However, as will be shown in the triangle fire case, the value of β for an unsteady fire falls within the range of values for steady fires of similar HRR magnitude. This suggests that a range of steady calibration cases can be run to determine a β versus HRR curve that can be used for calorimetry of an unknown HRR.

Table 4

Errors in the HRR reconstructions of experimental fires with FDS

Peak HRRBest βl2 error (kW)Mean absolute error (kW)Percent of mean HRR (%)
200 kW steady0.1012613.36.65
400 kW steady0.0529626.46.60
400 kW triangular0.0515014.16.41
Peak HRRBest βl2 error (kW)Mean absolute error (kW)Percent of mean HRR (%)
200 kW steady0.1012613.36.65
400 kW steady0.0529626.46.60
400 kW triangular0.0515014.16.41

The best FDS solution for the 200 kW steady fire (β=0.1) still increases over time but to a lesser extent as when only the ceiling DFTs were used in the cost function (Fig. 12). The mean absolute error for this reconstruction is 13.3 kW or 6.65% of the mean experimental HRR which is a great improvement over inverting with CFAST and is of the order of the low error limit established by the synthetic inversion.

Fig. 12
Inversion for a 200 kW steady fire on experimental data using FDS with β= 0.1
Fig. 12
Inversion for a 200 kW steady fire on experimental data using FDS with β= 0.1
Close modal

The inverse algorithm reconstructed the 400 kW steady experimental HRR with a similar level of accuracy to the 200 kW case (6.6% of the mean experimental HRR) as seen in Fig. 13. We note that there is a bias toward over-prediction at the 200 kW HRR and under-prediction at the 400 kW HRR. Yeager conducted a study on the uncertainty in HRR measurements using room-scale oxygen consumption calorimetry and methane with a known flowrate [24]. Four test fires were run that burned at steady HRRs for 3 min intervals and a sampling rate of one sample every 5 s, the same rate as the inverse reconstruction in this work. Two of Yeager's tests contained intervals of steady burning at 400 kW for 3 min separated by 3  min intervals of steady burning at approximately 57 kW. Three of the measured HRR profiles compared to the expected HRR from the supplied methane flowrate are presented in Fig. 13.

Fig. 13
Inversion for a 400 kW steady fire on experimental data using FDS with β= 0.05 (a) and room-scale oxygen consumption HRR measurements for a 400 kW steady fire from Yeager (b)
Fig. 13
Inversion for a 400 kW steady fire on experimental data using FDS with β= 0.05 (a) and room-scale oxygen consumption HRR measurements for a 400 kW steady fire from Yeager (b)
Close modal

The mean absolute errors as a percentage of the peak HRR are 6.2%, 7.0%, and 10.2% for test 2 part A, test 2 part B, and test 3 part A, respectively. These errors are comparable to the 6.6% error in the 400 kW steady HRR reconstructed solution using FDS. The main difference between the inverse reconstructions and room-scale oxygen consumption calorimetry is the response time. Yeager reports a response time of 30 s to reach 95% of the target HRR while the inverse reconstruction reaches 400 kW within 10 s after ignition. After the reconstructed and calorimetry HRR profiles reach 400 kW, both solutions fluctuate ±50 kW from the true solution. This indicates that the inverse HRR reconstruction framework could be a viable alternative to the more expensive room-scale oxygen consumption calorimetry.

For the 400 kW triangular fire, the FDS reconstruction under-predicts the experimental HRR the greatest during initial 50 s, but otherwise tracks the experimental HRR closely as seen in Fig. 14 with under-prediction occurring to a lesser extent during the remainder of the fire. The mean absolute error is 14.1 kW or 6.41% of the mean experimental HRR and the best β is the same as the 400 kW steady HRR case at 0.05. Additionally, we note that the 400 kW triangle fire data were not used in the calibration of the Nusselt number fit for the ceiling DFTs. This indicates flexibility of the correlation to be applied to transient HRR profiles within the range of calibrated steady profiles (i.e., 200–400 kW in this case). Future use of this methodology could entail a series of steady calibration experiments with known HRRs over the range of interest where both the Nusselt number and β versus HRR functions could be calibrated before testing items of interest.

Fig. 14
Inversion for a 400 kW triangular fire on experimental data using FDS with β=0.05
Fig. 14
Inversion for a 400 kW triangular fire on experimental data using FDS with β=0.05
Close modal

The mean absolute error as a percentage of the mean HRR is on the same order as the best error achieved when inverting on synthetic data. We note that the synthetic data source is net heat flux on gypsum wall and ceiling panels while DFTs were used in the experiments. The steel plates of the DFTs respond faster to changes in the thermal environment than gypsum as the thermal diffusivity of steel is approximately 2.5× that of gypsum. The increased sensitivity of the DFTs to changes in gas temperature and incident radiation might yield improved accuracy over using gypsum as a heat flux sensor. The bias in the inverse reconstructions with FDS might be further alleviated by expanding the HRR range of the experiments used to generate the heat transfer coefficient at the ceiling and including a larger number of midlevel wall DFTs in the experiments.

5 Conclusions

An algorithm was built that inverted for the HRR of burning items in a compartment with the option of using two forward models, and the limits of employing each model in the inversion framework were explored. The models were initially tested on synthetic net heat flux data before moving to experimental net heat flux data measured by DFTs.

Generally, using FDS as the forward model resulted in a lower error at the expense of 40× the computational cost of inversion with CFAST. The only case where a lower error in the reconstructed HRR was achieved using CFAST was for a triangular fire with a peak HRR of 200 kW and synthetically generated data. A joint inversion model that switches between prediction with CFAST and FDS based on the HRR should be explored. Not surprisingly, the inversion process is less accurate when there are errors in the forward models, necessitating a calibration process for the space spanning the cases of interest. The Nusselt number correlation at the ceiling in FDS was calibrated to alleviate under-prediction of the heat transfer coefficient by the default model. Inverse problems typically require the specification of some set of hyper-parameters. In this case, the primary hyper-parameter investigated was the relative weighting of the heat flux data from the wall and ceiling sensors (β) in the cost function.

Comparison to room-scale oxygen consumption calorimetry data for a 400 kW steady fire showed that errors in the inverse reconstructed solution were on the same order as the calorimetry errors. These results indicate that with continued improvement to the characterization of the room and DFTs, the inverse HRR algorithm can produce comparable results to the oxygen consumption calorimeter. Inversion with FDS using data from a 400 kW triangular fire resulted in a mean absolute error of 3.53% of the peak HRR which is on the same order as the lowest error achieved by the inverse reconstruction framework on synthetic data. We note that the best value of β for the triangle case was the same as the steady case with the same HRR magnitude (400 kW), and the triangle fire data were not used in constructing the Nusselt number correlation. Future applications of this inversion framework will focus on more complex fires (e.g., burning liquid pools and real fuels) by leveraging steady calibration fires to construct the Nusselt number correlation and β versus HRR function for a range of HRRs.

Funding Data

  • U.S. National Science Foundation (Award No. 1707090; Funder ID: 10.13039/100000001).

Nomenclature

c =

ceiling

cp =

specific heat

CFAST =

consolidated model of fire and smoke transport

DFT =

directional flame thermometer

FDS =

fire dynamics simulator

HRR =

heat-release rate

k =

thermal conductivity

L =

convective length scale

lp =

plate thickness

Nc =

number of ceiling sensors

Nw =

number of wall sensors

Nu =

Nusselt number

obs =

observed value

pre =

predicted value

Pr =

Prandtl number

Q˙ =

heat-release rate

Re =

Reynolds number

T =

temperature

w =

wall

α =

absorptivity

β =

sensor weighting parameter

Δtre =

reconstruction interval

ϵ =

emissivity

ρ =

density

References

1.
Babrauskas
,
V.
, and
Peacock
,
R. D.
,
1992
, “
Heat Release Rate: The Single Most Important Variable in Fire Hazard
,”
Fire Saf. J.
,
18
(
3
), pp.
255
272
.10.1016/0379-7112(92)90019-9
2.
O'neill
,
M. J.
, and
Watson
,
E. S.
,
1966
, “
Differential Microcalorimeter
,” U.S. Patent No. 3,263,484.
3.
Quintiere
,
J.
,
2006
,
Fundamentals of Fire Phenomena
,
Wiley
,
Chichester, UK
.
4.
Lyon
,
R. E.
, and
Walters
,
R.
,
2002
, “
A Microscale Combustion Calorimeter
,” Federal Aviation Administration Office of Aviation Research, Washington DC, Report No. DOT/FAA/AR-01/117.
5.
Babrauskas
,
V.
,
1984
, “
Development of the Cone Calorimeter'a Bench-Scale Heat Release Rate Apparatus Based on Oxygen Consumption
,”
Fire Mater.
,
8
(
2
), pp.
81
95
.10.1002/fam.810080206
6.
Huggett
,
C.
,
1980
, “
Estimation of Rate of Heat Release by Means of Oxygen Consumption Measurements
,”
Fire Mater.
,
4
(
2
), pp.
61
65
.10.1002/fam.810040202
7.
Parker
,
W. J.
,
1984
, “
Calculations of the Heat Release Rate by Oxygen Consumption for Various Applications
,”
J. Fire Sci.
,
2
(
5
), pp.
380
395
.10.1177/073490418400200505
8.
ASTM
,
2007
, “
Standard Test Method for Heat and Visible Smoke Release Rates for Materials and Products Using an Oxygen Combustion Calorimeter
,” American Society for Testing and Materials, West Conshohocken, PA, Standard No. ASTM E 1354-04a.
9.
ISO
,
2016
, “
Reaction to Fire Tests—Room Corner Test for Wall and Ceiling Lining Products—Part 1: Test Method for a Small Room Configuration
,” Standard, International Organization for Standardization, Geneva, CH, Standard No. ISO 9705-1:2016.
10.
Babrauskas
,
V.
,
2002
,
SFPE Handbook of Fire Protection Engineering
, 3rd ed.,
National Fire Protection Association
,
Quincy, MA
.
11.
Richards
,
R.
,
Munk
,
B.
, and
Plumb
,
O.
,
1997
, “
Fire Detection, Location and Heat Release Rate Through Inverse Problem Solution—Part I: Theory
,”
Fire Saf. J.
,
28
(
4
), pp.
323
350
.10.1016/S0379-7112(97)00005-2
12.
Richards
,
R.
,
Ribail
,
R.
,
Bakkom
,
A.
, and
Plumb
,
O.
,
1997
, “
Fire Detection, Location and Heat Release Rate Through Inverse Problem Solution—Part II: Experiment
,”
Fire Saf. J.
,
28
(
4
), pp.
351
378
.10.1016/S0379-7112(97)00006-4
13.
Neviackas
,
A.
,
2007
, “
Inverse Fire Modeling to Estimate the Heat Release Rate of Compartment Fires
,” Master's thesis, University of Maryland, College Park, MD.
14.
Price
,
M.
,
Marshall
,
A.
, and
Trouvé
,
A.
,
2016
, “
A Multi-Observable Approach to Address the Ill-Posed Nature of Inverse Fire Modeling Problems
,”
Fire Technol.
,
52
(
6
), pp.
1779
1797
.10.1007/s10694-015-0541-7
15.
Overholt
,
K. J.
, and
Ezekoye
,
O. A.
,
2012
, “
Characterizing Heat Release Rates Using an Inverse Fire Modeling Technique
,”
Fire Technol.
,
48
(
4
), pp.
893
909
.10.1007/s10694-011-0250-9
16.
Kurzawski
,
A. J.
,
2018
, “
Inverse Modeling and Characterization of an Experimental Testbed to Advance Fire Scene Reconstruction
,” Ph.D. thesis, The University of Texas, Austin, TX.
17.
Kokel
,
P.
,
Weinschenk
,
C.
, and
Ezekoye
,
O.
,
2010
, “
Evaluation of Directional Flame Thermometer for Real-Time Inversion of Heat Flux
,”
ASME
Paper No. IHTC14-22917.10.1115/IHTC14-22917
18.
McGrattan
,
K.
,
McDermott
,
R.
,
Hostikka
,
S.
, and
Floyd
,
J.
,
2012
, Fire Dynamics Simulator (Version 6) User's Guide. National Institute of Standards and Technology, Gaithersburg, MD.
19.
Kokel
,
P. M.
,
2008
, “
Modeling of Directional Flame Thermometer for Real-Time Incident Radiation Measurements in Room Fire Testing
,” Master's thesis, University of Texas at Austin, Austin, TX.
20.
Alpert
,
R.
,
2008
,
SFPE Handbook of Fire Protection Engineering
, 4th ed.,
National Fire Protection Association
,
Quincy, MA
.
21.
Incropera
,
F. P.
,
De Witt
,
D. P.
, and
Lavine
,
A. S.
,
2007
,
Fundamentals of Heat and Mass Transfer
, 6th ed.,
Wiley
,
New York
.
22.
Brent
,
R. P.
,
1971
, “
An Algorithm With Guaranteed Convergence for Finding a Zero of a Function
,”
Comput. J.
,
14
(
4
), pp.
422
425
.10.1093/comjnl/14.4.422
23.
Virtanen
,
P.
,
Gommers
,
R.
,
Oliphant
,
T. E.
,
Haberland
,
M.
,
Reddy
,
T.
,
Cournapeau
,
D.
,
Burovski
,
E.
,
Peterson
,
P.
,
Weckesser
,
W.
,
Bright
,
J.
,
van der Walt
,
S. J.
,
Brett
,
M.
,
Wilson
,
J.
,
Jarrod Millman
,
K.
,
Mayorov
,
N.
,
Nelson
,
A. R. J.
,
Jones
,
E.
,
Kern
,
R.
,
Larson
,
E.
,
Carey
,
C.
,
Polat
,
İ.
,
Feng
,
Y.
,
Moore
,
E. W.
,
Vand erPlas
,
J.
,
Laxalde
,
D.
,
Perktold
,
J.
,
Cimrman
,
R.
,
Henriksen
,
I.
,
Quintero
,
E. A.
,
Harris
,
C. R.
,
Archibald
,
A. M.
,
Ribeiro
,
A. H.
,
Pedregosa
,
F.
,
van Mulbregt
,
P.
, and
Contributors
,
S.
,
2019
, “
SciPy 1.0–Fundamental Algorithms for Scientific Computing in Python
,”
Nat. Methods
, epub.https://www.nature.com/articles/s41592-019-0686-2
24.
Yeager
,
R. W.
,
1986
, “
Uncertainty Analysis of Energy Release Rate Measurement for Room Fires
,”
J. Fire Sci.
,
4
(
4
), pp.
276
296
.10.1177/073490418600400404