1464
Ind. Eng. Chem. Res. 2009, 48, 1464–1467
Determination of Distributed Activation Energy Model Kinetic Parameters Using Simulated Annealing Optimization Method for Nonisothermal Pyrolysis of Lignin Thilakavathi Mani, Pulikesi Murugan, and Nader Mahinpey* Faculty of Engineering, UniVersity of Regina, Saskatchewan, Canada S4S 0A2
The distributed activation energy model (DAEM) has been shown to be more descriptive of the pyrolysis reaction than other applicable models. In this study, the temperature dependency of the preexponential factor has been included in the DAEM. The model equation has been solved using Simpson’s 1/3 rule, and the kinetic parameters were determined using an optimization method. Simulated annealing method has been used to determine the DAEM kinetic parameters for the nonisothermal pyrolysis of lignin using thermogravimetric analysis (TGA) data. The nonisothermal pyrolysis of lignin was conducted at three different heating rates of 5, 10, and 15 °C/min under nitrogen atmosphere. Predicted results from the optimum kinetic parameters have been compared with the experimental data. The DAEM equation predicts the experimental data very well for different heating rates. Introduction Due to the scarcity of natural fuel resources and the abundance of biomass, interest in the conversion of biomass to useful energy products has evolved. Lignin is one of the major components present in biomass that can be converted into biooil through thermochemical processes.1,2 Lignin is the waste material obtained from the pretreatment of biomass for the fermentation process of producing biofuel. Since the presence of lignin affects the growth of microorganisms, lignin has to be removed from biomass before fermentation. The disposal of this waste lignin is one of the major issues in the biochemical conversion of biomass. Pyrolysis is the primary thermochemical process to convert biomass into valuable products, namely, solid char, liquid, and gas.3 Kinetics of pyrolysis reactions can be used for design and control of thermal decomposition processes. Application of the model-fitting method to nonisothermal data results in highly uncertain kinetic parameters.4 The alternative to the model-fitting method is the consideration of variation of the activation energy and the preexponential factor (frequency factor) with respect to temperature in the model equation.5 Therefore, for kinetic studies of complex reactions such as pyrolysis of biomass, the distributed activation energy model (DAEM) has been widely utilized.6,7 The model has been applied to represent the change in the activation energy of biomass with respect to temperature and overall conversion. The DAEM has some mathematical difficulties in determining kinetic parameters because of the complex structure of the model equation. In previous studies, these parameters were established by different methods, such as Marquardt nonlinear regression method, flexible simplex method, direct search method, and pattern search method.8 In the present work, the simulated annealing procedure is used for determining the kinetic parameters in the DAEM. Distributed Activation Energy Model. The distributed activation energy model is represented as follows when it is applied to represent the change in total volatiles:9 1-X)
∫
∞
0
(
exp -
∫
T
0
)
k0 -E⁄RT e dT f(E) dE β
(1)
where X is the total conversion of volatiles at temperature T, f(E) is the distribution curve of the activation energy to represent the differences in the activation energies of many first-order irreversible reactions, β is the linear heating rate of the pyrolysis reaction, and k0 is the preexponential factor. Usually, the preexponential factor is considered to be a constant over the temperature range studied. However, in some solid-state reactions, the preexponential factor is connected with the temperature through the following relationship:10 k0 ) k0 ′ Tm
(2)
where k0′ is a constant and the exponent m ranges from -1.5 to +2.5.11 In this study, the dependence of the frequency factor on temperature has been considered for the DAEM, which is represented as follows: 1-X)
∫
∞
0
(
exp -
k0 ′ β
∫
T
0
)
Tme-E⁄RT dT f(E) dE
(3)
The inner integral ∫0T Tme- E/RT dT in eq 3 does not have an exact analytical solution. Therefore, many authors have proposed approximations for the general temperature integral.11,12 In this work, the approximation proposed by Cai and Liu,12 given by eq 4, has been used for the general temperature integral.
∫
T
0
Tme-E⁄RT dT ) RTm+2 -E⁄RT 0.99954E + (0.58058 - 0.044967m)RT e (4) E E + (2.54 + 0.94057m)RT
Usually, f(E) is taken to be a Gaussian distribution with a mean activation energy E0 and a standard deviation σ, 1 -(E - E0)2⁄2σ2 e (5) √ σ 2π Simulated Annealing Optimization Method. Equation 3 was solved by numerical methods using Simpson’s 1/3 rule. The parameters k0′, E0, σ, and m are calculated using an optimization technique with initial values given in Table 1. The objective function (OF) to optimize kinetic parameters of a nonisothermal DAEM equation is f(E) )
n
* To whom correspondence should be addressed. Tel.: (306) 5854490. Fax: (306) 585-4855. E-mail:
[email protected].
OF )
∑ (X
2 i,cal - Xi,exp)
i)1
10.1021/ie8013605 CCC: $40.75 2009 American Chemical Society Published on Web 12/22/2008
(6)
Ind. Eng. Chem. Res., Vol. 48, No. 3, 2009 1465 Table 1. DAEM Kinetic Parameters for the Lignin Pyrolysis k0′ (1/(min.K))
E0 (kJ/mol)
σ (kJ/mol)
m
objective function value
10000 92670
100 135
10 33
1 2.8
16.9122 0.0537
initial value final value
where Xi,cal and Xi,exp are calculated and experimental values of the lignin conversion, respectively, and n is the data number. To determine k0′ and m values by minimizing the OF value, a certain optimization technique is needed. It is difficult to obtain information about gradient or higher derivatives of the objective function without explicit expression. Therefore, the optimization algorithm should be derivative-free and robust with respect to local optima and should require the least function evaluations for finding the optimum.8 In this work, we used the simulated annealing method for determining optimum kinetic parameters in the DAEM. The objective function was to minimize the sum of the square of errors between the experimental data and estimated values at all times. The parameters are varied one at a time randomly to obtain a new set of parameters. With each set of parameters, the objective function was determined and the difference in the objective function (∆f) with the old and new sets of parameters was calculated. If the new set of values improved the objective function, the move was accepted. Otherwise, the move was accepted with a probability of exp(∆f/T), where T is the simulated annealing temperature, a dummy variable that is used to control the acceptance of uphill moves. Initially, T was fixed at a higher value and periodically annealed by a proportional cooling schedule in the outer loop. At any specific temperature, the parameters were randomly varied a number of times in the inner loop. Thus, the optimum parameter values were obtained after T reached a desired lower value. The algorithm for the simulated annealing procedure can be found elsewhere.13 The simulated annealing method is a simple procedure in which no derivative of the function is required. No need for solving the matrix as in the case of the Marquardt nonlinear regression method. In other methods, such as pattern search methods, calculation steps are more for finding the search directions, etc., whereas, in the simulated annealing method, the variables are changed randomly and there are only two conditions: whether to accept the probability or not, depending on the error value. The program codes are given in the Appendix from which it is clear that this procedure has a lower number of calculation steps to obtain optimum values. Experiments Alcel lignin obtained from hardwood species from New Brunswick, Canada, has been used for this study. The results of the elemental analysis of the lignin were given in the literature.14 The experimental apparatus used for the thermogravimetric analysis (TGA) studies consisted of a DuPont 951 thermogravimetric analyzer and a DuPont 2100 data analysis system. Lignin pyrolysis was performed under nitrogen atmosphere at a flow rate of 45 cm3/min. The nonisothermal pyrolysis of lignin was conducted at three different heating rates of 5, 10, and 15 °C/min. All experiments were performed in triplicate. Results and Discussion Miura and Maki7 proposed a simple method for solving the distributed activation energy model. Using that method, the lower (83 kJ/mol) and upper (195 kJ/mol) values of the activation energy for lignin pyrolysis have been determined.
Figure 1. Objective function value versus iterative.
Figure 2. Comparison of experimental and predicted values of lignin pyrolysis at 5 °C/min.
The data were similar to those for the Kraft lignin, which has an activation energy in the range of 80-158 kJ/mol.15 The model equation has been solved with Simpson’s 1/3 rule using the limits for activation energy. The optimum kinetic parameters were determined using the simulated annealing method. The optimization program was written in MATLAB codes, which were used to compute the optimum values of kinetic parameters. The simulated annealing algorithm can handle optimization problems that are difficult to solve. Figure 1 shows the decrease in the objective function value with an increase in iteration number. Initially the function value decreased smoothly in a linear fashion and then sharply reduced. The fluctuation occurred because of the acceptance of the move with probability when the error was higher than the previous value. It took approximately 500 iterations to get the minimum objective function, and the program running time was about 10 min. This shows that the simulated annealing algorithm works well for complex equations such as the DAEM. The optimum values of the kinetic parameters, along with the function values, are given in Table 1. Figure 2 shows the comparison of experimental and predicted values of lignin conversion at 5 °C/min. Similar curves were obtained for other heating rates for lignin pyrolysis. The correlation coefficient between calculated and experimental data for the heating rates of 5, 10, and 15 °C/min are 0.9998, 0.9996, and 0.9998, respectively. The DAEM equation fits the experimental data very well. The comparison of activation energy distribution between calculated distributed function7 and Gaussian function is given in Figure 3. It is evident that the distribution of the activation energy calculated on the basis of the experimental runs of lignin pyrolysis was somewhat different, showing maximum activation energy in the range of 158-170 kJ/mol. Values of m ranging from 0 to 2.5 have been proposed for the case of reactions of desorption of gases from the surface of
1466 Ind. Eng. Chem. Res., Vol. 48, No. 3, 2009
Figure 3. Distributed function f(E) versus activation energy.
Figure 4. Variation of the frequency factor with activation energy.
solids, and the range -3/2 to 0 have been proposed for shrinkage processes depending on the sintering mechanism.16 During pyrolysis, various processes such as evaporation, visbreaking, and thermal cracking, etc., occur. This may be the reason for the optimum value of m ) 2.8, for lignin pyrolysis. From the results, it is clear that frequency factor k0 is not constant and varies with temperature. Kinetic analysis performed with this assumption made good agreement of model curve with experimental data. The variation of the frequency factor with activation energy is given in Figure 4. Therefore, it is concluded that the frequency factor for complex reactions such as lignin pyrolysis varies with temperature.
Nomenclature
The distributed activation energy model equation fits well with the pyrolysis of lignin using the Gaussian distribution function f(E). The temperature dependency of the preexponential factor was included in the DAEM. Kinetic parameters of the DAEM were determined from nonisothermal TGA data of lignin pyrolysis through the simulated annealing optimization method. The simulated annealing algorithm works well for determining parameters in complex equations such as the DAEM.
E ) activation energy (kJ/mol) k0 ) frequency factor (1/min) k0′ ) constant (1/(min Km)) m ) exponent of T OF ) objective function R ) universal gas constant T ) temperature (K) X ) lignin conversion β ) heating rate (K/min) σ ) standard deviation (kJ/mol)
Acknowledgment
Literature Cited
The authors wish to acknowledge the financial support of the Natural Science and Engineering Research Council of Canada (NSERC).
(1) Yoshida, T.; Matsumura, Y. Gasification of Cellulose, Xylan, and Lignin Mixtures in Supercritical Water. Ind. Eng. Chem. Res. 2001, 40, 5469. (2) López Pasquali, C. E.; Herrera, H. Pyrolysis of Lignin and IR Analysis of Residues. Thermochim. Acta 1997, 293, 39. (3) Wang, G.; Li, W.; Li, B.; Chen, H. TG Study on Pyrolysis of Biomass and its Three Components under Syngas. Fuel 2008, 87, 552. (4) Vyazovkin, S.; Wight, C. A. Model-Free and Model-Fitting Approaches to Kinetic Analysis of Isothermal and Nonisothermal Data. Thermochim. Acta 1999, 340–341, 53.
Conclusions
Appendix Program codes for the simulated annealing optimization procedure are given as follows:
Ind. Eng. Chem. Res., Vol. 48, No. 3, 2009 1467 (5) Vyazovkin, S. Computational Aspects of Kinetic Analysis, Part C. The ICTAC Kinetics Project—The Light at the End of the Tunnel? Thermochim. Acta 2000, 355, 155. (6) Cai, J.; He, F.; Yao, F. Nonisothermal nth-Order DAEM Equation and Its Parametric StudysUse in the Kinetic Analysis of Biomass Pyrolysis. J. Math. Chem. 2007, 42, 949. (7) Miura, K.; Maki, T. A Simple Method for Estimating f(E) and k0(E) in the Distributed Activation Energy Model. Energy Fuels 1998, 12, 864. (8) Cai, J.; Ji, L. Pattern Search Method for Determination of DAEM Kinetic Parameters from Nonisothermal TGA Data of Biomass. J. Math. Chem. 2007, 42, 547. (9) Miura, K. A New and Simple Method to Estimate f(E) and k0(E) in the Distributed Activation Energy Model from Three Sets of Experimental Data. Energy Fuels 1995, 9, 302. (10) Cai, J.; Liu, R. New Distributed Activation Energy Model: Numerical Solution and Application to Pyrolysis Kinetics of Some Types of Biomass. Bioresour. Technol. 2008, 99, 2795. (11) Wanjun, T.; Yuwen, L.; Xil, Y.; Zhiyong, W.; Cunxin, W. Approximate Formulae for Calculation of the Integral ∫T0 Tme- E/RT dT. J. Therm. Anal. Calorim. 2005, 81, 347.
(12) Cai, J.; Liu, R. New Approximation for the General Temperature Integral. J. Therm. Anal. Calorim. 2007, 90, 469. (13) Thilakavathi, M.; Basak, T.; Panda, T. Modeling of Esterase Production from Saccharomyces cereVisiae. J. Microbiol. Biotechnol. 2008, 18, 889. (14) Murugan, P.; Mahinpey, N.; Johnson, K.; Wilson, M. Kinetics of the Pyrolysis of Lignin Using Thermogravimetric and Differential Scanning Calorimetry Methods. Energy Fuels 2008, 22, 2720. (15) Ferdous, D.; Dalai, A. K.; Bej, S. K.; Thring, R. W. Pyrolysis of Lignins: Experimental and Kinetics Studies. Energy Fuels 2002, 16, 1405. (16) Criado, J. M.; Pe´rez-Maqueda, L. A.; Sa´nchez-Jime´nez, P. E. Dependence of the Preexponential Factor on Temperature. J. Therm. Anal. Calorim. 2005, 82, 671.
ReceiVed for reView September 9, 2008 ReVised manuscript receiVed November 4, 2008 Accepted November 6, 2008 IE8013605
