Kissinger method: the sequential approach and DAEM for kinetic study of rubber and gliricidia wood

This article is published under the Creative Commons CC-BY-ND License (http://creativecommons.org/licenses/by-nd/4.0/). This license permits use, distribution and reproduction, commercial and non-commercial, provided that the original work is properly cited and is not changed anyway. Abstract: Rubber and gliricidia are commonly used fuel wood types for industrial applications in Sri Lanka. The government has promoted the cultivation of gliricidia as an energy crop and rubber is generated as a forestry residue. A percentage of the thermal energy requirement in the industrial sector is fulfilled by fuel wood. Although these fuels are widely used in industry, scientific data related to their thermal decomposition is limited. These kinetic data are essential for the optimum usage of firewood and the design of biomass conversion systems. Kinetic parameters of pyrolysis are not available for rubber and gliricidia. Pyrolysis is a major step involved in thermal decomposition of biomass, and is an important process of deriving liquid and gaseous fuels for different applications. This study focused on the analysis of kinetic parameters of pyrolysis for rubber and gliricidia wood by sequential approach to Kissinger method and distributed activation energy model (DAEM).


INTRODUCTION
Biomass has gained importance as a renewable energy source, which is carbon neutral. In Sri Lanka the largest share of the primary energy supply, which is 42.5 % is given by biomass (Sustainable Energy Authority of Sri Lanka, 2014). Biomass fulfills 73.1 % of the energy requirement of the industrial sector in Sri Lanka (Sustainable Energy Authority of Sri Lanka, 2014). Gliricidia and rubber are two widely used wood types in industrial sector in Sri Lanka. Rubber is produced as an agricultural residue in rubber plantations, and gliricidia is planted as an energy crop and a supporting crop for other plants in tea plantations and minor export crops. Although the usage of biomass is high in industrial sector, technological interventions are limited in this field. Therefore, scientific studies related to fuel wood types are also limited. Scientific data related to thermal decomposition of these local fuel sources are important for designing and optimising the performance of thermal conversion systems.
Pyrolysis is the thermal degradation of solids in the absence of air or oxygen. It produces gases, solids and other liquid products. It is a process used to produce other fuels or chemicals and is a major step occurring in gasification and combustion. Biomass pyrolysis kinetics analysis is an important tool for the description of the effect of the process parameters on the feedstock conversion process (Zhang et al., 2006). Pyrolysis characteristics of different types of biomass vary with the biochemical composition. Many studies have been conducted to find the pyrolysis kinetics for fuel wood types such as Birch, Pine, Poplar and Eucalyptus (Van de Velden et al., 2010; Tapasvi et al., 2013;Hu et al., 2016). The kinetic data related to wood types such as gliricidia and rubber cannot be found in literature, and kinetic studies have not been done to evaluate the kinetic parameters of these wood types. Therefore, it is necessary to evaluate the kinetics of pyrolysis of these wood types for future developments in thermal conversion technology of gliricidia and rubber.
Thermal decomposition kinetics has been researched widely, and the kinetics obtained not only depends on the sample size, shape, flow rate and heating rate but also on the mathematics used to describe the experimental results. The kinetic parameters can be evaluated under isothermal or non-isothermal environments (Khawam & Flanagan, 2005). Either model free (iso-conversional methods) or model fitting methods can be used to predict the pyrolysis kinetics of biomass under isothermal or non-isothermal conditions (Martí-rosselló et al., 2016). A number of limitations exists in both model fitting and model free methods (Vyazovkin & Wight, 1999).
In model fitting methods kinetic parameters are defined according to previously defined reaction mechanisms. The requirement of assumptions on the reaction mechanism has reduced the popularity of model fitting methods (Khawam & Flanagan, 2005). Differential, Freeman Carroll (Freeman & Carroll, 1958) and Coats Redfern (Coats & Redfern, 1964) methods are some model fitting methods used for analysing nonisothermal pyrolysis kinetics. Jayaraman et al. (2017) used the Coats Redfern method to analyse the kinetics of coal biomass blends and observed an increase in activation energy and pre-exponential constant when compared to coal itself. This was described by the higher activation energy value and pre-exponential factor of biomass. The results obtained were in agreement with the values reported in literature. Coats Redfern method was used to evaluate the catalytic pyrolysis kinetics of biomass with the Friedman model free method, and Coats Redfern analysis showed that the reaction behaviour of catalytic biomass deviate from the simple first order reaction mechanism of virgin biomass to multistep reaction mechanism (Lu et al., 2009). Huang et al. (2011) presented a sequential method for the evaluation of non-isothermal kinetic parameters based on the method proposed by Kissinger (1957) on the maximum reaction rate. The method proposed in their research was used to evaluate kinetic parameters of seven kinds of biomass feedstocks. The proposed method was able to give results for the kinetic parameters, which are in the same range with the values mentioned in literature.
Iso-conversional methods are model free methods, which are used to describe the non-isothermal decomposition of solids (White et al., 2011). Isoconversional methods avoid the uncertainties added by assuming a specific thermal degradation model before evaluating the kinetic parameters. Flynn-Wall-Ozawa (FWO) (Ozawa, 1965;Flynn & Wall, 1966) method, modified Coats-Redfern method (Burnham & Braun, 1999) and Kissinger-Akahira-Sunose (KAS) method are some widely used iso-conversional methods (White et al., 2011). There are also other iso-conversional methods reported in literature (Vyazovkin, 2001;Hu et al., 2016). FWO, KAS and modified Coats-Redfern techniques have been applied to find the pyrolysis kinetics of soybean straw (Huang et al., 2016), while FWO and KAS methods have been used to find the pyrolysis kinetics of poplar wood (Slopiecka et al., 2011). FWO, KAS and Vyazovkin (2001) iso-conversional methods and pseudo component model fitting method have been used to evaluate pyrolysis kinetics of elephant grass (Collazzo et al., 2017). It was concluded that isoconversional methods are not suitable for evaluation of pyrolysis kinetics of elephant grass. Distributed activation energy model (DAEM) is another method, which can be used to evaluate kinetic parameters of biomass pyrolysis (White et al., 2011). The chemical complexity and variation of pyrolysis products have been taken into consideration in DAEM (Diblasi, 2008). The DAEM assumes that reaction proceeds through infinite number of independent parallel reactions, which have different activation energies. The activation energy variation is presented by a continuous distribution function. Distribution free or distribution fitting methods are used to evaluate the kinetic parameters in DAEM (Cai et al., 2014). Gaussian (Gašparoviè et al., 2012), Weibull (Kuo-Chao et al., 2009) or Gamma distribution (Burnham & Braun, 1999) functions have been used to describe the distribution of activation energy in distribution fitting method. Gauss distribution function was the widely used distribution function in DAEM (Gašparoviè et al., 2012). Pyrolysis kinetics of rice husk, bamboo and pine wood were studied using distribution fitting methods (Hu et al., 2016). The study assumed a Gaussian distribution for the activation energy and was able to obtain good match with experimental data. Gašparoviè et al. (2012) used Gaussian distribution to find pyrolysis kientics of wood (Gašparoviè et al., 2012). Miura and Maki (1998) approach is a distribution free approach of DAEM. This method was successfully used to evaluate pyrolysis kinetics of different biomass types Maki approach (Miura & Maki, 1998) of DAEM was used to analyse pyrolysis kinetics using the data provided by the thermo gravimetric analysis (TGA) of two widely used industrial fuel wood types, rubber and gliricidia. The pyrolysis kinetics obtained is useful for the determination of thermal decomposition characteristics of these wood fuels. Hence, these values can be used in designing and optimisation of thermal conversion systems and optimum utilisation of these local fuel woods in energy production.

Experimental study
Gliricidia and rubber wood was used in this experimental series. Gliricidia wood supplied for thermal applications and rubber wood, which was available in the market were used in the experiment. Proximate and ultimate analyses of gliricidia and rubber wood are shown in Table 1. Thermo gravimetric analyses were carried out in TGA analyser SDT Q600. Approximately 10 mg was used for the experiments. Three experiments were carried out for each wood type at 10 o Cmin -1 , 20 o Cmin -1 and 30 o Cmin -1 heating rate. These samples were heated from room temperature to the target temperature of 800 o C. The inert gas was nitrogen, which was purged at a flow rate of 100 mLmin -1 for all the experiments.

Kissinger method: the sequential approach
The procedure presented in Huang et al. (2011) for evaluation of pyrolysis kinetics was used in this study.
The reaction rate of pyrolysis is expressed by the following expression:

Kissinger method the sequential approach
The procedure presented in (Huang et al., 2011) for evaluation of pyrolysis kinetics was used in this study.
The reaction rate of the pyrolysis is expressed by the following expression The x is the mass conversion, t is time, k is the reaction rate constant and f(x) is the conversion function. The x is denoted as conversion hereafter throughout the article.
x in the above equation is expressed as m 0 is the initial mass of the sample, m f is the final mass of the sample and m is the mass of the sample at time t. The reaction rate constant is expressed by the following formula. ...(01) x is the mass conversion, t is time, k is the reaction rate constant and f(x) is the conversion function. The x is denoted as conversion hereafter throughout the article. x in the above equation is expressed as The x is the mass conversion, t is time, k is the reaction ra function. The x is denoted as conversion hereafter throughout The reaction rate constant is expressed by the following formula; The x is the mass conversion, t is time, k is the reaction rate co function. The x is denoted as conversion hereafter throughout the a x in the above equation is expressed as m 0 is the initial mass of the sample, m f is the final mass of the s sample at time t. The reaction rate constant is expressed by the foll If it is assumed that the decomposition of solids follows n th order o where A is the frequency factor, E is the activation energy, R is the universal gas constant and T is the temperature.
If it is assumed that the decomposition of solids follows n th order of conversion   =  The x is the mass conversion, t is time, k is the reaction rate function. The x is denoted as conversion hereafter throughout the x in the above equation is expressed as m 0 is the initial mass of the sample, m f is the final mass of the sample at time t. The reaction rate constant is expressed by the fo If it is assumed that the decomposition of solids follows n th order With the approximations given in (Otero et al., 2008) the integral written as Differential of equation (5) at the maximum reaction rate will give The heating rate is

 =  
Substituting equation (5) and (8) in equation (7)   At the maximum conversion rate Differential of equation (5) at the maximum reaction rate will give t The heating rate is Substituting equation (5) and (8) in equation (7)   Differential of equation (5) at the maximum reaction rate will give the following relationship.
With the approximations given in (Otero et al., 2008) the integral for written as Differential of equation (5) at the maximum reaction rate will give th The heating rate is Substituting equation (5) and (8) in equation (7)   With the approximations given in (Otero et al., 2008) the integral form of the above equation can be written as Differential of equation (5) at the maximum reaction rate will give the following relationship.
The heating rate is Substituting equation (5) and (8) in equation (7)   The heating rate β is The heating rate is Substituting equation (5) and (8) in equation (7)   At the maximum conversion rate By rearranging equations (6) and (10) the order of the reaction can easily be found using TGA data.
... (08) Substituting equations (5) and (8) in equation (7) The heating rate is Substituting equation (5) and (8) in equation (7)   At the maximum conversion rate By rearranging equations (6) and (10) the order of the reaction can easily be found using TGA data.
1 −    =    (11) ... (09) At the maximum conversion rate The heating rate is Substituting equation (5) and (8) in equation (7)   At the maximum conversion rate By rearranging equations (6) and (10) the order of the reaction can easily be found using TGA data.
The heating rate is Substituting equation (5) and (8) in equation (7)   At the maximum conversion rate By rearranging equations (6) and (10) the order of the reaction can easily be found using TGA data.
Substituting equation (5) and (8) in equation (7)   At the maximum conversion rate By rearranging equations (6) and (10) the order of the reaction can easily be found using TGA data.
1 −    =    (11) ... (11) Once the order of the equation is obtained the activation energy E and the pre-exponential factor can be calculated from the following relationships.
Once the order of the equation is obtained the activation energy E and the pre-exponential factor can be calculated from the following relationships.
For a first order reaction integration of equation (5) gives By differentiating equation (5) At the maximum reaction rate      ⁄ = 0 ... (12) Once the order of the equation is obtained the activation energy E and the pre-exponential factor can be calculated from the following relationships.
For a first order reaction integration of equation (5) gives By differentiating equation (5) At the maximum reaction rate      ⁄ = 0 ... (13) For a first order reaction integration of equation (5) gives Once the order of the equation is obtained the activation energy E and the pre-exponential factor can be calculated from the following relationships.
For a first order reaction integration of equation (5) gives By differentiating equation (5) At the maximum reaction rate      ⁄ = 0 Since ... (14) By differentiating equation (5) Once the order of the equation is obtained the activation energy E and the pre-exponential factor can be calculated from the following relationships.
For a first order reaction integration of equation (5) gives By differentiating equation (5) At the maximum reaction rate      ⁄ = 0 Since ... (15) At the maximum reaction rate Once the order of the equation is obtained the activation energy E and the pre-exponential factor can be calculated from the following relationships.
For a first order reaction integration of equation (5) gives By differentiating equation (5) At the maximum reaction rate      ⁄ = 0 Equation (16)

leads to
Once the order of the equation is obtained the activation energy E and the pre-exponential factor can be calculated from the following relationships.

The DAEM method
Fraction of solid converted to volatile matter is used to determine the py method (23).

The DAEM method
Fraction of solid converted to volatile matter is used to determine the py method (23).

The DAEM method
Fraction of solid converted to volatile matter is used to determine the pyro method (23).

...(21)
First order or n th order reaction kinetics can be obtained by the above method of calculation. The best fitting parameters can be calculated by the nonlinear square fit method using the below formula.

The DAEM method
Fraction of solid converted to volatile matter is used to determin method (23).
... (22) In this study it was assumed that the pyrolysis kinetics obey first order reaction. Therefore, equations which describe the first order decay were used to find the kinetic parameters. TGA data at a heating rate of 20 o Cmin -1 was used to calculate model based nonisothermal pyrolysis kinetics. The numerical calculations were accomplished by minimising the objective function (equation 22) using inbuilt solver fmin-search in Matlab software. This solver has the ability to optimise multiple parameters. The calculated values of E (equation 19) and A (equation 21) were given as the initial guess values for the solver.

The DAEM method
The fraction of solid converted to volatile matter is used to determine the pyrolysis kinetics in DAEM method (equation 23). In the above equation 12 In the above equation The inner dT integral and the outer dE integral makes it difficult to find exact analytical solution to DAEM. Therefore, many mathematical approximations have been presented in literature. Due to the complex structure of the DAEM the parameter estimation also has difficulties. The parameter estimation of DAEM is classified as either distribution free or distribution fitting method. Miura and Maki proposed a revised DAEM to estimate the distribution of f(E) without prior assumptions on f(E) (Miura, 1995;Miura & Maki, 1998). Therefore, Miura and Maki (1998) method was used in this study for the estimation of activation energy distribution.
Using mathematical simplifications and approximations proposed by Miura and Maki, equation (23) is presented in the below form  The value of E is calculated from the gradient and the value of A is calculated from the intercept from the relationship using equation (25) ... (24) The inner dT integral and the outer dE integral makes it difficult to find an exact analytical solution to DAEM. Therefore, many mathematical approximations have been presented in literature. Due to the complex structure of the DAEM, parameter estimation also has difficulties. The parameter estimation of DAEM is classified as either a distribution free or distribution fitting method. Miura and Maki proposed a revised DAEM to estimate the distribution of f(E) without prior assumptions on f(E) (Miura, 1995;Miura & Maki, 1998). Therefore, Miura and Maki (1998) method was used in this study for the estimation of activation energy distribution.
Using mathematical simplifications and approximations proposed by Miura and Maki (1998), equation (23) is presented in the below form 12 In the above equation The inner dT integral and the outer dE integral makes it difficult to find exact analytical solution to DAEM. Therefore, many mathematical approximations have been presented in literature. Due to the complex structure of the DAEM the parameter estimation also has difficulties. The parameter estimation of DAEM is classified as either distribution free or distribution fitting method. Miura and Maki proposed a revised DAEM to estimate the distribution of f(E) without prior assumptions on f(E) (Miura, 1995;Miura & Maki, 1998). Therefore, Miura and Maki (1998) method was used in this study for the estimation of activation energy distribution.
Using mathematical simplifications and approximations proposed by Miura and Maki, equation (23) is presented in the below form  The value of E is calculated from the gradient and the value of A is calculated from the intercept from the relationship using equation (25) ... (25) Variation of ln(β/T 2 ) versus 1/T is a straight line where E can be calculated. The following procedure is used for the calculation of E and f(E)  The variation of ln(β/T 2 ) versus 1/T is evaluated at heating rates of 10, 20 and 30 o Cmin -1 for the same levels of conversion(1-x).  The value of E is calculated from the gradient and the value of A is calculated from the intercept from the relationship using equation (25).

Pyrolysis kinetics by Kissinger method: the sequential approach
Kinetic parameters of pyrolysis of gliricidia and rubber were calculated using TGA data. Table 2 shows the kinetic parameters of gliricidia and rubber wood calculated by the sequential approach based Kissinger method for a heating rate of 20 o Cmin -1 . Figure 1 shows the reaction rate of gliricidia obtained by the experiments, and reaction rate obtained by the kinetic parameters calculated by the TGA data. For gliricidia, the activation energy is 107.19 kJmol -1 and the pre-exponential factor is 8.88 × 10 8 s -1 . Figure 2 shows the reaction rate profiles for rubber wood. Rubber has an activation energy of 83.44 kJmol -1 and a pre-exponential factor of 9.5 × 10 4 s -1 . The rubber wood decomposition curve clearly shows two separate peaks, whereas the gliricidia decomposition curve shows a shoulder at low temperature. This low temperature peak in rubber wood reaction rate and the low temperature shoulder in gliricidia reaction rate profile can be attributed to the

Pyrolysis kinetics by Miura and Maki method
Plot of ln(β/T 2 ) versus (1/T) for rubber wood is shown in Figure 3 and gliricidia is shown in Figure 4 for conversion levels range from 0.05 to 0.95. the whole range of temperatures of pyrolysis. This effect of decomposition of pseudo components was not visible from the predicted kinetics.

Pyrolysis kinetics by Miura and Maki method
Plot of ln(β/T 2 ) versus (1/T) for gliricidia and rubber wood are shown in Figure 3 and Figure 4, respectively for conversion levels ranging from 0.05 to 0.95. 16 Figure 3. ln(β/T 2 ) versus1/T at different conversion levels for gliricidia Figure 3: ln(β/T 2 ) versus1/T at different conversion levels for gliricidia hemicellulose decomposition (Diblasi, 2008). The highest reaction rate occurs at the cellulose devolatilisation peak. This corresponds to the 2 nd peak in Figure 2 and 1 st peak in Figure 1. The temperatures corresponding to these peaks are 635 K for gliricidia and 634 K for rubber wood. These temperatures agree with the temperature range reported in literature for occurrence of maximum reaction rate corresponding to cellulose pyrolysis (Yang et al., 2007). The decomposition of lignin spreads over A good correlation can be observed for conversion range between 0.15 and 0.8 for gliricidia. Correlations are high for conversions from 0.05 to 0.85 for rubber (Table 3). Both wood types showed an increase in activation energy at 0.85 conversion and a sharp decrease above 0.85 conversion (Figures 5 and 6). A similar observation was made for different wood types (Shen et al., 2011) and other biomass types (Chen et al., 2013) at higher conversions. Further decomposition of charcoal formed during the main pyrolysis stage is the reason for this increase in activation energy.
High correlation factor emphasises the linear and parallel relationships at conversion levels between 0.15 and 0.80 for gliricidia, and 0.05 and 0.85 for rubber wood. This implies that single or a set of parallel first order reactions occur during the pyrolysis of gliricidia and rubber woods (Miura & Maki, 1998). Therefore, DAEM method can be used to describe the pyrolysis kinetics of gliricidia and rubber woods.

CONCLUSION
The activation energy of rubber evaluated by the sequential approach to Kissinger method is 83.44 kJmol -1 , and for gliricidia it is 107.19 kJmol -1 . Pre-exponential factor is 8.88 × 10 8 s -1 and 9.5 × 10 4 s -1 for gliricidia and rubber, respectively.
Model free DAEM kinetics was obtained using Miura and Maki (1998) approach. At conversion below 0.15 the correlation factors are low for gliricidia wood. Therefore, the kinetic parameters obtained by the Miura   (2) and Maki (1998) approximation of DAEM cannot be used to describe pyrolysis kinetics of gliricidia at low conversions. At high conversions the relationship between (1/T) and (β/T 2 ) is not linear for both wood types. Therefore, the kinetic parameters obtained for conversions higher than 0.85 are not reliable for both gliricidia and rubber. High correlation factors demonstrate that Miura and Maki (1998)   The pre-exponential factor varies between 3.45x10 14 and 6.78x10 18 s -1 for gliricidia. Rubber wood has pre-exponential factor which varies between 6.45x10 08 and 2.61x10 12 s -1 .
High correlation factor emphasizes the linear and parallel relationships at conversion levels between 0.15 and 0.80 for gliricidia and 0.05 and 0.85 for rubber wood. This implies that single or a set of parallel first order reactions occur during the pyrolysis of gliricidia and rubber woods (Miura & Figure 6: Variation of activation energy with conversion for rubber E (kJmol -1 ) 19 gliricidia and considered in the following discussion. The kinetic parameters calculated between 0.05 and 0.85 conversion values were discussed in the following section for rubber wood. The activation energy values of gliricidia wood ranged from 190.58 kJmol -1 to 230.57 kJmol -1 . The activation energy value varied from 111.52 kJmol -1 to 179.07 kJmol -1 for rubber wood.