Application of Particle Swarm Optimization (PSO) Algorithm in Determining Thermodynamics of Solid Combustibles

: The thermodynamics of a solid are crucial in predicting thermal responses and ﬁre behaviors, and they are commonly determined by inverse modeling and optimization algorithms at constant heat ﬂux. However, in practical scenarios, the imposed heat ﬂux frequently varies with time, and related thermodynamics determination methods are rarely reported. In this study, the particle swarm optimization (PSO) algorithm and a 1D numerical model were utilized to determine temperature-dependent thermal conductivity and speciﬁc heat of beech wood and polymethyl methacrylate (PMMA). Surface, 3 and 6 mm in-depth temperatures were measured in three sets of ignition tests where constant and time-dependent heat ﬂuxes (HFs) were applied. In each set, PSO was implemented at individual HFs, and the average value was deemed as the ﬁnal outcome. Reliability of the optimized thermodynamics was veriﬁed by comparing with the reported values in the literature and predicting the experimental measurements that were not employed during parameterization. The results showed that wood thermodynamics attained under constant and time-dependent HFs in agreement with previously reported ones. Similar optimization procedures were conducted for PMMA, and good agreement with literature values was found. Using the obtained thermodynamics of wood under constant HF, the numerical model successfully captured the surface temperature at time-dependent HFs. Meanwhile, comparisons using wood temperatures at constant HFs and PMMA temperatures at linear HFs also veriﬁed the feasibility of PSO.


Introduction
Pyrolysis of solid combustibles frequently endangers personnel and environment due to its potential to cause fires. Pyrolysis of solids is closely related to their thermodynamics, such as temperature-dependent density (ρ), thermal conductivity (k) and specific heat (C p ). Therefore, accurate determination of these parameters is of great importance in assessing heat transfer and predicting pyrolysis and combustion behaviors [1]. For this reason, many existing studies have investigated the material pyrolysis, thermodynamics, ignition criteria and combustion process using experimental, analytical and numerical approaches.
Direct measurement of solid thermodynamics can be achieved by experimental methodologies. For instance, k can be measured by steady-state and transient techniques. Steadystate techniques require the measurement to be implemented in complete thermal equilibrium. These techniques generally need a very long time to reach required equilibrium and cannot be used for moist materials due to moisture migration. Transient techniques perform the measurement during a process of small temperature change and can be carried out relatively quickly. The transient plane source (TPS) method and non-steady-state hot wire method are the most used transient techniques. All these techniques postulate that no physical or chemical change occurs during measurement, and the temperature needs to be maintained at relatively low value ranges. Nevertheless, in fire research, solid fuels frequently suffer very high temperature, implying k of the substances cannot be determined by these techniques directly. As a result, indirect methods by inversely modeling benchscale experimental results have gained considerable popularity in recent years. With rapid development of computer technology, several heuristic optimization algorithms, such as genetic algorithms (GAs) [2,3], particle swarm optimization (PSO) [4], shuffling complex evolution (SCE) [5][6][7], hill climbing (HC) [8] and multi-algorithmic genetic adaptive multiobjective methods (AMALGAMs) [9], have been applied to automate the inverse modeling process with coupled numerical models.
In Gong's [10] work, GA was used to parameterize the hybrid pyrolysis mechanism. Ding [11] used GA and PSO to estimate the kinetics of a three-component parallel reaction mechanism. The accuracy and efficiency of these two algorithms were compared. A significant advantage of PSO over GA was the high convergence efficiency and accuracy compared to GA, but PSO may obtain locally optimal solutions for complex problems and blind search ranges [12]. Prathana [13] used a GA, PSO and the gradient-based method (GBM) to evaluate and quantify the performance of four selected biomasses and found PSO was better than the others. Xu [14] investigated pyrolysis of rape straw and proposed a multi-component parallel reaction scheme to determine the kinetics using PSO. Saman [15] compared the performance of a decision tree [16] and PSO in determining the optimal solar power source location. Mortaza [17] and Zhao [18] found the combination of an adaptive network fuzzy inference system (ANFIS) [19], artificial neural network (ANN) [20], PSO and optimization of their initial weights and thresholds could significantly improve the prediction accuracy. Other scholars [21][22][23][24] also conducted similar studies using PSO in combination with other optimization algorithms. Among these algorithms, PSO is most commonly used, and it has the advantages of fewer parameters, high accuracy, fast convergence and the ability to search collaboratively [25]. Most of existing studies using optimization methods and inverse modeling to determine thermodynamics of solids employed constant heat flux to achieve simplicity, such as a cone calorimeter [26] and CAPA II [27]. However, this experimental strategy is oversimplified since, in practical scenarios, the heat flux received by the solid changes with time due to fire growth and flame spread. Whether the thermodynamic parameters determined under constant heat fluxes are applicable to time-dependent heat flux conditions has not been revealed. Furthermore, relevant studies aimed at deriving thermodynamics of solids by inversely modeling experimental measurements under time-dependent heat fluxes are rarely reported.
To solve this issue, this work utilized PSO and a previously developed 1D numerical model [28] to determine temperature-dependent k and C p of beech wood and PMMA based on experimentally measured surface and in-depth temperatures. Both constant and time-dependent radiative heat fluxes were applied in tests to achieve more comprehensive heating conditions. The obtained thermodynamic parameters were validated by blind predictions of the experimental results which were not used in model parameterization and the reported values in the literature. The thermodynamics determined in this study are more applicable to predicting thermal response and fire behaviors of solids under time-dependent heating conditions, which are more frequently encountered in practical applications than the oversimplified constant heat flux scenarios.

Experiment Description
Experimentally measured surface and in-depth temperatures of beech wood and PMMA under constant and time-dependent radiative heat fluxes (HFs) were used for PSO combined with a 1D numerical model. All the heating and ignition tests were elaborated in our previous studies, and therefore only the main procedures are briefly introduced in this section. More detailed information about the experiments can be found in our recent publications. Ignition tests of wood were conducted under constant [29] and power-law time-dependent HFs [30], while the experiments involving PMMA were carried out under linear HFs [31].

Wood Ignition Tests under Constant Heat Fluxes
Beech wood circular samples with 15 mm thickness and 60 mm diameter were prepared with grain parallel to the radiation. The measured average density of dried samples was 654 ± 20 kg m −3 . A hole with 0.5 mm diameter and 12 mm depth was drilled on the rear surface of the sample for 3 mm in-depth temperature measurements. Two 0.5 mm diameter type K thermocouples were used to measure the surface and in-depth temperatures. Detailed information of sample dimensions and the locations of thermocouples are demonstrated in Figure 1. The ignition experiments were performed using a heating system consisting of a control box and a cone heater to achieve four constant HFs, 25,30,40 and 50 kW/m 2 . The sample was embedded in a 25 mm thick and 90 cm square ceramic fiber plate that served as insulation material. When the prescribed HF was calibrated and stabilized before each test, the shutter of the heater was closed, and the equipped sample was loaded 3 mm below the heater. Subsequently, the shutter opened, and the test commenced. The test was terminated when a sustained visible flame was observed. The relevant input parameters of beech wood used in the numerical model during optimization and the optimized k and C p , which are given in the following sections, are listed in Table 1.

Wood Ignition Tests under Constant Heat Fluxes
Beech wood circular samples with 15 mm thickness and 60 mm diameter were prepared with grain parallel to the radiation. The measured average density of dried samples was 654 ± 20 kg m −3 . A hole with 0.5 mm diameter and 12 mm depth was drilled on the rear surface of the sample for 3 mm in-depth temperature measurements. Two 0.5 mm diameter type K thermocouples were used to measure the surface and in-depth temperatures. Detailed information of sample dimensions and the locations of thermocouples are demonstrated in Figure 1. The ignition experiments were performed using a heating system consisting of a control box and a cone heater to achieve four constant HFs, 25,30,40 and 50 kW/m 2 . The sample was embedded in a 25 mm thick and 90 cm square ceramic fiber plate that served as insulation material. When the prescribed HF was calibrated and stabilized before each test, the shutter of the heater was closed, and the equipped sample was loaded 3 mm below the heater. Subsequently, the shutter opened, and the test commenced. The test was terminated when a sustained visible flame was observed. The relevant input parameters of beech wood used in the numerical model during optimization and the optimized k and Cp, which are given in the following sections, are listed in Table 1.

Wood Ignition Tests under Power-Law Heat Fluxes
The experimental apparatus and procedures of wood ignition tests under power-law HFs are similar to those under constant HFs. The difference is that a proportion-integration-differentiation (PID) temperature controller, Delta DT320, was used to adjust the output power of the heater to achieve the power-law radiation. Meanwhile, an additional 0.5 mm diameter type K thermocouple was embedded inside the interior of sample to record the 6 mm in-depth temperature. In other words, surface, 3 and 6 mm in-depth temperatures were simultaneously collected in each test. Detailed information of sample dimensions and the locations of thermocouples are demonstrated in Figure 2. Originally, the impact of moisture content (based on dry weight) on ignition behaviors was examined using four moisture contents, 0%, 15%, 25% and 38%. However, only the measured temperatures of oven-dried samples were used for optimization to ensure consistency with other studies in the literature. Meanwhile, dry wood minimizes the effect of cracks during heating, which is intensified for moist wood. The calibrated power-law HFs can be expressed as:

Wood Ignition Tests under Power-Law Heat Fluxes
The experimental apparatus and procedures of wood ignition tests under power-law HFs are similar to those under constant HFs. The difference is that a proportion-integrationdifferentiation (PID) temperature controller, Delta DT320, was used to adjust the output power of the heater to achieve the power-law radiation. Meanwhile, an additional 0.5 mm diameter type K thermocouple was embedded inside the interior of sample to record the 6 mm in-depth temperature. In other words, surface, 3 and 6 mm in-depth temperatures were simultaneously collected in each test. Detailed information of sample dimensions and the locations of thermocouples are demonstrated in Figure 2. Originally, the impact of moisture content (based on dry weight) on ignition behaviors was examined using four moisture contents, 0%, 15%, 25% and 38%. However, only the measured temperatures of oven-dried samples were used for optimization to ensure consistency with other studies in the literature. Meanwhile, dry wood minimizes the effect of cracks during heating, which is intensified for moist wood. The calibrated power-law HFs can be expressed as:  Table 2.

PMMA Ignition Tests under Linear Heat Fluxes
The PMMA ignition experiments were similar to those of wood under time-dependent HFs, and the difference was that finite thick samples, 1, 1.5, 3, 6 and 10 mm thick, and six linear HFs were used. Considering the fact that 1-6 mm thick PMMA samples suffered deformation during heating, especially for 1-3 mm ones, which significantly undermines the 1D premise of the numerical model, only the surface and in-depth temperatures of 10 mm thick samples were used for optimization. Detailed information of sample dimensions and the locations of thermocouples are demonstrated in Figure 3. The calibrated linear HFs can be expressed as: where ' ' 0 q  is an arbitrary initial HF, taken as 6.5 kW/m 2 in our tests, and a is a constant determining the increasing rate of HF. The values of a of the six linear HFs are listed in Table 3. The relevant parameters of PMMA used in the numerical model and the optimized k and Cp, which are given in the following sections, are listed in Table 4.

PMMA Ignition Tests under Linear Heat Fluxes
The PMMA ignition experiments were similar to those of wood under time-dependent HFs, and the difference was that finite thick samples, 1, 1.5, 3, 6 and 10 mm thick, and six linear HFs were used. Considering the fact that 1-6 mm thick PMMA samples suffered deformation during heating, especially for 1-3 mm ones, which significantly undermines the 1D premise of the numerical model, only the surface and in-depth temperatures of 10 mm thick samples were used for optimization. Detailed information of sample dimensions and the locations of thermocouples are demonstrated in Figure 3. The calibrated linear HFs can be expressed as: where . q 0 is an arbitrary initial HF, taken as 6.5 kW/m 2 in our tests, and a is a constant determining the increasing rate of HF. The values of a of the six linear HFs are listed in Table 3. The relevant parameters of PMMA used in the numerical model and the optimized k and C p , which are given in the following sections, are listed in Table 4.

PMMA Ignition Tests under Linear Heat Fluxes
The PMMA ignition experiments were similar to those of wood under time-dependent HFs, and the difference was that finite thick samples, 1, 1.5, 3, 6 and 10 mm thick, and six linear HFs were used. Considering the fact that 1-6 mm thick PMMA samples suffered deformation during heating, especially for 1-3 mm ones, which significantly undermines the 1D premise of the numerical model, only the surface and in-depth temperatures of 10 mm thick samples were used for optimization. Detailed information of sample dimensions and the locations of thermocouples are demonstrated in Figure 3. The calibrated linear HFs can be expressed as: q  is an arbitrary initial HF, taken as 6.5 kW/m 2 in our tests, and a is a constant determining the increasing rate of HF. The values of a of the six linear HFs are listed in Table 3. The relevant parameters of PMMA used in the numerical model and the optimized k and Cp, which are given in the following sections, are listed in Table 4.     Table 4. Key parameters of PMMA and optimized k and C p .

Numerical Model
A previously established 1D numerical model [28] involving solid pyrolysis and heat transfer in a multilayer composite was employed to simulate the measured surface and in-depth temperatures. Utilization of a 1D numerical model was determined by the experimental condition simulated. In all ignition tests, the four sides and bottom surface of samples were wrapped with thermal insulation material to minimize the heat loss on these surfaces, leaving only the top surface exposed to radiation. Such sample assembly follows the protocols in standard ignition/flammability tests, such as the cone calorimeter [26] and fire propagation apparatus (FPA) [32], to achieve a 1D heating condition. Consequently, using a 1D numerical model is appropriate. In the original model, a multilayer structure composed of different combustibles can be considered. However, in this study, only a thick layer, at least 10 mm thickness, was involved. Consequently, only the governing equation in this layer was introduced. Only surface absorption of radiation was considered, and in-depth absorption was ignored even for PMMA. Meanwhile, it was hypothesized that the evolved volatiles by pyrolysis are released instantaneously upon generation, and therefore no mass transfer is involved. Only the condensed phase is simulated, and the external heat flux from the heater is treated as the boundary condition. The energy and mass conservation equations are expressed as: where ρ is the density, T is the temperature of solid, x is the spatial variable, r is the reflectivity of the solid surface, S v is the pyrolysis rate, ∆H v is the reaction heat, T ∞ is the environmental temperature, C g is the specific heat of gas and . m is the mass flux in the controlled volume. The pyrolysis reaction rate can be written as: where A is the pre-exponential factor, E a is the activation energy and R is the ideal gas constant. The transient mass flux in a controlled volume at fixed location x and the total transient mass flux can be calculated as: where l is the thickness of sample. The initial and boundary conditions are: where ε, σ and h c are the absorptivity/emissivity of the top surface, Stefan-Boltzmann constant and convection coefficient of the top surface, respectively. The grid size and time step used were 0.05 mm and 0.5 s, respectively, and the duration time of simulation was determined by the relevant experimental results. More detailed information of the numerical model as well as the computation methodologies can be found in Ref. [28]. Since k and C p are the target parameters in the following optimization process and pyrolysis of the solid exerts limited effect on surface and in-depth temperatures prior to ignition, thermal decomposition in the condensed phase is neglected to accelerate the convergency efficiency in optimization.

PSO Algorithm
The particle swarm optimization (PSO) algorithm was first proposed based on the process of finding the best feeding area for a flock of birds [33]. As an intelligent algorithm, PSO simulates the process of optimal decision making. It is assumed that at the beginning of foraging, the movement of the flock is chaotic, and over time the birds in random positions learn from each other and share foraging information within the flock. Each bird combines its own experience and the information transmitted by its peers to estimate the value of the current position to find food. Compared with other optimization algorithms, the PSO algorithm has higher convergence efficiency and accuracy [11]. In many optimization problems, the PSO algorithm converges to the optimal solution faster than GA. However, it also may fall into the local optimal solution, resulting in low convergence accuracy. Considering this drawback of PSO, the search range selection is carried out manually in current study to improve its convergence efficiency. PSO is applied to estimate temperature-dependent k and C p of wood and PMMA combining the numerical model and experimentally measured temperatures at constant and time-dependent HFs.
In PSO, the input parameters of the numerical model excluding k and C p were obtained from the literature, our previous publications, and experimental measurements. Then, several sets of temperature-dependent k and C p of wood and PMMA were collected from the literature, and their average values of k and C p were calculated. For each parameter, the lower bound of the search range was set to 0.1 times the average value, while the upper bound was determined to ensure the average value is also the mean value of the search range. Four constants need to be estimated in k and C p expressions: where k 0 , k T , C p,0 and C p,T are constants. The objective function necessary to quantify the discrepancy between experimental and numerical results in each iteration is defined as: where F obj is the objective function; the subscript h, exp and num denote the depth of the monitor point in tests (0, 3 and 6 mm), experimental and numerical results, respectively; and n exp is the total number of the experiment data points. The lower value of F obj indicates better fitness and vice versa. The particle swarm size and the iteration number were selected to be 80 and 50, respectively. Combining the numerical model, heuristic optimization algorithms and experimental measurements to determine unknown properties is a common routine. More specifically, it is commonly called the curve-fitting method [34]. In this method, the experimental curves are simulated by the numerical model, and the target unknown parameters serving as inputs are adjusted by the PSO algorithm in each iteration until a prescribed error tolerance between numerical and experimental curves is achieved, indicating convergence. All computation was implemented using in-house MATLAB (2018b) code.

Optimization Results of Wood
In order to obtain more reliable k and C p of wood, PSO is implemented at each constant HF, from 25 to 50 kW m −2 , using experimentally measured surface and 3 mm indepth temperatures, and the optimization results are compared with experimental data in Figure 4. In each experiment, at least three repeatable tests were conducted to guarantee the reproducibility and to estimate the experimental uncertainty. All the calculated uncertainty ranges reflected by the lengths of error bars can be found in Refs. [29][30][31]. Nevertheless, these error bars are not depicted in Figure 4 and the following figures to avoid potential location congestion. In all circumstances, the numerical profiles fit the experimental ones well despite some minor divergence since enforcing minimal discrepancy is essential for any optimization algorithm. In all subgraphs, the deviation of surface temperature is larger compared to 3 mm in-depth temperature, which is presumably due to the small cracks emerging over surface during heating and the generated thin char layer before ignition. These influential factors are not involved in the numerical model. Good fitness of 3 mm in-depth temperature is found in all subgraphs, which is particularly important for accurate prediction of k and C p since these two parameters are highly relevant to the heat transfer inside the condensed phase.
Energies 2023, 16, x FOR PEER REVIEW

Optimization Results of Wood
In order to obtain more reliable k and Cp of wood, PSO is implemented at ea stant HF, from 25 to 50 kW m −2 , using experimentally measured surface and 3 depth temperatures, and the optimization results are compared with experimen in Figure 4. In each experiment, at least three repeatable tests were conducted to gu the reproducibility and to estimate the experimental uncertainty. All the calculated tainty ranges reflected by the lengths of error bars can be found in Refs. [29][30][31]. theless, these error bars are not depicted in Figure 4 and the following figures t potential location congestion. In all circumstances, the numerical profiles fit the mental ones well despite some minor divergence since enforcing minimal discrep essential for any optimization algorithm. In all subgraphs, the deviation of surfa perature is larger compared to 3 mm in-depth temperature, which is presumably the small cracks emerging over surface during heating and the generated thin ch before ignition. These influential factors are not involved in the numerical mode fitness of 3 mm in-depth temperature is found in all subgraphs, which is particul portant for accurate prediction of k and Cp since these two parameters are highly r to the heat transfer inside the condensed phase.  Figure 5a,b, respectively, and their linear temperature dependen  The optimized k and C p of wood under each HF in Figure 4 are portrayed in Figure 5. Apparently, for each parameter, the optimized results at different HFs are very close to each, suggesting high accuracy of the optimization strategy. The average values of k and C p are given in Figure 5a,b, respectively, and their linear temperature dependencies are listed in Table 1  A similar optimization strategy was also implemented for wood under the f power-law HFs listed in Table 2, and the comparison between experimental and numer results are demonstrated in Figure 6. The difference is that an additional set of 6 mm depth temperatures were involved. Again, PSO was performed at each HF, and there four temperature-dependent expressions were obtained for k and Cp. Under all power-HFs, the numerical curves fit the experimental ones well except for the end segment surface temperature before ignition. Under power-law HFs, the heating times are m longer than those under constant HFs in Figure 5. As a result, the prolonged heating t facilitates the accumulation of yielded char layer over surface due to wood pyrolysis. B wood pyrolysis and char oxidation were not considered in the numerical model, therefore the measured surface temperatures were underpredicted. However, in the jority portion of surface temperature curve where the impact of pyrolysis can be igno the numerical model matches the experimental measurements well. A similar optimization strategy was also implemented for wood under the four powerlaw HFs listed in Table 2, and the comparison between experimental and numerical results are demonstrated in Figure 6. The difference is that an additional set of 6 mm in-depth temperatures were involved. Again, PSO was performed at each HF, and therefore four temperature-dependent expressions were obtained for k and C p . Under all power-law HFs, the numerical curves fit the experimental ones well except for the end segments of surface temperature before ignition. Under power-law HFs, the heating times are much longer than those under constant HFs in Figure 5. As a result, the prolonged heating time facilitates the accumulation of yielded char layer over surface due to wood pyrolysis. Both wood pyrolysis and char oxidation were not considered in the numerical model, and therefore the measured surface temperatures were underpredicted. However, in the majority portion of surface temperature curve where the impact of pyrolysis can be ignored, the numerical model matches the experimental measurements well.
Again, the optimized k and C p under the four power-law HFs and their average values are illustrated in Figure 7. A similar variation trend and moderate scatter are observed in each group. The attained average k and C p are listed in Table 1.
The optimized average k and C p of beech wood at constant and power-law HFs were compared with the reported ones in Refs. [29,35], as demonstrated in Figure 8. In Gong's [29] work, these two parameters were optimized using GA under different constant HFs, while in Ross's [35] study they were directly measured by experimental methods, and therefore constant k was obtained, and the linear temperature dependency of C p was only applicable in a narrow temperature range. In Figure 8a, both the optimized k under constant and power-law HFs in current study agree relatively well with the one in Ref. [29] but are slightly higher than the constant value in Ref. [35]. In Figure 8b, our optimized C p at constant HF agrees well with the one in Gong's [29] work where the same experimental data but also the GA algorithm were utilized. This agreement again confirms the high reliability of the optimization strategy employed in this study. Meanwhile, the optimized C p of wood at power-law HF fits Ross's [35] result well. The relatively large divergence between the two sets of C p is primarily attributed to the contained moisture content in wood samples of constant HF tests [29]. In normal environment conditions, the moisture content of wood typically varies in the range of 8-10% [35]. However, in power-law HF Energies 2023, 16, 5302 9 of 16 tests [30] and Ref. [35], the wood samples were oven-dried. A high value of C p of water caused this discrepancy.
four temperature-dependent expressions were obtained for k and Cp. Under all power-law HFs, the numerical curves fit the experimental ones well except for the end segments of surface temperature before ignition. Under power-law HFs, the heating times are much longer than those under constant HFs in Figure 5. As a result, the prolonged heating time facilitates the accumulation of yielded char layer over surface due to wood pyrolysis. Both wood pyrolysis and char oxidation were not considered in the numerical model, and therefore the measured surface temperatures were underpredicted. However, in the majority portion of surface temperature curve where the impact of pyrolysis can be ignored, the numerical model matches the experimental measurements well.   Table 1. The optimized average k and Cp of beech wood at constant and power-law HFs were compared with the reported ones in Refs. [29,35], as demonstrated in Figure 8. In Gong's [29] work, these two parameters were optimized using GA under different constant HFs, while in Ross's [35] study they were directly measured by experimental methods, and therefore constant k was obtained, and the linear temperature dependency of Cp was only applicable in a narrow temperature range. In Figure 8a, both the optimized k under constant and power-law HFs in current study agree relatively well with the one in Ref. [29] but are slightly higher than the constant value in Ref. [35]. In Figure 8b, our optimized Cp

Optimization Results of PMMA
Similar optimization procedures were also conducted for PMMA under the six linear HFs in Table 3, and the comparison between experimental and numerical results is depicted in Figure 9. Under all linear HFs, good agreement is observed despite some minor divergence. The small disagreements of in-depth temperatures in Figure 9a-c are presumably caused by the potential inaccurate locations of the thermocouple beads during tests, while the slight overprediction of surface temperature at the end segments is associated with the neglect of endothermic pyrolysis reaction of PMMA. Furthermore, the average values of k and Cp, listed in Table 4, based on the six optimization runs were plotted and compared with the reported ones in Stoliarov's [36] work, as shown in Figure 10. Both our optimized k and Cp show similar variation trends to those in Ref. [36] though piecewise functions are used in that work. The discrepancy between our and Stoliarov's results may be attributed to the different types and processing of PMMA used.

Optimization Results of PMMA
Similar optimization procedures were also conducted for PMMA under the six linear HFs in Table 3, and the comparison between experimental and numerical results is depicted in Figure 9. Under all linear HFs, good agreement is observed despite some minor divergence. The small disagreements of in-depth temperatures in Figure 9a-c are presumably caused by the potential inaccurate locations of the thermocouple beads during tests, while the slight overprediction of surface temperature at the end segments is associated with the neglect of endothermic pyrolysis reaction of PMMA. Furthermore, the average values of k and C p , listed in Table 4, based on the six optimization runs were plotted and compared with the reported ones in Stoliarov's [36] work, as shown in Figure 10. Both our optimized k and C p show similar variation trends to those in Ref. [36] though piecewise functions are used in that work. The discrepancy between our and Stoliarov's results may be attributed to the different types and processing of PMMA used.

Optimization Results of PMMA
Similar optimization procedures were also conducted for PMMA under the six linear HFs in Table 3, and the comparison between experimental and numerical results is depicted in Figure 9. Under all linear HFs, good agreement is observed despite some minor divergence. The small disagreements of in-depth temperatures in Figure 9a-c are presumably caused by the potential inaccurate locations of the thermocouple beads during tests, while the slight overprediction of surface temperature at the end segments is associated with the neglect of endothermic pyrolysis reaction of PMMA. Furthermore, the average values of k and Cp, listed in Table 4, based on the six optimization runs were plotted and compared with the reported ones in Stoliarov's [36] work, as shown in Figure 10. Both our optimized k and Cp show similar variation trends to those in Ref. [36] though piecewise functions are used in that work. The discrepancy between our and Stoliarov's results may be attributed to the different types and processing of PMMA used.

Parameter Validation
The reliability of the optimized k and Cp was verified by simulating the experimental temperatures that were not utilized in parameterization. Figure 11 shows experimental surface, 3 and 6 mm in-depth temperatures of wood under a power-law HF, HF1 in Table  2, and the corresponding simulation results using the optimized k and Cp at constant HF in Table 1. Apparently, the numerical model successfully captures the measured surface temperature except for the end segment. Nevertheless, moderate and unacceptable discrepancies are observed for 3 and 6 mm in-depth temperatures. These deviations are mainly attributable to the existence of moisture content in wood samples in constant HF tests as discussed in the preceding section. Higher k and Cp at constant HF lead to lower predicted in-depth temperatures as depicted in Figure 11. However, the moisture content at the surface is eliminated quickly after exposure of radiation and thus has little impact on surface temperature history.

Parameter Validation
The reliability of the optimized k and C p was verified by simulating the experimental temperatures that were not utilized in parameterization. Figure 11 shows experimental surface, 3 and 6 mm in-depth temperatures of wood under a power-law HF, HF1 in Table 2, and the corresponding simulation results using the optimized k and C p at constant HF in Table 1. Apparently, the numerical model successfully captures the measured surface temperature except for the end segment. Nevertheless, moderate and unacceptable discrepancies are observed for 3 and 6 mm in-depth temperatures. These deviations are mainly attributable to the existence of moisture content in wood samples in constant HF tests as discussed in the preceding section. Higher k and C p at constant HF lead to lower predicted in-depth temperatures as depicted in Figure 11. However, the moisture content at the surface is eliminated quickly after exposure of radiation and thus has little impact on surface temperature history.

Parameter Validation
The reliability of the optimized k and Cp was verified by simulating the experimental temperatures that were not utilized in parameterization. Figure 11 shows experimental surface, 3 and 6 mm in-depth temperatures of wood under a power-law HF, HF1 in Table  2, and the corresponding simulation results using the optimized k and Cp at constant HF in Table 1. Apparently, the numerical model successfully captures the measured surface temperature except for the end segment. Nevertheless, moderate and unacceptable discrepancies are observed for 3 and 6 mm in-depth temperatures. These deviations are mainly attributable to the existence of moisture content in wood samples in constant HF tests as discussed in the preceding section. Higher k and Cp at constant HF lead to lower predicted in-depth temperatures as depicted in Figure 11. However, the moisture content at the surface is eliminated quickly after exposure of radiation and thus has little impact on surface temperature history.  under each HF of 30, 40 and 50 kW m −2 can be validated using the same method. Meanwhile, the optimized k and Cp of wood under each power-law HF can also be verified by simulating the experimental measurements of the remaining power-law HFs. Good agreement is detected in all these validation cases, and they are not exhibited for concision.

Impact of Search Range on PSO
In the tests described in the preceding sections, the lower bound of the search range during PSO was set to 0.1 times the manually selected average value, while the upper bound was determined to ensure the selected average value was also the mean value of A similar validation methodology was further implemented for PMMA under linear HF, as demonstrated in Figure 13 where the optimized k and C p under linear HF of HF 4 are verified under the other five linear HFs. In all validation cases, the numerical model well captures the variation trends and magnitudes of all surface, 3 and 6 mm in-depth temperatures, which again implies high credibility and accuracy of the optimization method and numerical model. under each HF of 30, 40 and 50 kW m −2 can be validated using the same method. Meanwhile, the optimized k and Cp of wood under each power-law HF can also be verified by simulating the experimental measurements of the remaining power-law HFs. Good agreement is detected in all these validation cases, and they are not exhibited for concision. A similar validation methodology was further implemented for PMMA under linear HF, as demonstrated in Figure 13 where the optimized k and Cp under linear HF of HF 4 are verified under the other five linear HFs. In all validation cases, the numerical model well captures the variation trends and magnitudes of all surface, 3 and 6 mm in-depth temperatures, which again implies high credibility and accuracy of the optimization method and numerical model.

Impact of Search Range on PSO
In the tests described in the preceding sections, the lower bound of the search range during PSO was set to 0.1 times the manually selected average value, while the upper bound was determined to ensure the selected average value was also the mean value of

Impact of Search Range on PSO
In the tests described in the preceding sections, the lower bound of the search range during PSO was set to 0.1 times the manually selected average value, while the upper bound was determined to ensure the selected average value was also the mean value of the search range. In order to examine the impact of search range on optimization results, an additional expanded large search range was employed to implement the optimization. The lower bound of the large search range is 0.01 times the selected average value, while the upper bound was obtained through the same strategy as mentioned above. Wood experimental results under 25 kW m −2 constant HF were used to conduct the optimization, and the comparison is illustrated in Figure 14, where LSR and SSR denote large and small search ranges, respectively. When using SSR, good agreement is found for both surface and 3 mm in-depth temperatures, implying convergence and high accuracy. However, if LSR is used, appreciable divergence emerges even though the objective function, F obj in Equation (11), remains unchanged with increasing iterations during optimization, suggesting convergence. This unacceptable deviation is presumably caused by the inherent flaw of the PSO algorithm in that it may be trapped in a local optimal solution [12]. Therefore, it can be concluded that utilization of a reasonable narrow search range can effectively suppress the potential occurrence of local optimal solution. Meanwhile, a reasonable narrow search range can also effectively diminish the potential compensation effect, which is still an unsolved issue in determining unknown parameters using the optimization algorithm as interpreted by Vyazovkin [37]. the search range. In order to examine the impact of search range on optimization results, an additional expanded large search range was employed to implement the optimization. The lower bound of the large search range is 0.01 times the selected average value, while the upper bound was obtained through the same strategy as mentioned above. Wood experimental results under 25 kW m −2 constant HF were used to conduct the optimization, and the comparison is illustrated in Figure 14, where LSR and SSR denote large and small search ranges, respectively. When using SSR, good agreement is found for both surface and 3 mm in-depth temperatures, implying convergence and high accuracy. However, if LSR is used, appreciable divergence emerges even though the objective function, Fobj in Equation (11), remains unchanged with increasing iterations during optimization, suggesting convergence. This unacceptable deviation is presumably caused by the inherent flaw of the PSO algorithm in that it may be trapped in a local optimal solution [12]. Therefore, it can be concluded that utilization of a reasonable narrow search range can effectively suppress the potential occurrence of local optimal solution. Meanwhile, a reasonable narrow search range can also effectively diminish the potential compensation effect, which is still an unsolved issue in determining unknown parameters using the optimization algorithm as interpreted by Vyazovkin [37].
Particle swarm size (PSS) and iteration number (IN) in PSO may also pose a nonnegligible impact on optimization results. To investigate their effects, an additional optimization run where both PSS and IN were doubled was implemented using the same experimental data and the same SSR, and the optimization results are depicted in Figure 14. Apparently, little discrepancy exists between the original case and the modified case, implying the original PSS and IN are sufficiently high for accurate optimization and convergence.

Conclusions
Two important thermodynamic parameters, k and Cp, of two representative solid combustibles, beech wood and PMMA, were quantitatively determined by the PSO algorithm, a 1D numerical model and the related measured surface and in-depth temperatures in ignition tests under constant, power-law and linear HFs. Linearized temperature-dependent k and Cp were obtained through the employed approach. In each set of HFs, the optimization was carried out under each HF, and the average value was deemed the final result. For wood at constant HFs, good agreement was found between experimental and optimization results, and the scatter among the optimized thermodynamics is relatively Particle swarm size (PSS) and iteration number (IN) in PSO may also pose a nonnegligible impact on optimization results. To investigate their effects, an additional optimization run where both PSS and IN were doubled was implemented using the same experimental data and the same SSR, and the optimization results are depicted in Figure 14. Apparently, little discrepancy exists between the original case and the modified case, implying the original PSS and IN are sufficiently high for accurate optimization and convergence.

Conclusions
Two important thermodynamic parameters, k and C p , of two representative solid combustibles, beech wood and PMMA, were quantitatively determined by the PSO algorithm, a 1D numerical model and the related measured surface and in-depth temperatures in ignition tests under constant, power-law and linear HFs. Linearized temperature-dependent k and C p were obtained through the employed approach. In each set of HFs, the optimization was carried out under each HF, and the average value was deemed the final result. For wood at constant HFs, good agreement was found between experimental and optimization results, and the scatter among the optimized thermodynamics is relatively small, implying high accuracy of the utilized method, while under power-law HFs, good fitness was also observed for surface and in-depth temperatures of wood. Larger scatter of the optimized thermodynamics was detected, and an alternative pair of k and C p was obtained. All these thermodynamic parameters fell into the same ranges of the reported ones in the literature. The deviation of these two sets of k and C p was attributed to the moisture content difference of the utilized samples. Subsequently, the same optimization strategy was further implemented for PMMA under six linear HFs. Good agreement was again found for surface and in-depth temperatures under all HFs, and the attained k and C p matched the correlations reported in the literature well. Credibility of the optimized wood thermodynamics under constant HFs was verified by emulating the measured temperatures of wood under a power-law HF, and good agreement existed for surface temperature. The divergences for 3 and 6 mm in-depth temperatures were caused by moisture content difference. Additionally, the optimized thermodynamics of wood at 25 kW m −2 were validated by predicting the experimental temperatures under the remaining constant HFs. Meanwhile, the optimized thermodynamics of PMMA under HF4 were verified by simulating the measured temperatures under the remaining linear HFs. In all these scenarios, the numerical curves matched the experimental data well, suggesting high reliability of the obtained k and C p . Furthermore, the impact of search range of PSO on optimization results was investigated by using an additional expanded search range. The outcomes showed that larger search range featured a higher possibility of being trapped in local optimal solution, and therefore a narrow search range is strongly suggested. Using doubled particle swarm size and iteration number had little effect on optimization results, indicating the originally adopted values were sufficiently high for convergence.
A potential weakness of current study is the fact that pyrolysis in the solid before ignition was ignored. This assumption is reasonable and would introduce negligible influence for piloted ignition wherein pyrolysis exerts little impact [34]. Nevertheless, this hypothesis is unacceptable for autoignition, where intense thermal decomposition occurs in the vicinity of surface before ignition, implying pyrolysis cannot be neglected. Accurate predictions of surface temperature and mass loss rate necessitate full knowledge of pyrolysis kinetics. In other words, kinetic parameters, such as preexponential factor and activation energy in the Arrhenius equation can also be quantitatively estimated by inversely modeling the measured temperatures and mass loss rate during autoignition or combustion using proper optimization algorithms. All these deserve further in-depth exploration in future work.