Global Kinetic Analysis of Complex Materials - Energy & Fuels (ACS

VOLUME 13, NUMBER 1

JANUARY/FEBRUARY 1999

© Copyright 1999 American Chemical Society

Reviews Global Kinetic Analysis of Complex Materials Alan K. Burnham and Robert L. Braun Lawrence Livermore National Laboratory, Livermore, California 94551 Received April 8, 1998

A variety of global kinetic models are reviewed, including first-order, nth-order, nucleation, and sequential models as well as models having Gaussian, Weibull, and discrete activation-energy distributions. The important characteristics of these various models are outlined, with guidance in how to select the correct model. Some of the models have similar characteristics, and the parameter relationships among similar models are discussed. The comparison includes the relationship between conversion-dependent parameters determined by modified Friedman and Coats-Redfern isoconversion methods and reactivity distribution parameters determined by nonlinear regression of rate or fraction-reacted profiles. A new method for deriving discrete activation-energy distribution parameters having ln(A) ) a + bE is also presented. Data accuracy requirements are discussed briefly. Kinetic analyses are given for a variety of materials, including synthetic polymers (polyethylene, polystyrene, polydimethylenenaphthalene, polysulfone, and polyvinyl acetate), petroleum sources rocks (including well-preserved algal kerogens and the Bakken and Monterey shales), oil shales (including kukersite), and the Illinois and Pittsburgh Premium coal samples.

Introduction The derivation of activation energies and frequency factors from chemical reaction data has permeated much of the scientific literature since the original empirical work by Arrhenius in 1889. While the traditional methods for analyzing first- and second-order gaseous and solution reactions are commonly taught in undergraduate courses, these methods are often not applicable to reactions of practical interest involving polymers, minerals, fossil fuels, and biochemical materials. Reactions of practical interest often involve a complex set of sequential and parallel unimolecular and bimolecular reactions that are often impossible to characterize at any significant level of detail. In the case of oil exploration, a further challenge is to reliably extrapolate kinetics derived from experiments taking minutes to weeks in a laboratory to millions of years in nature. The

kinetic analysis must capture the essence of the exceedingly complex reaction set in a tractable mathematical way. Consequently, a variety of other mathematical techniques, often called global kinetic analysis, have been developed to characterize the kinetic behavior of these types of reactions. In our experience, complex materials commonly deviate from first-order kinetic behavior in two principal ways: one involving nucleation kinetic behavior and the other involving distributed reactivity. Nucleation kinetic behavior can be considered random initiation of the reaction at a variety of sites, each of which then expand in space. The reaction rate is proportional to the reaction interface area, which increases during the initial phases of the reaction and then decreases as the reacted regions start to overlap. In contrast, distributed reactivity reactions involve a mixture of labile and refractory

10.1021/ef9800765 CCC: $18.00 © 1999 American Chemical Society Published on Web 12/10/1998

2

Energy & Fuels, Vol. 13, No. 1, 1999

Reviews

reactions, and the reactivity per unit of unreacted material decreases as the reaction progresses. The objective of this paper is to review the variety of kinetic methods that have been developed for analyzing complex reactions. The review attempts to teach how to easily determine which method is most appropriate for the situation at hand. Important aspects include an assessment of how errors in temperature measurement influence the derived kinetic parameters and the prediction of reactivity far from the regime of measurement. Types of Global Kinetic Models This paper treats only some of all conceivable reaction systems and specifically excludes diffusion-controlled reactions. It assumes that the true reaction system is too complex to be characterized in any fundamental way, so the reaction is described in terms of pseudo (or lumped) species, which are themselves complex materials or mixtures. Absolute concentration is not important, as all species are characterized in terms of the fraction of their initial or final value. The basic building block for all reactions is a pseudo-unimolecular reaction

Pseudo-unimolecular: -dx/dt ) ai dyi/dt ) k f(x) (1)

∑

where x is the fraction of the initial material unreacted, f(x) is a mathematical function of the unreacted initial material, yi is the ith product of the reaction, and ∑ai ) 1. The simplest case is that of a pseudo-first-order reaction, for which f(x) ) x. Other more complex functions will be discussed later. The yi values represent, for example, a partitioning into gaseous, liquid, and solid products. The notation used here is not the same as that often used in the kinetic literature. There, R often denotes the fraction reacted, so 1 - x ) R. The present notation is consistent with Braun and Burnham1 and Burnham et al.,2 but Braun et al.3 used x to denote the fraction reacted, so x ) R in that case. However, the present notation is more convenient for describing sequential reactions, as x can be considered as a normalized concentration. The pseudo-unimolecular reactions can be combined in both parallel and serial (sequential) fashion

important extension of the serial reaction is the alternate pathway network

Alternate Pathway:

-dx/dt ) k2 f2(x) + k1 f1(x) (4a) y/dt ) k2 f2(x) - k3 f3(y)

(4b)

dz/dt ) k3 f3(y) + k1 f1(x)

(4c)

The new reaction of the alternate pathway, k1f1, represents the direct conversion of the initial material, x, to the final material, z. Obviously, if k1 equals zero, eq 4 reduces to eq 3. If f2 and f1 are the same functional form, the fraction of material going directly to z relative to that going via the intermediate, y, is the ratio of k1/k2. If k2 is zero, eq 4 reduces to eq 1. The f ’s can take several forms. The most common is a simple first-order reaction, in which case

f(x) ) x

(5)

This case results in the simple solution for constant k in which the reactant undergoes an exponential decay. Another common form is

f(x) ) xn

(6)

commonly known in global kinetics as a pseudo-nthorder reaction, because it does not depend on the absolute concentration of reactant as in ordinary gas or solution kinetics. If n < 1 in a condensed-phase reaction, it is usually considered a shrinking core reaction. A shrinking sphere has n ) 2/3, a shrinking cylinder has n ) 1/2, and a receding plane has n ) 0. In contrast, n > 1 can be interpreted as a γ distribution of reactivity. The final form to be considered is a generalized nucleation-type reaction,4 also called a sigmoidal or autocatalytic reaction, for which

f(x) ) xn(1 - qx)m

(7)

where each kifi can itself be a parallel reaction set. One

where n is still the reaction order, m is a nucleation parameter, and q is an initiation parameter. If m is zero, eq 7 reduces to eq 6. If n, m, and q all equal 1, f(x) reduces to x(1 - x), which is the Prout-Tompkins model.5 The initial rate for such an equation is identically zero, but if and when the reaction does start, the reaction accelerates with time until the consumption of reactant causes the reaction rate to reach a maximum and then decay to zero. The Prout-Tompkins equation was originally chosen empirically because the reaction history of potassium permanganate decomposition had precisely that characteristic. An alternate logarithmic form for a nucleation model, the Avrami-Erofeyev model,6,7 has been shown by series expansion and regrouping to be mathematically equivalent to the generalized nucleation model equation

(1) Braun, R. L.; Burnham, A. K. Energy Fuels 1987, 1, 153-161. (2) Burnham, A. K.; Braun, R. L.; Coburn, T. T.; Sandvik, E. I.; Curry, D. J.; Schmidt, B. J.; Noble, R. A. Energy Fuels 1996, 10, 4959. (3) Braun, R. L.; Burnham, A. K.; Reynolds, J. G.; Clarkson, J. E. Energy Fuels 1991, 5, 192-204.

(4) Sestak, J.; Berggren, G. Thermochim. Acta 1971, 3, 1-12. (5) Prout, E. G.; Tompkins, F. C. Trans Faraday Soc. 1944, 40, 488498. (6) Erofeev, B. V. C. R. Dokl. Akad. Sci. USSR 1946, 52, 511-514. (7) Avrami, M. J. Chem. Phys. 1939, 7, 1103-1112; 1940, 8, 212224; 1941, 9, 177-184.

Parallel: -dx/dt )

∑i∑j aij dyij/dt ) ∑ ajkj f(xj)

(2)

where j represents the jth component of x, ∑ ajxj ) x, yij is the ith product of reaction component j, ∑i∑jaij ) 1, and ∑iaij ) aj.

Serial: -dx/dt ) k2 f2(x)

(3a)

dy/dt ) k2 f2(x) - k3 f3(y)

(3b)

dz/dt ) k3 f3(y)

(3c)

Reviews

Energy & Fuels, Vol. 13, No. 1, 1999 3

with q ) 1 and m and n having certain ordered pairs.8 The justification for the q < 1 comes from two considerations. First, the reaction cannot otherwise start if 1 - x equals one, so q < 1 might represent flaws in the unreacted starting material. Second, one detailed model for polymer decomposition9,10 can be rearranged into a similar function form as eq 7, except there is a nonzero initial reaction rate.2 The derivation of eq 7 is admittedly qualitative rather than rigorous, but it will be seen in the following that it has great utility, and utility is the primary objective of global kinetics. The Arrhenius Equation Most chemical kinetic models assume that the temperature dependence of the rate constant, k, is given by the Arrhenius equation

k ) A exp(-E/RT)

Integrating the Reaction Rate To this point, the reaction history laws have been expressed in differential form. What is normally desired, however, is the amount of any given reactant or product as a function of some thermal history. The rate laws can be integrated either numerically or for simple cases analytically. The general integrated solutions for first- and nthorder reactions, respectively, are

(8) Erofeev, B. V. In Reactivity of Solids, Proceedings of the 4th International Symposium (Amsterdam, 1960); Elsevier: Amsterdam, 1961; pp 273-282. (9) Simha, R. J. Chem. Phys. 1956, 24, 796-802. (10) Ozawa, T. J. Therm. Anal. 1970, 2, 301-324. (11) Glasstone, S.; Laidler, K. J.; Eyring, H. The Theory of Rate Processes; McGraw-Hill: New York, 1941; pp 25-27, 191-193. (12) Benson, S. W. Thermochemical Kinetics; Wiley: New York, 1968; pp 13, 56-66.

(10)

∫0tk(T) dt]1/(1-n)

(11)

and

x ) [1 - (1 - n)

For isothermal conditions, k is constant, so the integral equals kt

x ) exp[-kt]

(12)

x ) [1 - (1 - n)kt]1/(1-n)

(13)

and

(9)

where kB is Boltzmann’s constant and h is Planck’s constant.11,12 Unlike for the Arrhenius equation, however, the preexponential factor depends on temperature. While the appearance of T in the preexponential factor shows that the Arrhenius equation is indeed an empirical equation with only a qualitative theoretical basis, the distinction is not significant from a practical sense. If we assume Es/R ) 25 000 K, the transitionstate rate constants at 427 and 477 °C are 0.004502 and 0.05217 s-1, respectively. If the transition-state rate constants are then fitted to the Arrhenius equation, one obtains A ) 4.1E+13 s-1 and Ea/R ) 25 725 K (in fact, Ea ) Es + RT). For comparison, kT/h is 1.51E+13 s-1 at 452 °C. The temperature dependence of the preexponential factor has been absorbed into the activation energy, with a compensating change in A, so the Arrhenius equation works over the 50 °C interval with an accuracy exceeding one part in 104. More significant is how well the Arrhenius equation works far from the region of the fit. Simple calculations find that the Arrhenius rate constant is low by 24% and 7%, respectively, at temperatures of 100 and 800 °C. Alternatively, one can ask how much difference in temperature would be required to make the rate constants equal at these extremes. The answer is 1.5 °C at 100 °C and 3 °C at 800 °C. These differences are

∫0tk(T) dt]

x ) exp[-

(8)

where A is the preexponential or frequency factor, E is the activation energy, R is the gas constant, and T is the absolute temperature. Although eq 8 is widely used, many fail to recognize that it, also, is an empirical equation with only qualitative justification. Transitionstate theory, which considers the case of reactants and products of a simple reaction in equilibrium with an intermediate activated complex, gives an equation very similar to the Arrhenius equation

k ) (kBT/h)exp(-Es/RT)

unimportant for any practical situation, so the Arrhenius equation can be considered an adequate approximation for any situation in which the true temperature dependence of the preexponential factor is some small power of temperature.

Although the general nucleation model cannot be integrated analytically, the Prout-Tompkins limit can. For isothermal conditions

x ) 1/(e-kt(1 - x0)/x0 + 1) ≈ 1/(e-kt/x0 + 1) (14) where x0 is the initial fraction reacted and the second equality holds if x0 , 1. A value of x0 ) 0.01 is equivalent to using q ) 0.99 in eq 7. Many reactions are conducted with temperature increasing linearly with time, i.e., a constant heating rate dT/dt ) Hr. The integration is more difficult in this case and cannot be done exactly. The most common solution for the kinetic integral, assuming an Arrhenius dependence for k(T), is

∫0tk(T) dt ) (ART2/HrE)(1 - 2RT/E + ...)exp(-E/RT) (15) The 2RT/E term is small relative to 1 and often ignored, but it does make a few degree difference in Tmax and is easy to include. (More accurate expressions than 1 2RT/E have been proposed,13-15 e.g., (1 + 4RT/E)-0.5 and (1 + 2RT/E)-1, but the differences among the various options are 10 times less important than the error from not using any of them as long as E/RT > 20, which is usually true.) Substituting eq 15 into eqs 10 and 11 gives (13) Rogers, F. E.; Ohlemiller, T. J. J. Macromol. Sci. Chem. 1981, 169-185. (14) Popescu, C.; Segal, E. Thermochim. Acta 1984, 253-257. (15) Flynn, J. H. Thermochim. Acta 1997, 300, 83-92.

4


Reviews

∑p∆Up]

xp ) exp[-

x ) exp[-(1 - 2RT/E)(ART2/HrE)exp(-E/RT)] (16)

(22a)

xp ) [1 - (1 - n)Σp∆Up]1/(1-n)

and

x ) [1 - (1 - n)(ART2/HrE)(1 1/(1-n)

2RT/E)exp(-E/RT)]

(17)

Remembering that dx/dt ) -kxn, where n ) 1 for a firstorder reaction, the reaction rate is easily calculated from eqs 16 and 17

dx/dt ) -A exp[-E/RT (1 - 2RT/E)(ART2/HrE)exp(-E/RT)] (18) and

dx/dt ) -A exp(-E/RT)[1 (1 - n)(ART2/HrE)(1 - 2RT/E)exp(-E/RT)]n/(1-n) (19) for n ) 1 and n * 1, respectively. Alternatively, one can use n ) 1.001 in either eq 17 or 19 and obtain results insignificantly different from eqs 16 and 18. Most reactions of practical importance occur at neither isothermal nor constant heating rate conditions. Even for reactions at nominally isothermal conditions, there is usually a finite heat-up time in which some fraction of the material has reacted. Consequently, numerical integration techniques are important not only for complex rate equations that cannot be integrated analytically but also for mathematically simple rate equations undergoing a complex thermal history. One class of numerical integration methods breaks the thermal history into a series of isothermal segments. The simplest of these is an explicit one-step method

xp ) xp-1 - kp-1f(xp-1)(tp - tp-1)

(20)

where xp is the fraction remaining at time p; kp-1 is the reaction rate at the previous time; f(xp-1) is taken from eq 5, 6, or 7; and (tp - tp-1) is the time interval. For a single reaction, this method gives surprisingly small errors. For example, a maximum error of ∼0.005x occurs for a first-order reaction using a time step of 60 s and a heating rate of 1 °C/min. Higher-order numerical methods can attain any level of desired accuracy. For a first-order reaction, another simple method is available. Since the fraction remaining is a simple exponential, it is equal to a product of exponentials.16 For increased accuracy, the average rate constant over the pth time interval in question is used, giving

xp) Πp exp[-(tp - tp-1)(kp + kp-1)/2] ) exp[-

∑p(tp - tp-1)(kp + kp-1)/2]

(21)

This method is several times more accurate than eq 20, with the factor varying with the thermal history. Increased accuracy can be obtained more efficiently than by reducing the isothermal time step by approximating the thermal history as a series of constantheating-rate segments. First, we rewrite eqs 10 and 11 in a generalized form of eq 21 (16) Wood, D. A. Am. Assoc. Pet. Geol. Bull. 1988, 72, 115-134.

(22b)

for first-order reactions and nth-order reactions, respectively, and where

∆Up )

∫tt

p+1

p

k(T)dt

(23)

The usual constant heating rate solution, eq 15, cannot be used in this case because of the nonzero lower integration limit. Instead, we ordinarily use the threeterm rational approximation for the exponential integral17

∆Up )

{

[

]

(E/RTp+1)2 + a1(E/RTp+1) + a2 1 Tp+1kp+1 1 Hr (E/RT )2 + b (E/RT ) + b

[

Tpkp 1 -

p+1

1

p+1

2

]}

(E/RTp) + a1(E/RTp) + a2 (E/RTp)2 + b1(E/RTp) + b2

2

(24)

For a constant heating rate, this equation is more accurate than eq 16. Even greater accuracy (1 part in 107) can be obtained by using the five-term rational approximation given in the same article, but that accuracy is rarely required. Instead, the key to this method is to make the time steps small enough when d2T/dt2 is large that the heating-rate segments faithfully represent the true thermal history. For isothermal segments, one can use either an infinitesimally small heating rate or the analytic isothermal solution. None of these simple integration methods work well for networks of kinetic equations that create a stiff system of differential equations. In this case, advanced numerical integration techniques must be used. We have typically used the numerical integration package LSODE18 for this purpose, and we also use it for integrating the nucleation rate expression (eq 7) and the serial and alternate pathway reaction networks (eq 3 and 4), except for scoping calculations in spreadsheets where the simple one-step explicit or a two-step trapezoidal-rule running sum is sufficiently accurate. Conventional Methods of Kinetic Analysis The ordinary method of analyzing a first-order reaction is to determine a rate constant at each temperature from the slope of either ln(-dx/dt) vs time or ln(x) vs time. The activation energy and frequency factor are then determined from a fit of ln(k) vs 1/T to a straight line. Alternatively, one can assume that each data point represents a measurement of the rate constant and determine A and E directly from a plot of ln(ln(x)) vs 1/T. If either ln(-dx/dt) or ln(x) vs time are not linear, an nth-order reaction is ordinarily considered. If the curve is concave upward, n would be greater than one, and if (17) Gautschi, W.; Cahill, W. F. In Handbook of Mathematical Functions; Abramowitz, M., Stegun, I. A., Eds.; Report AMS 55; National Bureau of Standards: Washington, DC, 1964; pp 228-231. (18) Hindmarsh, A. C. ACM SIGNUM Newsletter 1980, 15 (4), 1011. (19) Friedman, H. L. J. Polym. Sci., Part C 1963, 6, 183-195.

Reviews


the curve is concave downward, n would be less than one. Rearranging eq 13, one can then plot x1-n vs t for various values of n to optimize linearity, obtain (1 n)k from the slope, then obtain A and E from an Arrhenius plot. A more general method of kinetic analysis valid for any thermal history is due to Friedman.19 Combining eqs 1, 6, and 8 gives

-dx/dt ) A exp(-E/RT)xn

(25a) n

ln(-dx/dt) ) -E/RT + ln(Ax )

(25b)

A plot of ln(-dx/dt) at a given fraction reacted versus the 1/T value at which that conversion is reached for several different thermal histories will be linear with a slope equal to -E/R and an intercept of ln(Axn). In his original analysis, Friedman considered only the possibility that there was a single nth-order reaction present. That would imply the same value of -E/R at each level of conversion, and a plot of the intercepts (i.e., ln(Axn)) vs ln(x) would be linear with a slope of n and an intercept of ln(A). In fact, E/R is frequently not the same for all levels of conversion for complex materials, and the plot of ln(Axn) vs ln(x) is rarely linear in our experience, so a different interpretation should be applied to the method. This interpretation will be discussed in the section on distributed reactivity. A second method originates with Coats and Redfern.20 First, eqs 16 and 17 are rearranged into

ln[ln(x)/T 2] ) -E/RT + ln[-AR(1 - 2RT/E)/HrE] (26) and

ln[(1 - x1-n)/((1 - n)T 2)] ) -E/RT + ln[-AR(1 - 2RT/E)/HrE] (27) Data from a single heating rate can be used to derive A and E assuming either a first- or nth-order reaction. Plotting the left-hand side of either equation versus 1/T gives a slope of E/R. In both cases, 1 - 2RT/E is relatively constant over the range of most data. For eq 27, various values of n must be assumed and n is chosen on the basis of which line is the straightest. Aspects of the Friedman and Coats-Redfern analyses can be combined for analyzing data at multiple constant heating rates. First, eq 26 is rearranged into

ln[Hr/(T2(1 - 2RT/E))] ) -E/RT + ln[-AR/(E ln(x))] (28) At fixed conversions for each of the heating rates, the left-hand side for each heating rate is plotted versus 1/T at that heating rate, giving a family of straight liness one for each conversion level. The slope of each line gives E/R, which is substituted into the intercept to obtain A. Since the left-hand side is a weak function of E, the process must occur interatively by first assuming a value of E, then recalculating the left-hand side until convergence occurs. This is a slightly more accurate process than moving (1 - 2RT/E) into the intercept and assuming it is a constant. As with the Friedman (20) Coats, A. W.; Redfern, J. P. Nature 1964, 201, 68-69.

method, E and A are obtained for each level of conversion examined. This method was first used, to our knowledge by Braun et al.,3 although the ln(fraction remaining) term appeared in the ordinate rather than in the intercept and although similar integral isoconversion methods have been reported earlier.21,22 Unfortunately, a programming error in code versions prior to December 1997 replaced the fraction remaining with the fraction reacted in one statement, so the values of A reported in several publications3,23-25 are in error away from midconversion. The interpretation of a variation of A and E with conversion is discussed in the next section. A final conventional method originates with Kissinger.26,27 It uses the temperature of the maximum reaction rate at multiple heating rates to derive A and E. For a first-order reaction, the relationship is derived by first differentiating the rate equation with respect to T and setting it equal to zero, as it is at its maximum value:

d(dx/dt)/dT ) 0 ) -x dk/dT - kdx/dT

(29)

The Arrhenius expression is substituted into the temperature derivative of k, and the right-hand temperature derivative is expanded into (d/dx)(dx/dT), giving

0 ) -kx(E/RTmax2 - Ae-E/RT/Hr)

(30)

This can only be true when the expression inside the right-hand parentheses is zero. Setting it so and taking natural logarithms yields

ln(Hr/RTmax2) ) - E/RTmax + ln(A/E)

(31)

Equation 31 is very similar to eqs 25b and 28 in that the Tmax-shift method (what we call eq 31) is really a special case of determining A and E at a fixed conversion. Equation 31 is exact for both first-order26 and nthorder reactions27 and is an excellent approximation for nucleation28 and distributed reactivity models.1 Distributed Reactivity Models A common approximation for the decomposition kinetics of complex materials is the parallel reaction model. This model assumes that the distribution of reactivity caused by the reaction complexity can be represented by a set of independent, parallel reactions, each with their own frequency factor and activation energy, i.e., eq 2. Usually, it is further assumed that all reactions share the same frequency factor, so the reactivity distribution is represented by a continuous distribution of activation energies. The validity of this approximation will be investigated in a later section, but for the moment, the basic integrated kinetic equa(21) Ozawa, T. Bull. Chem. Soc. Jpn. 1965, 38, 1881-1886. (22) Flynn, J. H. J. Therm. Anal. 1983, 27, 95-102. (23) Reynolds, J. G.; Burnham, A. K. Org. Geochem. 1995, 23, 1119. (24) Reynolds, J. G.; Burnham, A. K.; Mitchell, T. O. Org. Geochem. 1995, 23, 109-120. (25) Reynolds, J. G.; Burnham, A. K. Energy Fuels 1997, 11, 8897. (26) Kissinger, H. E. J. Res. Natl. Bur. Stand. 1956, 57, 217-221. (27) Kissinger, H. E. Anal. Chem. 1957, 29, 1702-1706. (28) Chen, D.; Gao, X.; Dollimore, D. Thermochim. Acta 1993, 215, 109-117.

6


Reviews

tion for this simplified case is

x)

∫0∞exp[-∫0tk(T) dt]D(E)dE

(32)

where D(E) is the distribution of activation energies and

∫0

∞

D(E)dE ) 1

(33)

Specific mathematical forms appearing in the literature are the Gaussian,1,29-32 Weibull,33 and Gamma distributions.34-40 The distribution can also be a finite discrete distribution of arbitrary form, in which case the integral in eq 32 would be replaced with a summation.41-44 In practice, there is less difference between continuous and discrete distributions than one might first imagine, because the most efficient way of computing reaction histories of reactions having continuous distributions is to approximate the continuous distribution by a discrete distribution. To our knowledge, the first application of an activation-energy distribution kinetic model is due to Vand.45 Vand considered the case of irreversible electrical resistance changes upon annealing of evaporated metallic films at both constant temperature and for a constant heating rate. Due to the lack of digital computing at that time, Vand needed to make certain mathematical approximations that are unnecessary today. First, he characterized individual reactions by a time related to the temperature (or heating rate) and their activation energy and frequency factor. He next assumed that, at a given time, only those reactions with characteristic times longer than the given time contribute to the reaction. An approximate energy distribution could then be calculated from the instantaneous reaction rate. Pitt46 used the same technique to calculate approximate activation-energy distributions for coal pyrolysis. Hashimoto et al.47 and Miura48 combined aspects of the Friedman approach with an approximation similar to Vand’s in order to derive approximate activation-energy distributions for a variety of coals. (29) Anthony, D. B.; Howard, J. B. AICHE J. 1976, 22, 625-656. (30) Campbell, J. H.; Gallegos, G.; Gregg, M. L. Fuel 1980, 59, 727732. (31) Lewellen, P. C.; Peters, W. A.; Howard, J. B. Symp (Int.) Combust., [Proc.] 1977, 16, 1471-1480. (32) LePoutre, Etude de Modelisation de la Genese des Hydrocarbures dans les Bassins Sedimentaires; Report ENSPM-IFP, 1981. (33) Lakshmanan, C. C.; White, N. Energy Fuels 1994, 8, 11581167. (34) Boudreau, B. P.; Ruddick, B. R. Am. J. Sci. 1991, 291, 507538. (35) Taurutis, W. J., Jr. Geochim. Cosmochim. Acta 1993, 13491350. (36) Ho, T. C.; Aris, R. AICHE J. 1987, 33, 1050-1051. (37) Aris, R. AICHE J. 1989, 35, 539-548. (38) Crickmore, P. J. Can. J. Chem. Eng. 1989, 67, 392-396. (39) Astarita, G. AICHE J. 1989, 35, 529-532. (40) Kemp, R. R. D.; Wojciechowski, B. W. Ind. Eng. Chem. Fundam. 1974, 13, 332-336. (41) Tissot, B. P.; Espitalie, J. Rev. Inst. Pet. 1975, 30, 743-777. (42) Ungerer, P. In Thermal Phenomena in Sedimentary Basins; Durand, B., Ed.; Technip: Paris, 1986; pp 235-246. (43) Ungerer, P.; Pelet, R. Nature 1987, 327, 52-54. (44) Burnham, A. K.; Braun, R. L.; Gregg, H. R.; Samoun, A. M. Energy Fuels 1987, 1, 452-458. (45) Vand, V. Proc. Phys. Soc. (London) 1943, 222-246. (46) Pitt, G. J. Fuel 1962, 41, 267-274. (47) Hashimoto, K.; Miura, K.; Watanabe, T. AICHE J. 1982, 28, 737-746. (48) Miura, K. Energy Fuels 1995, 9, 302-307.

Anthony and Howard29 popularized the Gaussian distribution of activation energies for characterizing coal pyrolysis kinetics. In this case,

D(E) ) (2π)-1/2σ-1exp[-(E - E0)2/2σ2]

(34)

where σ is the distribution parameter, is substituted into eq 32. Increasing the Gaussian distribution parameter, σ, makes the reaction profile at a constant heating rate broader and more symmetric. The Gaussian distribution was applied shortly thereafter to oil shale pyrolysis,30 biomass pyrolysis,31 and petroleum formation (Le Poutre,32 as quoted by Ungerer42). Some early coal researchers (e.g., Wiser49 and Sklyar et al.50) had success with pseudo-nth-order rate laws. As for Gaussian distributions, nth-order reactions have a long time tail for isothermal conditions. For a constant heating rate, increasing the reaction order also makes the reaction profile broader, but unlike the Gaussian distribution, increased n makes the reaction profile more skewed to high temperature. We have routinely analyzed coal and petroleum source rock data with a pseudo-nth-order model as well as other distribution models and have found similar mean activation energies for the nth-order and distribution models.3 There was a strong correlation between reaction order and the Gaussian σ, suggesting that the two parameters accomplished much the same purpose. Until recently, we thought that activation-energy distribution models are more rigorously correct, because the distribution models clearly relate to observed variations in product composition and because the shape of the reaction profile for petroleum source rocks was more closely matched by the Gaussian distribution.3 However, researchers considering the early diagenesis of sedimentary organic matter34,35 and others considering general chemical engineering36-40 have shown that a Gamma distribution of reactivity is identical to both a power law (in time) reaction and a pseudo-nth-order reaction. The Gamma distribution was first proposed as a method to describe chemical reactivity distributions at constant temperature by Aris.51 In this case

x)

∫0∞D(k)e-kt dk

(35)

which represents an integration of first-order exponential decays over the reactivity distribution. For the Gamma distribution, D(k) is given by

D(k) ) aνkν-1e-ak/Γ(ν)

(36)

where ν is a distribution shape parameter, a is a measure of the average lifetime of the mixture, and Γ() is the Gamma function. If ν is greater than one, there is no contribution in the distribution at k ) 0 and the distribution peaks at some finite value of k. If 0 < ν < 1, the distribution is weakly singular at k ) 0 but exponential-like for large k. Integrating this equation over k yields (49) Wiser, W. H.; Hill, G. R.; Kertamus, N. J. Ind. Eng. Chem. Process Des. Devel. 1967, 6, 133-138. (50) Sklyar, M. G.; Shustikov, V. I.; Virozub, I. V. Int. Chem. Eng. 1969, 9, 595-602. (51) Aris, R. Arch. Ration. Mech. Anal. 1968, 27, 356-64.

Reviews


x ) aν/(a + t)ν

(37)

Differentiating with respect to time gives the reaction rate

dx/dt ) νaν/(a + t)ν+1 ) (ν/a)x1+1/ν

(38)

Now it is apparent that ν/a is the effective rate constant and 1 + 1/ν is the apparent reaction order. Consequently, it is also apparent that fitting to an nth-order rate law is equivalent to fitting to a Gamma distribution of reactivity. What is less clear from this discussion is how to include temperature dependence with the Gamma distribution approach. It will be shown in the next section that the most appropriate way is with a distribution of frequency factors, but a distribution of activation energies as in eq 32 is an excellent approximation. Another distribution form used in the literature is the Weibull distribution30,52

D(E) ) (β/η)[(E - γ)/η]β-1 exp{-[(E - γ)/η]β}

(39)

in which η is a width parameter, β is a shape parameter, and γ is the activation-energy threshold. The mean activation energy, E0, is also defined for later use

E0 ) γ + ηΓ(1/β + 1)

(40)

where Γ( ) again is the Gamma function. The Weibull distribution was originally derived to model the distribution of mechanical failure strengths.52 It was first used, to our knowledge, as a distributed chemical reaction law by Dorko53 (as cited by Brown et al.54). Burnham et al.55 suggested it as an alternative to a combination of a Gaussian distribution with an nthorder reaction to account for the width and asymmetry of coal pyrolysis profiles. Subsequently, Lakshmanan et al.33 explored the use of the Weibull distribution in much greater detail for petroleum source rocks. The final distribution form to be considered is the discrete distribution. In its most general form, there is no relationship between A and E for the various discrete reactions. This can lead to non-Arrhenius behavior as the various reactions change their relative reactivity ordering with large temperature changes.56 More useful are simplifications in which (i) ln(A) is a linear function of E,42 (ii) A has the same value for all reactions and the reactivity distribution is reflected in an activationenergy distribution,43,44,57 and (iii) E has a single value and the distribution in reactivity is reflected in a distribution of frequency factors.58 The second simplification is most common and is used by petroleum exploration geochemists to analyze hundreds of petroleum source rocks annually. The keys to its widespread, (52) Weibull, W. J. J. Appl. Mech. 1951, 18, 293-296. (53) Dorko, E. A.; Bryant, W.; Regulinski, T. L. Anal. Calorim. 1974, 3, 505. (54) Brown, W. E.; Dollimore, D.; Galwey, A. K. Reactions in the Solid State. Comprehensive Chemical Kinetics; Bamford, C. H., Tipper, C. F., Eds.; Elsevier: Amsterdam, 1980; Vol. 22, pp 55, 83. (55) Burnham, A. K.; Oh, M. S.; Crawford, R. W.; Samoun, A. M. Energy Fuels 1989, 3, 42-55. (56) Golikeri, S. V.; Luss, D. AICHE J. 1972, 18, 277-282. (57) Hanbaba, P. Reaktionskinetische Untersuchungen zur Kohlenwasserstoffentbindung aus Steinkohlen bie niedrigen Aufheizgeschwindigkeiten. Dissertation, University of Aachen, 1967. (58) Forbes, P. L.; Ungerer, P. M.; Kuhfuss, A. B.; Riis, F.; Eggen, S. AAPG Bull. 1991, 75, 873-893.

routine usage are its flexibility to fit almost any reactivity distribution and an efficient nested nonlinearconstrained linear regression routine first used by Burnham et al.44 and later described in more detail by Sundararaman et al.59 Comparison and Selection of the Various Kinetic Models The key to selecting the correct kinetic model is to understand the characteristics of the various models. The objective of this section is to graphically compare the models. This comparison is done within the context of KINETICS for Windows 95/98/NT, our kinetics analysis program with a graphic user interface.60 Many times at least part of the kinetic data to be analyzed has either a constant temperature or a constant heating rate. If so, visual inspection of that data gives the first indication of what model to use. First, however, it is necessary to determine whether the measurement device reports x, 1 - x, or dx/dt. A thermogravimetric analyzer (TGA) measures x (fraction remaining), an integrating product collection device such as a graduated test tube measures 1 - x, and an evolved gas analyzer ordinarily measures dx/dt, which is proportional to product concentration in a flowing gas. Having one, of course, allows the others to be obtained by simple mathematical transformation. One warning is that a TGA often includes an unreactive residue that should not be included as unreacted material. This can be treated either by including a constant term in the analysis or, preferably, by differentiating the data (converting it to reaction rates) so that any constant term drops out. For isothermal data, two characteristics determine the appropriate model. First, is a plot of ln(x) or ln(dx/dt) versus time linear, concave upward, or concave downward? Second, does the reaction rate have its maximum value at initial time? The first plot determines whether the reaction is first-order or whether a more complicated model is required. If the plot is concave upward, the reaction model could be either an nth-order reaction with n > 1 or a multiexponential decay having anywhere from two to a continuous distribution of rate constants. If the plot is concave downward, either an nth-order model with n < 1, a serial model, or a nucleation model is required. If the initial rate does not have its maximum value at time zero, one must next determine whether that characteristic is due to the intrinsic chemical kinetics or to measurement artifacts. If the former, it can be used to distinguish an nth-order reaction with n < 1 from a serial or nucleation reaction. On the other hand, the initial rise could be due either to a finite heatup time of the sample in the apparatus or to dispersion during product transit to the detector. If due to a finite heatup time, one must either discard the early part of the reaction data or use a time-dependent temperature and numerically integrate the calculated rate over the exact thermal history. KINETICS is capable of the latter. If due to a dispersion effect, one must either deconvolve (59) Sundararaman, P.; Merz, P. A.; Mann, R. G. Energy Fuels 1992, 6, 793-803. (60) Available for U.S. Government work at http://www.doe.gov/ html/osti/estsc/estsc.html and for commercial purposes at http:// www.humble-inc.com/indexhis.htm.

8


Figure 1. Comparison of three types of kinetic models, all having A ) 1.0E+13 s-1 and E ) 51 kcal/mol, at a contant temperature of 500 °C.

the true data from the transfer function prior to kinetic analysis or convolve the transfer function with the calculated reaction data prior to comparison to the measured data during nonlinear regression. In practice, we have found the latter course to be more stable, so KINETICS takes that approach. The transfer function can be either a broadened delta function or a broadened step function. Figure 1 shows a selection of calculated isothermal reaction rate curves on a semilog plot for several kinetic models. Limited space prevents showing all possible parameter ranges, but it should be evident from those shown that while it is easily possible to reduce the number of plausible kinetic models to a few, it is usually not possible to select only one, in part because nature hardly ever follows any simple rate law exactly. Particular difficulty arises in distinguishing between nthorder models with n > 1 and distributed activation energy models. Sometimes the distinction between two potential kinetic models can be enhanced by examining an alternative heating schedule. A constant heating rate is very commonly used, in part because it can easily sample the entire reaction with equivalent detail in a single experiment. This is not to say that a single experiment is capable of deriving reliable kinetic parameters, but it is very useful as a means to select a kinetic model for further investigation. Further, reac-

Reviews

Figure 2. Comparison of three types of kinetic models, all having A ) 1.0E+13 s-1 and E ) 51 kcal/mol, at a heating rate of 10 °C/min.

tion-rate data is usually more sensitive to details of the reaction mechanism, since any inflection due to partial overlap of separate reactions is amplified by the derivative. Figure 2 shows a selection of calculated reactionrate data for the same kinetic models shown in Figure 1. (In this case, neither nonlinear sample heatup due to, for example, the heat load from an endothermic reaction nor dispersion from the reactor transfer function are considered.) The reaction profiles differ in both width and shape. It is important to realize that profile width alone for a single heating rate is not a good way to distinguish kinetic models. If A and E are changed in concert to maintain a constant Tmax, the reaction profile width is inversely proportional to A and E. Consequently, if a single reaction profile consisting of overlapping parallel reactions is fitted to a single first-order reaction, the derived A and E will be smaller than any in the reaction set. This effect was first demonstrated for a discrete distribution by Hanbaba57 and studied in more detail for reactions having a Gaussian energy distribution by Braun and Burnham.1 For example, fitting a first-order reaction to a reaction profile having a Gaussian energy distribution of 4% of the mean value, within the range for natural materials, will result in an activation energy only one-half the true mean value. Likewise, the profile shape of a nucleation reaction having n ) 1 in eq 7 is independent of m for 0 < m 1. This is shown more clearly in Figure 10 for polystyrene, which clearly has a shape similar to a second-order reaction. Consequently, nonlinear regression analysis allowing both n and m to float in eq 7 returns reaction orders near 2 for PS, PSF, and PVAC-1. That m is less than 1 for these three reactions is important. An nth-order reaction alone is too broad to fit both the reaction profile and the shift in Tmax with heating rate. Second, this particular polystyrene sample has been pyrolyzed at constant temperature3 and found to have an acceleratory phase at the beginning of the reaction. This characteristic is incompatible with a simple nth-order reaction. A slightly better fit (12% reduction in residuals) is obtained also for PDMN if all both m and n are allowed to float. The derived value of n ) 1.2 is only slightly greater than 1. The high-temperature peak of PVAC is different from the others in that it is broader than a first-order reaction. It is also more symmetric, suggesting a reaction order near 1.5. That is precisely what is obtained by nonlinear regression, as shown in Table 4. Comparison of the measured and calculated reaction curves (Figure 11, bottom) indicates the model fits well, even though slight changes in profile shape may be occurring with heating rate. It is also noteworthy that not all eight regression parameters (2 A’s, 2 E’s, 2 n’s, 1 m, and relative fraction) were varied simultaneously. Instead, four successive regression analyses were conducted with 3-6 parameters floating at a time, leading to a rapid overall convergence. These polymers give a good test of the significance of the modified Friedman and Coats-Redfern analysis. Only the modified Coats-Redfern method is presented in Table 5, due to its superior match of Tmax in Figure 8, but similar qualitative conclusions can be derived from the modified Friedman method. PE and PSF have quite uniform A’s and E’s over the reaction profile, which agree well with the values in Table 2. The parameters for PS and PDMN tend to drift with conversion, in opposite directions, but it is difficult to say if the drift is significant. The values at midconversion agree well with the results in Table 4. PVAC provides the most interesting result. The values at 0.5

Reviews


Table 3. Kinetic Parameters Derived Using the Relationship ln(A) ) a + bE from Two Model Reactions: A Discrete Distribution Having ln(A) ) 13.77 + (3.466E-04)E and a Third-Order Reaction a b

discrete distribution

third-order reaction

13.63 3.491E-04

51.02 -4.113E-04

fraction

A (s-1)

E (cal/mol)

fraction

A (s-1)

E(cal/mol)

0.1 1.2 2.1 5.3 9.9 15.5 21.1 15.8 9.9 7.1 5.0 3.4 1.9 1.5 0.1

5.481E+11 1.073E+12 2.103E+12 4.119E+12 8.067E+12 1.58E+13 3.096E+13 6.062E+13 1.188E+14 2.326E+14 4.556E+14 8.925E+14 1.748E+15 3.424E+15 6.705E+15

38 392 40 317 42 243 44 169 46 094 48 020 49 946 51 871 53 797 55 722 57 648 59 574 61 499 63 425 65 350

1.7 19.5 27.3 19.9 12.6 7.4 4.5 2.7 1.7 1.0 0.7 0.4 0.3 0.2 0.2 0.2

2.338E+13 1.549E+13 1.027E+13 6.807E+12 4.512E+12 2.991E+12 1.982E+12 1.314E+12 8.708E+11 5.772E+11 3.826E+11 2.536E+11 1.681E+11 1.114E+11 4.894E+10 2.15E+10

49 203 50 203 51 203 52 203 53 203 54 203 55 203 56 203 57 203 58 203 59 203 60 203 61 203 62 203 64 203 66 203

Table 4. Approximate and Nonlinear Regression Analysis of Kinetic Data for Several Simple Synthetic Polymers property

PE

PDMN

PS

first-order A, s-1 E, cal/mol std error, cal/mol avg exp/calc width avg asymmetry n from width n from asymmetry m from width m from width and asym A corrected for m, s-1

1.80E+16 62 722 2719 0.66 0.60 0.44 0.95 0.42 0.40 2.86E+16

Approximate Kinetic Analysis 3.25E+16 7.90E+14 61 843 53 683 1172 2120 0.61 0.71 0.78 1.10 0.39 0.51 1.17 1.53 0.46 0.36 0.54 0.55 6.37E+16 1.57E+15

A, s-1 E, cal/mol n m

2.25E+16 62 224 (1.00) 0.42

Nonlinear Regression Analysis 3.22E+16 1.25E+16 60 647 56 154 1.21 1.81 0.59 0.63

Figure 10. Normalized reaction-rate curve used to compare profile shapes with nth-order reactions.

and 0.6 conversion, during the transition between the two peaks, are discarded because of standard deviations greater than 10 kcal/mol. The values at 0.3 and 0.8 fraction reacted agree well with those in Table 4 from the other analysis methods. Two qualifications are important. First, the exact kinetic behavior of a given polymer depends on its molecular weight distribution and degree of branching. The reaction profile of a high-density polyethylene

PSF

PVAC-1

PVAC-2

1.26E+15 61 740 7270 0.65 1.10 0.43 1.53 0.43 0.66 3.02E+15

1.08E+11 35 787 773 0.62 1.26 0.40 1.70 0.46 0.78 3.20E+11

2.03E13 50 433 2666 1.21 0.96 1.45 1.37 (0.00) (0.00) 2.03E+13

3.49E+15 61 367 2.19 0.78

1.44E+12 37 526 1.70 0.74

1.14E+15 56 139 1.66 (0.00)

examined in our laboratory by thermogravimetric analysis at 5 °C/min followed the 0.5 reaction order curve fairly closely. Also, the temperature at its maximum decomposition rate was about 10 °C higher than for the low-density film. Second, these models are not complete mechanistic representations of the reaction rate. For example, the shape of the polystyrene decomposition reaction profile changes heating rate. This may suggest the need for a more complicated model having two pathways whose relative importance changes with temperature, but such an analysis is beyond the scope of this article. Cellulose and Biomass. The kinetics of cellulose, biomass, and wood products are important for energy, fire safety, and other considerations. Two materials of particular interest to us are cellulose itself and waste paper characteristic of municipal solid waste. Reynolds and Burnham25 found cellulose and paper dunnage to be characteristically different in their kinetic properties: cellulose has a very narrow reaction profile characteristic of a nucleation reaction, and paper has a multicomponent reaction profile with at least three contributing reactions. They used the modified ProutTompkins model (three-parameter nucleation model) for cellulose and the discrete activation energy model for paper, both as-manufactured and after various levels of hydrothermal pretreatment. The principal paper

16 Energy & Fuels, Vol. 13, No. 1, 1999

Figure 11. Observed and calculated reactions rates for PDMN, PSF, and PVAC. For clarity, only the highest and lowest heating rate is shown for PSF. The lines are the fits.

pyrolysis peak is at about the same temperature as for cellulose. Alternatively, Chang et al.66 (1996) used parallel first- and nth-order reactions to fit similar reaction profiles for waste paper. (66) Chang, C. Y.; Wu, C. H.; Hwang, J. Y.; Lin, J. P.; Yang, W. F.; Shih, S. M.; Chen, L. W.; Chang, F. W. J. Environ. Eng. ASCE 1996, 299-305.

Reviews

Antal et al.67 appear to challenge the conclusion of Reynolds and Burnham25 that a nucleation (or autocatalytic) reaction is required for cellulose by reasserting, in essence, that a first-order reaction fits all primary-stage pyrolysis data as long as the inability of the same kinetic parameters to fit all heating rates is interpreted as a temperature measurement error that can be absorbed into the frequency factor. While it is true that a uniform temperature shift in the data will be reflected primarily in the frequency factor,44 the ability to fit a single reaction profile does not prove that a first-order model is correct. LLNL data from both TGA and Pyromat also have the same shape as a first-order reaction when plotted as for PE in Figure 10, but that correspondence does not address the value of the nucleation parameter, m. Deriving 50% conversion temperatures of 317.65 and 338.25 °C at 2 and 10 °C/min from their kinetic parameters, we obtain E ) 53 683 cal/mol from eq 31 compared to 57 361 cal/mol they obtain from fitting the reaction profile. Given that a 1 °C error over their 5-fold heating rate range causes a 2500 cal/mol error, these two reaction profiles cannot realistically constrain the activation energy to any better than 49-58 kcal/mol. The low end of the range is more consistent with the literature as a whole and supports our notion that the profile is narrower than a simple first-order reaction. Aggarwal et al.68 recently reported an activation energy of ca. 58 kcal/mol for cellulose from fitting a firstorder reaction to a single heating rate, despite an earlier paper69concluding that the Avrami-Erofeyev model (similar to the Prout-Tompkins model) is most appropriate. More seriously, they report an activation energy of 113 kcal/mol (greater than virtually all chemical bonds) for corn starch, which decomposes at a slightly lower temperature than cellulose with a sharper nonisothermal profile. We consider these values as additional support that a first-order fit to a single heating rate cannot yield a reliable activation energy because of the insensitivity of profile shape to deviations from first-order character for nucleation and autocatalytic reactions. Additional support for a nucleation reaction comes from both groups. Isothermal data of Ainscough et al.70 appears to show acceleratory behavior in the early stages, even though the later stages are described well by a first-order reaction. Also, analysis of low-temperature isothermal data of Varhegyi et al.71 yields kinetic parameters similar to our three-parameter nucleation analysis of high-temperature data. First, conversions were calculated for the first six experiments in their Table 1 from the fraction of weight remaining by assuming 100% conversion at the weight remaining (typically 3-9%) after 30 min at 500 °C. (Consistent with their model, the residual reactivity is more prop(67) Antal, M. J., Jr.; Varhegyi, G. Energy Fuels 1997, 11, 13091310. (68) Aggarwal, P.; Dollimore, D.; Heon, K. J. Therm. Anal. 1997, 50, 7-17. (69) Dollimore, D.; Holt, B. J. Polym. Sci., Polym. Phys. 1973, 11, 1703-1711. (70) Ainscough, A. N.; Dollimore, D.; Holt, B.; Kirkham, W.; Martin, D. In Reactivity of Solids, Proceedings of the 7th Intl Symp.; Chapman and Hall: New York, 1972; pp 543-552. (71) Varhegyi, G.; Jakob, E.; Antal, M. J., Jr. Energy Fuels 1994, 8, 1345-1352.

Reviews

Energy & Fuels, Vol. 13, No. 1, 1999 17 Table 5. Modified Coats-Redfern Analysis of Five Synthetic Polymersa PE

PDMN

PS

PSF

PVAC

fraction reacted

A

E

A

E

A

E

A

E

A

E

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

2.80E+14 7.72E+14 1.39E+15 2.83E+15 4.74E+15 6.84E+15 6.65E+15 6.90E+15 5.85E+15

57 807 58 886 59 501 60 347 60 934 61 321 61 153 61 062 60 635

1.28E+17 1.03E+17 5.11E+16 3.26E+16 2.15E+16 9.10E+15 5.02E+15 2.20E+15 5.70E+14

65 217 64 437 63 180 62 337 61 579 60 230 59 244 57 942 55 909

8.49E+14 1.20E+15 2.11E+15 4.18E+15 7.02E+15 2.48E+16 9.42E+16 5.52E+17 4.73E+18

55 261 55 264 55 764 56 492 57 061 58 680 60 452 62 890 66 055

2.94E+13 1.31E+14 1.64E+14 1.29E+14 1.34E+14 1.46E+14 1.24E+14 2.38E+13 5.80E+14

57 959 59 620 59 633 59 051 58 934 58 940 58 626 56 158 61 970

5.27E+11 3.40E+11 4.35E+11 1.28E+12

39 528 38 656 38 852 40 278

7.13E+17 3.11E+16 9.95E+15

64 129 60 490 59 615

a

A is in s-1; E is in cal/mol.

erly considered as secondary pyrolysis of a retrograde char rather than cellulose pyrolysis.) Fitting these experiments simultaneously to a nucleation reaction using the exact thermal history and conversions up to 370 °C yielded Arrhenius parameters of A ) 2.57E+14 s-1, E ) 46 969 cal/mol, and a nucleation parameter, m, of 0.15. The latter result shows that the activation parameters reported by Varhegyi et al.70 do not extrapolate to their low-temperature experiments as they claim.66 Further, it reinforces our conclusion that the mean activation energy is in the mid-to-high 40 kcal/mol range and that the narrow nonisothermal profile causing the higher apparent first-order activation energies is due to some acceleratory process that is still not yet adequately understood. Type I and II Kerogens. Type I and II kerogens are the organic matter in oil shales and the source of most of the world’s crude oil. Understanding their kinetics is important both for oil shale processing and petroleum exploration. Type I and II kerogens have somewhat elusive definitions, sometimes given in terms of their elemental composition and position on a Van Krevelan diagram, sometimes given in terms of their pyrolysis yield, and sometimes given in terms of the composition of the pyrolysis products. There is a rough correlation between type I kerogens and algal kerogens formed in lacustrine environments and between type II kerogens and more diverse kerogens formed in a marine environment. Regardless, the variation in chemical composition of these kerogens is reflected in a corresponding variation in their chemical kinetics. Although small in relative abundance, some wellpreserved kerogens required a nucleation kinetic model.2 Again, the characteristic of a reaction profile that is narrower than that calculated from first-order parameters derived from eq 31 correlated with the presence of an acceleratory phase during isothermal pyrolysis. The five such samples we have investigated have the following ordered pairs of relative width and asymmetry: Uinta Basin core sample Brotherson 1-23B4 (0.73, 0.52); Ravnefjeld marine shale sample GGU24 (0.67, 0.70); Indonesian lacustrine source rock sample OLS-1 (0.76, 0.58); and boghead coal samples Frejus (0.72, 0.56) and Autun (0.91, 0.59). The Brotherson sample appears to have slightly greater asymmetry (deviation from 1.0), while the GGU24 sample is more symmetric. A plot of normalized rate versus conversion is shown for three of the samples in Figure 12. The Brotherson and Frejus samples follow the first-order curve fairly well, while the GGU24 sample migrates

Figure 12. Normalized reaction profile shapes for three kerogens having narrow reaction profiles. GGU24 deviates the most from a first-order shape, and an improved fit is obtained by allowing reaction order to be optimized in the nucleation model.

Figure 13. Fit of sample GGU24 to a four-parameter nucleation kinetic model.

from the first-order to the second-order curve after the maximum reaction rate. As a result, it is not surprising that a slightly better fit to the GGU24 data is obtained by allowing n to float along with m, resulting in a reduction of the residuals from 0.454 to 0.372. The derived parameters are A ) 5.38E+14 s-1, E ) 56021 cal/mol, m ) 0.43, and n ) 1.28, only slightly different than reported earlier2 for n ) 1. The improved fit is shown in Figure 13. Although oil shale never met the optimistic prospects of the 1970s in the United States, it still produces significant quantities of oil in China and Estonia. Chemical kinetic expressions have been published earlier for Anvil Points Green River shale in the United States and for the Maoming and Fushun shales of China.3 Pyrolysis data was also acquired in our labora-


Reviews

Table 6. Modified Friedman and Coats-Redfern Analyses of the Kukersite Shale modified Fr

modified CR

1-x

A

E

A

E

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1.17E+14 2.36E+14 2.72E+14 3.03E+14 2.77E+14 1.99E+14 1.22E+14 5.51E+13 7.46E+12

53 539 54 606 54 887 55 095 55 013 54 588 53 941 52 891 50 321

2.65E+13 8.43E+13 1.46E+14 2.22E+14 2.73E+14 2.66E+14 2.19E+14 1.35E+14 6.39E+13

51 360 53 015 53 841 54 477 54 818 54 831 54 609 53 986 53 044

tory a few years ago for the Kukersite of Estonia but never published. The reproducibility of these data is among the best ever at LLNL, so it is a useful standard material. Approximate kinetic analysis of 13 Kukersite pyrolysis profiles from 0.9 to 51 °C/min gives A ) 7.0E+14 s-1 and E ) 56 136 cal/mol, with a standard error of 480 cal/mol. The profile averaged 1.10 times as broad as the corresponding first-order profile and had an asymmetry of 0.62, compared to the first-order value of 0.628. Because the profile is slightly broader without becoming significantly more symmetric, a Gaussian distribution is more appropriate than an nth-order reaction to account for the broadness. The approximate value of σ from the width is 1.2% of Eo. The modified Friedman and Coats-Redfern results are given in Table 6. Standard errors of the energies were typically 400 cal/mol for the Friedman method and 300 cal/mol for the Coats-Redfern method. Even though the variation in A and E is small, both sets of AE pairs fall on the same tight line when ln A is plotted versus A. Even though Kukersite is close to a first-order reaction, Friedman’s plot of ln(Axn) vs ln(x) is not linear, confirming that such plots are rarely useful. Nonlinear regression of three files (one at each heating rate) to the nth-order, Gaussian, and Weibull distribution models all reduced the residual sum of squares to less than one-third of the first-order model. All three obtained mean activation energies (∼54.6 kcal/mol) closer to that from the shift in Tmax (55.8 kcal/mol for these data subset), while the first-order fit compromised the profile shift for the profile width, yielding E ) 52.2 kcal/mol. The discrete model reduced the residuals another 5-fold, giving the result in Figure 14. The principal activation energy is within the range from eq 31 and the other distribution models. A final point of interest is that although the average decomposition kinetics of kukersite are slightly faster than for typical Green River oil shale (e.g., Tmax at 25 °C is 469 °C rather than 475480 °C), the slightly broader reaction profile for the kukersite means that the reaction is 99% complete at about the same time. A common mechanism quoted in the oil shale literature is the sequential conversion of kerogen to bitumen to oil and gas (i.e., eq 3), where bitumen is defined by being soluble in organic solvents but still contained within the oil shale and oil is defined as condensable liquid that has been distilled from the shale. Miknis et al.72 report bitumen and oil yields as a function of time (72) Miknis, F. P.; Turner, T. F.; Berdan, G. L.; Conn, P. J. Energy Fuels 1987, 1, 477-483.

Figure 14. Discrete activation-energy model fit to Estonian Kukersite, giving A ) 3.12E+14 s-1 and the following energy distribution, E(%of total): 47(0.08), 48(0.08), 49(0.42), 50(0.05), 52(2.87), 53(3.68), 55(83.56), 57(5.58), 58(2.77), 59(0.3), 61(0.62).

for both a Colorado and Kentucky oil shale, as well as review earlier literature. From the change in the maximum bitumen yield with temperature, they conclude that the activation energy for bitumen formation must be higher than for oil and gas for a sequential reaction mechanism. However, they also note that earlier data could be reconciled with theirs if most of the kerogen is converted directly to oil and gas and only part of it is converted indirectly via bitumen. Furthermore, the maximum bitumen yield argument is compromised by the fact that more coke and less oil is formed at low temperatures due to enhanced liquidphase retrograde reactions. More recently, Burnham et al.2 noted that a sequential reaction for oil formation is incompatible with isothermal pyrolysis results, because the maximum oil and gas generation rate is at initial time. Instead, they proposed that the alternate pathway model of eq 4 is more appropriate, consistent with the conclusion of Ziegel and Gorman.73 Only in cases for which there is an acceleratory oil generation phase should eq 3 be considered. In that case, however, Burnham et al.2 found the nucleation model to be a better and more robust fit. The Bakken shale is a major petroleum source rock in North Dakota and Saskatchewan. This sample was analyzed both at LLNL using a Pyromat II instrument (heating rates of 1, 7, and 50 °C/min) as well as at Advanced Fuel Research using a TGA-Plus instrument (weight loss plus tar by FTIR, heating rates of 3 and 30 °C/min).74 The Tmax at 25 °C/min agrees well for the FTIR-tar and total weight-loss measurements, indicating that most of the weight loss for this material is an oil-like substance. The Pyromat Tmax is about 5 °C lower, not bad agreement considering the completely independent methods of calibration and measurement. The pyrolysis profile width and asymmetry are similar for the three methods and suggest a relatively symmetric activation-energy distribution. The approximate kinetic parameters are also similar, with the activation energy from the Pyromat being intermediate. The complete regression analysis also gives similar distributions, (73) Ziegel, E. R.; Gorman, J. W. Technometrics 1980, 22, 139-151. (74) Burnham, A. K. Pyrolysis Kinetics for the Bakken Shale; Lawrence Livermore National Laboratory Report UCRL-ID-109622; 1992.

Reviews


Figure 15. Comparison of kinetic predictions for three sets of kinetic parameters for the Bakken shale derived from different measurement techniques. Table 7. Partial Summary of Kinetic Parameters Derived for the Bakken Shalea Tmax(25) exp/calc asym Aapprox Eapprox σapprox Adisc 40 000 41 000 42 000 43 000 44 000 45 000 46 000 47 000 48 000 49 000 50 000 51 000 52 000 53 000 54 000 55 000 56 000 57 000 58 000 59 000 a

pyromat

FTIR-tar

TGA

457.2 1.25 0.86 1.43E+13 49 931 2.1 3.84E+13 0 0 0 0 0.1 0.9 0.8 0.4 0 6.0 11.5 40.6 20.4 10.0 4.5 0 4.8 0 0 0

462.4 1.53 0.9 9.02E+13 53 045 3.2 1.71E+13 1.9 0 0 2.0 0 0 0 0 9.6 0 41.8 4.0 31.3 0.0 1.9 4.2 0 0 3.3 0

462.9 1.29 0.89 7.67E+11 46 142 2.5 1.30E+14 0 0 0 2.6 0 0 2.6 0 0.6 1.8 0 11.2 8.6 28.3 16.8 9.1 7.2 0 0 11.2

A in s-1, E in cal/mol, and σ in % of E.

except that the FTIR-tar activation energies are slightly lower than the TGA energies. This is because of the relative sparseness of the FTIR-tar data and the noise in the TGA data, suggesting there is more uncertainty in these kinetic parameters. Figure 15 compares the predictions of the three sets of kinetic parameters at three heating rates: (1) at 10 °C/min, in the midrange of their measurement, (2) at a heating rate of 1 °C/h, which gives generation temperatures similar to 72 h at the specified temperature, as is typical for hydrous pyrolysis, and (3) at a typical geological heating rate of 3 °C/million years. The TGAplus kinetic parameters are slightly wider, especially at the extremes of the conversion curve. The midgeneration temperatures at 10 °C/min mirror the differences in the Tmax at 25 °C/min given in Table 7. The extrapolation to hydrous-pyrolysis-like time frames produces oil in the expected temperature range of 300350 °C. The geologic predictions are similar, with oil being generated primarily between 100 and 150 °C. It is particularly interesting that the differences between midpyrolysis temperatures at 10 °C/min and the activa-

Table 8. Conversion-Dependent Arrhenius Parameters for the Ventura Basin Shale Sample by Modified Friedman (Rates) and Coats-Redfern (Integrals) Methods modified Friedman fraction energy A (s-1) reacted cal/mol for n ) 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

55 557 57 617 57 778 59 254 61 321 63 943 67 941 73 312 86 718

1.49E+15 5.53E+15 5.10E+15 1.18E+16 3.96E+16 1.85E+17 1.99E+18 4.44E+19 1.20E+23

A (s-1) Miura 1.04E+16 1.80E+16 1.06E+16 1.87E+16 5.21E+16 2.11E+17 2.07E+18 4.31E+19 1.29E+23

modified Coats-Redfern energy A (s-1) for cal/mol n)1 53 082 55 342 56 379 57 217 58 357 60 008 62 716 67 035 76 443

3.03E+14 1.36E+15 2.52E+15 4.04E+15 7.85E+15 2.12E+16 1.12E+17 1.57E+18 4.76E+20

A (s-1) Miura 1.57E+15 3.31E+15 3.85E+15 4.30E+15 6.17E+15 1.26E+16 5.07E+16 5.32E+17 1.13E+20

tion energies have canceled for the Pyromat and Bakkentar cases, giving virtually identical predictions over the main stage of oil generation. Although most kerogens have a relatively small variation in their frequency factor as a function of the extent of conversion, some kerogens do. One example is a medium-sulfur kerogen from the Ventura Basin, reported earlier by Reynolds et al.24 Results from the modified Friedman and Coats-Redfern analyses are given in Table 8, along with the associated frequency factors calculated from the equations of Miura. Keep in mind that the modified Coats-Redfern frequency factors given in earlier papers are in error by the factor ln(x/(1 - x)). The activation energy is broad enough for this sample that Miura’s frequency factors are probably more accurate than those determined by merely substituting n ) 1 into the respective intercepts. When the frequency factor deviates by this much from a constant value, a noticeable error can be made upon extrapolation to time-temperature scales much shorter or longer than those used to calibrate the parameters. As a result, improved extrapolation accuracy can be obtained by allowing ln(A) to be a linear function of E. The activation-energy distributions determined from the conventional and improved methods are given in Table 9. (The distributions have been simplified a little by combining some minor peaks contributing less than 1% of the reaction.) The activation energy and frequency factor are about the same at the peak of the distribution as well as can be interpolated. The sum of squared residuals decreases from 0.048 for the conventional discrete model to 0.031 for the extended discrete model. The physical significance of the highest energies and frequency factors in the extended model is questionable. One might even argue that much of the decrease in


Reviews

Table 9. Comparison of Kinetic Parameters for the Ventura Basin Shale Sample Determined by the Conventional and Extended Discrete Activation-Energy Distribution Model discrete model % of rxn

A (s-1)

0.19 0.2 0.45 1.12 1.83 9.77 19.06 20.41 15.72 12.88 6.97 4.65 3.23 0.77 1.19 1.56


extended discrete model

E (cal/mol) % of rxn 46 000 48 000 49 000 51 000 53 000 55 000 56 000 57 000 58 000 59 000 60 000 61 000 63 000 64 000 67 000 69 000

0.22 0.84 1.16 5.63 36.58 28.22 14.64 5.57 1.99 1.91 0.59 0.71 0.42 0.38 0.52 0.61

A (s-1)

E (cal/mol)


34 603 39 566 44 529 49 493 54 456 59 419 64 383 69 346 74 309 79 272 84 236 89 199 99 125 104 089 120 887 148 758

reactivity at high conversion should be due to decreased concentration of condensable groups, which have increasing difficulty finding each other. This effect should formally be reflected by a decreased frequency factor. However, the highest energy channels account for only a few percent of the reaction, so there is little practical consequence. Coal. The coal pyrolysis literature has many kinetic analyses using both nth-order and distributed reactivity models. In one of our publications on this topic, we reported extensive analyses using Gaussian and discrete activation-energy models and imply that they are superior to nth-order models. While the discrete model is clearly superior to either the Gaussian or nth-order models, eq 36 and Figure 6, along with the experience of Boudreau and Ruddick,34 suggests that the nth-order reaction has some advantages over the Gaussian model. This might explain the early success of the nth-order approach.49,50 Maki et al.75 recently criticized the use of a Gaussian distribution model for analyzing coal pyrolysis data. We agree that a single Gaussian is of marginal value. While a single Gaussian model fits the evolution of some individual species and the main part of the pyrolysis profile, it cannot fit the entire pyrolysis profile. An example of such a poor fit to the fraction reacted is given in Figure 16 (top) for their TGA data for Pittsburgh No. 8 coal from the Argonne Premium Coal Bank. In addition to the poor visual fit, the residual sums of squares vary by only 1% as A varies by 7 orders of magnitude, so the kinetic parameters derived by fitting fraction-reacted data are essentially meaningless. Maki et al.75 say further that this indicates a distributed activation-energy model with a constant frequency factor is a very poor approximation and cannot extrapolate well to flash heating conditions. While a variable frequency factor model does have some advantages for some samples, that advantage needs to be put in proper perspective. Table 10 gives the residual sum of squares of several more appropriate activation-energy distribution models having a single frequency factor. The best simple model is an nth-order reaction. A model with two (75) Maki, T.; Takatsuno, A.; Miura, K. Energy Fuels 1997, 11, 972977.

Figure 16. Comparison of the ability of three distributed activation-energy models to fit fractional weight-loss data for Pittsburgh No. 8 coal. The top figure shows the Gaussian model fit, the middle figure shows an nth-order fit, and the bottom figure shows a discrete-model fit.

parallel Gaussian distributions having the same frequency factor fits significantly better, but the best fit is to the discrete distribution model. Figure 16 shows the nth-order model gets the profile shape about right. The

Reviews


Figure 17. Activation-energy distribution for Pittsburgh No. 8 TGA data of Maki et al.74 Table 10. Summary of the Residual Sums of Squares for Fitting Various Distributed Activation-Energy Models to Fractional Weight-Loss Data from Maki et al.74 for Illinois No. 6 and Pittsburgh No. 8 Coal Gaussian Weibull Gaussian-nth nth-order 2 Gaussians discrete

ILL No. 6

PIT No. 8

0.188 0.069 0.020 0.019 0.005 0.001

0.799 0.270 0.197 0.195 0.013 0.011

discrete distribution fits the data almost perfectly. The parameters for this fit are given in Figure 17. Maki et al. also point out that the single Gaussian model does not extrapolate well to heating rates of 3000 °C/min. However, calculations with the discrete model parameters show that this failure is primarily due to the incorrect shape of the Gaussian distribution, not the constant A approximation. The fractions reacted calculated for heating at 3000 °C/min to 485 and 590 °C, respectively, and then holding there for 10 s, are 0.37 and 0.73, which are close to those estimated from their Figure 7. Given the suggestion48,76 that the single-frequency factor approximation may not be accurate for coal samples, because of their very wide energy distribution, we explored the utility of the extended discrete model on TGA data for Pittsburgh No. 8 coal kindly provided by Dr. K. Miura. First, conversion-dependent kinetic parameters are presented in Table 11. Kinetic parameters derived by the extended discrete distribution model are shown in Figure 18. The energies are evenly spaced, so the E and A values for the unlabeled channels are easily calculated. Allowing the frequency factor to increase with activation energy reduced the sum of squared residuals from 1.08E-02 to 0.74E-03. The fit is excellent, though not visually better than for the regular discrete model in Figure 16 for a figure that size. The nature of the minimization surface is shown in Figure 19. The modulation of the surface is due to the finite number of energies, but the resulting affect on b is negligible compared to typical experimental error. The physical significance of the highest energies and frequency factors is again questionable, and perhaps some of the decreased reactivity at the highest conversion could be due to a decreased frequency factor. (76) Burnham, A. K.; Schmidt, B. J.; Braun, R. L. Org. Geochem. 1995, 23, 931-939.

Figure 18. Activation-energy distribution for Pittsburgh No. 8 data fitted to the extended discrete model.

Figure 19. Minimization surface for Pittsburgh No. 8 data fitted to the extended discrete distribution model in which ln(A) ) a + bE.

Figure 20. Comparison of the dependence of A on E analyzing TGA data for Pittsburgh No. 8 coal by the Miura integral method and the extended discrete model.

A more careful comparison of the various methods provides important insights. First, notice that the A factor from Miura’s method decreases by nearly 10 times over the course of 20-40% conversion while the activation energy is approximately constant. This decrease is shown more clearly in Figure 20, which compares the finer-grained analysis provided by Miura with the LLNL analysis using his method. The decrease here is unmistakable and can be related to the breakdown of the broad distribution approximation of Miura’s method, as shown in Figure 7. The reaction-rate profile of the main decomposition stage is only 60% broader than a first-


Reviews

Table 11. Conversion-Dependent Kinetic Parameters for Pittsburgh No. 8 TGA Data of Maki et al.74 modified Friedman

modified Coats-Redfern

1-x

A (n ) 1)

A (Miura)

E (cal/mol)

std err

A (n ) 1)

A (Miura)

E (cal/mol)

std err

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1.59E+14 4.80E+14 1.81E+14 2.02E+14 6.23E+14 2.93E+16 5.87E+20 4.06E+25 1.03E+14

2.84E+15 2.59E+15 5.41E+14 4.37E+14 1.20E+15 6.23E+16 2.01E+21 2.58E+26 3.78E+14

53 829 56 596 55 824 56 557 58 999 65 972 83 645 108 197 72 066

1324 2626 1114 80 275 2909 4001 9740 98203

1.83E+13 9.75E+14 4.49E+14 3.26E+14 4.13E+14 1.80E+15 1.02E+18 3.51E+26 1.80E+20

9.48E+13 2.38E+15 6.85E+14 3.48E+14 3.25E+14 1.07E+15 4.62E+17 1.19E+26 4.26E+19

49 385 56 412 56 079 56 160 57 050 59 933 70 892 106 895 93 092

4442 1007 1219 663 429 1244 3118 12 970 79 698

Figure 21. Comparison of activation-energy distributions derived from the Miura integral method (equivalent to the modified Coats-Redfern method) and the extended discrete model. Both allow A to be a function of E.

order reaction, corresponding to an approximate Gaussian energy distribution of only 3-4% of Eo. As mentioned earlier, Miura’s approximation tends to break down for Gaussian distributions less than 5%. Perhaps part of the decrease in A is also due to some nth-order character associated with decreasing concentrations of groups capable of retrograde reactions. The frequency-factor function from the extended discrete model traverses the middle of that oscillation, thereby producing the same mean reactivity with a weaker dependence of A on E. The resulting activationenergy distributions from the Miura method and the extended discrete model are similar, as shown in Figure 21. Summary and Conclusions Well-chosen kinetic models fit the thermal decomposition data for synthetic polymers, oil-prone kerogens, and coal samples over a wide range of times and temperatures. The key to finding a model that will extrapolate well outside its calibration temperature range is to thoroughly decouple the effects of time, temperature, and extent of reaction. A variety of simple models with different characteristics are available from the literature, including ones that produce product over narrower and wider time-temperature intervals than first-order models. Linear polyolefin materials are usually described best by nucleation models, which have an initial acceleratory

phase for isothermal conditions and a narrow pyrolysis profile at a constant heating rate. Cellulose and wellpreserved alkal kerogens have the same kinetic characteristics. Other more complex materials, such as type II kerogens and coals, have isothermal profiles that are concave upward and are broader than a first-order reaction at a constant heating rate. These characteristics require a distributed reactivity model. A variety of distributed reactivity models are available. A pseudo-nth-order model is equivalent to having an exponential-like distribution of frequency factors, which can be mimicked accurately over a fairly broad temperature by an exponential-like distribution of activation energies. A pseudo-nth-order model usually provides a better fit than a Gaussian activation-energy distribution model for natural materials such as coals because of the asymmetry of the reactivity distribution. A Weibull energy distribution and a Gaussian distribution of nth-order models provide additional flexibility. The most versatile distributed reactivity models have a discrete energy distribution that is able to conform to the subtleties of the pyrolysis profile. The conventional discrete distribution assumes the same frequency factor for all parallel reactions. However, frequency factors derived from modified Friedman and Coats-Redfern analyses indicate that the uniform frequency factor approximation is not always valid. However, these methods are not as accurate for determining reactivity distribution parameters. Consequently, a new method has been derived in which the discrete activation-energy distribution is derived by assuming a linear relationship between the logarithm of the frequency factor and the activation energy. This model provides improved accuracy for the initial and final stages of the reaction for some samples when the kinetics are extrapolated far outside their range of calibration. This extension is most important for flash coal pyrolysis and natural gas generation. Acknowledgment. This work was performed under the auspices of the U.S. Department of Energy by the Lawrence Livermore National Laboratory under Contract W-7405-Eng-48. We acknowledge the helpful correspondence with Dr. Kouichi Miura that prompted us to discover our programming error for the modified Coats-Redfern method and to develop the extended discrete distribution model. We also appreciate the experimental contributions of Dr. John G. Reynolds. EF9800765

Global Kinetic Analysis of Complex Materials - Energy & Fuels (ACS

Recommend Documents