Kinetic Study on the Pyrolysis of Medium Density Fiberboard : Effects of Secondary Charring Reactions

The reaction models employed in the kinetic studies of biomass pyrolysis generally do not include the secondary charring reactions. The aim of this work is to propose an applicable kinetic model to characterize the pyrolysis mechanism of medium density fiberboard (MDF) and to evaluate the effects of secondary charring reactions on estimated products yields. The kinetic study for pyrolysis of MDF was performed by a thermogravimetric analyzer over a heating rate range from 10 to 40 ◦C/min in a nitrogen atmosphere. Four stages related to the degradation of resin, hemicellulose, cellulose, and lignin could be distinguished from the thermogravimetric analyses (TGA). Based on the four components and multi-component parallel reaction scheme, a kinetic model considering secondary charring reactions was proposed. A comparison model was also provided. An efficient optimization algorithm, differential evolution (DE), was coupled with the two models to determine the kinetic parameters. Comparisons of the results of the two models to experiment showed that the mass fraction (TG) and mass loss rate (DTG) calculated by the model considering secondary charring reactions were in better agreement with the experimental data. Furthermore, higher product yields than the experimental values will be obtained if secondary charring reactions were not considered in the kinetic study of MDF pyrolysis. On the contrary, with the consideration of secondary charring reactions, the estimated product yield had little error with the experimental data.


Introduction
Medium density fiberboard (MDF) is a kind of wood composite material which is widely used in the construction and housing furniture.According to Food and Agriculture Organization (FAO) statistics [1], the total MDF production was more than 100 million m 3 over the world in 2017.Specially, in China, the number was more than 59 million m 3 .A large amount of waste MDF was generated from the need for regular replacement of furniture due to aging, mechanical damage, fire burning, etc. [2].At the beginning, most of the waste MDF was sent to burn for cooking or space heating, while a small portion was disposed in landfills [3].With the rapid rising of MDF waste, land occupation and environmental problems are becoming more serious.From the perspective of protecting the natural environment and the sustainable use of resources, it is a good choice to recycle this large amount of waste as energy fuel.Pyrolysis has been proved to be a promising thermochemical conversion technology for resource and energy recovery [4].On one hand, biomass can be converted into low molecular weight products during the pyrolysis process with low environmental damage.On the other hand, the products generated can be recycled as useful chemicals or gaseous fuels.Hence, it is important to investigate the pyrolysis behavior and pyrolysis kinetics of MDF.
Thermogravimetric analysis (TGA) is widely used to study the pyrolysis behavior of biomass as it is a high-precision method to obtain necessary information to reveal the kinetic mechanism [5][6][7].
Energies 2018, 11, 2481 2 of 17 Hashim et al. employed TGA to analyze thermal properties of MDF [8].Additionally, some quantitative methods are usually coupled with TGA curves to determine kinetic parameters, such as model-free methods and model-fitting methods.For example, Yorulmaz et al. calculated the thermal kinetic constants for waste wood samples (including MDF) by using the Coats-Redfern model-fitting method [9].Han et al. estimated the activation energy of waste MDF by the Ozawa model-free methods [10].Although these studies may indicate the existence of independent reaction of biomass component, the underlying reaction mechanism is not involved.
The single biomass component competitive scheme, which was first proposed by Shafizadeh and Chin [11], is commonly used for wood pyrolysis modelling in previous studies.For example, Bastiaans et al. assumed a single reaction to extract Arrhenius parameters and to estimate CO release from MDF [12].Additionally, competitive schemes for biomass multi-components (cellulose, hemicellulose, and lignin) have also been proposed by Koufopanos et al. [13] and Miller and Bellan [14].As depicted in Figure 1, for each biomass component, formation of char, tar, and gas complete with each other in the primary pyrolysis.The process of tar to further react to produce more gases and chars is called secondary charring [15].It is worth noting that this multi-component scheme was never applied in kinetic studies, only the simplified mechanisms of gas and char formation (the black part in Figure 1) were commonly employed in previous studies [16][17][18][19].Among them, Li et al. determined the kinetic parameters of MDF using the genetic algorithm coupled with Kissinger's method [19].However, recently reviews indicate that activation energy for each component calculated using the simplified mechanism varies widely [20,21].For example, the differences between the reported activation energies of hemicellulose are huge, from about 75 to 225 kJ/mol [15].The uncertainty estimation of kinetic parameters, which will result in a difference between the estimated product yields and the practical values, may be caused by heating dynamics in the experiment or the absence of secondary charring reactions in the kinetic studies [22].Lately, some scholars modified the kinetic scheme of biomass pyrolysis to include secondary charring reactions [22][23][24].However, the kinetic studies on the secondary charring reactions in the multi-component competitive scheme are rare.
In kinetic studies of biomass pyrolysis, the model-free methods and model-fitting methods were usually used to determine kinetic parameters.But, there are several problems with the model-fitting methods, such as the inability to choose the reaction order and the reaction model.The model-free method can only estimate the activation energy at specific conversion rate range for an independent model [25].Therefore, when analyzing thermogravimetric data coupled with secondary charring reactions, optimization algorithms must be employed to obtain the unknown parameters which cannot be obtained by the above methods.Differential evolution (DE) is an effective method to solve optimization problems in various research fields but rarely used in biomass pyrolysis analysis.Therefore, due to the convergence speed and ease of operation, DE was chosen for data optimization in this study.
The aim of this work is to propose an applicable kinetic model considering secondary charring reactions to better describing the thermal decomposition process of MDF.The differential evolution (DE) algorithm coupled with the proposed model was used to obtain the reaction kinetics of MDF pyrolysis and the effects of secondary charring reactions on products yields were also evaluated.MDF used in this work came from a wood factory in Hefei City, Anhui Province, China.According to the manufacturer, the main component of the MDF is pine fiber.At first, the MDF was milled, and then the particles with a size less than 125 µm were obtained through a series of Tyler sieves.The samples were dehydrated for at least 48 hours in an XGQ-2000 (BAIHUI, Dongguan, Guangdong, China) electric drying oven at 105 • C, and then sealed in plastics bags in preparation for following tests.According to ASTM standards, proximate analysis was performed using a LECO TGA 701thermogravimetric analyzer.Ultimate analysis was performed using a LECO TruSpec CHN to detect carbon, hydrogen, and nitrogen.Sulfur was detected using a LECO TruSpec S. The proximate and ultimate analysis results of MDF sample are shown in Table 1.

TGA Tests
Thermogravimetric analysis was carried out using a NETZSCH STA 449 F3 Jupiter (NETZSCH, Selb, Germany) thermal analyzer.The samples of 6 ± 0.5 mg were heated with a gas flow rate of 80 mL/min under a nitrogen atmosphere.The program temperature was increased from room temperature to about 800 • C at four selected heating rates of 10, 20, 30, and 40 • C /min.The real-time mass change and heat flux change were monitored using software in the whole process.TGA experiments were repeated three times and the results showed a good repeatability (mean uncertainly is less than 1%).

Kinetic Reaction Model
As illustrated in Section 1, lignocellulosic biomass can be considered to be composed of n components.In this work, the MDF sample is made by pine and less than 10% resin.Therefore, the independent reactions of four components: resin, hemicellulose, cellulose, and lignin, are applied to explain the decomposition of MDF.Table 2 lists the two models applied in this study detailed (Model II and Model I represent the models that considering the secondary charring reactions or not, respectively).For the convenience of writing, some parameters are listed here.f (α) is the conversion function that describes the reaction mechanism, g(α) is the integrated forms of f (α); m 0 , m t , and m ∞ represent the initial, instantaneous, and the final solid mass of samples, respectively; the kinetic triplets A is the pre-exponential factor, E is the activation energy and n is the reaction order; R is the gas constant; for component i, W i is the instantaneous mass fraction, W i,0 is the initial mass fraction, γ i , σ i , and ϕ i are the stoichiometry coefficient of tar and char, respectively.
It is known that the kinetic parameters of thermogravimetry can be effectively calculated based on the decomposition rate dα/dt which is assumed to be a function of temperature T and the degree of conversion α [26,27]: The conversion rate can be described as α = (m 0 − m t )/(m 0 − m ∞ ), and the kinetic constant of temperature k(T) is expressed by Arrhenius equation as k(T) = A exp(−E/RT).
The associated kinetic equations of models can be expressed as follows: • Model I The reaction rate for reactions in Model I can be calculated by: The total mass loss rate (dW/dt) tot and the total solid mass fraction W tot can be calculated as: • Model II Step 1 (primary pyrolysis): Similar to Model I, the mass loss rate (dW/dt) 1 and the solid mass fraction W 1 of step 1 can be calculated as: Step 2 (secondary charring reaction): The reaction rate for secondary charring reaction can be calculated by: Energies 2018, 11, 2481 5 of 17 where W tar,0 is the initial mass fraction of tar which is calculated from step 1: Therein, W tar , A tar , E tar , n tar denote the instantaneous mass fraction the kinetic triplets for tar generated by step 1, respectively.
Thus, the mass of char generated by step 2 W char,2 is calculated as: where θ is the coefficient of char in step 2. Therefore, the total solid mass fraction in the whole process can be expressed as: The total mass loss rate of whole process is then deduced as:

.2. Optimization Method-Differential Evolution Algorithm
Differential evolution (DE) is originally proposed by Rainer Storn and Kenneth Price, in 1995 [28].It is a kind of population-based adaptive global optimization algorithm.Assuming there are M individuals in a random population.An individual consists of n-dimensional vectors, which can be expressed as X i (g) = (x i,1 (g), x i,2 (g), . . .x i,n (g)), i = 1, 2, . . ., M; g = 1, 2, . . ., N. N is the maximum number of iterations.In the present study, one unknown parameter, such as the activation energy E, is an individual.The implementation of DE is as follows [29]:

• Initialization
The value of the jth dimension of the ith individual is calculated as: where x L ij , x U ij denote the lower and upper bands of the ith individual, respectively.rand (0, 1) denotes a random number distributed in [0, 1].

• Variation
A variation individual is calculated via three random individuals: where F is the scaling factor.

• Crossover
The new individual can be calculated as: where Cr is the cross rate distributed in [0, 1].

• Selection
In this step, the individual in the next iteration is selected by: where f (x) is the fitness of the individuals.In the present study, the fitness is evaluated by the objective function obj, which can be calculated as: where W tot and (dW/dt) tot are the calculated solid mass and mass loss rate in the pyrolysis process, respectively.exp means the corresponding experimental data.n j is the number of data points considered in each single heating rate and k is the number of heating rates.c is the weight coefficient, and c = 0.5 is set for equally weighted both mass and mass loss rate.
From the objective function, there are 28 unknown parameters in Model II: kinetic triplets (E i , A i , n i ), initial mass fractions (W i,0 ), tar yields (γ i ), and char yields (σ i , θ) for each component and tar.Since the initial mass fraction of the sum of the four components is equal to 1, the initial mass fraction of lignin can be calculated by W l,0 = 1 − W r,0 − W h,0 − W c,0 .Then the number of unknown variables is 27.So, there are 27 individuals in Model II and 19 in Model I.In order to obtain the optimal solution, appropriate algorithm parameters should be set.In this work, the population size, scaling factor, crossover probability and maximum iteration number of the optimization are set as 60, 0.5, 0.2, and 1000, respectively.The optimization process is implemented in MATLAB.

Estimation of Initial Guesses of Model Parameters for DE
If an appropriate initial guess of unknown variable is provided, the DE optimization efficiency can be greatly improved.Therefore, the kinetic triplets are firstly estimated prior to the optimization process.
The complementary use of model-free and model-fitting methods were suggested by Khawam and Flanagan [30] to determine solid state reaction kinetic triplets from experimental data.In this study, three model-free kinetic methods (Equations ( 21)-( 23)) [31][32][33] and one model-fitting method (Equation ( 24)) [34,35] were used to estimate the kinetic triplets.
The kinetic models used in this study were presented in Section 2.2.THE optimization methodology (DE) was introduced for the estimation of kinetic parameters.Moreover, to improve the DE optimization efficiency, three model-free kinetic methods and one model-fitting method were presented for the initial guesses of the kinetic triplets.

Thermogravimetric Analysis
The mass fraction (TG) and mass loss rate (DTG) profiles in the cases of 10, 20, 30, and 40 • C/min are illustrated in Figure 2. As shown in Figure 2a, a similar phenomenon occurs in any thermogravimetric analysis can be observed: a shift of the TG curves at higher temperature happens by the increase of the heating rate.In Figure 2b, there are less identifiable inflection points in the |DTG| curves.Li et al. [19] estimated the possible components and sub-reactions by comparing DTG and DDTG (the second derivatives of TG curves).Figure 3  It is noted that model-free methods can calculate the kinetic parameters without equation expression of g(α).
The kinetic models used in this study were presented in Section 2.2.THE optimization methodology (DE) was introduced for the estimation of kinetic parameters.Moreover, to improve the DE optimization efficiency, three model-free kinetic methods and one model-fitting method were presented for the initial guesses of the kinetic triplets.

Thermogravimetric Analysis
The mass fraction (TG) and mass loss rate (DTG) profiles in the cases of 10, 20, 30, and 40 °C/min are illustrated in Figure 2. As shown in Figure 2a, a similar phenomenon occurs in any thermogravimetric analysis can be observed: a shift of the TG curves at higher temperature happens by the increase of the heating rate.In Figure 2b, there are less identifiable inflection points in the |DTG| curves.Li et al. [19] estimated the possible components and sub-reactions by comparing DTG and DDTG (the second derivatives of TG curves).Figure 3    It is noted that model-free methods can calculate the kinetic parameters without equation expression of g(α).
The kinetic models used in this study were presented in Section 2.2.THE optimization methodology (DE) was introduced for the estimation of kinetic parameters.Moreover, to improve the DE optimization efficiency, three model-free kinetic methods and one model-fitting method were presented for the initial guesses of the kinetic triplets.

Thermogravimetric Analysis
The mass fraction (TG) and mass loss rate (DTG) profiles in the cases of 10, 20, 30, and 40 °C/min are illustrated in Figure 2. As shown in Figure 2a, a similar phenomenon occurs in any thermogravimetric analysis can be observed: a shift of the TG curves at higher temperature happens by the increase of the heating rate.In Figure 2b, there are less identifiable inflection points in the |DTG| curves.Li et al. [19] estimated the possible components and sub-reactions by comparing DTG and DDTG (the second derivatives of TG curves).Figure 3    Based upon Section 2.2.3., the Arrhenius plots for the conversion from 0.1 to 0.9 by model-free models FWO, Starink and KAS are shown in Figure 4a-c, respectively.Differences appear at α = 0.1, α = 0.2, and α = 0.8.The activation energies at different conversion obtained from the slopes of plots are presented in Figure 4d and Table 3. High correlation coefficients R 2 of fitting curves at different conversions are demonstrated.The values of E at each conversion determined by different model-free methods are almost the same.In addition, the values of E increase gradually at 0 ≤ α ≤ 0.1, followed by a sharp increase at 0.1 < α ≤ 0.2.Then E remains almost constant at 0.2 < α ≤ 0.8.Finally, a second sharp increase of E occurs at 0.8 < α ≤ 1.0.These ranges are coincidence with the position where the differences observed in Figure 4a-c.Therefore, four stages which related to four sub-reactions of four components can be divided in the pyrolysis process by α = 0.1, α = 0.2 and α = 0.8.Based upon Section 2.2.3., the Arrhenius plots for the conversion from 0.1 to 0.9 by model-free models FWO, Starink and KAS are shown in Figure 4a-c, respectively.Differences appear at α = 0.1, α = 0.2, and α = 0.8.The activation energies at different conversion obtained from the slopes of plots are presented in Figure 4d and Table 3. High correlation coefficients R 2 of fitting curves at different conversions are demonstrated.The values of E at each conversion determined by different model-free methods are almost the same.In addition, the values of E increase gradually at 0 ≤ α ≤ 0.1, followed by a sharp increase at 0.1 < α ≤ 0.2.Then E remains almost constant at 0.2 < α ≤ 0.8.Finally, a second sharp increase of E occurs at 0.8 < α ≤ 1.0.These ranges are coincidence with the position where the differences observed in Figure 4a-c    In Figure 4d, the temperatures corresponding to the four stages can be observed.Table 4 lists the pyrolysis process based upon conversion and temperature at different heating rates.In the case of 10 • C/min, the first stage in the range of 25-214 • C refers mainly to the degradation of resin [19,36].The second stage of 214-274 • C corresponds to the degradation of hemicellulose [36].The third stage in the range of 274-363 • C is attributed to the decomposition of cellulose [37].Finally, the last stage, temperatures above 363 • C, is mainly dominated by the decomposition of lignin [38].These stages are comparable to the analysis of resin, hemicellulose, cellulose, and lignin found in [19] and [39].

Estimated Initial Guess of Model Parameters
It is known that the activation energy which is related to the reaction rate represents the minimum energy required to start a reaction.Therefore, a large change of activation energy against conversion corresponds to the occurrence of several reactions [31,40].As shown in Figure 4d, E remains almost constant at 0.2 < α ≤ 0.8 whereas it increases significant in other ranges.Thus, the pyrolysis process at 0.2 < α ≤ 0.8 may be expressed by a single reaction model g(α), as listed in Table 5.The activation energies calculated by the slopes of CR plots coupled with 17 types of reaction models are listed in Table 5.
According to Chen et al. [31], if the mean value of activation energies for different heating rates at 0.2 < α ≤ 0.8 based upon one certain model is most close to the value by the model-free methods, then this model may dominate the reactions in this stage.As a result, the average value of E at 0.2 < α ≤ 0.8 based upon the 3rd-order reaction model g(α) = [1 − (1 − α) −2 )]/(−2) (183.88 kJ/mol) is closest to the value by the model-free methods (159.70 kJ/mol, as shown in Table 3).Thus, the 3rd-order reaction model may be the most appropriate model to express the pyrolysis process in the range of 0.2 < α ≤ 0.8.
Since the higher correlation coefficients R 2 , the E α obtained by FWO are used to determine the pre-exponential factor.The lnA α values evaluated by FWO with the 3rd-order reaction model are plotted in Figure 5.It can be seen that lnA α varies with E α , and a fitting line which shows the so-called kinetic compensation effect [41] can be obtained: ln A α = 0.17058E α + 11.797, R 2 = 0.993.According to Yao et al. [42], if a compensation effect is observed, the choice of 3rd-order reaction model is appropriate.At the same time, the compensation effect parameters, a = 0.17058 and b = 11.797, can describe the whole pyrolysis process well for that they are not affected by experimental conditions [43].Thus, the evaluations of lnA α based upon the compensation effect are listed in Table 3.The mean values of E and lnA of the four stages will be used as the initial guesses for the optimization process below.

Model Parameter Optimization by DE
The optimized parameters of the two models are listed in Table 6.As discussed above, the initial values of E and lnA for four components are the average values of E and lnA for the corresponding four stages which are listed in Table 3.According to Santaniello et al. [22], the initial values of E and lnA for tar are 124 kJ/mol and 17.9 ln s −1 , respectively.The initial values of reaction orders n are assumed to be 1.According to the manufacture and Munir et al. [44], the initial values of mass fractions of resin, hemicellulose and cellulose are assumed to be 0.1, 0.18, and 0.48, respectively.The initial values of φ, σ, γ, and θ are all assumed to be 0.5.The lower limit values of search range for parameters are set to 10% of the respective initial values.The upper limit values of search range for ni are set to 8 [19].The upper limit values for other parameters are set to 200% of the respective initial values.In Equation (18), it may cost a lot of computational time when data from all heating rates are used.Since 20 °C/min is a typical heating rate for slow pyrolysis [15], the experimental data of 20 °C/min and 30 °C/min are selected to obtain optimized parameters by DE.

Model Parameter Optimization by DE
The optimized parameters of the two models are listed in Table 6.As discussed above, the initial values of E and lnA for four components are the average values of E and lnA for the corresponding four stages which are listed in Table 3.According to Santaniello et al. [22], the initial values of E and lnA for tar are 124 kJ/mol and 17.9 ln s −1 , respectively.The initial values of reaction orders n are assumed to be 1.According to the manufacture and Munir et al. [44], the initial values of mass fractions of resin, hemicellulose and cellulose are assumed to be 0.1, 0.18, and 0.48, respectively.The initial values of ϕ, σ, γ, and θ are all assumed to be 0.5.The lower limit values of search range for parameters are set to 10% of the respective initial values.The upper limit values of search range for n i are set to 8 [19].The upper limit values for other parameters are set to 200% of the respective initial values.In Equation (18), it may cost a lot of computational time when data from all heating rates are used.Since 20 • C/min is a typical heating rate for slow pyrolysis [15], the experimental data of 20 • C/min and 30 • C/min are selected to obtain optimized parameters by DE.The optimized parameter values are listed in Table 6.It can be found that: (1) the activation energies for resin, hemicellulose, cellulose, and lignin obtained from Model I and Model II are in conformity with previous study by Li et al. [19].Li et al. reported the activation energies of 130-164 kJ/mol, 120-194 kJ/mol, 128-192 kJ/mol, and 144-240 kJ/mol for resin, hemicellulose, cellulose, and lignin, respectively.Meanwhile the pre-exponential factors also lie in the presented range; (2) According to the study of Munir et al. [44], the components contents in lignocellulosic biomass are: 12-24% hemicellulose, 43-54% cellulose, and 17-29% lignin.The estimated values of hemicellulose, cellulose and lignin obtained by Model I are 28.67%,33.88%, and 23.62%, respectively.The hemicellulose content is greater than the reported value whereas the cellulose content is less than the reported value.On the contrary, the values of 20.04% hemicellulose, 43.48% cellulose, and 26.78% lignin obtained by Model II are precisely in the reported ranges.
Comparisons between the experimental TG/DTG data and calculated values are presented in Figure 6.Good fitting qualities are observed at heating rates from 10 to 40 • C/min for both Model I and Model II.These optimized parameters can be found to satisfy not only the heating rates where they are calculated from (20 and 30 • C/min) but also other heating rates (10 and 40 • C/min).However, for Model I, predicted TG values are slightly larger than experimental data when the temperature is higher than 350 The optimized parameter values are listed in Table 6.It can be found that: (1) the activation energies for resin, hemicellulose, cellulose, and lignin obtained from Model Ι and Model II are in conformity with previous study by Li et al. [19].Li et al. reported the activation energies of 130-164 kJ/mol, 120-194 kJ/mol, 128-192 kJ/mol, and 144-240 kJ/mol for resin, hemicellulose, cellulose, and lignin, respectively.Meanwhile the pre-exponential factors also lie in the presented range; (2) According to the study of Munir et al. [44], the components contents in lignocellulosic biomass are: 12-24% hemicellulose, 43-54% cellulose, and 17-29% lignin.The estimated values of hemicellulose, cellulose and lignin obtained by Model Ι are 28.67%,33.88%, and 23.62%, respectively.The hemicellulose content is greater than the reported value whereas the cellulose content is less than the reported value.On the contrary, the values of 20.04% hemicellulose, 43.48% cellulose, and 26.78% lignin obtained by Model II are precisely in the reported ranges.
Comparisons between the experimental TG/DTG data and calculated values are presented in Figure 6.Good fitting qualities are observed at heating rates from 10 to 40 °C/min for both Model Ι and Model II.These optimized parameters can be found to satisfy not only the heating rates where they are calculated from (20 and 30 °C/min) but also other heating rates (10 and 40 °C/min).However, for Model Ι, predicted TG values are slightly larger than experimental data when the temperature is higher than 350 °C in all cases.Additionally, the peak values of the predicted |DTG| curves are smaller than the experimental data although the predicted results can capture the characteristic regions of the experimental |DTG| curves.In summary, the total TG-DTG curves calculated from both Model Ι and Model II fit the experimental data well, while components DTG curves obtained by models are different.Therefore, the effects of the consideration of secondary charring reactions in reaction models on components and products will be discussed in the following section.

Effects of Secondary Charring Reactions on the Components and Products
For each heating rate, the predicted TG curves of each component of the two models are plotted in Figure 8.It shows that the initial temperatures at which degradation for resin, hemicellulose, cellulose, lignin begin are around 225, 280, 360, and 400 °C for Model Ι and Model II, which are consistent with the temperatures corresponding to the components peak DTG.It indicates that the secondary charring reactions have little effect on the characteristic regions of pyrolysis process.However, components contents are very different when secondary charring reactions are considered or not.Differences in components contents will affect the estimated yields of products.Figure 9 presents gas, tar, and char yields at different heating rates calculated by Model Ι and Model II.It shows tiny differences in the yield of each product at different heating rates for two models.In the case of 10 °C/min, the gas and char yields increase with temperature, and remain almost constant above 700 °C.Additionally, higher gas and char yields are calculated by Model Ι.For tar yield, it increases when T < 350 °C, and then there will be a reduction.The most possible reason is that the tar cracking which may be one of the reactions involved in secondary charring occurs at around 400 °C [23].In summary, the total TG-DTG curves calculated from both Model I and Model II fit the experimental data well, while components DTG curves obtained by models are different.Therefore, the effects of the consideration of secondary charring reactions in reaction models on components and products will be discussed in the following section.

Effects of Secondary Charring Reactions on the Components and Products
For each heating rate, the predicted TG curves of each component of the two models are plotted in Figure 8.It shows that the initial temperatures at which degradation for resin, hemicellulose, cellulose, and lignin begin are around 225, 280, 360, and 400 • C for Model I and Model II, which are consistent with the temperatures corresponding to the components peak DTG.It indicates that the secondary charring reactions have little effect on the characteristic regions of pyrolysis process.However, components contents are very different when secondary charring reactions are considered or not.Differences in components contents will affect the estimated yields of products.Figure 9 presents gas, tar, and char yields at different heating rates calculated by Model I and Model II.It shows tiny differences in the yield of each product at different heating rates for two models.In the case of 10 • C/min, the gas and char yields increase with temperature, and remain almost constant above 700 • C. Additionally, higher gas and char yields are calculated by Model I.For tar yield, it increases when T < 350 • C, and then there will be a reduction.The most possible reason is that the tar cracking which may be one of the reactions involved in secondary charring occurs at around 400 • C [23].The char yield generated by MDF can be determined by weighing the residue amount after the experiment [45,46].In Figure 6, the mass loss (TG) was slow and its variation (DTG) could be ignored above 700 °C.Therefore, the solid residue amount at 800 °C was regard as the experimental char yield, which was marked by pentagram in Figure 9. Obviously, the char yield calculated by Model Ι is higher than experimental data.The percentage difference in the char yields between models and experiment is determined by: Diff.% = |  −   |   × 100% (23) where wm and wexp represent the char yield obtained from model and experiment, respectively.Table 7 lists the percentage differences in the char yields calculated by two models.The average error of char yield for Model Ι in four heating rates is 23.6% while the average error is 2% for Model II.Therefore, the components' contents and products' yields can be better estimated by Model II.In conclusion, it has been validated that the proposed model considering the secondary charring reactions coupled with the DE method can be applied to obtain the proper kinetic parameters.Additionally, the percentage difference in the char yields can be decreased from 23.6% to 2% when secondary charring reactions are considered.The char yield generated by MDF can be determined by weighing the residue amount after the experiment [45,46].In Figure 6, the mass loss (TG) was very slow and its variation (DTG) could be ignored above 700 • C. Therefore, the solid residue amount at 800 • C was regard as the experimental char yield, which was marked by pentagram in Figure 9. Obviously, the char yield calculated by Model I is higher than experimental data.The percentage difference in the char yields between models and experiment is determined by: Diff.% = w m − w exp w exp × 100% (23) where w m and w exp represent the char yield obtained from model and experiment, respectively.Table 7 lists the percentage differences in the char yields calculated by two models.The average error of char yield for Model I in four heating rates is 23.6% while the average error is 2% for Model II.Therefore, the components' contents and products' yields can be better estimated by Model II.In conclusion, it has been validated that the proposed model considering the secondary charring reactions coupled with the DE method can be applied to obtain the proper kinetic parameters.Additionally, the percentage difference in the char yields can be decreased from 23.6% to 2% when secondary charring reactions are considered.
illustrates the corresponding |DTG| and |DDTG| curves at 10 • C/min.There are four local minimum values in the |DDTG| curve, which correspond to the possible peak DTG of each sub-reaction.This phenomenon can also be observed for other heating rates.Therefore, four components are expected to react in the pyrolysis process of MDF.The temperatures corresponding to the peak DTG of each component may be 228, 303, 368, and 400 • C. Energies 2018, 11, x FOR PEER REVIEW 7 of 17

Figure 5 .
Figure 5. lnA versus E using FWO and the third order reaction model at 0.2 < α ≤ 0.8.

Figure 5 .
Figure 5. lnA versus E using FWO and the third order reaction model at 0.2 < α ≤ 0.8.

Figure 6 .
Figure 6.Comparisons between experimental TG/DTG data and predictions based on optimized parameters: (a) β = 10 • C/min; (b) β = 20 • C/min; (c) β = 30 • C/min; and (d) β = 40 • C/min.To further verify the reliability of the optimized parameters, the experimental curves of |DDTG| and |DTG| of MDF at 10 • C/min with the mass loss rates of four components calculated by Model I and Model II are presented in Figure 7.It can be found that the temperatures corresponding to the peak DTG of four components matches the temperatures corresponding to the four local minimum

Figure 7 .
Figure 7. Experimental curves of |DDTG| and |DTG| of MDF at 10 °C/min with four components by: Model Ι is shown in dashed lines; Model II is shown in solid lines.

Figure 7 .
Figure 7. Experimental curves of |DDTG| and |DTG| of MDF at 10 • C/min with four components by: Model I is shown in dashed lines; Model II is shown in solid lines.

Figure 8 .
Figure 8. Comparisons of components mass between Model Ι and Model II at four heating rates: (a) 10 °C/min; (b) 20 °C/min; (c) 30 °C/min; and (d) 40 °C/min.Model Ι is shown in dashed lines, and Model II is shown in solid lines.

Figure 8 .
Figure 8. Comparisons of components mass between Model I and Model II at four heating rates: (a) 10 • C/min; (b) 20 • C/min; (c) 30 • C/min; and (d) 40 • C/min.Model I is shown in dashed lines, and Model II is shown in solid lines.

Table 1 .
Proximate and ultimate analysis of MDF sample.
a By difference.

Table 3 .
Calculation results of E by FWO, KAS and Starink methods and evaluation values of lnA.

Table 3 .
Calculation results of E by FWO, KAS and Starink methods and evaluation values of lnA.

Table 4 .
Pyrolysis process based upon conversion and temperature at different heating rates.

Table 5 .
Calculation values of E based on the CR method for 0.2 ≤ α ≤ 0.8.

Table 6 .
Optimized values of kinetic parameters for Model I and Model II by DE.

Table 6 .
Optimized values of kinetic parameters for Model I and Model II by DE.
• C in all cases.Additionally, the peak values of the predicted |DTG| curves are smaller than the experimental data although the predicted results can capture the characteristic regions of the experimental |DTG| curves.

Table 7 .
The percentage differences in the char yields between models and experiment.

Table 7 .
The percentage differences in the char yields between models and experiment.