Energy & Fuels 1994,8, 1158-1167
1158
Art i d e s A New Distributed Activation Energy Model Using Weibull Distribution for the Representation of Complex Kinetics Chittur Chandrasekharan Lakshmanan* and Noam White BHP Research-Melbourne Laboratories, P.O. Box 264, Rosebank MDC, Clayton 31 69, Victoria, Australia Received March 21, 1994. Revised Manuscript Received June 10, 1994@
A new distributed activation energy model (DAEM) is developed using Weibull distribution to mathematically represent the kinetics of complex thermal transformations. Also developed is a grid search procedure t o estimate the kinetic parameters represented by the DAEM. This procedure uses approximate analytical solutions to the time-temperature integrals of DAEM. These approximations reduce the computational efforts involved in parameter estimation substantially and allow rapid analysis of the three-dimensional error surface of the objective function over a wide range of kinetic parameters while retaining the true nonlinearity of the model. The existence of multiplicity-a number of sets of kinetic parameters which describe the observed data, generally within the limits of experimental accuracy-is also demonstrated using experimental data on kerogen decomposition. Several sets of isothermal and nonisothermal experimental data are used to validate the DAEM. This work seems to be the first published one where Weibull distribution is used to describe the kinetics of complex thermal transformations.
Introduction Kinetic models t o quantify the effect of time and temperature on thermal decomposition of organic matter (e.g., pyrolysis) are often needed for the purposes of design and scale-up of chemical reactors in which such transformations take place. It may not be always possible to develop rigorous kinetic models using mechanisms based on fundamental structural chemistry. Therefore, properly validated, pragmatic engineering approaches are frequently used to develop these models. A number of models, based on finite and infinite numbers of first order reactions have been These models use continuous and discrete distributions of activation energies to account for the effect of temperature on the progress of the reactions. These kinds of representations reduce the number of unknowns (parameters of the model) t o a reasonable number which can then be estimated using parameter estimation techniques and experimental data. Often these estimation techniques, except for simple linear models, need numerical methods. Anthony and Howard1 used a normal distribution for the representation of activation energies. Lakshmanan et have published * Corresponding author. E-mail address:
[email protected]. @
Abstract uublished in Advance ACS Abstracts. Julv 15. 1994.
(1)Anthony, D.E.; Howard, J. B. AIChE J. 1976,22*(4),625-656.
(2)Agarwal, P.K.;Genetti, W. E.; Lee, Y. Y. Chem. Eng. Sci. 1986, 41 (9). 2373-2383. , ~~
I!
~~~
~
~~
(3)Lakshmanan, C.; White, N. Proc. 15th. Aust. Chem. Engg. Conf., CHEMECA-87, 1987, 44.1-44.7. ( 4 )Ungerer, P.; Pelet, R. Nature 1987, 327 (71, 52-54. (5)Sezen, Y.Int. J.Heat Mass Transfer 1988,31 (61,1319-1322. (6)Lakshmanan, C.C.;Bennett, M. L.;White, N. Energy Fuels 1991, 5,110-117.
a paper on the application of such models t o explain generation of petroleum from an exploration point of view. Included in that paper are the discussions on problems associated with multiplicity and correlation between kinetic parameters, namely, compensation effect. During the development of the kinetic model for petroleum generation, in geological time scale, certain problems were encountered by the authors which prompted them to develop a new distributed activation energy model (DAEM)using Weibull distribution. This statistical distribution has some interesting properties and generates a variety of distributions. Application of this distribution to a number of experimental data is the subject of this paper. It is anticipated that the model developed in this paper will be useful for a number of process chemical engineering applications (e.g., combustion of coal and oil shale). Before proceeding to the development of such a model, it is necessary to consider some fundamental properties of this distribution.
Properties of Weibull Distribution In 1951 Weibul17 derived a three-parameter probability density function known as the Weibull distribution. This was first applied t o study fatigue of materials. It is well suited to represent many empirical distributions. The probability density function RE)(as applied to describe the distribution of activation energies) for this distribution is given by (7) Weibull, W. J.Appl. Mech. 1951, 18, 293-296.
0887-0624/94/2508-1158$04.50/0 0 1994 American Chemical Society
Distributed Activation Energy Model
Energy &Fuels, Vol. 8, No. 6, 1994 1159
0, andB > 0 (1) for all other values of E , 7 , ,i3
=0
where E is the activation energy expressed in kcal/mol. There are three parameters, namely, q, the scale parameter; p, the shape parameter; and y , the threshold or location parameter characterizing the distribution. A number of different distributions can be generated by a suitable choice of these parameters. For /3 = 1, the Weibull distribution coincides with exponential distribution. For values of p >1, the distribution becomes “bell shaped”, but becomes positively skewed for values of p > 1. As p increases, the Weibull distribution approaches the normal distribution more and more closely. In fact, for p = 4,the Weibull and normal distributions become almost indistinguishable. Choosing a threshold value of activation energy ( y ) , implies that reactions with activation energies less than this value do not occur. This is a useful feature to overcome one of the problems associated with the parameter estimation step, namely, the estimation of negative activation energies resulting from numerical function minimisation methods. Negative activation energies have no physical significance. The mean (p)of the distribution, first moment, equals the mean activation energy and is given by p = y -t qF(;+
1)= E o
Eo is expressed in kcdmol. The variance of the distribution, second moment, is given by
(3) where the general form of the Gamma function with argument ( x 1)is given by
+
Note CJ and y are also expressed in kcal/mol.
Distributed Activation Energy Model (DAEM) Using Weibull Distribution Consider the weight loss of organic matter such as in the thermal decomposition of shale or coal, as given by the following first-order kinetic model:
Mt -- -kR, dt
At the initial time (t = O),Rt = 100%. k is the rate constant, with units of reciprocal time. This initial value problem can be integrated to give t
Rt = expJ-k
dt
(6)
0
Equation 6 gives the residual weight of organic matter, a measure of chemical conversion, as a function of time. Let the thermal decomposition be assumed to consist
of a large number of independent, parallel chemical reactions of first order. The effect of temperature on the rate constants is given by the Arrhenius law. It is now assumed that the activation energies (E)of these parallel reactions follow a Weibull distribution and the preexponential factors are equal to KO. With these assumptions, the residual weight of organic matter at any time t during its decomposition (isothermal and/or nonisothermal) is given by ca
[
1
R, = Jexp - J k o exp Y
{
-
R(To+Pt*)
Most laboratory experiments involving thermal decomposition of organic matter are conducted at temperatures linearly varying with time from a starting temperature. This linear heating rate is accounted by P in eq 7. If nonlinear heating rates are used, eq 7 and the equations derived from it will have to be modified accordingly. Equation 7 considers all reactions with positive values of activation energy above a threshold value y. RE)represents the probability density function given by eq 1. At this stage a distribution in KO is not considered.
Parameter Estimation Using Grid Search Technique: Outline of the Methodology A number of kinetic models using normal distribution have been published. These assign a fixed value of preexponential factor and estimate the mean and standard deviation of the distribution. For complex reactions such as pyrolysis and combustion, it is often difficult if not impossible to provide a valid justification for fixing the preexponential factor at an arbitrary value. Further, it has been shown by Lakshmanan et ala6that the preexponential factors in such models have a strong correlation, some times referred to as the compensation effect, with mean activation energy. In view of these findings, preexponential factor will also be estimated as one of the parameters of the kinetic model in this paper. Numerical values of the kinetic parameters (mean activation energy, standard deviation, and the preexponential factor) are required to use DAEM for a given thermal history. Estimation of mean activation energy and standard deviation corresponds to the estimation of the parameters characterising the Weibull distribution, namely, 7 , the scale parameter; p, the shape parameter; and y , the threshold or location parameter. Once these are known, eqs 2-4 can be used to obtain mean activation energy and standard deviation. As there are only three parameters to be estimated, a threedimensional ( y , p, and KO) grid search procedure has been adopted. TI, will be internally adjusted as will be shown later. This approach retains the true nonlinearity of the model in the estimation of its parameters. A region for grid search is defined by specifying the lower and upper limits of the values of y , p, and KO. The number of grid points in each direction is specified and the predicted values of residual weights are obtained using eq 7. The sum of the squares of residual, between predicted and experimental residual weights, (SSR) is obtained using the parameters which characterize the grid point. This process is repeated until SSR is
1160 Energy & Fuels, Vol. 8, No. 6, 1994
Lakshmanan and White
Table 1. Emerimental Data Used To Validate DAEM data source of data thermal history number and type of experiments Fission Track dieitised from Gleadow et al..1° isothermal 77 thermal annealing experiments Green et al.," with one observation at the end of each experiment Laslett et a1.,12 and Duddy et al.13 Fission Track digitised from Gleadow et al.,l0 nonisothermal 49 thermal annealing experiments Green et a1.,l1 with one observation at the Laslett et a1.,12and end of each experiment Duddy et aL13 Mackenzie digitised from Mackenzie and 31 observations from one 25 "C/min and Quigley Quigley14 thermogravimetric experiment (Rock-Eva1 dataa) at one heating rate Yallourn data obtained by the authors 5,15, and 25 Wmin 3 thermogravimetric experimenb with 16,28, and 19 observations, respectively Glen Davis data obtained by the authors 5,10, and 15 Wmin 3 thermogravimetric experiments with 25,21, and 20 observations, respectively Kimmeridge data obtained by the authors 10,15, and 20 Wmin 3 thermogravimetric experiments with 50,28, and 34 observations, respectively Ishiwatari et al. digitised from Ishiwatari et al.15 isothermal at 150,200, 3 sealed tube experiments 11, 11, and 310 "C and 11observations, respectively hydrous pyrolysis data obtained by the authors isothermal at 290 and 2 sealed tube experiments with 370 "C 7 and 9 observations, respectively anhydrous pyrolysis data obtained by the authors isothermal a t 290, 330, 3 sealed tube experiments with and 370 "C 10, 10, and 9 observations, respectively a Rock Eva1 Analysis: A special experimental procedure, widely used by petroleum exploration companies, to evaluate the hydrocarbon generation potential of kerogens. For the purposes of modelling, the data from a Rock-Eval Analysis can be considered as equivalent to a thermogravimetric data. See ref 16: Petroleum Formation and Occurrence; Tissot, B. P., Welte, D. H., Eds.; Springer-Verlag: New York, 1984; p 509. I
obtained for each grid point. Then the error surface (a three dimensional representation of the sum of the squares of deviation) is analyzed for its shape and the associated valley. The best set of parameters is then determined as the coordinates of that grid point which produced the lowest SSR. In this process, a number of sets of parameters may be found for which the SSR is either equal or approximately equal to the SSR corresponding to the best set of kinetic parameters. It must be noted that a detailed knowledge of the sensitivities of the parameters and of the errors associated with the experimental measurements are usually required to make this subjective consideration to define multiple sets. Computing predicted values of residual weight requires the evaluation of the time-temperature integral (inner integral) of eq 7. For complicated thermal histories, numerical evaluation of this integral may be the only choice. The time-temperature integral in D m M can be analytically evaluated only for isothermal conditions. But it is possible to obtain analytical values for the approximate forms of integral for nonisothermal conditions involving linear heating rates.
Evaluation of the Time-Temperature Integral Using Approximate Forms of the Integrand Consider eq 7 describing the thermal decomposition of organic matter. Using the approximation, valid for situations where RTIE
