Assessing the influence of the input variables employed by fire dynamics simulator (FDS) software to model numerically solid-phase pyrolysis of cardboard

Understanding a material’s fire behaviour implies to know the thermal decomposition processes. Thermal analysis techniques are widely employed to study thermal decomposition processes, especially to calculate the kinetic and thermal properties. Cardboard boxes are widely employed as rack-storage commodities in industrial buildings. Hence, the characterization of the cardboard is considered a key factor for fire safety engineering, because it enables the determination of its thermal behaviour at high temperatures. The employment of mathematical or computational models for modelling the thermal decomposition processes is commonly used in fire safety engineering (FSE). The fire dynamics simulator (FDS) software is one of the most commonly used computational fluid dynamics softwares in FSE to address thermal analysis. To properly set up FDS and obtain accurate results, the numerical values of the thermal and kinetic properties are needed as input data. Owing to the large number of variables to be determined, a preliminary study is bound to be helpful, which can well assess the influence of each variable over the pyrolysis model, discarding or restricting their influence. This study, based on the Monte Carlo method, presents a sensitivity analysis for the variables utilized as input data by the FDS software. The results show the conversion factor α, i.e. the mass involved in each reaction, and the triplet kinetics have a major impact on the reproduction of the thermal decomposition process in fire computer modelling.


List of symbols A
Pre-exponential factor (s -1 ) E a Activation energy (kJ kmol -1 ) n j Reaction order of the reaction j (-) H r Heat of reaction (kJ kg -1 ) r Reaction rate at temperature T s T s Surface temperature (°C) q Density (kg m -3 ) C p Specific heat (kJ kg -1 K -1 ) k Conductivity (W m -1 K -1 ) e Emissivity (-) g Absorption coefficient (m -1 ) Y s;i Quotient between the density of solid material i produced by the reaction (q s;i ) at surface temperature T s divided by the initial density of the solid material (q s ð0Þ) (-) r ij Reaction rate (kg s -1 ) v si 0 j Yield produced of the material i by the reaction j r i 0 j Residue produced of the material i by the reaction j a Coefficient the conversion factor of reactant (-) _ q 000 s;c Heat release rate per unit volume produced by the chemical reaction of the sample (kJ m -3 ) m Sample mass (kg) M Reacting material P Submaterial generated as product of the reaction

Introduction
To understand the fire behaviour of materials, it is necessary to know their decomposition processes. The thermal analysis techniques, such as thermogravimetric (TG) analysis [1] or its derivative (DTG) and differential scanning calorimetry (DSC) [2], have been widely used for years [3] to obtain thermal and kinetic properties of a material that undergoes thermal decomposition under external heat flux.
Thermal decomposition processes occur through a series of multiple reactions, which are provided by the decomposition scheme describing these reactions (e.g. products generated, energy released, temperatures, and whether they are simultaneous or non-competitive).
In a simultaneous thermal analysis (STA) test, which integrates TG and DSC experiments, the samples have the milligram scale to avoid the thermal diffusion occurring inside them and the influence of specific heat or conductivity. While the effects of different initial sample masses were analysed in [4], Comesaña et al. [5] studied the evolution of internal heating of the sample and thermal lag during an STA test. In that study, despite the small size of the samples, variables such as conductivity (k), specific heat (C p ), emissivity (e), and absorption coefficient (g) determined the internal heating of the sample.
Pyrolysis models describe the decomposition of a material that undergoes several reactions under external heat flux. Computational fluid dynamic (CFD) software simulation packages employ manifold pyrolysis models, such as GPyro [6], ThermaKin [7], FireFOAM [8], Pyropolis [9], and fire dynamics simulator (FDS) [10]. We used the FDS software in this study. FDS is a free opencode CFD software developed by the National Institute of Standards and Technology (NIST), which has been widely employed and validated by several research centres and universities to model phenomena involving fire, such as pyrolysis processes [11,12] or full-scale fire tests [13,14]. Its open access, the wide range of customizations of the reactions and material properties, and its tested reliability make FDS a suitable software for conducting the solidphase pyrolysis analysis in this study. Typically, FDS and the other CFD software packages require as input data the terms that compose the Arrhenius equation [15][16][17] (A, E a , and n, also known as triplet kinetics) and the thermal properties (C p , k, e, etc.) to represent the thermal decomposition processes. Because some thermal and kinetic properties cannot be achieved experimentally, other methods such as model-fitting methods [18], isoconversional methods [19], or numerical approach [20] are necessary.
FDS allows the users to customize the material properties and the thermal decomposition process. To properly reproduce the decomposition processes with FDS, not only the definition of the reactions and their characteristics, but also the values of the thermal properties are needed. Meanwhile, owing to the high degree of freedom required by the user to define the solid-phase pyrolysis, it is highly likely to require several properties, which should be introduced as input data. The authors in [8] estimated that at least nine properties are required for each reaction. The more the complex scheme used to model the pyrolysis, the more the variables are required to be introduced in the numerical approaching method, which slows down the process. Consequently, it may be helpful to conduct a preliminary analysis that assesses the influence of each variable over the FDS pyrolysis model and discards or restricts the less influential ones.
Previous studies have analysed the sensitivity of the variables by applying different methodologies to estimate their influence. There are mainly two approaches for addressing a sensitivity analysis: local sensitivity analysis and global sensitivity analysis. The nature of the analysed system or mathematical problem is the capital factor for selecting a method. Local sensitivity analysis methods [21] (e.g. analytical methods) are useful for analysing systems with few degrees of freedom and have well-known inputs and boundary conditions. These methods allow obtaining accurate results rapidly. By contrast, global sensitivity analysis methods [22] (e.g. Monte Carlo) are suitable for studying complex systems. The Monte Carlo method is a stochastic method adequate for analysing systems that have several degrees of freedom, have inputs with a significant degree of uncertainty, or have boundary conditions that cannot be determined easily. Finally, the Monte Carlo method is completely useful when no information is available about the internal distribution of the variables. However, this method requires a computer for its execution. In addition, owing to its stochastic nature, calculations could take longer than the analytical methods and the results obtained are not as accurate as deterministic methods. Moreover, the accuracy of the results depends on the number of analysed cases; that is, the more the cases studied, the better the results obtained.
In Atherton et al. [23], an analytical sensitivity analysis of the differential equations, which model chemical reactors, was developed; the paper presented the results of influence of each parameter. The study [24] presented an analysis of kinetic mechanisms based on the Fourier amplitude sensitivity test (FAST) method. This analysis allowed assessing the influence of kinetic variables on the predicted concentrations. Reference [25] presented a review of parametric sensitivities for several numerical techniques, including those with a stochastic character, such as curve fitting and derivative extraction; direct differential methods; green function method; the aim method; or finite differences. A recommendation of each analysed method was presented depending on the nature of the system. One of the conclusions of this study was that a sensitivity analysis should be conducted habitually when a kinetic model is used, because it is useful to measure the error propagation as it provides information about the mechanistic structure. The study [26] analysed the effect of the different variables used in temperature-programmed reduction (TPR) tests via deterministic methods, and subsequently, an estimation of the kinetic variables was conducted based on the results. Turányi [27] presented a review of manifold theoretical and numerical tools for sensitivity analysis (including deterministic and stochastic) and their application to the chemical kinetic area. Several sensitivity methods for chemical models, such as local sensitivity analysis, Monte Carlo methods, FAST, and Bayesian sensitivity analysis, were presented in [28]. Besides, it presented a methodology to accelerate the computation of the sensitivity indices employed in these methods. In Santos et al. [29], a parametric sensitivity analysis of the independent parallel reaction model was applied to study the pyrolysis of sugarcane bagasse. The study was dealt with by applying a software to solve and execute the sensitivity analysis of differential algebraic equations. In Batiot et al. [30], the influence of different parameters in the differential equation (Arrhenius law) used to model the mass loss of a solid material was studied, using the statistical Sobol approach.
The FDS software has been used to conduct several sensitivity analysis studies. Capote et al. [31] presented the results of the sensitivity analysis carried out. They evaluated the uncertainness of the input variables in using two different pyrolysis models: FDS and GPyro. In Suard et al. [32], a sensitivity analysis was carried out with several fire models, among which FDS was used to evaluate the physical input parameters such as the fuel, ventilation, or compartment as responses for fire safety studies. In the study, three alternative methods, Monte Carlo, fractionalfactorial, and full fractional, were applied for obtaining similar results between them. In Zhao et al. [33], the influence of parameters such as conductivity, specific heat, or emissivity on the pyrolysis behaviour of a mediumdensity fibreboard was investigated through a sensitivity analysis using FDS.
All these previous studies highlighted the necessity to know the effects of variation in the variables utilized to represent the thermal decomposition processes. This study aims to present the results of the sensitivity analysis carried out for assessing the influence of each input variable of FDS to model the decomposition process. Because many researchers have focused on modelling mass loss processes, this process will be considered as the prediction target. The material employed to carry out the study is cardboard.
This study is divided into the following three steps: The first step involves obtaining the experimental DTG curve. The second step involves the creation of the reference case, which is achieved using a combination of FDS software and a numerical approaching method. More details are provided in the corresponding section. Finally, the third step involves the sensitivity analysis, where each variable employed to model the thermal decomposition process is modified and compared with the reference case.

Experimental tests
To begin with, we need to obtain experimental data and ensure the repeatability; therefore, three STA tests were executed, following the ASTM-E1131-08 TG standard [1], under the same atmosphere and heating rate conditions. The samples of the solid cardboard, with a density of 526 kg m -3 and initial masses of 18.923, 14.380, and 13.174 mg, were previously conditioned in an environment with 50% of relative humidity (± 5%) at 23°C (± 3°C). After being chopped to obtain small pieces of milligrams, they were introduced into a 6-mm-diameter crucible made of Al 2 O 3 without lid, to allow decomposition under air atmosphere. The heating rate was 10 K m -1 from 30 to 800°C. The oxygen concentration during the test was 21% with a density of 1.293 (mg mL -1 ) and pressure of 1 atm. The inlet flow of the furnace and purge gas was 40 mL min -1 for all cases. Figure 1 shows the DTG (1st derivate of TGA) and DSC curves for the three tests.
The three tests had similar DTG curves. Usually cardboard is mainly composed of cellulose, hemicellulose, and lignin [34]. The decomposition of lignin occurs gradually between 160 and 900°C, the decomposition of hemicellulose occurs between 180 and 300°C, and the decomposition of cellulose occurs between 300 and 400°C [35]. These decomposition processes are shown in Fig. 1, where the first peak agrees with the decomposition of hemicellulose and the second one with the cellulose. Hereafter, the experimental curve refers exclusively to test 1. The experimentally obtained curves were similar to those obtained in other studies [36].

Creation of the reference case
Once the experimental data were available, the next step was to define a set of variables, used as input data, by using the pyrolysis model. This set was considered the reference case (Case 0). Pyrolysis and thermal oxidative decomposition processes occurred through a series of multiple reactions, which arose from the variation in the internal properties and structures of the material. Therefore, the reactions and the manner in which they occurred, i.e. the reaction scheme, should be developed.
To define the reactions, we employed the pyrolysis model used by the FDS software [10], which is based on the Arrhenius law [15]. Equation 1 represents the reaction rate, which is the rate of change of the mass fraction of the solid material i as a function of time in reaction j: where the term Y s;i is the quotient between the density of solid material i (q s;i ) produced by the reaction at temperature T s divided by the initial density of the material (prior to the reaction, q s 0 ð Þ). The term r ij is the reaction rate at temperature T s . Finally, the terms v si 0 j r i 0 j express the residue produced by the reaction r i 0 j with a yield of v si 0 j .
The value of the reaction rate r ij is calculated by the Arrhenius equation: Equation 2 contains the triplet kinetic, the pre-exponential factor A ij , activation energy Ea ij , the reaction order n s;ij , and the heterogeneous reaction order represented by term n O 2 ij . The term X n O 2 ij O 2 is used to simulate the effect of the oxidative reaction.
The term Y s;i can also be written as a function of a, the conversion factor of the reactant; in other words, the amount of mass involved in the reaction is expressed as: where m 0 is the mass of the sample at the beginning of the process, m is the mass of the sample at temperature T s , and m f represents the final mass of the sample, i.e. when the final reaction of the process is complete. Equation 2 can be represented as a function of a: The term f a ð Þ is the reaction mechanism [37]. The FDS software employs a reaction mechanism based on the reaction order: where n j is the reaction order n of the reaction j.
Once the mass loss rate is modelled, it is necessary to assess the energy released by the sample, which can be measured as: where _ q 000 s;c x ð Þ is the energy released or absorbed by reaction j and represents the heat of reaction (H r ) from reaction j and material i multiplied by the reaction rate and initial density. Further details about how FDS assesses the energy released in a decomposition process can be found in [10].
Having defined the reactions, the reaction scheme is needed to be determined, which depends on the tested material. The cardboard or any of its main components (lignin, cellulose, and hemicellulose) have been widely analysed. In the literature, different reaction schemes have been proposed to model the pyrolysis of cellulosic materials [38][39][40][41][42][43]. Given the results of [44], where the reaction scheme is concluded to be as simple as possible, but taking Temperature/°C Mass loss rate/% °C -1 Fig. 1 DTG-experimental curves into account a certain level of complexity that allows the model to represent the decomposition process, this study employed an approach similar to that proposed in [38], where the cardboard undergoes a four-reaction decomposition process (non-competitive reactions), creating three intermediate fictitious materials and a final residue. Each reaction produces a residue and releases burning gases, following the same idea as that proposed in [20,45,46]. In FDS, the generic reaction i can be defined as where M i represents the reacting material; Q represents the heat necessary to trigger the reaction (usually it is supplied by an external source, e.g. the STA furnace); O 2 represents the oxygen consumed by the reaction (in case of a reaction that consumes oxygen); P i represents the fictitious submaterial produced, which reacts if there is a subsequent reaction; F i represents the gas fuel released by the reaction; G i represents the non-burning gas released by the reaction (e.g. water vapour); and R i represents the residue produced by the reaction. The coefficients v i p , v i f , v i g , and v i r indicate the amounts of each product created by reaction i and are related to the conversion factor a. There is only nonburning gas in cardboard decomposition, where the water vapour is released in the first reaction at 100°C. Because a very small amount of water is present, v i g is considered to have a null value. It can also be considered that the residue R is produced only in the last reaction, as the final residue of the thermal decomposition process. To summarize, for the three first reactions, v p þ v f ¼ 1, and for the last reaction, v r þ v f ¼ 1. For example, in reaction 1, if the value m 1 p is equal to 0.6, 60% of material M 1 does not react, and hence, 40% is fuel gas (m 1 f ¼ 0:4). Consequently, in reaction 2, P 1 is 60% of M 1 .
Considering Eq. 7 and the reaction scheme proposed in [38], Eqs. 8-11 describe the reaction scheme used: Equations 4 and 6 summarize the variables needed to be introduced as input data in an FDS input file to define each reaction i (Eqs. [8][9][10][11]. Besides the reaction kinetics, the FDS input file requires to set the thermal properties of the material and intermediate fictitious materials. These thermal properties are conductivity (k), specific heat (C p ), emissivity (e), absorption coefficient (g), and heat of reaction (H r ). The density (q) is also included as the property of each material. The total number of variables required to model the pyrolysis process by employing the proposed reaction scheme, and including the thermal properties, is 41. Next, Fig. 2 shows the reaction scheme and variables employed.
In this study, once the reaction scheme is established, it remains invariable for all cases, because we do not aim to assess the accuracy of the reaction scheme.
In the FDS pyrolysis model, the volume of the simulated sample remains constant, that is, phenomena such as swelling and shirking are not taken into account by the model. Because of this feature, the mass is directly related to the density. In the FDS model, any variation in the mass loss is determined by the term a, as Eq. 4 shows. In other words, the mass of the sample at each moment of the thermal processes is proportional to the initial density multiplied by coefficients m i p , m i f , m i g , and m i r of each reaction i, i.e. the conversional factor a i .
Having identified the variables, it is necessary to achieve their values. Because most of the input thermal and kinetic properties of fire computer models cannot be obtained directly from thermal analysis tests, a numerical approaching method combined with a pyrolysis model will enable us to obtain the numerical values of all variables. This methodology has been extensively applied in the recent years to estimate the kinetic properties in an STA test in [20,47,48] and in bench-scale tests such as fire propagation apparatus (FPA) and cone calorimetric [8,49,50]. The most widely used mathematical methods include genetic algorithms [6,51,52], stochastic hill climber method [50,53], and shuffle complex evolution (SCE) [8,20,54]. In the present study, we selected SCE [55], according to [56], which is recommended for material pyrolysis property estimation. Next, Fig. 3 shows schematically how the approaching method was executed.  The numerical process requires the initial range for each variable. The numerical method selects one value of each variable, within the initial range. The iterative process reduces the error between the experimental and simulated curves in each loop. In this case, the variables are the kinetic and thermal properties and the error is the mean squared error (MSE-Eq. 12) between both DTG simulated and experimental curves. The process stops when, after five consecutive loops, the error does not decrease below 0.001%. When the numerical process stops, the DTG simulated curve obtained can be considered as the best approached one, obtaining the values of the variables. Next, Eq. 12 shows the MSE employed by the SCE method as a function to minimize its value: where k represents the total number of points measured, in this case 155 (from 30 to 800°C, each 5°C), and x p indicates the values of the simulated and experimental DTG at temperature l. After 22,400 iterations, the SCE method obtains the DTG simulated curve shown in Fig. 5, in Results section. The values that generate the curve are summarized in Tables 1 and 2. These values are employed as reference in the sensitivity analysis (Case 0).

Sensitivity analysis
Once the reference values of the variables are obtained (Tables 1 and 2), the influence of each can be analysed by a sensitivity analysis using a Monte Carlo methodology [57]. The Monte Carlo method is a non-deterministic method, i.e. the greater the number of random cases studied, the more accurate are the results. Monte Carlo is combined with the FDS simulation software, i.e. the results obtained owing to the modified values are produced after a simulation process. Therefore, it is necessary to establish a balance between the number of modified cases and the time required to process the results. To ensure that the results obtained by the mutated variables are feasible, we limit the range of variations up to ± 40% of its corresponding    Tables 1 and 2). Once the unfeasible values are discarded, each variable is modified or mutated 50 times, which seems to be enough to determine the influence of each variable. Two limits have particular properties because the values outside the range have no sense: emissivity has values between zero and one and the lower limit of the density is zero.
The sensitivity study has two phases: First, it carries out the simulation of each modified case serially. FDS takes less than a second to execute each simulation, producing as output file the DTG mutated curves. Following the fourreaction scheme, only one variable is mutated in each case and the remaining maintain a reference value of Case 0. These modified cases are known as mutated ones. In total, 2050 mutated cases are created and simulated, and the simulations are carried out using a Pentium 4 2 GB Ram. The whole simulation process took * 7 h. Finally, in the processing phase, the mutated DTG curves were compared with the reference DTG curve. This phase took * 8 h to process all cases and compare them. To compare the results and determine the influence of the mutated variable on the DTG curve, four parameters were compared: (1) the area below the DTG curves, (2) MSE, (3) the value of the largest DTG peaks, and (4) the temperature of the largest DTG peak. Next, Fig. 4 shows the sensitivity analysis process.
To measure the influence of each variable k in parameter s, we employed the normalized term b in s k of Eq. 13, which is a coefficient between the relative variation in the variable k and that in parameter s: in s k (between 0 and 1) is obtained by dividing the influence of each variable k over parameter s by the most influential parameter. The normalization of this term is helpful in establishing a comparison between variables because the most influential variables can mask the effect of the less influential ones.

Results
First, the results obtained by the SCE method used to achieve the reference values are shown. Figure 5 includes the DTG experimental curve used as a goal and the DTG simulated curve. The simulated DTG curve reproduces the DTG experimental curve from Test 1. According to Eq. 12, the MSE has a value of 0.0086. The results achieved by SCE in combination with the FDS pyrolysis model confirm the usefulness of the methodology and reaction scheme utilized in this study, as proved previously [6,20,47].
Next Tables 1 and 2 collect the values of the 41 variables achieved by SCE. Table 1 gathers the values of the variables related to the kinetics for all reactions (reaction 1 to reaction 4), according to the reaction scheme defined in Fig. 2. Table 2 shows the thermal properties of each material (material 1 to residue). With these input values, FDS can model the DTG curve (case 0) from Fig. 5.
Below, the results of the sensitivity analysis are presented. To show them properly and avoid enlarging the figures with the variables with limited influence, the figures do not show variables with an effect less than 0.05. Next, Fig. 6 shows the results of the four parameters analysed.
A priori, according to Eqs. 4 and 6, only A, E a , n, and H r could have influence in the modelling of thermal processes. Figure 6a, which shows the sensitivity to changes in DTG area, reveals not only their influence in thermal processes but also the participation of variables such as e and C p in the processes. In total, 33 of 41 variables influence the DTG area, which shows that this parameter is the most influenced by most of the variables. The reaction order n is the most important and has special influence in reactions 1, 3, and 4 described by the reaction scheme. Figure 6b presents the influence of the variables in DTG error, which was measured using the MSE shown in Eq. 12. Each mutated DTG curve was compared with the DTG curve of case 0. Unlike the DTG area, only six variables modify this parameter with an effect higher than 0.05: kinetic variables such as A, E a , and n. Figure 6c indicates the modification of the value of the largest DTG peak, i.e. the maximum value of the mass loss rate. Seven variables influence it, especially A and E a .  Fig. 6 Influence of variables over the parameters related with DTG curve: a area below DTG curve, b DTG MSE, c value of largest DTG peak, and d temperature of the largest DTG peak Figure 6d shows the effect of the variables within the temperature range at which the maximum peak of DTG occurs. This is the less easily influenced parameter by the variables. While the DTG area is affected by 33 variables, only 3 variables have influence in the temperature range of the DTG peak, i.e. A and E a of reaction 1 and E a of Reaction 2. Table 3 presents a general outlook of the influence of each variable in each parameter, where x indicates variables with an influence greater than 0.05.
In light of the results, only the triplet kinetics, i.e. A, E a , and n, has a considerable impact in the DTG curve, especially for the first reaction. The thermal properties have a unique influence in the DTG area, which means that by modifying k, C p , e, g, and H r , it is not possible to change the shape of the DTG curve or displace it; they only modify the size of the area enclosed by it.
As far as density (q) is concerned, because the mass loss rate involved in each reaction not only depends on the density but also the conversional factor (a), i.e. coefficients m i p , m i f , m i g , and m i r , its influence is limited. It only has effect in the DTG area and peak. Because this study aims to analyse the effects of the variables, not studying the reaction scheme, the reaction scheme of the mutated cases is identical to case 0. This characteristic limits the influence of density on DTG curve.

Conclusions
This study aims to analyse the role of the variables of the FDS software that allow the user to model cardboard thermal decomposition processes. To fulfil the objective, a sensitivity study based on a Monte Carlo method was carried out. All input variables for which FDS allows the user to configure the pyrolysis process were assessed by analysing how the DTG curve was modified.
In summary, by selecting a proper reaction scheme for the thermal decomposition of the cardboard and achieving the values of the triplet kinetics (A, E a , and n), a suitable fitting of DTG can be obtained using the FDS software. The less influential variables can be set up by employing reference values from the bibliography.
Owing to the high number of variables allowed by the FDS software to model the pyrolysis, it is highly recommended to carry out a sensitivity analysis prior to their estimation. The results presented in this study are helpful to discard, for pyrolysis processes that employed a similar reaction scheme (consecutive reactions), the variables without any influence. In case the reaction scheme involves parallel reactions, or a combination of parallel and consecutive reactions, it is recommended to conduct a similar analysis. The results obtained are only valid for the pyrolysis model utilized by FDS.