Weibull Mixture Model for Modeling Nonisothermal Kinetics of

Aug 17, 2007 - been obtained by a mixture of Weibull distribution functions. ... proposed the nth-degree polynomial logistic regression model for fitt...
0 downloads 0 Views 92KB Size
Download PDF
J. Phys. Chem. B 2007, 111, 10681-10686

10681

Weibull Mixture Model for Modeling Nonisothermal Kinetics of Thermally Stimulated Solid-State Reactions: Application to Simulated and Real Kinetic Conversion Data Junmeng Cai* and Ronghou Liu* Biomass Energy Engineering Research Center, School of Agriculture & Biology, Shanghai Jiao Tong UniVersity, 800 Dongchuan Road, Shanghai 200240, P. R. China ReceiVed: May 15, 2007; In Final Form: June 29, 2007

The possibility of applying Weibull mixture model for the fitting of the nonisothermal kinetic conversion data has been investigated. It has been found that the kinetic conversion data at different heating rates can be successfully described by one or the linear combination of few Weibull distribution functions. Several simulated and real kinetic conversion traces have been analyzed. An optimal fitting of the kinetic conversion data has been obtained by a mixture of Weibull distribution functions. The results obtained have shown that the obtained conversion curves calculated by the model proposed in this paper are in agreement with the raw kinetic conversion data.

Introduction In the solid-sate kinetic analysis, it is convenient to express the reaction sequence in terms of the kinetic conversion degree, R, defined as

R)

m0 - m m0 - mf

(1)

where m is the mass or mass fraction of a sample at a certain time (t) or temperature (T), m0 and mf are the initial and final values of the reaction. The most common heating profile used for studying thermally stimulated solid-state reactions is the linear heating program.1 Therefore, this paper focuses on the nonisothermal kinetics in solids. The kinetic conversion (R - T) curves can be obtained through nonisothermal analysis. Different researchers have stated that the use of the differential conversion (dR/dT - T) curves instead of the kinetic conversion curves makes easier the identification of the kinetics,2-4 mainly because small changes in the conversion curves are magnified in the corresponding differential conversion curves. However, small errors in the conversion curves are also magnified in the differential conversion curves. There are some high fluctuations in the differential conversion curves, and some filtrations of the data are needed for smooth curves. A large number of smoothing algorithms have been developed, ranging from a polynomial algorithm to the Gaussian smoothing algorithm.5-8 The kinetic analysis of complex reactions, such as overlapping reactions, has received the attention of a great number of authors.9-11 Overlapping reactions could be detected looking at the dependence of the activation energy on the conversion.12 To separate the individual steps of overlapping processes, many methods have been proposed.13-18 In the paper of Wagner et al.,15 the sum of two Gaussian functions were used for fitting exothermal peaks obtained through differential scanning calorimetry of amorphous semiconducting alloys. Naya et al.16 proposed the nth-degree polynomial logistic regression model * Corresponding auithors. E-mail: Liu R., [email protected]; Cai J., [email protected] or [email protected].

for fitting the TGA curves. In the papers of Cao et al.17 and Barbadillo et al.,18 the logistic mixture function has been applied for description of thermogravimetric analysis (TGA) curve. Concerned about overcoming the smoothing problem when analyzing the nonisothermal data and deconvolution overlapping reactions, in this paper we aim to show a new method to model the nonisothermal kinetic conversion data by one or the linear combination of few Weibull distribution functions. Weibull Mixture Model The model proposes decomposition of the kinetic conversion curve into several Weibull distribution functions, assuming that decomposition kinetics of each component of the sample is represented by one or the sum of few functions. Even in the case of homogeneous materials, it is supposed that several different structures may exist, each one following its specific kinetics that may be different from the others.17,18 In this model, it is assumed that a kinetic conversional curve may be fitted by one or the linear combination of few Weibull distribution functions: k

R(T) )

wi(1 - e-[(T-T )/η ] i) ∑ i)1 0

i

β

(2)

where i ) 1, 2, ..., k represent different components from the kinetic conversion process, not necessarily different chemical compounds, T is the temperature (T g T0, T0 is the starting temperature of the process and acts as the location parameter of Weibull distribution function), wi stands for the weight of the ith component of Weibull mixture model, ηi is the scale parameter of the ith Weibull distribution function, ηi > 0, and βi is the shape parameter of the ith Weibull distribution function, βi > 0. According to the dimensional analysis, βi is expressed in dimensionless, and ηi is expressed in K. When T f ∞, the value of the R(T) function has to tend to 1. Therefore, from eq 2, the following relation can be obtained: k

wi ) 1 ∑ i)1

10.1021/jp0737092 CCC: $37.00 © 2007 American Chemical Society Published on Web 08/17/2007

(3)

10682 J. Phys. Chem. B, Vol. 111, No. 36, 2007

Cai and Liu

TABLE 1: η and β of the Weibull Distribution Function Obtained from the Simulated First-Order Model Process at Different Heating Rates with Corresponding Values of R2 heating rate (K min-1)

η (K)

β

R2

8 16 32 64

124.51812 141.63583 159.67474 178.70715

4.91478 5.33373 5.72947 6.10023

0.99952 0.99961 0.99967 0.99971

Estimation of the Parameters For fitting the nonisothermal kinetic conversion data with Weibull mixture model, some estimation of the parameters in eq 2 is needed. To evaluate the parameters, the residual sum of squares of deviations of actual and calculated model prediction values is commonly used. It is defined by the following relation:19 nd

RSS )

(Ra,j - Rc,j)2 ∑ j)1

(4)

where Ra,j is the actual value of the kinetic conversion degree in the data point j, Rc,j is the value calculated from eq 2 in the corresponding point, and nd is the number of data points. There are a number of algorithms for minimizing eq 4 enabling us to find the best values of the parameters. The fundamentals of the optimization methods were described in detail by Jongen et al.20 In this study, the Levenberg-Marquardt method has been used for the calculation of the values of the parameters that minimize eq 4. The Levenberg-Marquardt method is a method of nonlinear optimization, which minimizes the function that is a sum of squares of nonlinear functions. Many methods use the gradient of the function to be minimized. The Levenberg-Marquardt method uses a Jacobian instead of gradient. More detailed information about the LevenbergMarquardt method can be found in the literature.21 For performing the Levenberg-Marquardt method, either general purposed mathematical software or a computer program developed in any programming language is used. In this study,

the Origin software has been employed for performing the optimization procedure. One of the problems that appear when using the method is to choose some initial values for the parameters to estimate.17 To do this, one possibility consists of trying to estimate the inflection point by observing the kinetic conversion curve. Because this method is not easy and requires previous expertise, we have proposed a method based on the idea of assuming that the data follow one Weibull distribution function, eq 2. We use the linear form of eq 2, given by the relation

ln[-ln(1 - R(T))] ) β ln(T - T0) - β ln η

(5)

The values of β and η would be obtained from the slope and intercept of the plot of ln[-ln(1 - R(T))] versus ln[T - T0]. Application to a First-Order Model Process For nonisothermal conditions, when the temperature varies with the time with a constant heating rate, q ) dT/dt, the reaction rate of a theoretical first-order model process can be kinetically described by the following expression:

dR A -E/RT (1 - R) ) e dT q

(6)

where A and E are the Arrhenius parameters (frequency factor and activation energy, respectively) and R is the gas constant. The integrated form of eq 6 is represented as

A -ln(1 - R) ) ( q

∫0Te-E/RT dT - ∫0T e-E/RT dT) 0

(7)

where T0 is the starting temperature of linear heating program. Let us introduce the following function

I(E,T) )

∫0Te-E/RT dT

(8)

The above function is generally termed the temperature integral or Arrhenius integral,22-24 which does not have an exact

Figure 1. Simulated data obtained from a first-order model process and the Weibull distribution function model fitting: (a) 8 K min-1; (b) 16 K min-1; (c) 32 K min-1; (d) 64 K min-1; (9) simulated data; (s) model prediction.

Weibull Mixture Model for Nonisothermal Kinetics

J. Phys. Chem. B, Vol. 111, No. 36, 2007 10683

Figure 2. Simulated dR/dT data and its corresponding model prediction: (a) 8 K min-1; (b) 16 K min-1; (c) 32 K min-1; (d) 64 K min-1; (9) simulated data; (s) model prediction.

Figure 3. Comparison between the simulated data and model prediction of the distributed activation energy model process: (a) 5 K min-1; (b) 10 K min-1; (c) 20 K min-1; (d) 40 K min-1; (9) simulated data; (s) model prediction.

analytical solution but can be evaluated by numerical integration. With regard to eq 8, eq 7 takes the form

R ) 1 - e-(A/q)[I(E,T)-I(E,T0)]

(9)

The nonisothermal simulated data, which unlike experimental data are not affected by noises, were analyzed. The simulated data for A ) 1010 min-1, E ) 125.4 kJ mol-1, T0 ) 500 K, and the linear heating rates 8, 16, 32, and 64 K min-1 were calculated through the numerical calculation of eq 9 using the Mathematica software, which is a powerful general computation system. Applying one Weibull distribution function to model the data, the values of the parameters for all heating rates were established and included in Table 1. The overlay of the simulated data and the corresponding model prediction is shown in Figure 1, where it can be seen that the model matches the data very well.

In the kinetic analysis, the differential kinetic conversion data are usually used. It would be of interest to compare the accuracy of the corresponding Weibull distribution function model for the estimation of differential kinetic conversion data. Substituting eq 9 into eq 6, the theoretical expression of dR/ dT can be obtained

dR A -E/RT-(A/q)[I(E,T)-I(E,T0)] ) e dT q

(10)

Additionally, the reaction rate of the first-order model process can be given by first derivative of one Weibull distribution function with respect to temperature. Thus, we can get the fitting dR/dT:

(

)

dR β T - T0 ) dT η η

β-1

e-[(T-T0)/η]

β

(11)

10684 J. Phys. Chem. B, Vol. 111, No. 36, 2007

Cai and Liu Application to a Distributed Activation Energy Model Process The distributed activation energy model has been proved successful in describing the decomposition kinetics of various complex materials, such as coal,25 biomass,26 and oil shale.27 The model assumes that the distribution of reactivity by a set of independent and parallel reactions, each with its own frequency factor and activation energy. 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.28,29 Assumptions and restrictions of the distributed activation energy model, and the derivation of its equation can be found in the literature.30

Figure 4. Overlay of the raw data, the Weibull mixture model prediction, and the single model components obtained by the fitting in the case of an epoxy-hexahydrophthalic anhydride.

heating rate (K min-1)

η

β

R2

5 10 20 40

171.93509 186.18762 201.11778 216.77885

3.58287 3.7926 3.99961 4.2035

0.99997 0.99995 0.99992 0.99989

ηi (K)

βi

wi

1 2 3

182.4589 94.53706 268.02999

12.00298 4.17204 14.76312

0.34021 0.10734 0.55245

T -E/RT T0e

dT

g(E) dE

(12)

In the above equation, g(E) is the distribution of activation energies. Usually, g(E) is taken to be a Gaussian distribution with a mean activation energy of E0 and a standard deviation σ, so

g(E) )

TABLE 3: Values of the Parameters of the Weibull Mixture Model Obtained in the Case of an Epoxy-hexahydrophthalic Anhydride components

∫0∞e-(A/q) ∫

R(T) ) 1 -

TABLE 2: Weibull Mixture Model Parameters Obtained in the Case of the Distributed Activation Energy Model Process

1 -0.5(E-E0)2/σ2 e σx2π

(13)

Using eq 13, one may express eq 12 for the yield as

R(T) ) 1 -

The simulated dR/dT data and its corresponding model prediction are shown in Figure 2. It can be observed that the model prediction fits the simulated dR/dT data very well from Figure 2.

1 σx2π

∫0∞e-(A/q) ∫

T -E/RT T0e

dT-[(E-E0)2/2σ2]

dE (14)

The simulation data for the distributed activation energy model process are used. The parameters of eq 14 have been set as follows: E0 ) 125 kJ mol-1, σ ) 10 kJ mol-1, A ) 9 × 1010 min-1, T0 ) 400 K. The process has been simulated at

Figure 5. Overlay of the conversion data obtained from the foamed polyurethane decomposition in air, the Weibull mixture model fitting, and its components.

TABLE 4: Values of the Parameters of the Weibull Mixture Models Obtained in the Case of the Foamed Polyurethane Decomposition in Air heating rate (K min-1)

η1 (K)

β1

w1

η2 (K)

β2

w2

R2

3 6 9

130.96746 137.09918 140.54804

3.89954 4.06908 3.95711

0.28107 0.27203 0.28263

305.12054 349.57246 378.17207

3.98555 3.88266 3.91852

0.71893 0.72797 0.71737

0.99991 0.9999 0.99989

Weibull Mixture Model for Nonisothermal Kinetics four linear heating rates of 5, 10, 20, and 40 K min-1. A number of mathematical approaches have been pursued to deal with the double integration of eq 14.31,32 In this study, Simpson’s onethird rule33 for integration has been used for the numerical solution of eq 14. For more details of the numerical calculation of the distributed activation energy model equation, the readers are referred to the literature.34,35 Applying one Weibull distribution function to model the simulated data, the values of the parameters for all heating rates were established and included in Table 2. Figure 3 shows the overlay of the simulated data and the model fitting. It can be obtained that the proposed method allows for optimizing the fitting of the conversion curves of the distributed activation energy model process from Figure 3. Application to a Kinetic Conversion Curve Obtained from an Epoxy-Hexahydrophthalic Anhydride The kinetic conversion data of the degradation of an epoxyhexahydrophthalic anhydride was obtained from the literature.17 The experiment was performed using a single heating ramp from 375 to 720 K at a heating rate of 10 K min-1 and a purge of argon at a rate of 50 mL min-1. In this experiment a curing reaction and hardener volatilization occur before epoxy resin degradation, giving a complex conversion curve where threestage processes, two of them clearly overlapped, are evident.16 The best fitting was obtained with three Weibull distribution functions, and considering the residual standard error, a higher number does not improve significantly the quality of the fitting. The estimated parameters are listed in Table 3. Figure 4 shows the overlay of the raw data, the Weibull mixture model prediction, and the single model components obtained by the fitting. The first process corresponds to reactant volatilization, and the others correspond to different decomposition processes. The fitting matches the kinetic conversion data from Figure 4. Application to Foamed Polyurethane Decomposition in Air The experimental data on the thermal decomposition of foamed polyurethane in air have been analyzed.36 The process was studied by thermogravimetry under synthetic air at different heating rates (3, 6, and 9 K min-1). An optimal fitting was obtained with two Weibull mixture model components. Table 4 lists values of the parameters of Weibull mixture models corresponding to three heating rates. The fitting of the conversion data by the mixture of Weibull distribution functions are shown in Figure 5, where a very good match between the experimental data and the model prediction can be obtained. Conclusions In this paper, we have developed a new method to optimize the fitting of the kinetic conversion data by one or the linear combination of few Weibull distribution functions. It overcomes the noise problem when analyzing the kinetic conversion data and the differential conversion data. The separation into single functions provides an approach to separate some overlapping processes, because some processes can be modeled by one or the sum of few single functions. The computations have been successfully performed for (1) a first-order model process, (2) a distributed activation energy model process, (3) an epoxy-hexahydrophthalic anhydride degradation, and (4) a thermal decomposition of foamed polyurethane in air.

J. Phys. Chem. B, Vol. 111, No. 36, 2007 10685 References and Notes (1) Criado, J. M.; Pe´rez-Maqueda, L. A.; Sa´nchez-Jime´nez, P. E. Dependence of the preexponential factor on temperature: Errors in the activation energies calculated by assuming that A is constant. J. Thermal Anal. Calorim. 2005, 82, 671-675. (2) Friedman, H. L. Kinetics of thermal degradation of char-foaming plastics from thermogravimetry - Application to a phenolic resin. Polym. Sci. 1963, 6C, 183-195. (3) Budrugeac, P. Differential non-linear isoconversional procedure for evaluating the activation energy of non-isothermal reactions. J. Thermal Anal. Calorim. 2002, 68, 131-139. (4) Caballero, J. A.; Conesa, J. A. Mathematical considerations for nonisothermal kinetics in thermal decomposition. J. Anal. Appl. Pyrol. 2005, 73, 85-100. (5) Enke, C. G. Signal-to-noise ratio enhancement by least-squares polynomial smoothing. Anal. Chem. 1976, 48, 705-712. (6) Whittem, R. N.; Stuart, W. I.; Levy, J. H. Smoothing and differentiation of thermogravimetric data by digital filters. Thermochim. Acta 1982, 57 (2), 235-239. (7) Va´rhegyi, G.; Till, F. Computer processing of thermogravimetricmass spectrometric and high pressure thermogravimetry data. Part 2. Smoothing and differentiation. Thermochim. Acta 1999, 329, 141-145. (8) Liu, N.; Chen, H.; Shu, L.; Zong, R.; Yao, B.; Statheropoulos, M. Gaussian smoothing strategy of thermogravimetric data of biomass materials in an air atmosphere. Ind. Eng. Chem. Res. 2004, 4087-4096. (9) Criado, J. M.; Gonza´lez, M.; Ortega, A.; Real, C. Discrimination of the kinetic model of overlapping solid-state reactions from non-isothermal data. J. Thermal Anal. 1988, 34, 1387-1396. (10) Ross, J.; Vlad, M. O. Nonlinear kinetics and new approaches to complex reaction mechanisms. Annu. ReV. Phys. Chem. 1999, 50, 5179. (11) Wilburn, F. W. Kinetics of overlapping reactions. Thermochim. Acta 2000, 354, 99-105. (12) Flynn, J. H. The effect of heating rate upon the coupling of complex reactions. I. Independent and competitive reactions. Thermochim. Acta 1980, 37 (2), 225-238. (13) Kurajica, S.; Bezjak, A.; Tkalcˇec, E. Resolution of overlapping peaks and the determination of kinetic parameters for the crystallization of multicomponent system from DTA or DSC curves: I. Non-isothermal kinetics. Thermochim. Acta 1996, 288, 123-135. (14) Kurajica, S.; Bezjak, A.; Tkalcˇec, E. Resolution of overlapping peaks and the determination of kinetic parameters for the crystallization of multicomponent system from DTA or DSC curves. II. Isothermal kinetics. Thermochim. Acta 2000, 360, 63-70. (15) Wagner, C.; Va´zquez, J.; Villares, P.; Jime´nez-Garay, R. A theoretical method for resolving overlapping peaks in differential scanning calorimetry. Mater. Lett. 1994, 18 (5-6), 280-285. (16) Naya, S.; Cao, R.; Artiaga, R. Local polynomial estimation of TGA derivatives using logistic regression for pilot bandwidth selection. Thermochim. Acta 2003, 406, 177-183. (17) Cao, R.; Naya, S.; Artiaga, R.; Garcı´a, A.; Varela, A. Logistic approach to polymer degradation in dynamic TGA. Polym. Degrad. Stab. 2004, 85, 667-674. (18) Barbadillo, F.; Fuentes, A.; Naya, S.; Cao, R.; Mier, J. L.; Artiaga, R. Evaluating the logistic mixture model on real and simulated TG curves. J. Thermal Anal. Calorim. 2007, 87, 223-227. (19) Slova´k, V. Determination of kinetic parameters by direct non-linear regression from TG curves. Thermochim. Acta 2001, 372, 175-182. (20) Jongen, H. T.; Meer, K.; Triesch, E. Optimization Theorey; Kluwer Academic Publisher: Boston, 2004. (21) Chong, E. K. P.; Z˙ ak, S. H. An Introduction to Optimization; John Wiley & Sons, Inc.: New York, 2001. (22) Cai, J.; Yao, F.; Yi, W.; He, F. New temperature integral approximation for nonisothermal kinetics. AIChE J. 2006, 52 (4), 15541557. (23) Cai, J.; He, F.; Yi, W.; Yao, F. A new formula approximating the Arrhenius integral to perform the nonisothermal kinetics. Chem. Eng. J. 2006, 124, 15-18. (24) Straink, M. J. Activation energy determination for linear heating experiments: deviations due to neglecting the low temperature end of the temperature integral. J. Mater. Sci. 2006, doi: 10.1007/s10853-006-10677. (25) Gu¨nes¸ , M.; Gu¨nes¸ , S. A study on thermal decomposition kinetics of some Turkish coals. Energy Source, Part A: RecoVery, Utilization, and EnVironmental Effects 2005, 27 (8), 749-759. (26) Cai, J.; Ji, L. Pattern search method for determination of DAEM kinetic parameters from nonisothermal TGA data of biomass. J. Math. Chem. 2006, doi: 10.1007/s10910-006-9130-9. (27) Schenk, H. J.; Dieckmann, V. Prediction of petroleum formation: the influence of laboratory heating rates on kinetic parameters and geological extrapolations. Marine Petroleum Geol. 2004, 21 (1), 79-95.

10686 J. Phys. Chem. B, Vol. 111, No. 36, 2007 (28) Burnham, A. K.; Braun, R. L. Global kinetic analysis of complex materials. Energy Fuels 1999, 13 (1), 1-22. (29) Rostami, A. A.; Hajaligol, M. R.; Wrenn, S. E. A biomass pyrolysis sub-model for CFD applications. Fuel 2004, 83, 1519-1525. (30) Please, C. P.; McGuinness, M. J.; McElwain, D. L. S. Approximations to the distributed activation energy model for the pyrolysis of coal. Combust. Flame 2003, 133 (1-2), 107-117. (31) Donskoi, E.; McElwain, D. L. S. Optimization of coal pyrolysis modeling. Combust. Flame 2000, 122, 359-367. (32) Va´rhegyi, G.; Szabo´, P. Kinetics of charcoal devolatilization. Energy Fuels 2002, 16, 724-731.

Cai and Liu (33) Su¨li, E.; Mayers, D. F. An Introduction to Numerical Analysis; New York Cambridge University Press: Cambridge, U.K., 2003. (34) Gu¨nes¸ , M.; Gu¨nes¸ , S. The influences of various parameters on the numerical solution of nonisothermal DAEM equation. Thermochim. Acta 1999, 336, 93-96. (35) Cai, J.; He, F.; Yao, F. Nonisothermal nth-order DAEM equation and its parametric study - use in the kinetic analysis of biomass pyrolysis. J. Math. Chem. 2006, doi: 10.1007/s10910-006-9151-4. (36) Mamleev, V.; Bourbigot, S.; Bras, M. L.; Duquesne, S.; Sˇ esta´k, J. Thermogravimetric analysis of multistage decomposition of materials. Phys. Chem. Chem. Phys. 2000, 2, 4796-4803.