On the importance and modeling of in-depth spectral radiation absorption in the pyrolysis of

In-depth radiation absorption in pyrolysis models is commonly modeled with Beer’s law using a constant absorption coefficient. This approximation neglects the effects of spectral variation of absorption coefficient, material emissions, and directional dependency of penetrated thermal radiation. This research aims to address the importance of these factors in flammability analysis of black poly(methyl methacrylate) (PMMA). Accordingly, a pyrolysis model including the separated full spectrum correlated k-distribution (SFSCK) method for spectral modeling, the ordinate weighting method (OWM) with finite volume method (FVM) for directional distribution of the thermal radiation is introduced as the most detailed model. By gradually simplifying this model, five other pyrolysis models are developed with different levels of details in radiation calculations. Simulations are first done using the optimized thermal conductivity and specific heat capacity for the most detailed model, then by optimizing these parameters for each pyrolysis model separately. The results show that the pyrolysis model with a depth and source temperature-dependent absorption coefficient, two-flux method, and diffuse boundary with a scenario-dependent effective reflectivity can predict the measured mass loss rate and ignition time with accuracy similar to the most detailed model, and can be recommended for future pyrolysis simulations.


Introduction
Pyrolysis models are used for simulating the high-temperature thermal decomposition of condensed materials by solving the numerical equations describing the heat and mass transfer with simultaneous chemical reactions.The commonly included mechanisms in most of the presented pyrolysis models for non-charring materials are heat conduction, decomposition kinetics, material shrinkage with decomposition [1], and in-depth radiation absorption [2,3].In reality, other mechanisms, such as the heat exchange between gas flow and solid material, bubble formation, source irradiation absorption by pyrolysis gas and bubbles, temperature-dependent evolution of decomposition kinetics and material absorption, and interaction with oxygen may also be involved in the thermal decomposition.Including these mechanisms in the formulation by considering e.g.,permeability, porosity, and pressure field of the generated gas, more complex pyrolysis models have been reported in literature [4,5].Although more complex models are conventionally assumed to provide higher accuracy, their uncertainty increases with their complexity [6].Therefore, there should be a compromise between the accuracy and complexity of pyrolysis models.
The simplifications of pyrolysis models are due to the lack of measured material properties, incapability in measuring the actual material property independently, the anticipated small effect of some mechanisms compared to the uncertainty in experimental data [7,8], and lowering the computational time.To assess and propose appropriate simplified models, sensitivity analyses for the input parameters [9,10], and the considered mechanisms [11] are usually performed.To validate the pyrolysis models, their results are compared with experimental data of mass loss rate or heat release rates, sometimes temperatures.Due to the compensation effect in parameters' estimation, even the simple pyrolysis models may capture the essential features of the experimental data.To avoid significant compensation effect, the pyrolysis models should solve the important physics of the processes without adding new unknown parameters.
In-depth radiation absorption is one of the mechanisms that may significantly affect the pyrolysis of non-charring polymers [2,3].In most of the presented pyrolysis models, in-depth radiation absorption is simplified to Beer's law with constant absorption coefficient [8,12].This method ignores at least three possibly important features of thermal radiation: directional and spectral dependencies, and medium emissions.The directional distribution of the penetrated thermal radiation has not been addressed in literature yet.Even though the Beer's law has served well in many pyrolysis models, it is theoretically valid only https://doi.org/10.1016/j.firesaf.2022 for single beam of radiation, not for heat flux, which is always threedimensional phenomenon.Besides, at the interface of the material with air, due to the difference in refractive indexes, redistribution of thermal radiation may occur [13].
The spectral nature of thermal radiation is commonly taken into account in gas phase simulations [14,15], and various models have been presented, such as weighted sum of gray gases (WSGG) [16,17], spectral line based weighted sum of gray gases (SLW) [18,19], and full spectrum correlated-k distribution (FSCK) [20,21].In case of condensed-phase pyrolysis, the spectral modeling has only been reported in [22][23][24]  in which spectral radiative transfer was modeled to extract the preignition temperature profile within the sample [25,26].Using the measured absorption coefficient of clear PMMA [26] and applying the two-flux method, Kacem et al. [23] showed that spectral modeling of in-depth radiation leads to more accurate predictions of the mass loss rate and temperature profile, compared to the gray modeling.Applied radiative heat fluxes in their experiments were 14 and 18 kW∕m 2 , which are near the critical heat flux of ignition for clear PMMA.At higher external heat fluxes, the in-depth radiation absorption is anticipated to be even more important for the ignition time and mass loss rate.The optimal level of simplifications for thermal radiation details in pyrolysis modeling has not yet been determined.To accurately estimate the in-depth radiation absorption, the spectral and directional aspects of penetrating thermal radiation through condensed materials and the interface effect should be precisely taken into account.In the current research, thermal radiation modeling is applied to the flammability prediction of condensed materials.Therefore, the optimal level of details may be a compromise between the computational cost and the obtained, additional accuracy of mass loss rate and ignition predictions.However, adding new details may introduce new unknown parameters in modeling and therefore may increase the compensation effect [6].We therefore aim to improve the accuracy of estimated in-depth radiation absorption without adding uncertainty in terms of unknown parameters.
In the present work, in-depth radiation is modeled with great detail, and its effect on the prediction of pyrolysis models is investigated to suggest the most efficient simplified model.Six different pyrolysis models are developed, starting from the most detailed model and gradually simplifying the models.The most detailed model is based on the recently developed FSCK method for the in-depth radiation in liquid fuels [27].We later extended this model for solid materials, where material temperatures can by much higher, by separating the solution of radiation intensity originating from the medium emission from that of the boundary irradiation [28].The SFSCK model was shown to provide higher accuracy with a modest increment in calculation time [28].The simplified, gray model for the spectral absorption coefficient is based on the depth and source temperature-dependent effective absorption coefficient that we have developed and tabulated using the measured spectral absorption coefficient of black PMMA in the wavelength region of 0.25 to 25 μm [29].To include the interface effect, caused by the different refraction indexes of two phases, Alinejad et al. [13] presented the ordinate weighting method (OWM), which accurately and efficiently includes directional alternation of thermal radiation at the interface.To evaluate the performance of the models, they are used for predicting anaerobic thermal decomposition of black PMMA and the ignition time.The parameter optimization using the cone calorimeter data [30,31] is only carried out for thermal conductivity and specific heat capacity.First, these parameters are estimated only for the most detailed model, and the predictions of different models with this set of data are compared to reveal the effects of the simplifications.Then, thermal conductivity and specific heat capacity are separately estimated for each pyrolysis model to determine the ability of different models in predicting the experimental data through compensation effects.Finally, the most efficient model for the indepth radiation transfer is proposed based on simplicity, accuracy, and computational cost.

Model descriptions
To study the flammability of black PMMA by including various levels of details in thermal radiation modeling, six different pyrolysis models are presented in Table 1.
Due to applying non-constant absorption coefficient for gray method, model 6 in Table 1 is not identical with the most widely used pyrolysis models in literature with a constant value for absorption coefficient.Except for the thermal radiation related sub-models listed in Table 1, models 1 to 6 share the other equations of the pyrolysis model given in Sections 2.1 and 2.2.The simulations are done to study the thermal decomposition of black PMMA in cone calorimeter experiments with the structure shown in Fig. 1.
The indixes 0, 1, and 2 in Fig. 1 are applied for air, black PMMA, and backing layer, respectively.

Energy conservation equation
By considering the in-depth thermal radiation absorption and energy consumption due to the decomposition reactions as source terms and neglecting the heat transfer between gas products and the condensed material, the energy conservation equation for the black PMMA layer is written as: where  , ,   , and   are specific heat capacity, density, and thermal conductivity of black PMMA, respectively and  , q′′′  ,   ,  and ℎ , are temperature, radiative heat source, reaction rate of component , considered number of reactions (components), and heat of pyrolysis for the component  of black PMMA, respectively.For the upper boundary condition of Eq. ( 1) convective heat transfer is applied for models 1 to 6.For model 6, due to the lack of applying RTE, the thermal radiation emission from the surface is added to the boundary condition.For the bottom boundary, the equality between conducted heat from PMMA to the interface and from the interface to the backing layer is applied.
The assumed material of the backing layer is Insulfrax ® LTX TM blanket [32].Neglecting the in-depth radiation absorption and decomposition of this material, the energy conservation equation for the backing layer is written as: where  , ,   , and   are specific heat capacity, density, and thermal conductivity of backing layer, respectively.The adiabatic boundary condition is applied for the bottom boundary of the backing layer.Doing a polynomial fitting for the given data in [32], the thermal conductivity of the backing layer is given as:

Decomposition kinetics
Using the TGA data of black PMMA at heating rate of 10 K∕min, three components of humidity (water),  − PMMA, and  − PMMA with mass fractions of  1 ,  2 , and  3 , respectively are considered for the raw black PMMA.The assumed reaction mechanisms are: Using the Arrhenius law, the following equation gives the mass loss rate per unit volume or reaction rate (  ) for each component: where   ,   , , and  ′′′  are pre-exponential factor, activation energy of reaction , universal gas constant, and instantaneous mass per unit volume of the component , respectively.
Summing up the reaction rate of the components, the total mass loss rate per unit volume is given as: The mass loss rate per unit area is obtained by integrating Eq. ( 8) over the sample thickness.It is worth mentioning that using a single heating method may impact the estimated kinetic parameters [33].Consequently, applying a fixed set of values for the kinetic parameters leads to a compensation effect from the optimized parameters.However, the investigation of this issue and showing its effect in pyrolysis modeling is not in the scope of the current article.

Spectral modeling of thermal radiation
To calculate the radiative heat source in Eq. ( 1), spectral and directional aspects of thermal radiation should be taken into account.The strength of materials in absorbing and emitting radiation varies with wavelength and this makes radiative heat transfer a spectrally dependent phenomenon.Therefore, for obtaining an exact solution of thermal radiation, the RTE should be solved for very narrow wavelength intervals in which the radiation can be assumed monochromic, and the absorption coefficient could be considered constant.Solving the RTE using this method at each time step, makes the simulations computationally demanding.For the spectral aspect, two methods of gray and SFSCK [28] are applied.

Gray method
One of the simplest approaches for the high computational cost of spectral calculations, which has been popularly used in fire simulations, is applying an averaged value for the spectral absorption coefficient.This approach is called gray method in which the RTE needs to be solved only once: where ,   , , and  ef f are the total intensity, blackbody total intensity, direction vector, and averaged (or effective) value of absorption coefficient, respectively.The measured mean value of absorption coefficient showed a dependency to the heater type and the sample thickness in experiments [34,35].To consider the mentioned effects, Alinejad et al. [29] proposed a new spectrally averaged value of absorption coefficient (i.e.,  ef f ) varying with source (heater) temperature ( Sc ) and the distance from the surface of the sample ().Accordingly, using the presented look-up table in [36],  ef f in Eq. ( 9) is determined for each computational cell based on the heater (or flame) temperature and the distance from the interface (depth).We determined the values of  ef f for models 3 to 6 of Table 1 applying the provided table for the averaged absorption coefficient of black PMMA in [36].

SFSCK method
To efficiently include the spectral nature of absorption coefficient in solution of thermal radiation in condensed materials, Alinejad et al. [27] adopted the FSCK method previously presented for the gas phase [37].In this method, a highly fluctuating spectrum of the absorption coefficient is converted to a monotonically increasing diagram that is easy to work with.Picking a few points from the new diagram of absorption coefficient (k values in g space) and applying a quadrature method, the radiative heat source is determined with good accuracy and lower computational time.Later investigations showed that this method cannot preserve the medium emission accurately [38].Therefore, the SFSCK method was developed by the same authors [28] for the accurate estimation of both medium emission and absorption.It is based on separating medium and incident radiation originated intensities.According to this method, the intensity is first divided into two components originating from the domain boundary (  ) and medium emission (  ).The RTE to obtain   is: (10) and the following equation is for   : where  is the reordered wavelength and  and  are the parameters of SFSCK method and are determined from the look-up table given as the supplementary material to this article.

Solution methods of RTE
Three different methods are used for solving the directional dependence of intensities arising from the RTE of Eqs. ( 9), ( 10) and (11).

Beer's law
The Beer's law normally is applied for the intensity but in fire research it is applied for determining the heat flux distribution [39] inside the condensed materials.This method gives the following equation for the radiative heat flux ( q′′ rad ): where q′′ in , ρ01 , and  ef f are the incident heat flux at the sample surface, the averaged reflectivity at the surface of the sample as shown in Fig. 1, and the effective absorption coefficient [29] as described in Section 2.3.1.The radiative heat source for the energy conservation equation is then determined by a derivation of Eq. ( 12):

Two-flux method
This method considers only one ordinate for each of forward and backward thermal radiation transports.The RTE for the forward direction is: For the backward direction, RTE is written in a similar way [28].
Eq. ( 14) needs a boundary condition at the sample surface that is given in Section 2.5.For the backward direction, the boundary condition is: where ρ12 is the averaged reflectivity at the interface of the black PMMA and the backing layer as shown in Fig. 1.After calculating  + and  − , radiative heat flux is given for the forward and backward directions by ( q′′ rad ) + =  + and ( q′′ rad ) − =  − , respectively.Finally, the radiative heat source is determined as:

FVM
In this method, forward and backward directions are discretized to some control angles.To solve the RTE for one-dimensional geometry, the FVM gives the following discrete form: where  , ,   1, , and   2, are the intensities at the node , the upper and the bottom faces of the control volume , respectively.The boundary condition for Eq. ( 17) is given in Section 2.5.To determine the face intensities, the STEP scheme is applied in the present research.The subscript  is the index of control angle as shown in Fig. 2. The parameter  is the grid size and   is the size of control angle .
Considering  divisions for the polar angle and  divisions for the azimuthal angle (i.e.,  =  × ), the control solid angle   is calculated as: where  shows the index of control angle and defined as  =  × .In Eq. ( 17),   , is defined for the control volume faces ( = 1, 2) as: where n is the unit normal vector of face j of the control volume i.The detailed derivation of this method can be found in [40].To apply this direction method with SFSCK method, radiative heat source is calculated as: To calculate the radiative heat flux from Eq. ( 20) using  control angles to cover the whole 4 solid angles and applying Gauss-Chebyshev quadrature method, the following equation is obtained: where  shows the number of applied quadrature points that is 8 in the present research and   is the weight parameter of th quadrature in Gauss-Chebyshev quadrature method.More details about this method can be found in [41].
In case of gray method, radiative heat source is calculated as:

Radiation treatment of the upper boundary
To solve the RTE using FVM, a proper boundary condition for the upper surface of the black PMMA should be given.In the present research, two types of boundaries are applied: Diffuse and Fresnel boundaries.Fig. 3 shows the difference of these two boundary types.

Diffuse boundary
In a diffuse boundary, incident intensity and internal reflection is distributed equally among the different control angles that is described as: where  +  ,  −  , ρ01 , and ρ10 are forward radiation intensity within control angle , backward radiation intensity within control angle , the averaged reflectivity at the interface from the air side to the black PMMA side, and the averaged reflectivity from the black PMMA side to the air side (see Fig. 1), respectively.Reflectivity is a function of the polar angle, but for the diffused boundary, an average value of reflectivity is applied (see Eq. ( 36)).To apply a realistic value for ρ01 , the following equation is used: 0 cos() sin() (24) where  01 () is the reflectivity function and is calculated using the Fresnel relation as: In Eq. ( 25),  0 and  1 are related together through the Snell's law: Using the Eq. ( 24), the calculated values of ρ01 for gasification and flaming conditions in cone calorimeter are 0.05 and 0.089, respectively.Similar equations are used to calculate ρ10 .

Fresnel boundary
In the Fresnel boundary, the incident intensity is not distributed equally among the control angles, and redistribution of radiation intensity is done at the black PMMA side of the interface.The difference refractive indexes of PMMA and the gas-phase alter the distribution of the intensity passing through the interface.To calculate intensity within each control angle at the PMMA side of the interface, the OWM [13] is applied that gives the intensity within each control angle ( +  ) as: where  is the number of considered divisions for polar angle at one side of the interface and   is the weighting parameter for the incident intensity coming from control angle  at the first side to the control angle  at the second side of the interface and is defined as:  (28) where  , and  ,+1 are the corresponding angles of   and  +1 after transmission and are given using Snell's law.The averaged reflectivities ρ01, and ρ01, in Eq. ( 28) are defined as:  (30) where  01 is calculated using Eq. ( 25).Similar calculations are done for ρ10, in Eq. ( 27).In Eq. ( 27),  , is the incident intensity within the control angle  and is calculated as: where  and  ℎ are the emissivity and the temperature of the heater, respectively.The parameter   is the weight of the incident intensity for the control angle  that is calculated as: cos() sin() ∫  +1   cos() sin() (32) where  ℎ,1 and  ℎ,2 are the boundary angle of the incident intensity and   and  +1 are the boundaries of the control angle  as shown in Fig. 4.

Experimental data
In this research, the experimental data of the MaCFP project [31] is applied for the parameter estimation and validation of the simulation results.These data are the output of an international collaboration providing measurement data of black PMMA in the 2021 IAFSS-MaCFP Condensed Phase Workshop [30].The applied experimental data of MaCFP in current research are TGA and DSC for the heating rate of 10 K/min, mass loss rate, and ignition time at heat fluxes of 25, 50, and 65 kW∕m 2 .For TGA, the measured data by the University of Lille and for DSC, the measured data by the University of Maryland are applied.Table 2 summarizes the mass loss rate and heat release data used in the present research.The applied black PMMA samples used by different research groups are from the same manufacturer [30].

Model calibration by parameter optimization
To simulate pyrolysis using the presented models of Table 1, the material properties and kinetic parameter should be known.These parameters are estimated by performing an optimization with the where  is either the mass fraction of TGA data, or the normalized heat flow of DSC data, or the mass loss rate data of anaerobic experiment with cone calorimeter.In Eq. (33),  shows the number of selected points from the experimental data.
To model the chemical kinetic part, three reactions are assumed for the thermal degradation of black PMMA that gives 8 independent unknown parameters to estimate: , and  2 .For inert environment, using the TGA data for the heating rate of 10 K/min mentioned in Section 3.1, these 8 parameters were estimated and to estimate the heats of reactions in Eq. ( 1) (ℎ ,1 , ℎ ,2 , and ℎ ,3 ), the DSC data for the heating rate of 10 K/min mentioned in Section 3.1 has been applied.Finally, for the pyrolysis simulations, we first determined a set of optimized data for  , and   using the most detailed model (i.e., model 1).This parameter estimation aims to reveal the effect of different simplifications in estimation of in-depth radiation absorption.Then, optimized  , and   are determined for each pyrolysis model separately, aiming to address the compensation effect in capturing the experimental data and the comparison of estimated parameters.For the parameter estimation, the mass loss rate of black PMMA exposed to a constant heat flux in a cone calorimeter device with an anaerobic environment is applied as the target of the optimization.By using the data of anaerobic cone calorimeter experiments, the complexities related to the combustion in the gas phase are avoided.In-depth radiation absorption shows its effect substantially at high heat   [36] ρ01 (Cone) 0.050 (Eq.( 24)) ρ10 0.584 (Eq.( 24)) ρ01 (Flame) 0.089 (Eq.( 24)) ρ12 1.0 (Assumed) a The superscript  shows the inert environment and the superscript  shows the air environment.fluxes [8].Hence, the given experimental results for the heat flux of 65 kW∕m 2 in Table 2 is selected as the target test case for the optimization of pyrolysis models.The applied parameters except  , and   in different pyrolysis models are summarized in Table 3.

Model testing
To check the performance of presented models in Table 1, all the listed experimental data in Table 2 and the measured ignition times are applied.Two simulations sets have been done to predict the mass loss rate of black PMMA for the gasification experiments at the heat fluxes of 25, 50, and 65 kW∕m 2 .First, the simulations of the six models have been done by applying the estimated  , and   of the most detailed pyrolysis model (i.e., Model 1), and the results are marked with the letter ''".Second, the simulations of six models have been conducted using the optimized  , and   for each pyrolysis model separately, and the results are marked with the letter ''".Through these two sets of simulations, we address the efficient level of detailed thermal radiation modeling needed for pyrolysis modeling of black PMMA.To predict the ignition time using the presented models with their estimated  , Fig. 5. Mesh-independence study using model 1 and heat flux of 65 kW∕m 2 .and   , the critical mass flux of 0.003 kg∕m 2 ∕s has been applied as the ignition criterion of black PMMA.The critical mass flux is the most widely used criterion in literature [45].The predicted ignition times are compared with experimental data [31] for the heat fluxes of 65, 50, and 25 kW∕m 2 .

Results and discussion
To do the simulation for different pyrolysis models, we prepared MATLAB codes.The finite-difference method was used to solve the energy conservation equation.To solve the system of the discretized equations, we applied the tridiagonal matrix algorithm.The grids size can vary according to Eq. ( 8).The threshold to remove a grid is reaching 0.1% of its initial size.The applied time step is 1 s for all the pyrolysis methods.The RTE and energy conservation equations are coupled and, therefore, an iterative method is applied to solve both energy conservation equation and RTE.The criterion for the convergence of the energy conservation equation is: F. Alinejad et al.where the superscripts show two consecutive iterations, and the subscripts show the grid index.The convergence criterion for the RTE is given as [46]: where  is the irradiation, and is defined as:

Mesh independence study
To check the mesh-independence of the results, three simulations are done for the pyrolysis model 1 using different number of grids for anaerobic decomposition of a black PMMA layer in cone calorimeter at a heat flux of 65 kW∕m 2 .In the pyrolysis models, the thermal radiation equation is solved only within the black PMMA due to its strong radiation absorption.Therefore, the mesh-independence study is done only for the black PMMA layer.For the backing layer, a computational grid with 100 cells is applied in all the simulations.The difference between the results of different grids, as shown in Fig. 5, is less than 2% which is mostly related to the last stage of decomposition when most of the grids are removed during the simulations.Accordingly, all the simulations have been done using 60 grids for the black PMMA layer.

Anaerobic decomposition of black PMMA
To simulate the black PMMA decomposition at various heat flux levels in gasification experiments, estimating  ,BP and   has been done using the experimental data for the heat flux of 65 kW∕m 2 .The estimated properties have been then used in the simulations for the heat fluxes of 25 and 50 kW∕m 2 .The experimental data used in the present research have been provided by different research groups participated in IAFSS-MaCFP experimental campaign and their uncertainties are given by one standard deviation (i.e., ṁ′′  ± ).Model 1 represents the most detailed and complicated radiation modeling in pyrolysis calculations.To investigate the effects of simplifying thermal radiation in pyrolysis modeling, some simulations of all the models have been done using the estimated parameters of model 1 as shown in Fig. 6.The obtained results showed that using the same material property, all the pyrolysis models except model 6 produce fairly similar results.Model 6 overestimates the experimental data  Comparing the results of model 1 and model 3 shows that using the proposed gray method with the mean absorption coefficient, which depends on the depth and the heater temperature, could produce the mass loss rate similar to application of SFSCK method.Also, comparing the results of model 3 and model 4 shows that the proposed way of calculating the averaged reflectivity (i.e., Eq. ( 24)) well approximates the Fresnel effect at the boundary.Finally, comparing the results of model 4 and model 5 shows that using two-flux method as the RTE solver could produce fairly similar predictions of mass loss rate to FVM.It suggests the sufficiency of two-flux method for solving radiative heat transfer in pyrolysis modeling.
The obtained heat source distributions inside the material by different pyrolysis models are shown in Fig. 7 for the heat flux of 65 kW∕m 2 at the beginning of experiment and the time of 100 s.The heat source distribution at the beginning of the experiment shows the estimation of the models for the medium absorption of source irradiation.The heat source distribution at the time of 100 s shows the estimation of the models for both medium absorption and medium emission.Compared to the model 1, model 6 overestimates the heat source.Therefore, applying the  ,BP and  BP of model 1, models 6 overestimates the mass loss rate of the black PMMA as shown in Fig. 6.For the models 2 to 5, the obtained heat source distributions are close to that of model 1 as they overestimates the heat source at some depths and underestimates the heat source at other depths.Therefore, using the same values for  ,BP and  BP , their results are close together.
For the next set of simulations, parameter estimations were done for each pyrolysis model separately using the experimental data of 65 kW∕m 2 .Model specific  ,BP and  BP were then used to simulate gasification at heat fluxes 25 and 50 kW∕m 2 .The results in Fig. 8 show that all the pyrolysis models could capture the experimental data with their own optimized  ,BP and  BP .The models could not capture the last phase of the gasification due to simplifying the 3D experimental sample to a 1D model.In the experiments, different parts of the sample are decomposed at different rate compared to the sample's center.Therefore, applying the heat flux of the center for whole surface could not produce the exact total mass loss rate from the sample even with optimizing the  ,BP and  BP .Thermal radiation absorption by the pyrolysis gases reduces the radiative heat flux at the surface of the material [47], and also bubbles layer absorption changes the in-depth radiation absorption for the black PMMA [48].We neglected these details in our modeling that could be a reason for the inaccurate prediction of mass loss rate at this part.Another potential reason is neglecting the permeability of the porous backing layer that could absorb the molten PMMA.It provides slower pyrolysis compared to the ideal 1D pyrolysis.
The heat source distributions inside materials obtained by different pyrolysis models are shown in Fig. 7 for the heat flux of 65 kW∕m 2 at the time of 100 s.For the beginning of the experiment, the heat source distribution is the same as what is already shown in Fig. 7.As seen, model 6 overestimates the heat source compared to the model 1.
However, the differences in heat source predictions are compensated applying the different estimated  ,BP and  BP optimized for each pyrolysis model separately.
To assess the computational costs of different models, some simulations were done to predict anaerobic decomposition of black PMMA at heat flux of 65 kW∕m 2 using MATLAB codes with similar programming logic for different pyrolysis models on a same computer (Intel core i7 6500U processor).The obtained results are shown in Table 4.As model 1 considers the details of thermal radiation, it increases the simulation time considerably.As shown in Table 4, adding the spectral calculations to the pyrolysis models is computationally expensive.In conclusion, while providing more simplicity and lower computational cost, model 5 efficiently and accurately accounts for in-depth radiation absorption.
The estimated  ,BP and  BP optimized for different models are shown in Fig. 10.Assuming a linear relation for change of  ,BP and  BP with temperature, four fitting parameters should be obtained in each optimization.According to Fig. 10, the estimated values of  ,BP and  BP for the models 1 to 5 are close together, but the estimated values for model 6 are remarkably different, especially for  ,BP .It shows that simplifying the thermal radiation in model 6 has been compensated with quite different values of material thermal properties.The reported values of  ,BP using DSC measurements in [39] indicates that the estimated  ,BP for model 6 could not be realistic.In summary, using optimization for parameter estimation in pyrolysis models compensates the neglected or simplified mechanisms through providing different set of parameters which are dependent on the applied pyrolysis models.
The  ,BP and  BP optimized for different models are given in Table 5.

Ignition time prediction
To predict ignition time, a set of kinetic parameters were estimated using the TGA data of air environment for the heating rate of 10 K∕min [49], and were applied only for the first grid.The kinetic parameters for the air environment are given in Table 3.It is worth mentioning that the effect of oxidative pyrolysis for the black PMMA is negligible [50].The experimental data of ignition time for the heat fluxes of 25 and 65 kW∕m 2 show higher uncertainty even when the outliers data are excluded.For the heat flux of 50 kW/m 2 , there is only one reported experimental data.The predictions of different models   at different heat fluxes are shown in Fig. 11.As seen, models 1 to 5 could capture the experimental data for the heat fluxes of 50 and 65 kW∕m 2 , but they overestimate the ignition time at the heat flux of 25 kW∕m 2 .Model 6 exhibits different behavior in predicting the ignition times compared to other models.

Sensitivity analysis
While there is a good estimation of the reflectivity at the surface of the sample in literature [2], there is not any estimation for other reflectivities shown in Fig. 1.In the present research, the reflectivities ρ10 and ρ12 are used for the first time in pyrolysis modeling, and therefore, their effects in prediction of mass loss rate are investigated.The models 1 and 3 treats the reflectivity at interfaces accurately with considering Fresnel boundaries, but for models 2, 4, and 5, an averaged value of reflectivity is calculated using Eq. ( 24) in the diffuse boundary condition.Several simulations with model 4 for the heat flux of 65 kW∕m 2 were done applying different values for the reflectivities as shown in Fig. 12.The results show that the entire mass loss rate curve is quite sensitive to the value of ρ10 , but, ρ12 shows its effect only on the peak region of the mass loss rate curve.Due to the strong radiation absorption of black PMMA, thermal radiation is mostly absorbed within a very thin layer as shown in Figs.7 and 9. Therefore, the sensible effect of ρ12 is appeared when the layer is thin enough.In the current research, the value of 0.584 for ρ10 has been applied.However, for ρ12 due to the lack of radiation properties for the backing layer and similarity of the obtained peak region to the experimental data, the value of unity has been applied.

Conclusion and remarks
In this research, implementing the most up-to-date advances of thermal radiation modeling in condensed materials, we investigated the effect of detailed thermal radiation modeling on the prediction of black PMMA flammability characteristics by pyrolysis models.Six different pyrolysis models were developed, representing different levels of details in the thermal radiation calculations.The simulation results were compared with measured ignition times and mass loss rates in anaerobic decomposition at heat fluxes of 25, 50, and 65 kW∕m 2 .The following items summarize the main conclusions of this research: • The simplest model overestimated the heat source inside the material, compared to the most detailed radiation calculations.This can be due to the oversimplification of radiation heat transfer, including spectral variation of absorption coefficient, directional distribution of intensities, and skipping medium emission.• The most detailed spectral modeling of in-depth radiation was done by the SFSCK method [28].A simpler method to account for the spectral nature of thermal radiation was a gray method with non-constant mean absorption coefficient which depends on the source temperature and the depth.The results showed that this approximation works well for the pyrolysis models considering the scale of experiments uncertainty.• To account for the directional distribution of penetrated thermal radiation, FVM and two-flux method were used.The higher number of ordinates in FVM had only negligible effect on pyrolysis modeling.Hence, the two-flux method can be chosen as an optimal model for solving thermal radiation in condensed materials.• The effect of the air-material interface was included by OWM for the assumed Fresnel interface.To simplify, a diffuse boundary condition with averaged reflectivity was developed.Although the results of this simplified approach showed an overestimation of heat source close to the interface, the effect was weak in mass loss rate predictions.Therefore, this simplified approach is suggested as a proper approximation of Fresnel boundaries.• Results of the simulations and the sensibility of the estimated optimized values of  , and   for different models indicated that the pyrolysis model with gray medium and two-flux RTE solution (i.e., Model 5) produces results that are similar to the

F
.Alinejad et al.

Fig. 4 .
Fig. 4. Angular span and boundary angles for the cone heater over the black PMMA sample.

Fig. 6 .
Fig. 6. Results of anaerobic decomposition simulation for different heat fluxes using estimated parameters of Model 1.

Fig. 7 .
Fig. 7.The heat source distribution inside the material exposed to 65 kW∕m 2 by different pyrolysis models using the  ,BP and  BP optimized for model 1 at the times of left) 0 s and b) 100 s.

Fig.
Fig. Results of anaerobic decomposition simulation for different heat fluxes using each model specific estimated parameters.

Fig. 9 .
Fig.9.The heat source distribution inside the material exposed to 65 kW∕m 2 obtained by different pyrolysis models at 100 s using the  ,BP and  BP optimized for each model separately.

F
.Alinejad et al.

Fig. 11 .
Fig. 11.Comparison of predicted ignition times by different pyrolysis models and measured ignition times for different heat fluxes.

Table 1
Defining different pyrolysis models.

Table 3
The values of different parameters used in simulations.

Table 4
Simulation time of different pyrolysis models predicting anaerobic decomposition of black PMMA at heat flux of 65 kW∕m 2 .