Kinetic analysis of solid-state reactions: Precision of the activation energy calculated by integral methods

1 Kinetic analysis of solid-state reactions: Precision of the activation energy calculated by integral methods L.A. PÉREZ-MAQUEDA*, P.E. SÁNCHEZ-JIMÉNEZ AND J.M. CRIADO Instituto de Ciencia de Materiales de Sevilla. C.S.I.C.-Universidad de Sevilla. Avd. Americo Vespucio 41092 Sevilla. Spain Abstract The integral methods are extensively used for the kinetic analysis of solid-state reactions. The Arrhenius integral function [p(x)] does not have an exact analytical solution. Thus, different approaches, accomplishing the condition that ln g(α) is a linear function of either 1/T or a predetermined function of T, have been proposed for this integral to determine the activation energy from a linear plot of the logarithm of g(α) versus some function of T. The first approach was proposed by Van Krevelen and after that, a number of authors developed new approaches, very often with the scope of increasing the precision of the Arrhenius integral as checked from the standard deviation of the p(xa) function determined from these approximation with regards to the true value of the p(x) function. Besides this method, those proposed by Doyle, Horowitz and Metzger, Coats and Redfern, MacCallum and Tanner and Gyulay and Greenhow are very popular for determining activation energies. In fact, we have found more than 4500 citations (1300 in the last five years) for the papers were these methods were proposed. However, a systematic analysis of the errors involved in the determination of the activation energy from these methods is still missing. A comparative study of the precision of the activation energy as a function of x and T computed from the different integral methods has been carried out. Keywords: Arrhenius integral, integral methods, solid state reactions, errors in activation energy. 2 1. INTRODUCTION Thermally stimulated solid-state reactions, such as decompositions, solid-solid reactions, crystallizations, etc, are, in general, heterogeneous processes. The reaction rate of such processes can be kinetically described, when it takes place under conditions far from equilibrium, by the following expression:1 )()(   fTf dt d (1) where t is the time and α is the extent of reaction ranging from 0 before the process starts to 1 when it is over. Thus, the left hand side term in eq. (1) is the reaction rate. The right hand side term in eq. (1) consists of two terms, i.e. f(T) and f(α), being f(T) a function that describes the dependence of the reaction rate with the temperature (T). Usually, this dependence is described by the Arrhenius equation: RTE eATf / )(   (2), being A the preexponential factor of Arrhenius, E the activation energy and R the gas constant. Additionally, f(α) is a term that describes the dependence of the reaction rate with the mechanism of the process. Different functions have been proposed in literature for describing the kinetic mechanism of the solid-state reactions. These mechanisms are proposed considering different geometrical assumptions for the shape of the material particles (spherical, cylindrical, planar) and driving forces (interface growth, diffusion, nucleation and growth of nuclei). Some of the most common equations proposed for these reactions are included in Table 1. The most common heating profile used for studying solid-state reaction is the linear heating program. Under these experimental conditions, T changes in a wide range of values and a entire α-T curve is recorded in a single experiment. For linear heating rate conditions eq. (1) can be written 3 )( /    fe A dT dRTE  (3) being β the heating rate. Many of the experimental methods used to perform kinetic analysis of solid-state reactions are based in the measurement of the evolution of an integral magnitude, i.e. proportional to the extent of reaction, such as mass loss, released gas, amount of contraction, as a function of temperature. To perform the evaluation of such experimental data, it is necessary either to numerically differentiate the experimental data or to integrate eq. (3):    x x xp R AE dx x e R AE g)()( 2   (4), being x=E/RT. This expression can be written in the logarithmic form: )(lnln))(ln( xp R AE g   (5), Under linear heating rate program, eqs. 4 and 5 do not have an exact analytical solution to p(x) and, therefore, the solution cannot be expressed in a closed form.2 Although, other T-t profiles, such as parabolic or hyperbolic programs, yield to analytical solutions to the Arrhenius integral, they are very seldom used. Thus, several approximated equations have been proposed for p(x) under linear heating program. The approximations for p(x) most commonly used in the determination of the activation energy are those proposed by Coats and Redfern,3,4 Doyle,5-7 Horowitz and Metzger,8 MacCallum and Tanner,9,10 Gyulai and Greenhow,11,12 and Van Krevelen.13 All these approximations have been obtained either by simplifications of the series expressions or in an empirical way. For a given kinetic model, the resulting equations lead to a linear correlation where the kinetic activation energy is easy obtained from the slope. The number of publications where these integral methods are used for 4 determining activation energies is vast. Thus, about 4500 citations can be found in the literature for the original papers3-13 where these equations are proposed. Besides, the popularity of these integral methods has not decreased, as indicated by their more than 1300 citations just in the last years, i.e. 2000-2004. In these last five years, the approach with more citations has been that of Coats and Redfern3,4 with about 590 citations, followed by those of Horowitz and Metzner8 and Doyle5-7 with 230 and 102 citations, respectively (information on the number of citations have been obtained from ISI Web of Science data base). Nevertheless, independently of the approximation used, every g(  ) leads to a high linear correlation coefficient and, therefore, it is not possible to discriminate the kinetic model from a single experimental curve. Additionally, the resulting activation energy values are very much dependant on the g(  ) function assumed for the analysis (these limitations are extended not only to integral methods but also to any procedure that uses a single linear heating rate curve14,15 ). Thus, in principle, the integral methods should be only used under the two following circumstances: (i) when the kinetic model is already known for obtaining the activation energy or (ii) when the activation energy is known for determining the kinetic model. Nevertheless, a new question arises about the precision of the activation energy values determined by these popular integral methods because, as mentioned above, they are based in approximations to the p(x) function and their precision for the estimation of the kinetic parameters are still in doubt, thus some authors have claimed that these methods are imprecise.16-19 Some studies have estimated the errors in the approximated p(x) functions by comparing the resulting values with those calculated by numerical integration, concluding that the errors are quite large. Fig. 1 shows as a way of example the evolution of the relative error of the Coats and Redfern approximation for the estimation of the p(x) function versus x. This figure indicates that the error decreases 5 with x, being significantly large for values of x commonly found in literature for solidstate reactions. These findings have been used as an argument for invalidating these approximated equations in the estimation of the kinetic parameters. Nevertheless, the aim of the aforementioned approximations is the determination of the activation energy and not the accurate computation of p(x). Taking into account that the integral methods are so widely extended and that there is some controversy in their precision, it would be of interest to estimate the precision of such methods for the determination of the activation energy. The aim of the present paper is to perform a comparative study of the precision of the most extensively used approximations to p(x) in the determination of the activation energy. 2. ERRORS IN THE ACTIVATION ENERGY 2.1. Coats and Redfern method. The Coats and Redfern3,4 approach to the Arrhenius integral is the following:            aa x axx e xp a2 1)( (6) the subscript a stands for approximated. In general, the expression more commonly used is the simplified form: 2 )( a x ax e xp a   (7) This approach is named sometimes in literature as Fisher approach.20 By introducing eq. (7) into eq. (4), it follows RTE a aa eT E RA g/ 2 )(     (8) By taking natural logarithms, eq. (8) results 6 RT E E RTA ga a a   2 ln))(ln( (9) Thus, the activation energy could be easily obtained from the slope of the line resulting of plotting ln(g(  ))-2ln(T) versus 1/T. The relative error ε of the activation energy (Ea) calculated by the Coats and Redfern equation can be defined by the following equation: 1001100%                R E R E E EE a a  (10) By differentiating eq. (9): T R E T ga2 /1 )(ln       (11) and by differentiating eq. (5): x xp R E T xp T g        ))(ln( /1 ))(ln( /1 )(ln  (12) Thus, from eqs (11) and (12), it follows xx xp R E R Ea2))(ln(     (13) that substituting in eq. (10) leads to 1001 2))(ln( %           xx xp  (14) This equation indicates that the values of ε% depend on x=(E/RT), and, therefore, on the value of the activation energy and of the range of temperature of the process. The values of ε% have been computed by means of the Mathcad software by numerical integration of the p(x) function using a tolerance (precision in the calculus) of 10-5. The resulting ε% values as a function of the parameter x are included in Table 2. The values 7 included in Table 2 illustrate that there is a significant influence of x in the precision of the calculated activation energy values. Thus, ε% ranges from almost -20% for x=2 to less than -1% for x values larger than 20; in the limit, for x=  , the error is cero. 2.2 Doyle method. The Doyle approach to the Arrhenius integral is the following:5-7 aa xxp 4567.0315.2))(log(    (15) From eq. (15) and eq. (5), it follows RT E R EA ga4567.0315.2log))(log(    (16) Thus, the activation energy can be obtained from the slope of the line resulting from plotting the left hand side of eq. (16) as a function of 1/T: R E T ga 4567.0 /1 ))(log(     (17) The relative error ε% (eq. (10)) for the activation energy obtained by the Doyle method can be obtained from eqs. (12) and (17): 1001 ))(ln( 4567.0 4343.0 %                x xp  (18) The values of ε% have been computed by the same procedure as described in the latter section and the resulting error values are included in Table 2. 2.3. Horowitz and Metzger method. The integer equation after assuming the Horowitz and Metzger approach8 to the p(x) function is the following: 8 2 052.1052.1 33.5ln))(ln( s a s aa RT E TR EA g     (19) where  is a characteristic temperature such that  =T-Ts, being Ts an arbitrary reference temperature. From eq. (19), it is clear that the activation energy is obtained from the slope of the line resulting of plotting the left hand side of eq. (19) versus  , or versus T that yields the same slope: 2 ))(ln( s a RT E T g    (20) From eqs. (12) and (20), the relative error ε% (eq. (10)) in the activation energy obtained by the Horowitz and Metzger8 results: 1001 ))(ln( %        dx xp  (21) Table 2 includes the errors estimated by eq. (21) for the activation energy calculated by the Horowitz and Metzger approach. 8 2.6. Van Krevelen method Considering the Van Krevelen et al approximation13 to the exponential integral of Arrhenius, eq. (5) has the logarithmic form:  T TR E T TE A ga RT E a a a ln )1( 1 368.0 ln))(ln( max max max max                         (22), where Tmax is the temperature at the maximum thermogravimetric rate. The activation energy is determined from the slope of the line resulting from the plot of ln(g(α)) as a function of lnT: )1(ln ))(ln( max     TR E T ga  (23) 9 Thus, the Van Krevelen et al method,13 even though it is a integer equation, for the determination of the activation energy, it requires of the differential experimental curve to obtain the Tmax value for eq. (23). As in the previous sections, the relative error can be calculated from eq. (12) and (23), resulting: 1001 1))(ln( %                xx xp  (24) The resulting values for the error are included in Table 2. 2.4. MacCallum and Tanner method The decimal logarithmic form of the integer equation (eq. (5)) using the approach proposed by MacCallum and Tanner9,10 to the p(x) function results: T E E R EA ga a aa 217449 4828.0log))(log( 4351.0     (25) Thus, the activation energy can be calculated from the slope of the line resulting from the plot of log(g(α)) as a function of 1/T: a E T g217449 )/1( ))(log(     (26) The relative error ε% (eq. (10)) of the activation energy can be calculated from eqs. (12) y (26): 1001 217 449))(ln( · 217 4343.0 %        xRTx xp R  (27) In this case the error depends both on x and T. Table 3 includes the errors in the activation energy as estimated by means of eq. (27). 16 TABLE 1. f(  ) and g(  ) kinetic functions Mechanism Symbol f() g() Phase boundary controlled reaction (contracting area) R2 21 )1(   21 1 12 ()  Phase boundary controlled reaction (contracting volume) R3 32 )1(   31 1 13 ()  Random nucleation followed by an instantaneous growth of nuclei. (Avrami-Erofeev eqn. n =1) 1F1 )1(   )1ln(   Random nucleation and growth of nuclei through different nucleation and nucleus growth models. (Avrami- Erofeev eqn.) An   n n11 )1ln()1(     n/1 )1ln(   Two-dimensional diffusion D2 11   ln( )  ()ln()11    Three-dimensional diffusion (Jander equation) D3         3/1 3/2 112 )1(3    2 3/1 11        Three-dimensional diffusion (Ginstling-Brounshtein equation) D4 3 21 1 13 () /      12 3 1 23   () 1This equation represents an Avrami-Erofeev kinetic model with n=1 instead of a first order reaction. The symbol A1 would be more proper. 17 TABLE 2. Values of the relative error (ε%) for the activation energy calculated by means of the Coats and Redfern, Doyle, Horowitz and Metzger, and Van Krevelen et al. equations as a function of the parameter x (E/RT). x Coats and Redfern Doyle Horowitz and Metzger Van Krevelen 2 -19.72 71.43 80.28 30.28 5 -4.76 28.60 35.26 15.23 10 -1.47 12.72 18.53 8.53 20 -0.42 4.20 9.58 4.58 30 -0.20 1.25 6.47 3.13 50 -0.07 -1.17 3.92 1.92 100 -0.02 -3.02 1.98 0.98  0 -4.90 0 0 TABLE 3. Values of the relative error (ε%) for the activation energy calculated by means of the MacCallum and Tanner equation as a function of the parameter x (E/RT) and the temperature (T). x 400 600 800 1000 1200 2 -48.55 -5.14 16.56 29.59 38.27 5 -15.80 1.56 10.24 15.45 18.92 10 -6.59 2.09 6.43 9.03 10.77 20 -2.59 1.75 3.92 5.22 6.09 30 -1.38 1.51 2.96 3.83 4.40 50 -0.47 1.26 2.13 2.65 3.00 100 0.17 1.04 1.47 1.73 1.91  0.78 0.78 0.78 0.78 0.78 TABLE 4. Values of the relative error (ε%) for the activation energy calculated by means of the Gyulai and Greenhow equation as a function of the parameter x (E/RT) and the temperature (T). x 400 600 800 1000 1200 2 50.25 52.93 54.86 56.37 57.61 5 15.84 17.9 19.39 20.55 21.51 10 4.04 5.89 7.23 8.27 9.14 20 -1.21 0.55 1.82 2.81 3.63 30 -2.43 -0.69 0.56 1.54 2.35 50 -2.72 -0.99 0.26 1.24 2.04 100 -1.70 0.05 1.31 2.30 3.12  18.90 18.90 18.90 18.90 18.90 18 TABLE 5. Values of the activation energies (Ea) and errors (ε%) obtained of the analysis of the simulated curves included in Figs 2a and 2b by means of the different integral methods Simulated curve Fig. 2a (x  5)* Simulated curve Fig. 2b (x  20)* Ea (kJ mol-1)ε% Ea(kJ mol-1) ε% Coats and Redfern 33.5 -4.4 99.6 -0.352 Doyle 44.4 27.0 103.4 3.40 Horowitz and Metzger 48.5 38.7 114.2 14.20 MacCallum and Tanner 38.4 9.7 100.8 0.75 Gyulai and Greenhow 41.7 19.1 100.6 0.63 Van Krevelen et al. 41.9 19.8 106.9 6.93 *The average value of x has been obtained from E/RT  =0.5, where T  =0.5 is the temperature corresponding to  =0.5. TABLE 6. Values of the activation energies (Ea) and errors (ε%) obtained of the analysis of the experimental results for the thermal decomposition of BaCO3 obtained under high vacuum (Fig. 3) by means of the different integral methods. Method Correlation factor Ea  %* Coats and Redfern 0.999 212.5 -0.42 Doyle 0.999 218.6 2.39 MacCallum and Tanner 0.999 222.9 4.45 Gyulai Greenhow 0.999 220.1 3.14 Horowitz Metzger 0.999 230.3 7.92 Van Krevelen et al. 0.999 222.6 4.31 Friedman 0.999 213.4 - * The errors of the activation energies have been calculated using the activation energy determined by the Friedman equation as the accurate value (because no approximation is involved in the method). 19 Figure Captions Fig. 1. Evolution of the relative error of the Coats and Redfern approach for the estimation of the p(x) function versus the value of x. The relative error has been defined by the expression: (pa(x)-p(x)/ p(x))·100, being pa(x) the value obtained by the Coats and Redfern approximation and p(x) the value obtained by numerical integration. Fig. 2. Simulated curves (a) β= 10 K min-1, an A2 kinetic model and the following kinetic parameters: E=35 kJ mol-1 and A=10 min-1; and (b) β= 1 K min-1, an F1 kinetic model and the following kinetic parameters: E=100 kJ mol-1 and A=108 min-1 Fig. 3. Experimental TG curve obtained for the BaCO3 under high vacuum at a heating rate of 0.2 K min-1. 20 400 600 800 1000 0.0 0.2 0.4 0.6 0.8 1.0 b  T / K a Fig. 1 21 800 900 1000 1100 1200 0.0 0.2 0.4 0.6 0.8 1.0  T / K Fig. 2