J. Phys. Chem. 1994, 98, 10598-10605
10598
Parametrization of Chemically-Activated Reactions Involving Isomerization Andrei Kazakov, Hai Wang, and Michael Frenklach’ Department of Materials Science and Engineering, The Pennsylvania State University, University Park, Pennsylvania 16802-2303 Received: July 18, 1994@ Parametrization of theoretical rate coefficients for chemically-activated reactions involving isomerization of activated complexes was examined. The theoretical rate coefficients were computed for chosen test reaction systems by using the Rice-Ramsperger-Kassel-Marcus theory. The low- and high-pressure-limit rate coefficients were found to be proportional to [MI’, where [MI is the bath gas density and i identifies a specific reaction channel. Several parametrization approaches to the pressure dependence of the rate coefficients were tested. The formulas based on a falloff broadening factor and the generalized interpolation between the pressure limits produced generally large errors. Modifications were introduced that improve the accuracy of these previously suggested formulas. The most promising parametrization formula discovered in the present study contains only three temperature-dependent adjustable parameters and demonstrates a high level of accuracy, with maximum errors below 1% for most cases tested, including both chemically-activated and unimolecular reactions.
Introduction Detailed modeling of gas-phase chemical kinetics has become a powerful tool in studies of combustion, atmospheric, and materials chemistry. Such models require the knowledge of reaction rate coefficients over wide ranges of experimental conditions. Many reactions of practical interest proceed through formation of energized addusts and hence their rate coefficients depend not only on temperature but also on pressure. Theoretical evaluation of rate coefficients using existing statistical theories as RRKM,1-3 SACM,’-’ or QRRK839can provide the necessary information with reasonable accuracy, and the qualie of theoretical predictions improved remarkably in recent years. However, repetitive calculations of the rate coefficients are computationally demanding and thus cannot be used within the framework of kinetic modeling. The general practice, therefore, is to use empirical or semiempirical approximate formulas. Parametrization of simple unimolecular processes has been extensively studied.10-20 Less attention has been given to the fitting of rate coefficients for more complex processes, such as those involving chemical activation and isomerization of the energized adducts. The main unsolved problem is the theoretical evaluation of the rate coefficients themselves, as the rigorous solution of master equations for such systems is still under d e v e l ~ p m e n t . ~ Golden ~ - ~ ~ and c o - ~ o r k e r s ~suggested ~-~~ a simplified RRKM treatment of chemically-activated reactions using the pseudo-steady-state approximation (PSSA) for activated complex and Troe’s weak-collision model” for its collisional stabilization. The approach was extended to more complex chemically-activated systems, including multichannel dissociation and isomeri~ation.~~-*~ Similar formalism was used in QRRK calculations.8~30-34However, no analysis was undertaken of rate coefficient parametrization, except for a s u g g e s t i ~ to n ~use ~ ~Troe’s ~ ~ formalism of pressure broadening. In this paper, we extend further the PSSA approach and cast the results in a form amenable systematic parametrization for reaction systems of higher complexity, i.e., those including isomerization of energized adducts. The PSSA method and the resulting expressions are given in the next section. We then review the approaches to parametrization pertinent to the present @
Abstract published in Advance ACS Absrrucfs, September 1, 1994.
0022-3654/94/2098- 10598$04.50/0
Figure 1. Schematicdiagram of chemically-activated process involving isomerization of activated complexes, T i .
work and describe the reaction systems chosen as test cases. After that we report the results obtained with each of the principal parametrization forms and present a new formalism which offers the best accuracy in representation of the theoretical rate coefficients, for both unimolecular and chemically-activated bimolecular reactions. We conclude with a brief summary.
Pseudo-Steady-StateAnalysis A complex chemical activation process including multiple isomerizations of activated complexes can be presented schematically as shown in Figure 1. The process starts with initial reactants A and B, which form the first energized complex TI. It may undergo a sequence of isomerizations resulting in activated complexes Tz, ...,T,. Each activated complex Ti may, in turn, dissociate with the formation of products Di or be stabilized by collisions with the bath gas resulting in a stabilized adduct Si. Upon application of the pseudo-steady-state approximation to all activated complexes TI, ...,T,, we obtain rate expressions for all reaction channels. Thus, for overall reactions
and
i = 1,2, ..., n, we have the following expressions for the thermal
0 1994 American Chemical Society
J. Phys. Chem., Vol. 98, No. 41, 1994 10599
Parametrization of Chemically-Activated Reactions rate coefficients
high-pressure rate coefficients (see Appendix), i
(3)
(104
km(s,i)=
i
m
ko(s,i)= / ? i w . j m ( n Z j , o ) P ( EdE ) = ko(s,i) [MI ' Eo j=l
(4)
1
-Lo
1
where
,
I
( n k j ( E ) ) P ( E )dE
km(s,i)[MI'-' (lob)
J=1
j= 1
2, =
k,(d,i) =
j= 1
In eqs 10 and 11, [MI is the density of the bath gas and In these equations, Kq is the equilibrium constant of reaction A B =+ S i , P(E) the Boltzmann distribution function, ,and EO the energy barrier for the transformation of S I to TI. All the microscopic rate coefficients appearing in eqs 3-5 are defined in Figure 1. The above equations can be readily extended for some special cases of interest. For example, if an activated complex Ti has m dissociation channels, it can be taken into account by modifying the corresponding expression of Zi as
+
zi =
ki(E) m
2. J.0 = Z.(k J SJ. = 0, i = 1, ...,J)
(12)
Equations 10 and 11 yield the pressure dependence of the low- and high-pressure-limit rate coefficients as power-law functions. It is worthwhile to point out that the functional form of the pressure dependence is determined only by the sequential number i of the corresponding activated complex in the consecutive isomerization (Figure 1) and is not affected by the occurrence of multichannel dissociations or branched isomerizations.
(6) Parametrization Approaches
k-i(E) + kS,i(E)+ ckD,i,r(E)- Zi+lk-(i+l)(E)
The mathematical development presented above shows that each reaction channel has well-defined low- and high-pressure where k ~ , i , is , the microscopic rate coefficient of the rth limits, and thus one can attempt to apply interpolation formulas dissociation channel of Ti. Equations 3-5 represent an extenused for unimolecular reactions. The first approach we consider sion of the approach used by Golden and c o - ~ o r k e r s . * ~ ~ ~ ~ follows the formalism suggested by Troe,11-13 with the notion W e ~ t m o r e l a n dobtained ~~ similar equations for QRRK calculaof correcting the PSSA-derived rate coefficient by a broadening tions. factor, F, For the description of the collisional stabilization, we adopted the weak-collision model suggested by Troe' for unimolecular k(i) = kh(i)F (13) reactions; Le., rate coefficient ks,i is considered to be a product where of the collision frequency, wi, and collision efficiency, pi, F l
The collision efficiency was calculated by using the approximation suggested by Gilbert et al.13
where FE is the energy-dependent density of states
(9)
and kg is the Boltzmann constant. The asymptotic behavior of eqs 3-5, in the limits of small and large collision frequencies, defines the respective low- and
In eqs 13 and 14, h(i)and km(i)are the low- and high-pressure limits of k(i), while the latter is an approximation sought for k(s,i) or k(d,i) given by eq 3 or 4, respectively. For a single-activated-complex sequence (Le., for n = l), eq 14 gives the familiar Lindemann rate coefficient derived by assuming a steady state for the activated complex and neglecting the energy dependence of microscopic rate coefficients. However, a similar treatment for n > l results in expressions mathematically more involved than eq 14, with the complexity and number of parameters increasing with n but without exhibiting accuracy significantly higher than eq 14, as indicated by computational tests performed in the course of the present study. We thus continued with eq 14, which offers numerical simplicity, the number of parameters independent of n, and the correct asymptotic behavior for all chemically-activated reaction channels. We note that kh defined by eq 14 is in essence one-
Kazakov et al.
10600 J. Phys. Chem., Vol. 98, No. 41, 1994 half the harmonic mean35 of the two limits. Hence, we will refer to eq 14 as the harmonic-mean-of-limits formula (HML). The broadening factor F is usually present in the form
F = [Fc(T)]x'pr' where T is the temperature, P, the reduced pressure,
(15)
2kpr=1(i)
F, = -
(16) ko and k p , = l ( i ) the rate coefficient calculated at P, = 1. For unimolecular reactions, Pr is defined as the ratio of the lowand high-pressure rate coefficients. To be consistent with this, we generalize the definition of the reduced pressure to (17) which makes P,(i) proportional to [MI for all reaction channels
i. For the functional form of X(P,) appearing in eq 15, Troel1J2 suggested a simple relation, a Lorentzian, to describe the falloff broadening for unimolecular reactions at low temperatures
x = { 1 + (log P,)2}-1
(18)
At higher temperatures the falloff broadening exhibits asymmetry and Gilbert et al.13 proposed a modified form of eq 18
X = {1
+ [(log P, + c)/(N - d(l0g P, + c))I2}-' (19)
where parameters c, d, and N were fitted over a series of unimolecular reactions. Recently, Wang and Frenklachlg suggested a Gaussian form of pressure broadening,
x = exp{-[(log
P, - a>/a]*}
(20)
where parameters a and o were fitted as polynomial functions of temperature. Tested on several unimolecular reactions, eq 20 was found to give a better accuracy than eqs 18 and 19. The difference between the Gaussian and Lorentzian formalisms was found to be caused by a narrower peak and wider spreading tails of the latter. We note that eq 18 represents the functional form of the Cauchy distribution, known to possess no finite moments due to its long tail^.^^,^^ Gardiner14proposed a different approach for parametrization of a unimolecular rate coefficient. Instead of correcting an approximate rate coefficient, he suggested to interpolate between the low- and high-pressure limits, and k,, as kuni= (kt
+ k:)"'
+ kf(i)]""
+
+
can be viewed as an approximation to eq 22 and that eq 22 faithfully describes the falloff broadening, justifies the use of eq 14 in Troe's broadening formula, eq 13, as applied to chemically-activated reactions involving isomerization of activated complexes.
Test Cases Test calculations were performed at the Rice-RamspergerKassel-Marcus (RRKM) level of theory for two reaction systems: the addition of acetylene to vinyl, depicted schematically in Figure 2a, and the addition of acetylene to n-C4H3 (1buten-3-yn-l-y1), shown in Figure 2b. Following the notations introduced earlier, we identify the individual channels of the two reaction systems as follows: C2H3
+ C2H2 -nC4H.j k(s,1)
(23d
+ C2H2-C,H4 + H
C2H3
k(d2)
and
where parameter a is considered to be a function of temperature. This type of interpolation was previously shown by Churchill and to be applicable to a myriad of physical processes with asymptotic behavior. The mathematical form of eq 21 is that of the generalized mean35of the two limits, and hence we will refer to it as the generalized-mean-of-limits formula (GML). Applying it to the chemically-activated reaction system of interest to the present study, we obtain
k(i) = [k;f(i)
Figure 2. Potential energy diagrams of the test reaction systems: (a) CzH3 C2H2 and (b) n-CJ-Is C2H2. The energy units are kcal/mol.
n-C4H3 n-C,H3
(22)
We note that for a = - 1, eq 22 is reduced to eq 14, or stated differently, eq 14 is but one particular case of the more general interpolation formula. As will be shown below, eq 21 exhibits a very high degree of accuracy in representing the results of detailed calculations. Taking these facts together, that eq 14
-
+ C2H2
k(s,2)
phenyl
+ C2H2-benzyne + H li(W
(24c) (244
The rate coefficients of these reactions and the corresponding low- and high-pressure limits were calculated from eqs 3, 4, 10, and 11, respectively. Computations were performed with an in-house RRKM code2*were the sums of energy states for rovibrational degrees of freedom were calculated by using the direct-count algorithm of Beyer and Swinehart.@ Active
Parametrization of Chemically-Activated Reactions
J. Phys. Chem., Vol. 98, No. 41, 1994 10601
109
f
x Id
'@
I
100'
I -4
-2
0
2
4
Log p, Figure 3. Rate coefficients and corresponding low- and high-pressure limits (in cm3 mol-' s-l) calculated for the reaction channels of the C2H3 C ~ H system Z in argon at 1500 K.
+
rotations were taken into account following the method of Astholz et aL41 Densities of states were calculated by the Whitten-Rabinovitch a p p r o x i m a t i ~ n . ~The ~ . ~RRKM ~ parameters used in the present study were those reported in ref 29. All the calculations were carried out with argon as a bath gas over a temperature range of 500-2500 K.
Results and Discussion Falloff Behavior. The calculated rate coefficients of reactions 23 and 24 and the corresponding low- and high-pressure limits, each case depicted at a different temperature, are presented in Figures 3 and 4. For the C2H3 CzH2 reaction system, illustrated for a temperature of 1500 K, the rate coefficients of all the channels exhibit a smooth falloff behavior between the two pressure limits (Figure 3). A different situation was observed for the n-C4H3 CzHz reaction system, whose results are shown in Figure 4 for a low temperature, 500 K: the rate coefficients of reactions 24a and 24b demonstrate a pronounced bending, the rate coefficient of reaction 24c has a plateau at the center of the falloff region, and the rate coefficient of reaction 24d shows a distinct linear part in the falloff. Such irregular behavior is nonetheless consistent with the detailed physical model, as explained next. The true low-pressure regime, the one described by eqs 10a and 1lb, is attained when for every i (i = 1, ...,n>the collisional stabilization of the activated complex, T,, is slow compared to the rest of its reaction channels. Likewise, the true high-pressure regime, expressed by eqs 10b and l l b , is achieved when the collisional stabilization is relatively fast, again for all activated complexes. The transition from the low- to high-pressure regime is smooth when the kinetic dominance of the collisional stabilization at each isomerization step i occurs roughly at the same time. When, however, the individual high-pressure kinetic
+
+
Figure 4. Rate coefficients and corresponding low- and high-pressure limits (in cm3 mol-' s-l) calculated for the reaction channels of the n-C4H3 CzHz system in argon at 500 K.
+
regimes are attained at markedly separate conditions, the falloff transition may become nonsmooth. And this is exactly what happens in the case depicted in Figure 4. From the energy diagram given in Figure 2b, it can be seen that phenyl is located in a deep potential well. This causes its collisional stabilization to be practically irreversible and dominant at low temperatures. As a consequence, with an increase in pressure the stabilization of phenyl quickly enters the highpressure regime. On the other hand, the stabilization of n-C6H5 is still in the low-pressure regime; Le., Plol is low compared to k2 (the notations follow those given in Figure 1). This causes the appearance of the irregular intermediate regimes observed in Figure 4. With the further increase in pressure, the stabilization of n-CsH5 also enters the high-pressure regime, at which point the whole reaction system goes to the true high-pressure limit. At higher temperatures, the reverse of isomerization and the thermal dissociation of Tz,energized phenyl, compete with the collisional stabilization of Tz. As a result of this, the highpressure regime for the stabilization of phenyl is reached at higher pressures and the falloff becomes smooth. The complexity of the kinetic competition also results in the temperature dependence of the low- and high-pressure rate coefficients to exhibit strong non-Arrhenius behavior. Although the modified Arrhenius expression-AT exp(-BIT) with A , n, and 8 used as fitting parameters-provides adequate fits over a temperature range 500-2500 K in general, in some cases the maximum relative error exceeded 50% (e.g., for the highpressure limits of reactions 22b and 23b) or even 100% (e.g., for the high-pressure limit of reaction 22d). Thus, further modifications of the Arrhenius formula are needed to cover wide temperature ranges for such complex systems. One way of doing it is to introduce additional terms (B'lp,B"/p,etc.) into the exponential function.
Kazakov et al.
10602 J. Phys. Chem., Vol. 98, No. 41, 1994 1.o
2ot”r
0.8
10
0.6 0.4
0.2 0.0 1.o 0.8
0.6
.. . ...
0.4
..
... .....,
. .....
-20
0.2
4
-
2
0
2
4
4
-
2
0
2
4
0.0
Figure 5. Normalized broadening factor calculated for the reaction channels of the C a s + C2H2 system in argon at 1500 K: (a) reaction 23a; (b) reaction 23b; (c) reaction 23c.; (d) reaction 23d. The line markings are (0)RRKh4 calculations, (- - -) eq 18, (---) eq 19, and (-) eq 20. The falloff broadening center Fc was calculated from eq 16.
TABLE 1: Maximum Relative Errors Obtained in Fitting the Rate Coefficients of the C&IJ C& Reactions with Different Approximation Formulas maximum relative error, % broadening-factor generalized approach, eq 15, with interpolation betwn different X formulas pressure limits T,K eq18” eq 19” eq20” eq206 eq22‘ eqs22+26d
+
C2H3
500 1500 2500
3.3 9.3 12.7
1.9 6.6 5.8 C&I3
500 1500 2500
3.8 5.9 4.3
2.7 4.4 2.9 C&
1500 2500
21.6 20.4 23.6
11.1 15.5 13.2
500 1500 2500
19.7 21.5 22.2
11.7 12.5 10.3
500
C2H3
+ CZHZ-n-C& 0.5
1.2 2.8
1.3 2.7
+CzHz-Ca Kd. 1)
0.7 1.5 1.1
0.6 1.4 1.0
3-u -0.4
-0.8
I
0.1 ‘+--a
0.0 -0.1
2.4 6.2 8.4
+H
+ C2H2 -i-C&
0.0
8
Ir(S,l)
0.6
Log pr Figure 6. Percent deviations for the rate coefficients of the C2H3 CzHz reactions in argon at 1500 K obtained with different parametrization formulas. The marking of the lines are the same as in Figure 5 , with an additional line (- - -) corresponding to eq 22 with parameter a treated as a function of temperature only.
+
Log Pr
0.1 0.2 0.5
2.7 5.1 2.7
0.1 0.2 0.2
8.5 20.9 18.7
0.3 0.9 1.4
-0.2
0
1.2
k(s,2)
8.3 15.9 18.3
7.0 13.8 15.2
+ CzH2 E
L C&
8.2 12.0 9.3
7.5 10.0 7.7
+H
7.3 12.1 9.1
0.2 0.7 0.5
Fcis calculated by eq 16. Fc is fitted along with other adjustable parameters. Parameter a is a function of temperature only. Parameter a is a function of P,according to eq 26.
Quality of Fit. The normalized broadening factors, X, computed for the C2H3 C2H2 reactions from eqs 13- 16 with different parametrization formulas, are presented in Figure 5. All the curves have a characteristic bell shape and hence the use of empirical eqs 18-20 seems to be justifiable. The deviations in the predicted rate coefficient values obtained for this reaction system by using different parametrization formulas are summarized in Table 1 for three temperatures and plotted in Figure 6 for one of them, 1500 K. Examination of the results obtained indicates that all approximations tested for reactions 23a and 23b provide reasonable quality fits, similar to those obtained for unimolecular reactions,19with the Gaussian broadening, eq 20, giving the best accuracy. However, for reactions 23c and 23d the deviations between the RRKM and fitted X
+
Figure 7. Variations of Fc,a,and u as functions of temperature for the three-parameter Gaussian broadening, eqs 15 and 20: 0, reaction 23a; 0, reaction 23b; 0 , reaction 23c; W, reaction 23d. are substantially larger, with no approach appearing better than the other. This happens primarily because of lower Fc values, which causes the error in the predicted rate coefficient to be more sensitive to the error in X . This can be seen by the application of the propagation-of-error formula to eq 15, which gives AWk = Iln F,lAX. In addition, the broadening in the case of reactions 23c and 23d is quite asymmetric and becomes more pronounced as temperature increases. The values of F,,a,and 0 in eq 20 obtained for the C2H3 C2H2 reaction channels are shown in Figure 7. All these parameters appear to be fairly smooth functions of temperature and can be readily approximated by polynomials. Note also the much lower Fc values obtained for the isomerization channels, reactions 23c and 23d, than for those of the chemicalactivation channels, reactions 23a and 23b. The results for the n-C4H3 C2H2 reaction system are summarized in Table 2 and demonstrated in Figure 8 for a low temperature of 500 K. At the low temperatures, the quality of
+
+
J. Phys. Chem., Vol. 98, No. 41, 1994 10603
Parametrization of Chemically-Activated Reactions
TABLE 2: Maximum Relative Errors Obtained in Fitting the Rate Coefficients of the n-C4H3 C2H2 Reactions with Different Approximation Formulas maximum relative error, % broadening-factor generalized approach, eq 15, interpolation betwn pressure limits with different X formulas T, K eq 18" eq 19" eq 20" eq206 eq 22' eqs 22 + 26d
+
n-CdH3
500
376.0
1500
7.2
2500
6.2
1500 2500
244.1 4.6 3.7
45.7 3.7 2.4
76.5 19.6 19.1
500 1500 2500
61.7 20.6 19.9
36.3 12.5 12.4
178.4 5.9 5.4
6.7 0.1 0.2
130.7 4.0 2.9
5.6 0.1 0.1
k(d, 1)
36.7 1.1 0.8
6.2 1.0 0.7
+ CzH2 -phenyl k(s,2)
18.0 12.9 16.2
40.7 13.5 15.6 n-CdH3
8.2 1.2 1.3
+ C2H2 -C& + H
n-C4H3
500 1500 2500
k(s.1)
44.9 1.3 1.4
60.4 5.4 4.4 n-C4H3
500
+ C2H2 -n-C&
+ C2Hz
15.3 11.4 14.6
81.0 12.4 16.1
k(d 2)
18.7 12.1 12.2
5.4 0.5
0.8
+
benzyne H 15.3 62.9 10.3 11.5 10.6 11.0
5.2 0.5
0.9
F, is calculated by eq 16. Fcis fitted along with other adjustable parameters. Parameter a is a function of temperature only. Parameter d is a function of P, according to eq 26. (I
4
-
2
0
2
4
4
-
2
0
2
4
Log pr Figure 8. Normalized broadening factor calculated for the reaction channels of the n-C4H3 + C2H2 system in argon at 500 K: (a) reaction 24a; (b) reaction 24b; (c) reaction 24c; (d) reaction 24d. The line markings are the same as in Figure 5 . The falloff broadening center F, was calculated from eq 16. fit in this case is much worse than in the previous one. The presence of the intermediate falloff regime increases the width of the broadening bell and shifts it away from P, = 1 (Figure 8a,b). As a consequence of this shift, the values of Fc calculated from eq 16 do not represent the center of the falloff broadening and this, in turn, causes large errors in the predicted rate coefficients (see Table 2). As temperature increases and the falloff curves become smoother, the maximum deviations of the empirical fits drop. For reactions 24c and 24d, the computed values of Fc are again low, similar to the analogous channels of the C2H3 C2H2 system. The fitting error introduced by the shift in the broadening peak can be reduced if Fc is not computed from eq 16, but also fitted as a parameter, i.e., instead of fitting the normalized broadening factor X in eq 15, we fit F directly, as specified by
+
4.5
-4
-
2
0
2
4
4
-
2
0
2
4
Log P ,
Figure 9. Parameter a as a function of reduced pressure obtained by numerical solution of eq 22 for the CzH3 CzH2 reaction system at specified temperatures: (a) reaction 23a; (b) reaction 23b; (c) reaction 23c; (d) reaction 23d.
+
eq 13. In so doing, the quality of fit is improved, as demonstrated for the Gaussian form of broadening, eq 20, in Tables 1 and 2. However, even then, the absolute deviations are still rather large, up to 16% for the conditions tested. Searching for a more accurate representation of the pressure broadening, we discovered a formula which offers an impressive level of accuracy, which is discussed next. It is important to emphasize, however, that we use the absolute values of fitting errors as an indicator of the flexibility of parametrization formulas, as the deviation of 16% is in itself better than the typical quality of input data or sometimes the accuracy of the numerical model. New Interpolation Formula. We first notice that eq 22 with a single parameter, a, performs almost as good as eq 20 with three parameters (a, 0,and Fc). This may indicate that the interpolation-between-limits approach to parametrization of the pressure falloff region is more flexible than the use of the broadening factor. In the initial proposal of Gardiner,14 parameter a was assumed to be a function of temperature only, independent of pressure. We examined the behavior of a as a function of pressure by solving numerically eq 22 with the values of k(i), b(i), and k ( i )taken from the RRKM calculations. The results for the C2H3 C2H2 reaction system are shown in Figure 9. These results show that parameter a is clearly pressure dependent, with a characteristic bell shape whose peak position and width are affected by temperature. Asymptotic behavior of parameter a has finite limits (see Appendix),
+
ao(i)= a,(i) = -i-
1
(25)
where ao(i) and a=(i) are the low- and high-pressure limits, respectively, of a(i) in eq 22. We note that the limits are the same for the two pressure extremes and depend only on the sequential number of isomerization. Based on these considerations, we propose to modify the GML formula, eq 22, by introducing the dependence of a on P, in the following form
a = h exp{ -[(log P, - a)/a12} - i-'
(26)
where a , u, and h are fitted temperature-dependent parameters. Parameters a and u in eq 26 account for the width and shift of the falloff peak and parameter h serves the role of a switch function between the Gaussian form of the falloff region and the pressure limits, -i-'. The new approach was tested on the same two reaction systems as considered above. The maximum relative deviations
Kazakov et al.
10604 J. Phys. Chem., Vol. 98, No. 41, 1994
-
TABLE 3: Maximum Relative Errors Obtained in Fitting the Rate Coefficient of Reaction C a s C2H.4 -t H with Different Approximation Formulas maximum relative error, % RRKM calculations RRKM calculations with weak-collision model, eq 8 by solution of master eqs T, K eqs 15 20" eqs 22 + 26b eqs 15 20" eqs 22 26' 500 0.3 0.3 0.9 0.6 lo00 1.5 0.7 1.7 2.1 2000 2.4 0.6 2500 3.9 0.6 F, is fitted together with other adjustable parameters. Parameter a is a function of P, according to eq 26.
0.5 0.4 0.3
+
I
1.o"""""'"'"'
Ib
0.0
tj -1
.o
4.5
t3 3.5
1000
1500
2000
2500
T (K) Figure 10. Variations of h, a, and u in eq 26 as functions of temperature for the C2H3 C&Iz reaction system. The conditions and line markings are the same as in Figure 7.
+
obtained for the fitting with eqs 22 and 26 are reported in the last column of Tables 1 and 2. Inspection of these results indicates that eq 26 provides a significant improvement in the quality of fit by comparison to other approaches. The maximum deviations are typically below 1% and do not exceed 7% for the conditions tested. Figure 10 displays the parameters of eq 26 determined for the C2H3 CzH2 reaction system at different temperatures. These results demonstrate that parameters a, 0, and h are fairly smooth functions of temperature and can be adequately represented by polynomials. For example, for the test cases considered in the present study, fifth- to sixth-order polynomials produce results essentially indistinguishable from the original data and second-order polynomials maintain the relative deviations typically within 10-20%. Unimolecular Reactions. In light of the impressive performance of the new parametrization formula, eqs 22 and 26, as applied to chemically-activated systems, we decided to test it also on unimolecular reactions. For this purpose we chose the reaction
+
C,H,
-
C,H,
+H
+
the same conditions but using the @-approximation for the collisional efficiency, Le., those reported in the second and third columns of Table 3. These large errors appear to be the result of computational difficulties associated with the numerical solution of master equations. The results at 1000 K already indicate the dominance of these computational errors. At higher temperatures, the scatter in numerical results becomes very large, reaching 40% at 2500 K, making the parametrization analysis meaningless.
5.5
2.5 500
+
(27)
whose parametrization was examined in the previous study.lg The input RRKM data for this reaction were taken to be the same as in ref 19. Compared in Table 3 are the maximum relative errors obtained in fitting the RRKM values with two parametrization formulas: Gaussian broadening, eqs 15 and 20, and the modified GML formula, eqs 22 and 26. The second and third columns of Table 3 compare the results computed with the RRKM method described above, Le., using the approximate treatment of the collision efficiency, eqs 8 and 9. Inspection of these two sets of data indicates that the new parametrization approach, eqs 22 and 26, performs as well as in the case of chemical activation. The last two columns of Table 3 compare the results obtained from solution of master equations, for which purpose a modifiedu code of Gilbert and Smith45 was employed. Although the absolute errors reached in this case are still relatively small, their values are larger than those computed for exactly
Summary Parametrization of theoretical rate coefficients for chemicallyactivated reactions involving isomerization of activated complexes was examined. The pressure dependence of the lowand high-pressure-limit rate coefficients was found to be determined by the sequential number, i, of the corresponding activated complex in the chain of isomerizations and the type of reaction channel, stabilization or dissociation of the activated complex. Specifically, ko is proportional to DVl] and k,to [MI1-' for stabilization channels and ko is proportional to [MIo and k, to [MI' for dissociation channels of the ith activated complex. Several parametrization approaches were tested to describe the pressure falloff behavior for two selected test systems within a wide range of temperatures, 500-2500 K. The formulas based on a falloff broadening factor, eq 13, and the generalized interpolation between limits, eq 22, gave generally large errors. The accuracy could be improved by fitting entity Fc in the former approach along with other adjustable parameters, or by treating parameter a in eq 22 as a function of pressure. The latter approach appears to be particularly attractive for its low number of adjustable parameters and demonstrated a high level of accuracy, with maximum errors below 1% for most cases tested, including both chemically-activated and unimolecular reactions.
Acknowledgment. We thank Len Yu for help with the master equation solution. The computations were performed by using the facilities of the Pennsylvania State University Center for Academic Computing. The work was supported by Gas Research Institute, under Grant No. 5092-260-2454, and the Air Force Office of Scientific Research, Air Force Systems Command, USAF, under Grant No. F49620-94-1-0226. Appendix Asymptotic Behavior of Eqs 3 and 4 at the Pressure Extremes. We f i s t note that the pressure dependence of rate coefficients ks,, given by eq 7 can be factored out as (28) ks,' = @pi= P,w:[M] We substitute eq 28 into eq 5 and then consider the limits of Zj,j = 1, ..., n, at [MI 0. Upon substitution of the obtained
-
Parametrization of Chemically-Activated Reactions
J. Phys. Chem., Vol. 98, No. 41, 1994 10605
expressions of Zj,o,j = 1, ..., n into eqs 3 and 4, we obtain the low-pressure asymptotic expressions, eqs 10a and 1la, respectively. For the analysis of the high-pressure limit, we replace [MI in eq 5 by its inverse,
p = 1/[M] (29) and expand product Zj appearing in eqs 3 and 4 into a Taylor expansion with respect to p around p = 0. The first nonzero term of this expansion is
nJ=l
Substitution of eq 30 into eqs 3 and 4 gives the high-pressure asymptotic expressions, eqs 10b and 1lb, respectively. Asymptotic Behavior of Parameter a at the Pressure Extremes. Let us first consider the low-pressure limit. For this purpose, we rewrite eq 22 as
g
k,(i)
= [I
+
(G) ]
a ‘la
Parameter a is always negative for the Lindemann-type dependence and therefore the second term on the right vanishes as [MI 0. Hence, we can approximate eq 31 as
-
Invoking eqs 10, 11, and 17, we rewrite eq 32 as (33) Equation 33 describes the asymptotic behavior of k(i) defined by the interpolation formula, eq 20. This behavior can be also examined by expanding Zj in eqs 3 and 4 into a Taylor expansion with respect to [MI around [MI = 0. The first two terms of this expansion give
nj,l
(34) where function f has finite values and can be shown to approximate i{io(i)/k,(i)}l’i. Equating powers of [MI in eqs 33 and 34, we obtain am]-+^ = -i-1 as reported in eq 25. The asymptotic result at high pressures, a[M]-, = -i-l, is obtained in a similar manner: rewriting eq 22 as
%[
= 1
k,(i)
+
(G)] a
lla
(35)
expanding Zj into a Taylor expansion with respect to p around p = 0 and using the first two nonzero terms for evaluation of eqs 3 and 4, and then equating the powers of [MI of the two expressions.
References and Notes (1) Robinson, P. J.; Holbrook, K. A. Unimolecular Reactions; WileyInterscience: New York, 1972.
(2) Forst, W. Theory of UnimolecularReactions; Academic Press: New York, 1973. (3) Gilbert, R. G.; Smith, S. C. Theory of Unimolecular and Recombination Reactions; Blackwell Scientific: Oxford and Cambridge, Massachusetts, 1990. (4) Quack, M.; Troe, J. Ber. Bunsenges. Phys. Chem. 1974, 78, 240. ( 5 ) Quack, M.; Troe, J. Ber. Bunsenges. Phys. Chem. 1975, 79, 140. (6) Quack, M.; Troe, J. Ber. Bunsenges. Phys. Chem. 1975, 79, 469. (7) Quack, M.; Troe, J. Ber. Bunsenges. Phys. Chem. 1977, 81, 329. ( 8 ) Dean, A. M. J . Phys. Chem. 1985, 89, 4600. (9) Westmoreland, P. R.; Howard, J. B.; Longwell, J. P.; Dean, A. M. AIChE J . 1986, 32, 1971. (10) Troe, J. J . Chem. Phys. 1977, 66, 4745. (11) Troe, J. J. Phys. Chem. 1979, 83, 114. (12) Troe, J. Ber. Bunsenges. Phys. Chem. 1983, 87, 161. (13) Gilbert, R. G.; Luther, K.; Troe, J. Ber. Bunsenges. Phys. Chem. 1983, 87, 169. (14) Gardiner, W. C. In The 12th IMACS World Congress on Scientific Computations;Vichnevetsky, R., Borne, P., Vignes, J., Eds.; IMACS: Paris, 1988; p 58. (15) Oref, I. I. Phys. Chem. 1989, 93, 3465. (16) Pawlowska, Z.; Oref, 1. J . Phys. Chem. 1990, 94, 567. (17) Pawlowska, Z.; Gardiner, W. C.; Oref, I. J . Phys. Chem. 1993, 97, 5024. (18) Robertson, S. H.; Snider, N.; Wardlaw, D. M. J . Phys. Chem. 1993, 97, 7556. (19) Wang, H.; Frenklach, M. Chem. Phys. Lett. 1993, 205, 271. (20) Poole, J. S.; Gilbert, R. G. Int. J . Chem. Kinet. 1994, 26, 273. (21) Smith, S. C.; McEwan, M. J. J . Chem. Phys. 1989, 90, 4265. (22) Jackson, R. L. Chem. Phys. 1991, 157, 315. (23) Barker, J. R. J . Phys. Chem. 1992, 96, 7361. (24) Golden, D. M.; Larson, C. W. Twentieth Symposium (International) on Combustion; The Combustion Institute: Pittsburgh, PA, 1984; p 595. (25) Larson, C. W.; Stewart, P. H.; Golden, D. M. Int. J. Chem. Kinet. 1988, 20, 27. (26) Stewart, P. H.; Larson, C. W.; Golden, D. M. Combust. Flame 1989, 75, 25. (27) Diau, E. W.; Yu, T.; Wagner, M. A. G.; Lin, M. C. J . Phys. Chem. 1994, 98, 4034. (28) Wang, H. Detailed Kinetic Modeling of Soot Particle Formation in Laminar Premixed Hydrocarbon Flames. Ph.D. Thesis, The Pennsylvania State University, 1992. (29) Wang, H.; Frenklach, M. J . Phys. Chem., submitted. (30) Dean, A. M.; Westmoreland, P. R. Znt. J . Chem. Kinet. 1987, 19, 207. (31) Westmoreland, P. R.; Dean, A. M.; Howard, J. B.; Longwell, J. P. J . Phys. Chem. 1989, 93, 8171. (32) Ritter, E. R.; Bozzelli, J. W. J . Phys. Chem. 1990, 94, 2493. (33) Westmoreland, P. R. Combust. Sci. Technol. 1992, 82, 151. (34) Bozzelli, J. W.; Dean, A. M. I . Phys. Chem. 1993, 97, 4427. (35) Abramowitz, M.; Stegun, I. A., Eds. Handbook of Mathematical Functions; Dover: New York, 1972. (36) Kom, G. A,; Kom, T. M. Mathematical Handbook; McGraw-Hill: New York, 1968. (37) Johnson, N. L.; Leone, F. C. Statistics and Experimental Design in Engineering and the Physical Sciences; Wiley: New York, 1977. (38) Churchill, S. W.; Usagi, R. AIChE J. 1972, 18, 1121. (39) Churchill, S. W.; Usagi, R. Ind. Eng. Chem. Fundam. 1974, 13, 39. (40) Beyer, T.; Swinehart, D. F. Comm. Assoc. Comput. Machines 1973, 16, 379. (41) Astholz, D. C.; Troe, J.; Wieters, W. J . Chem. Phys. 1979, 70, 5107. (42) Whitten, G. 2.; Rabinovitch, B. S. J . Chem. Phys. 1963,38, 2466. (43) Whitten, G . 2.; Rabinovitch, B. S. J . Chem. Phys. 1964, 41, 1883. (44) Yu, C.-L., personal communication. (45) Gilbert, R. G.; Smith, S. C. Falloff calculation of falloff curves for unimolecular and recombination reactions through solution of the master equation; A Quantum Chemistry Program Exchange; Indiana University: Bloomington, IN; QCPE Program No. 460; 1986.
