Energy Fuels 2010, 24, 2841–2847 Published on Web 04/23/2010
: DOI:10.1021/ef1001265
Improved Method for the Determination of Kinetic Parameters from Non-isothermal Thermogravimetric Analysis (TGA) Data James Hillier, Trent Bezzant, and Thomas H. Fletcher* Chemical Engineering Department, 350 CB, Brigham Young University, Provo, Utah 84602 Received February 1, 2010. Revised Manuscript Received April 8, 2010
Thermogravimetric analysis (TGA) at a constant heating rate is a popular technique for characterizing reaction behavior. Determination of kinetic coefficients from such TGA data has given rise to a large range of reported activation energies for similar samples. In general, rate coefficients are determined using either a linearized integral method or a derivative method. A new method incorporating several well-known techniques is presented that avoids the problems with the linearized techniques as they are generally used. This paper compares results from the new method with results using either the integral or derivative method. Mathematically generated data were used with known kinetic parameters for the comparison. This paper demonstrates that the new method shows better agreement with the data than either the integral or derivative method, and in some cases, the error was an order of magnitude lower. This paper also shows that the integral method outperforms the derivative method for the generated data.
coefficient against the inverse of temperature. This reduction method linearizes the data, and the kinetic parameters can then be determined from the slope and intercept of the resulting linearized equation. Traditional non-isothermal experiments also involve linearizing an equation in some manner and fitting the resulting straight line to obtain the activation energy, pre-exponential factor, and reaction order.4,9-12 Isothermal and non-isothermal data require different methods to regress the kinetic parameters. The remainder of this paper will focus on two traditional and a newly proposed method for analyzing non-isothermal kinetic data. The starting equation in this paper for these three methods, shown in eq 1, is an assumed first-order global reaction with an Arrhenius form of the rate constant.
Introduction Thermogravimetric analyzers (TGAs) have been used by many investigators to measure mass release as a function of time and temperature to determine kinetic parameters.1-8 Isothermal techniques involve introducing the sample into a furnace held at constant temperature and observing the change of mass with time. Non-isothermal studies involve heating a sample in an apparatus, often linearly with time (i.e., a constant heating rate), and continually weighing it to develop a mass versus time curve. Isothermal and non-isothermal methods differ in the assumptions made and method of collecting data, and differences can arise from the imperfections in the data that are accentuated by a specific method. The non-isothermal method is often used for pyrolysis experiments. Once a mass (or conversion) versus time curve has been generated using the TGA with an experiment, the particular method of determining the activation energy can lead to differing results. Data from isothermal experiments are traditionally plotted in the form of the natural log of the reaction
dm ¼ - Ae- E=RT m dt
Traditionally, there have been two classes of methods using non-isothermal data to determine the kinetic parameters used in this equation.13 One method simply requires a linearization of this equation by dividing by the mass and then taking the natural log or log10 of both sides. The resulting equation is shown here as eq 2 and is referred to as the derivative method in this paper. ! dm E ð2Þ m ¼ lnðAÞ ln dt RT
*To whom correspondence should be addressed: Chemical Engineering Department, 350 CB, Brigham Young University, Provo, UT 84602. E-mail:
[email protected]. (1) Doyle, C. D. Kinetic analysis of thermogravimetric data. J. Appl. Polym. Sci. 1961, 5 (15), 285–292. (2) Freeman, E. S.; Carroll, B. The application of thermoanalytical techniques to reaction kinetics: The thermogravimetric evaluation of the kinetics of the decomposition of calcium oxalate monohydrate. J. Phys. Chem. 1958, 62 (4), 394–397. (3) Coats, A. W.; Redfern, J. P. Kinetic parameters from thermogravimetric data. Nature 1964, 201 (4914), 68–69. (4) Abu-Qudais, M.; Jaber, J. O.; Sawalha, S. Kinetics of pyrolysis of Attarat oil shale by thermogravimetry. Oil Shale 2005, 22 (1), 51–63. (5) Campbell, J. H.; Koskinas, G. H.; Stout, N. D. Kinetics of oil generation from Colorado oil shale. Fuel 1978, 57 (6), 372–376. (6) Burnham, A. K.; Braun, R. L. Global kinetic analysis of complex materials. Energy Fuels 1999, 13 (1), 1–22. (7) Li, S.; Yue, C. Study of pyrolysis kinetics of oil shale. Fuel 2003, 82 (3), 337–342. (8) Shao, J.; Yan, R.; Chen, H.; Wang, B.; Lee, D. H.; Liang, D. T. Pyrolysis characteristics and kinetics of sewage sludge by thermogravimetry Fourier transform infrared analysis. Energy Fuels 2008, 22 (1), 38–45. r 2010 American Chemical Society
ð1Þ
=
(9) Lee, S. Oil Shale Technology; CRC Press: Boca Raton, FL, 1991; p 261. (10) Torrente, M. C.; Galan, M. A. Kinetics of the thermal decomposition of oil shale from Puertollano (Spain). Fuel 2001, 80 (3), 327–334. (11) Rajeshwar, K. The kinetics of the thermal decomposition of green river oil shale kerogen by non-isothermal thermogravimetry. Thermochim. Acta 1981, 45 (3), 253–263. (12) Dogan, O. M.; Uysal, B. Z. Non-isothermal pyrolysis kinetics of three Turkish oil shales. Fuel 1996, 75 (12), 1424–1428. (13) Popescu, C.; Segal, E. Critical considerations on the methods for evaluating kinetic parameters from nonisothermal experiments. Int. J. Chem. Kinet. 1998, 30 (5), 313–327.
2841
pubs.acs.org/EF
Energy Fuels 2010, 24, 2841–2847
: DOI:10.1021/ef1001265
Hillier et al.
Equation 2 is in a linear form, so that when the left side of the equation is plotted on the ordinate and 1/T is plotted on the abscissa, the resulting graph will be a straight line when the reaction is first-order. The slope of the straight line yields the activation energy, and the pre-exponential factor can be determined from the intercept. This equation suffers from at least two problems. The first problem is an often numerically imprecise evaluation of the derivative. The second problem is that this equation naturally desensitizes the pre-exponential factor by taking the log of its value. The intercept is found by extrapolation of the data, and small errors in the intercept become large when the pre-exponential factor is determined. The other common method for determining kinetic parameters from non-isothermal data is to use an integral method. There are a few different versions of an integrated form that can be cast into a similar form.14 This method trades the simpler but imprecise derivative method for the solution to eq 1 by integration. However, the resulting integral of the righthand side of eq 1 has no analytical solution, and therefore, some approximations must be made to arrive at the solution. The Coats and Redfern3 derivation uses an approximation of the integral from an infinite sum. The resulting approximation is shown in eq 3 for a first-order reaction after neglecting some terms in the derivation. These assumptions have been shown to diminish the ability of the method to find the true set of kinetic parameters and are a source of the errors inherent in this method.15 The term C here represents the heating rate. ! - lnðmÞ AR 2RT E ¼ ln ln 1ð3Þ T2 CE E RT
Figure 1. (a) Derivative and (b) integral reduction of the same mathematically generated data (with noise).
Because eq 3 is linear, the activation energy can be determined from the slope by plotting the left-hand side against 1/T. The first term on the right-hand side is relatively constant and yields the pre-exponential factor. Other forms of the integral method are some variation of eq 3 with differing assumptions, leading to a variation of activation energies of about 10%.13,14 The Coats and Redfern method is not the only form of the integral method, but it is commonly found in the literature. For example, Popescu et al.13 compiled a table of some 16 techniques that can be classified as differential or integral. Subsequent references in this paper to the integral method will refer to the Coats and Redfern method. Previous “pencil and ruler” linearized methods allowed for only one set of data to be fit at a time, were developed to make up for the lack of computing power, and gave “the solution” without exploring the entire range of possible solutions. These methods also suffer from deciding which points to include in the analysis. The problem of which points to include can be seen in Figure 1. These simulated data were generated for a first-order rate with E = 175 kJ/mol, A = 1.59 1015 s-1, and noise of 0.3% of the normalized mass. The derivative was taken on a point-to-point basis. The natural log of a negative number was set to zero, which shows up as noise. The derivative curve in Figure 1a begins to curve upward at the highest temperature, because of the noise and the low amount of mass remaining. The integral curve in Figure 1a shows a region that is a little jagged when the change in mass is similar
to the noise level, particularly before much mass is reacted. This is the portion of the mass trace where the signal is less than the noise. If there were no noise, then jagged features from the integral method disappear but the curvature at the extremes of the derivative method does not. The curvature is especially pronounced in the derivative data at high temperatures (small values of 1/T) as the point to point derivative transitions from a high rate of change to zero. When fitting such data, one must make a judgment call regarding which points to include, because curvature or jagged features will affect the best fit line from which the kinetic parameters are determined. The traditional derivative method is more sensitive to the actual range of points included in the analysis (i.e., where the peak starts and ends) than the traditional integral method. This means that a traditional analysis may be somewhat arbitrary (as shown subsequently) based on which points are included. This is especially true for data that do not precisely follow the Arrhenius form of the rate constant. This paper proposes a new method to regress the kinetic parameters used in eq 1. The solution with the new method is not linearized; therefore, the pre-exponential factor is not determined by extrapolation. The solution does not neglect any terms, a problem found in the derivation of the integral methods. The problem can be further constrained by assuming that the heating rate has no effect on kinetic parameters (at low heating rates) and, hence, requires that the same kinetic parameters fit similar mass curves generated over a small range of heating rates.
(14) Popescu, C.; Segal, E.; Oprea, C. A comparison of some integral methods. J. Therm. Anal. Calorim. 1992, 38 (4), 929–933. (15) Cai, J.; Bi, L. Precision of the Coats and Redfern method for the determination of the activation energy without neglecting the lowtemperature end of the temperature integral. Energy Fuels 2008, 22 (4), 2172–2174.
Proposed New Method The new method proposed here fits the integral numerically and also requires that the mathematically calculated derivative 2842
Energy Fuels 2010, 24, 2841–2847
: DOI:10.1021/ef1001265
Hillier et al.
also fits the derivative of the data. Fitting the entire mass curve directly with the integrated form avoids both (a) the assumptions made by the analytical solutions and (b) the determination of which points to include in the analysis. This new method essentially weights most heavily on the portion of the curve where the rate of mass loss is the highest. Many of the techniques used in this new method are not unique, but the combination of the techniques is unique. The differential equation (eq 1) is manipulated into the form shown in eq 4 and then analytically integrated to the form shown in eq 5. The integral as shown in eq 5 can easily be evaluated with a numerical integration technique but requires a numerical optimization method to fit this to the given mass curves. The optimization method is required because, unlike the linear methods, these nonlinear methods do not have an easily determined slope and intercept from which the kinetic parameters can be determined. Iteration on the kinetic parameters is required to find the best fit to the data. Z t Z m dm ¼ - Ae- E=RT dt ð4Þ m0 m 0 Z m ¼ m0 expð
t
- Ae- E=RT dtÞ
from an arbitrary initial point chosen by the user to a point that is as least as good as the starting point. The end point is almost always a much better fit based on the sum-squared error. With each cycle, the system is “cooled”, meaning the probability of accepting a worse set of parameters decreases just like in the annealing process, where the molecules of a solid slowly lock into place. When the SA optimization is completed, a generalized reduced gradient algorithm finds the nearest local optimum. If a generalized reduced gradient numerical method was used alone then it would merely find the nearest local optimum without any hope of finding the global optimum. SA gives more hope of being in the region of the global optimum. Mathematically Generated Data Validation of this method was carried out by generating mass loss data from two parallel first-order equations with different activation energies. Randomly distributed noise of up to (0.2% was added to each data point. The noise is important because it is one of the many mechanisms that causes the linearized methods to mispredict the kinetic parameters. At each point, real data are some function of the kinetic parameters and error in the measurement. Equation 6 shows the general situation that the mass at the ith point is some function (represented with a first-order model) plus some error in measurement, ε. The noise levels chosen are similar to levels of noise observed by the authors in real systems. These simulated measurement errors are also the source of the curvature (seen in Figure 1) at the temperature extremes of the derivative curves and the jagged features at the low-temperature region of the integral curves. It is a combination of the noise and deviation from a first-order model that is the source of the difficulty in choosing which points to include in the analysis. Although not the purpose of this work, the authors would also like the reader to be aware of how linearizing affects a statistical treatment of the resulting parameters.
ð5Þ
0
Optimization The determination of a nonlinear set of parameters requires some form of numerical search. The optimization program OptdesX was used here,16 but any program with the necessary optimization routines can be used. OptdesX requires that the fitting and objective functions be provided through a C or FORTRAN subroutine and can be run on a Linux machine. Functions were programmed to fit and minimize a sumsquared error of the differences between the non-isothermal TGA traces and eqs 1 and 4. Both of these equations were derived from the Arrhenius equation for a first-order global reaction. The computational time to fit two replicates of each of the two heating rates (1 and 10 K/min) with their derivatives takes less than 5 min. In the analysis, each integral and derivative data point versus time was scaled to be of equal weight. The optimization method used here is called simulated annealing (SA). SA was so named because it mirrors the annealing process found in nature. A brief overview of SA is included here for convenience; an analogy between annealing and SA as well as further information about SA can be found in the work of Liu.17 In nature, heating a solid to a liquid state and then cooling it slowly to a solid generally allows the atoms to assume a minimum energy state. Likewise, a minimization system is started at a “high energy” state, and random perturbations are made to the parameters. The solution using the new parameters is compared to the current solution. If the new solution has a lower sum-squared error, it is automatically chosen; however, if the sum-squared error is larger, the program probabilistically decides, according to the Boltzmann factor, between moving to the state with a higher sum-squared error and staying in the current state. This is repeated for 1000 cycles until the program brings the system
mi ¼ f ðE, A, mi - 1 , Ti Þ þ εi
ð6Þ
Equation 6 shows the situation if the integral method is used. Once linearized, eq 6 leads to eq 7. Note how the error is now inside of the natural log term, and hence, instead of a normal error, it is now a form of log normal error. lnðmi Þ ¼ lnð f ðE, A, mi - 1 , Ti Þ þ ðεi ÞÞ
ð7Þ
This error means that a traditional statistical analysis based on the central limit theorem is no longer valid. The authors refer any interested readers to the body of knowledge available on the subject. Four validation experiments were performed, with kinetic parameters shown in Table 1. Each validation experiment consisted of four separate mass traces: simulated data were generated at heating rates of 1 and 10 K/min, with a replicate at each heating rate. The replicate data sets had random noise and, hence, were slightly different from each other. Kinetic parameters were hidden from the researchers until after the optimal kinetic parameters for that experiment were determined. A sample of the generated data is shown in Figure 2. The optimization procedure described in this paper was then applied to determine the kinetic parameters. The integral and derivative methods were also applied to each validation experiment. Numerical experiments were conducted on each data set to determine if the number of points included in the integral and derivative methods affected the predicted results. These experiments were conducted by selecting a set number of points to include, centering on the point of maximum mass loss in the first of the two reactions, and then including all points whose value of the derivative was at least 10% of the maximum rate of mass loss. The number of points is reported here as the number of points on each side of the maximum; for example, “10 points” means
(16) Parkinson, A. R.; Balling, R. J. The OptdesX design optimization software. Struct. Multidic. Optim. 2002, 23 (2), 127–139. (17) Liu, L. K. Discrete optimization of pipe networks systems by simulated annealing. Master’s Thesis, Brigham Young University, Provo, UT, 1990.
2843
Energy Fuels 2010, 24, 2841–2847
: DOI:10.1021/ef1001265
Hillier et al.
in Table 1. The integral and derivative methods used only the first of the four mass traces, which had a heating rate of 1 K/min. The results from the point selection experiments are shown in Figures 3-6. Figure 3 is a good example of why the integral method has persisted as a favorite method across the past few decades. Note how the activation energy values using the integral method agree quite well with the actual value. The poorest agreement for integral values of E in Figure 3 occurred when all of the points were included in the analysis. The derivative method obtained lower activation energies for the same experiments, with some compensation in the pre-exponential factor, especially when all of the points were included in the derivative analysis. Figure 4 shows a similar compensation in the pre-exponential factor. A simple calculation will reveal that at the activation energy values in Figure 3 (experiment 1) an increase of 10 kJ/mol will require an order of magnitude increase in the pre-exponential factor to result in the same Arrhenius rate coefficient (k). The comparisons shown in Figures 3-6 all support the concept that the choice of points can alter the predicted kinetic parameters. The derivative method is more sensitive to the point selection problem, although the integral method also varies at least 5 kJ/mol between that predicted for the 10 points on each side case to that predicted for all of the points. These figures also show that the integral method gives a better estimate than the derivative method. Therefore, if a researcher is deciding between the integral and derivative methods, the better choice is the integral method. However, of the three methods compared, the new method proposed in this paper most nearly matches the mathematically generated data in every case.
10 points on each side for a total of 20 points. The amount of “all of the points” varied somewhat between each case but was always above 100. The numbers of points were chosen to consider the sections of curve where the rates of mass loss were approximately 8, 15, 30, and 90% of the maximum mass loss. The actual choices are less important, but this shows that changes in the values of the kinetic parameter will result if an arbitrary decision is made as to which points to include. The pre-exponential factor for the integral method, as shown in eq 3, is dependent upon the temperature and will vary a few percentages depending upon the temperature selected for the analysis. This analysis used the temperature at the highest rate of mass loss to determine the pre-exponential factor in the integral method. The temperature range of 298-1223 K (25-950 °C) was the same in all experiments. This temperature range and the activation energy range for the first reaction are similar to pyrolysis reactions for organic material, such as coal, and the second reaction is similar to mineral decomposition, such as would be found in organic geochemistry reactions.
Results and Discussion The optimization program used all four mass traces to determine coefficients for each of the experiments listed Table 1. Kinetic Parameters Selected for the Validation Experiments first reaction -1
experiment
A (s )
1 2 3 4
1.59 10 1.32 109 2.35 1017 7.20 1010 15
second reaction
E (kJ/mol)
A (s-1)
E (kJ/mol)
175 175 225 150
1.45 1012 1.78 108 4.03 102 7.20 1010
250 219 85 225
Figure 2. Sample data at two heating rates with randomized noise. The kinetic coefficients used to generate data for the different experiments are listed in Table 1.
2844
Energy Fuels 2010, 24, 2841–2847
: DOI:10.1021/ef1001265
Hillier et al.
Figure 3. Experiment 1, kinetic parameter results from fitting of mathematically generated data. Corresponding pre-exponential factors are included for convenience.
Figure 4. Experiment 2, kinetic parameter results from fitting of mathematically generated data. Corresponding pre-exponential factors are included for convenience.
Because E and A are correlated, a comparison was made of the calculated values of k at temperatures representative of each peak. The temperatures selected were 500 K for the first reaction and 600 K for the second reaction. The values of E and A were taken from the curve fit using all of the data points (shown in Table 1). The comparison of k values is presented in Table 2. The error was calculated according to eq 8. error ¼
ðkmodel - kactual Þ kactual
involved generating data with an activation energy and pre-exponential factor initially hidden from the person performing the optimization. Initially, the first reaction of the third experiment initially gave an error of 32.52% using the new method. The source of this error was that the maximum allowable value for the pre-exponential factor in the optimization routine was not set high enough. Therefore, the optimization program was prevented from searching in the area of the correct answer. This problem was rectified by setting the bounds to include many orders of magnitude, using one of the other methods for a first guess. A final error of only 0.73% (shown in Table 2) was obtained when the regression was performed using the enlarged search area.
ð8Þ
It can be seen that the error generated from each type of regression method can change the predicted rate constant significantly, which in turn changes the predicted mass versus time curve. As indicated above, these numerical experiments 2845
Energy Fuels 2010, 24, 2841–2847
: DOI:10.1021/ef1001265
Hillier et al.
Figure 5. Experiment 3, kinetic parameter results from fitting of mathematically generated data. Corresponding pre-exponential factors are included for convenience.
Figure 6. Experiment 4, kinetic parameter results from fitting of mathematically generated data. Corresponding pre-exponential factors are included for convenience.
The new method for obtaining kinetic constants from nonisothermal data works because (1) it does not require extrapolation of the intercept to determine the pre-exponential factor, (2) it does not require a somewhat arbitrary determination of what points will be included in the regression, and (3) it does not require the same assumptions as the linear integral methods to be made in its derivation. The use of data from multiple heating rates also helps fit the kinetic parameters more closely.
A brief set of numerical experiments was performed to determine what pieces of the new method are most important. A set of numerical data was generated for one A-E pair with noise but for the same two heating rates (1 and 10 K/s) used above. The use of two heating rates with the linearized integral method gave relatively poor results, agreeing with the actual activation energy to within (35 kJ/mol. The use of only one heating rate with the new method generated an activation energy that was closer ((5 kJ/mol). The use of the nonlinear 2846
Energy Fuels 2010, 24, 2841–2847
: DOI:10.1021/ef1001265
Hillier et al.
Table 2. Comparison of k Values at (or near) the Temperature of the Maximum Reaction Rate for Each Reaction in a Two Parallel Reaction System first reaction at 500 K
second reaction at 600 K a
errora (%)
error (%) experiment 1
experiment 2
experiment 3
experiment 4
a
actual new method integral derivative method actual new method integral derivative method actual new method integral derivative method actual new method integral derivative method
0.00 -2.67 3.04 5960.29 0.00 -0.50 28.04 20095.95 0.00 -2.70 5.58 4607.10 0.00 -4.12 -95.10 10273.64
experiment 1
experiment 2
experiment 3
experiment 4
actual new method integral derivative method actual new method integral derivative method actual new method integral derivative method actual new method integral derivative method
0.00 -2.15 12.21 23530.43 0.00 -6.84 -65.48 17457.67 0.00 -0.73 -37.80 5986.69 0.00 4.77 19.97 19774.26
Error is calculated according to eq 6.
integral method with two heating rates but without fitting the derivative was even closer ((2 kJ/mol). The use of both the nonlinear integral and derivative methods came within (0.2 kJ/mol for this set of experiments. This new method is not limited to using a first-order equation or to linear heating rates. Any rate equation that describes a system and can be expressed mathematically could be used to determine best-fit parameters.
determination of what points are in the linear region of the graph after the data have been linearized. Fitting the entire mass curve and constraining the kinetic parameters to fit the derivative does not depend upon either extrapolation or a decision of where to start and stop the numerical analysis. A set of numerical experiments was performed on simulated data sets for two heating rates using the new method, a common integral method, and a common derivative method. This new method was shown to outperform the traditional derivative and integral methods in percent error, and between the two traditional methods, the integral method is shown to outperform the derivative method. The result is an improved method for determining the kinetic parameters. This method can in principle be used to fit as many mass traces as desired, but multiple curves at various heating rates are recommended.
Conclusions A new method is proposed to determine kinetic coefficients from non-isothermal data. The method involves using an optimization technique to fit the mass versus time curve as well as the derivative versus time curve. Previous methods have linearized the rate equation to make graphical analysis easier. The limitations of the linearized methods comes from both (a) the extrapolation to the intercept to determine the pre-exponential factor and (b) the somewhat arbitrary
Acknowledgment. The authors thank John Sowa and Stan Moore for demonstrating how to use OptdesX.
2847
