Volume 19 - Number 5 May 1980 Pages 333-41 0 International Edition in English Computer Simulation of the Kinetics of Complicated Gas Phase Reactions By Klaus H. Ebert, Hanns J. Ederer, and Gunther Isbarn**] Dedicated to Professor Matthias Seefelder on the occasion of his 60th birthday Modern digital methods and powerful computers make it possible to simulate the time behavior of chemical reactions. These calculations can be performed on systems containing an almost unlimited number of elementary reactions. Generally, however, the reaction models used should contain only those elementary reactions which describe the bulk of the conversion. Such a reaction model may be obtained by reduction of the complete set of elementary reactions. Another possibility is analysis of the chemical system starting from conditions ensuring a simple chemistry, which is generally the case at low temperatures and low conversions. The reaction model may then be extended into the range of the reaction variables (temperature, time) of interest. Mathematical simulations may be helpful during the development of the reaction model, and sometimes even decisive. These methods were applied to the pyrolysis of ethylbenzene and n-hexane, and to CO oxidation. They yield information on the reaction paths, the importance of particular elementary reactions, and reaction stability. Furthermore, quantitative data can be obtained concerning the influence of single elementary reactions on the product distribution. The sensitivity matrix shows, e.g., whether the determination of kinetic parameters of an elementary reaction from kinetic data of the overall reaction is possible in principle, and how high the accuracy of the rate constants should be for simulation of the reaction. Both results are important for modeling chemical reactions. 1. Introduction The rapid development of modem computers has markedly advanced the possibilities and the potency of chemical kinetics. Until recently the time behavior of only those chemical reactions could be determined exactly whose rate laws could be expressed as a system of differential equations r] Prof. Dr. K . H. Ebert, Dr. H. J. Ederer, DipLChem. G . Isbam Institut f i r Angewandte Physikalische Chemie and Sonderforschungshereich123 of the University Im Neuenheimer Feld 253, D-6900 Heidelberg (Germany) Angew. Chem. I n f . Ed. Engl. 19. 333-343 (1980) which can be integrated analytically. This is possible with elementary chemical reactions, which cannot be subdivided, and with a few simple reaction systems. It is not possible to obtain exact analytical expressions for the time dependence of most of the chemical systems of practical interest. Therefore the solution of the simultaneous differential equations requires other methods. For reactions involving radicals-and that is the majority of all homogeneous gas phase reactions-the best known approximation method for solving such systems of differential equations is the assumption of a steady state"'. This approxi- 0 Verlag Chemie. CmbH, 6940 Weinheim, 1980 0570-0833/80/0505-0333$02.50/0 333 mation implies that the reactive intermediates are consumed as fast as they are produced, resulting in a zero net formation rate. This greatly simplifies the mathematical treatment (see Section 3.2). Another method of solving complicated systems of differential equations is numerical integration, in which the time behavior is evaluated by stepwise computation. If there is a large difference in the rates of production i -d consumption of the various species, which is always the cL? in the initial period of radical reactions, very short time intervals have to be used for the calculations. In classical integration methods (e.g. Runge-Kutta this will entail considerable computational effort. Modern integration methods use variable time intervals, and are available as computer programs. Short time intervals are chosen in the area of large changes in the concentration of the species, and longer time intervals are used if these changes are smaller. Nevertheless, numerical integration of a large number of nonlinear simultaneous differential equations usually requires an enormous computational effort which can only be mastered with modern high capacity computers. Modern numerical methods tend to decrease computing time and storage size without reducing the quality of the results. Results are considered as “exact” if a reduction of the time intervals does not have any effect on the results. Before studying the differences in the results of kinetical calculations obtained by the “exact” method and by assuming a steady state, it is useful to check how well the laws of chemical kinetics generally describe the time behavior of a chemical reaction. The classical theory of the rates of thermal reaction^^^.^] is based upon two fundamental equationsf’]. The time law describes the rate of homogeneous gas phase reactions in dependence upon the partial pressures: The first term of the right side contains the components of first order reactions, and the second term those of second order reactions. The sign indicates the production (+) or the consumption (-) of the species n, in the reaction considered. The temperature dependence of the reaction rate constant can usually be expressed by the Arrhenius equation: k = A exp( - E J R T ) (2) The number of elementary reactions participating in a chemical reaction system can vary within wide limits. Examples of reactions consisting of only few elementary reactions [*] Both laws are based upon the Boltzmann distribution of molecular energies, which is satisfied within the usual reaction conditions of most chemical reactions. Chemical reactions are undergone predominantly by the energy-rich molecules of the ensemble; therefore, weak interactions, which are responsible for the numerous deviations of thermodynamic properties from “ideality”, are of little or no influence. In general, eqs. ( 1 ) and (2) are satisfied well. This means that reaction rates can be calculated with high accuracy in a wide range of the reaction variables if the values of the kinetic constants are known exactly. The accuracy of the calculations is usually very good relative to the experimental methods available for determining reaction rates. This is true without limitation for the elementary reactions which cannot be further subdivided and whose behavior is described by the stoichiometric equation. It should be noted that at low pressure energy transfer may influence the reaction rates, and that the assumption of the Boltzmann distribution may loose its validity in very fast reactions-at high temperatures or in reactions with a very low activation energy. Both possibilities cause deviations from eqs. ( I ) and (2) 334 (e.g. formation of hydrogen halides from the elements) can be found in any textbook of physical chemistry. Most of the gas phase reactions of practical interest consist of many-often several hundred-elementary reactions. Of course, they will all participate to some extent in the reaction, but the question to be answered is: which reactions determine the kinetic behavior and which exert such a slight influence that they may be neglected? Reaction mechanisms with a large number of elementary reactions almost exceed the capabilities of even the most powerful computers. As a rule of thumb, the computational effort increases with the square of the number of differential equations. It is therefore very important that complicated reaction mechanisms be reduced to a reaction model which describes the reaction with sufficient accuracy while still being suitable for kinetic calculations using either the “exact” method or the steady state assumption. 2. Reaction Models Reaction models should be made up of elementary reactions because of their general applicability. With a knowledge of the kinetic constants, absolute values of the reaction rates can be calculated. The number of elementary reactions increases with the number of atoms present in the reactants. For the fairly simple decomposition reactions of such small molecules as butane or pentane there exist more than 500 “possible” elementary reactions, even with the limiting assumption that no products are formed with a molecular weight exceeding twice that of the parent compound15]. On addition of other reactants, such as oxygen, water, eic., the number of elementary reactions is multiplied. However, a reaction model fulfills its purpose if it describes well enough the course of the reaction within the interesting range of the reaction variables (p, T ) , or if it provides answers to specific questions. It may happen that the production rate of a particular by-product is of more importance than an exact description of the formation of the major products. Of course, the lower the number of elementary reactions, the simpler becomes the mathematical treatment. However, an upper limit for this number cannot be given in general because the complexity of a reaction model also depends on the number of the participating species. However, if the number of elementary reactions approaches 50, a reduction should be considered. It is usually not easy to decide from the kinetic constants whether a particular elementary reaction belongs to the group of those which are involved in an “important” reaction path within the conversion, and therefore cannot be neglected. Furthermore, a reaction model has to be self-consistent and to contain reactions accounting for the consumption of all radicals. A useful approach to model reduction is the combination of certain species or elementary reactions; this is called lumping. However, the “kinetic constants” of the lumped reactions rarely remain constant, and the applicability of the corresponding rate expression is consequently limited. Lumping is favored if similar species (e.g. radicals from unsaturated hydrocarbons) form the same major product (e.g. soot)’6-91. In minimizing the number of elementary reactions for a reaction model the chemist usually adopts a different methAngew. Chem. Int. Ed. Engl. 19, 333-343 (1980) od. He tries to analyze the reaction mechanism using experimental data, and in so doing he first looks for the initial reactions and then proceeds to formulate those elementary reactions describing the fate of the reactive intermediates. In this way, the construction of a reaction model is continued by considering known chemical principles and attempting to formulate the formation of the products and intermediates with a minimum set of elementary reactions. Additional investigations about particular elementary are sometimes useful. Simultaneous kinetic calculations and quantitative comparisons with the experimentally determined product distribution may also be helpful in formulating the reaction model. It is convenient to start the development of a reaction model at low temperatures and low conversions, where the number of elementary reactions is small, and to extend it into the interesting regions of the reaction variables. This method is then always applicable if the model contains only few initial reactions, and if these form only few primary intermediates, as in the thermal decomposition of ethylbenzene (see Section 4.3). Furthermore, this method allows a clear description of reaction behavior and prediction of variations in the reaction behavior caused by variations of the reaction conditions or by the addition of various reactants. 3. Mathematical Methods 3.1. Numerical Treatment of Rate Equations For numerical treatment, the reaction model is transformed into a system of differential equations containing one differential equation for each species taking part in the reaction. In general this is not difficult. If the number of species is greater than ten, the system should be evaluated by computer. The transformed reaction model consists of linear and quadratic terms depending on the unimolecular or bimolecular nature of the corresponding elementary reactions. (In radical recombination reactions, third order terms may occur too.) To obtain a numerical solution of the differential equation system the values of the kinetic parameters (activation energy, pre-exponential factor) of each elementary reaction are needed. The number of elementary reactions whose kinetic constants are known with satisfactory accuracy is still small; such data have been compiled in tables”3 “I. If the kinetic constants are not known, their values should be estimated from chemical analogy or from thermodynamic considerations[‘‘ “1. The kinetic constants of reactions of the same type often vary within close limits. These “band widths” are compiled Table 1 Types of reactions in the pyrolysis of hydrocarbons and band widths of the kinetic parameters A [ s I: Is ‘ mol ‘1 (pre-exponential factor) and E., [kJ mol ‘1 (activation energy); M.M’: stable molecules; R.R’: radicals. Type no. Elementary reaction h o A E, -R+R’ M + R +R’+M’ R -M+R’ R+R’ + M(+M’) 15-17 7-1 1 12--15 7-1 1 10-12 200-360 30- 50 1CK--180 M R - R‘ Angew. Chem. Inr. Ed. Engl. 19, 333-343 (1980) 0 50- 80 in Table 1 for five types of reactions, which represent most of the elementary reactions occurring in the pyrolysis of hydrocarbons. For these estimates, it is important to have information about the sensitivity of the reaction model to variation of the values of respective reaction constants. More detailed data on the sensitivity of a particular model and its temperature dependence can be obtained only by simultaneous computations (cf. Section 5.2). Occasionally in simple systems quantitative data are available on the sensitivity of reaction constants, e.g. of the activation energy of unimolecular initial reactions. In complicated systems it is mostly very difficult, if at all possible, to predict the effect of the variation of a particular kinetic constant on a model. For some elementary reactions it is evident, however, that they do not affect the reaction model noticeably, even if the rates are varied considerably. Large differences in the sensitivities appear in all systems of differential equations in which the time constants differ by orders of magnitude; such systems are called ‘stiff. The well-known Runge-Kutta methods[21are not suitable for the numerical integration of stiff differential systems. Because of the huge difference in the time behavior of the individual elementary reactions the time step size has to be set extremely small in order to obtain correct results. Errors introduced by excessive step sizes in the initial period of the reactions are propagated also in regions of higher conversions. The computer time required to cover high conversions of a reaction is then extremely long. Fast methods for the solution of stiff differential equation systems have become an important field o f research in numerical mathematics. All modern “stiff solvers” contain an efficient self-adjusting optimum step size control, i. e. rapid variations require small step sizes, slow variations permit large step sizes[r91. A further important point is the algorithm for the numerical stability of a stiff solver. The algorithm has to prevent the accumulation of computational errors. The program of Gearizo],which is very often used by chemists, is based on a “multistep predictor-corrector’’ algorithm which also controls the step size. The program has been very successful; until now no case of a “wrong” solution of a chemical problem is known. Other methods for solving stiff differential equation systems are the implicit Runge-Kutta method of Wanner[Z’] and which requires a the semi-implicit algorithm of DeuflhardZZ] considerably smaller computer storage size than the program of Gear. 3.2. The Steady State Assumption The most important and most frequently applied approximative method in treating differential equation systems in chemical kinetics is based upon the quasi-steady state assumption (QSSA). In this approximation, the right-hand sides of the differential equations for the reactive intermediates are set equal to zero in order to transform these differential equations into algebraic equations. As a result, the remaining reduced system of differential equations (for the reactants and the stable products) is substantially less stiff than the original system. This may permit either an analytical solution of the differential equation system or the appli- 335 cation of considerably simpler numerical integration methods (e.g. Runge-Kutta). The QSSA does not mean, however, that the concentrations of all reactive intermediates are set constant during the reaction. According to the mathematical treatment with the singular perturbation t h e ~ r y f ~ ~ -which ~’I, will not be outlined in detail, there exist two types of radicals[281depending mainly on the type of reactions in which the radical is consumed: type rm, whose concentrations remain approximately constant, and type r,, whose concentrations decrease according to a function similar to that governing the disappearance of the reactant in the course of the reaction. The QSSA, as usually applied by chemists, neglects the time behavior of the initial period of the reaction in which the concentrations of all radicals increase to a limiting value corresponding to the “steady-state” valuef*].Generally, the smaller the steady-state concentrations are, the lower will be the conversion at which the steady state is attained. Therefore, the QSSA is better fulfilled at low temperatures and high pressures. The question then arises whether there exists a quantitative measure for the applicability of this approximation. The differential equations of the radicals can be transformed into a set of dimensionless differential equations. This procedure and a suitable rearrangement attributes a dimensionless constant E~ to each dimensionless differential quotient. The &,-valuesare composed of reaction rate constants and initial concentrations of the reactants, and it can be shown that they are characteristic values for the applicability of the QSSA. Although theory supplies exact solutions for calculating the G-values, in most cases more or less rough estimations are more readily available which give satisfactory information on the applicability of the QSSA. For thermal decomposition reactions (e. g. the pyrolysis of ethylbenzene; cf. Table 2) the numerical values of E , for the two groups of radicals rm and r, are approximated by equations (3) and (4), respectively: Em = &= rate constant of initial reaction (rate constant of the reaction (r, + M)) x Mo (3) rate constant of initial reaction rate constant of the reaction of r. (4) (Mo= initial concentration of the educt M). In cases in which there are more than one initiation or consecutive reactions of the radicals, the sum of the rate constants of these reactions has to be used. As the number of the reactions involved in consumption of a certain radical is usually small, calculations of E~ can be easily performed if the corresponding kinetic constants are known. The values E, and E. obtained have to be multiplied by a value T, which considers that part of the reactant M reacting via the initial reaction. Finally, the magnitude of the product E; x is the quantitative measure for the applicability of QSSA[Z8.291. Comparisons of QSSA calculations with “exact” calculations show that QSSA is a good approximation if E~ x % s l O - ’ . Rough estimations can be carried out up to E~ x %= lo-’. Theoretical studies confirm the general existence of these critical valu ~ s [ ~ O ~Another . question is whether the E; x % values of all [*I Quantitative information about the initial period can be obtained, however, by a suitable transformation of the time within the set of differential equations. 336 radicals participating in the reaction have to be calculated. For the model of hydrocarbon pyrolysis only the E~ x To-values of the “important” radicals have to be considered; “important” radicals are those formed in the reactions of the starting materials and those which react with the starting materials. As an example, the comparison of QSSA data with the “exact” calculations of the ethylbenzene pyrolysis are shown in Figure 1. The conversion of the reactant and of three “important” radicals as a function of time at 950 K are shownf281. The results show that the highest gi x %-values are obtained for the CH3 radical, and that the deviations are already considerable at &, x %= lo-’. Furthermore, it is shown that disagreement is higher at low pressures. 4. Applications 4.1. Oxidation of Isobutane and Isopentane The peroxide-initiated oxidation of isobutane and isopentane at relatively low temperatures was chosen as one of the first examples of the application of modem numerical methods of integration to complicated homogeneous gas phase reactions. The reaction models contained 20 and 32 elementary reactions, re~pectively[~’]. The product distribution of both reactions showed good agreement with experimental data. Another result of the calculations was that the concentrations of some of the radicals remained nearly constant over a longer period of reaction time, confirming the existence of a steady state under these conditions. 4.2. Pyrolysis of Alkanes A fundamental study of Allara and Edelson deals with the thermal decomposition of propane, n-butane, and n-pentanel5].Very extensive reaction models were used which contain almost all the possible elementary reactions (293, 515, and 586, respectively), assuming that the products have a molecular weight not greater than twice that of the reactant. The authors adopted the philosophy “that it is far better to err on the side having too many reactions than too few”. The calculations were compared with experimental data only for the pyrolysis of propane, and then only for the formation of methane. Nethertheless, interesting information was obtained by the calculations, e. g. on the contributions of particular reaction paths to the conversion, or on the concentration-time relation of selected radicals during the reaction, and on changes in the reaction behavior due to global variations of the reaction constants. Froment et al.[321studied a similar problem. A reaction scheme consisting of 133 elementary reactions was compiled, and it was assumed that this master set is applicable to the pyrolysis of saturated and unsaturated hydrocarbons generally. In order to set up a special reaction model, a suitable selection of elementary reactions (50 to 80) was taken from the master set by “trial and error” methods. The mathematical treatment of those reaction models was camed out using Gear’s program. Comparison of the calculated results with experimental data from a pilot reactor demonstrated that product distributions of the cracking reactions of ethane, propane, n-butane, and isobutane as well as of the correAngew. Chem. Inl. Ed. Engl. 19, 333-343 (1980) 005 01 fl. 005 01 t / S R05 01 t/. 005 01 tlS O3 i .,‘lo= 2 3.10‘7 aoi 003 ---.......---.-.-.-- ................f i........“.......... / 1 ...- tls - -...-.......... b 01 - 1 -3 I I 0.2 I -I I .,,,o.~-~~-s Fig. 1 . Variation of concentration with time (time-conversion diagram) of ethylbenzene (EB). methyl, a-methylbenzyl, and hydrogen radicals at 950 K and 3 x lo6 N m - 2 (30 bar) (left row) and at 3 x to’ Nm-’ (0.03 bar) (right row) according to QSSA (. . . . .) and “exact” calculations (-----). sponding olefins could be simulated satisfactorily. Distributions of the unstable products and the temperature profile along the flow reactor were further points of information obtained from the calculations. Simulation of the decomposition reactions of mixtures of the above mentioned hydrocarbons showed that the combination of reaction models by simple addition of the respective elementary reactions gave results in good agreement with experimental data. This is to be expected as long as the elementary reactions are independent and do not influence one another. Angew. Chem. Int. Ed. Engl. 19. 333-343 (1980) 43. Pyrolysis of Ethylbenzene The thermal decomposition of ethylbenzene was thoroughly investigated experimentally and comp~tationally[~~1. In order to establish the reaction model (Table 2) the stable and unstable species formed at low pressures (lo-* mbar) and 600-1400 K were determined by mass spectr~metry~~~l. Figure 2 shows the results, containing only those masses whose signal area contributed at least 1.5% to the total signal area. (Because of the different ionization propabilities of the 337 particular species this selection is somewhat arbitrary.) U p to 1100 K the following species were detected (the figures in parenthesis refer to the possible reaction(s) of formation in Table 2): methane (3), ethylene (8-1 I), benzene (4),benzyl (I), toluene (14), styrene (7), a-methylbenzyl (3-9, and bibenzyl (15), as well as undecomposed reactant. reactions (6, 10, 11). Reactions involving a hydrogen atom attack on the aromatic ring are well known in the literature[14.15,351. The experimental styrene-benzene ratio (Fig. 3) was used to adjust the competitive reactions (7) and (S), the ethylene-styrene ratio (Fig. 4) was used for the adjustment of the competitive reactions ( 5 ) and (6). With the model of Ta- Table 2. Reaction model of the pyrolysis of ethylbenzene (computer diagram). a: Kinetic parameters evaluated by analogy. Elementary reaction I ( 1 ) ( 2 ) C6H5CH2CH3 C6H5CH2CH3 3 4 5 6 C6H5CH2CH3 C6H5CH2CH3 C6H5CH2CH3 C6H5CH2CH3 ( ( ( ( ) ) ) ) ( 7 ( 8 ( 9 (10 + * + CH3 + C6H5 + + H + H + -+ * ) ) ) ) (11 ) (12 ) C6H5CHCH3 C6H5CHCH3 C2H5 C8H11 C6H7 C4H3 + (13 ) (14 ) CH3 + CH3 H + C6H5CH2 C6H5CH2 + C6H5CH2 * (15 ) -f + +. + +. +. + C6H5CH2 + CH3 C6H5 + C2H5 (15.3 /305.4 ) (16) (16.0 1347.3 ) a CH4 + C6H5CHCH3 C6H6 + C6H5CHCH3 H2 + C6H5CHCH3 C8H11 ( 7.821 (11.0 / (10.5 / (11.0 / C6H5CHCH2 + H C2H4 + C6H5 C2H4 + H C2H4 + C6H7 C2H4 + C4H3 C4H2 + H (14.85/125.5 (14.85/182.0 ( 1 3 . 0 /125.5 ( 1 3 . 0 /105.0 ( 1 3 . 0 /105.0 ( 1 3 . 0 /105.0 C2H6 C6H5CH3 C6H5CH2CHZC6H5 (10.34/ 29.3 31.4 21.0 46.0 ) (13) ) a ) a ) a (16) ) a (16) ) a ) a 1 a 0.0 ) I 0.0 (11.9 1 0.0 (11.0 Numerical integrations were performed with the most simple reaction models developed from these results (Table 2) using Gear’s program. The kinetic constants of the elementary reactions were taken from the literature. Values of constants not known in the literature were estimated from analogous reactions. The calculated results were compared with experimental data from a laboratory flow reactor using a gas chromatograph for analysis1331. The mass balance showed that part of the ethylene is formed from cleavage reactions of the aromatic ring already at relatively low temperatures. This was accounted for by (16) ) ) a a ble 2 the experimental results were in good agreement up to 1150 K. At higher temperatures the formation of benzene and toluene was markedly higher than the calculated values. The model could only be adapted by adding a second initial reaction (2) and the consecutive reaction (9). With these the reaction model of the pyrolysis of ethylbenzene was increased to 15 elementary reactions. Figure 5 shows that the model simulates the reaction very well at relatively short reaction times and within a wide temperature range. As Figure 6 shows, the agreement of the time-conversion plots is not as good. The simulations of the . . . . . . . . . . . . . . . * . . . . . . I . . . . . . . . . . . . . ’ . * . . . . . . . . . . . . . . . . . . .* . 6 . . . . . . * .. ... m 338 1 I I I 65-C5H, 78- Benzene m r n ~ m . a m m I .... . ...... ..... 1 79- C6H, - 89-CTH5 m m.... I - LL-Propane - 50- Butadiyne I . I... & . . I 39-C3H3 .......... m - . . . . . . m - , ........... S 15- Methyl 16-Methane 26-Acetylene - 28-Ethylene 90- Phenylcarbene 91- Benzyl - 32-Toluene - 102- Phenylacetylene - IOL-Styrene - 105-a-Methylbenzyl - 106-Ethylbenzene - 128-Naphthalene - 182-Bibenzyl I Anyew. Chem Int. Ed Engl 19. 3.11-341 (19Xnl benzene and ethylene formation agree reasonably well, but the calculated conversion values of ethylbenzene and styrene disagree strongly with the experiments. A possible reason is the existence of reactions of styrene, decompositions, as well 1.00 t =r 1.90 800 900 1000 1100 K Fig. 3. Styrene/benzene ratio in dependence of the temperature ( A 0 experimental values, ----- calculated values). 1.80 0 2 4 6 Ills 8 Fig. 6. Comparison of relative molar concentrations as function of time at 950 K (see Fig. 5). as reactions with radicals. Because these reactions are very similar to those of ethylbenzene, the calculations for benzene and ethylene are not so strongly affected. 4.4. Pyrolysis of n-Hexane Fig. 4. Ethylene/styrene ratio in dependence of the temperature (see Fig. 3) Fig. 5. Comparison of relative molar concentrations as function of temperature ( 0 ethylbenzene, A styrene; D ethylene: 0 benzene; ---- calculated curves). Angew. Chem. Inr Ed. EngI. 19. 333-343 (1980) In the last few years the pyrolysis of n-hexane has been in401. Hexane is a kind of key vestigated more substance. The molecule is large enough to be considered as a model for higher paraffins, but the product distribution is not too complicated. In order to describe the reaction at relatively low temperatures and conversions a reaction model consisting of 38 elementary reactionsI2'l is satisfactory (Table 3). Calculations showed that the reaction rate determined in a static apparatus by the overall pressure rise can be simulated well enough (Fig. 7). Table 4 shows that experimental and calculated product distributions also agree well. At higher conversions, however, product reactions have to be considered, the reaction model then extends to about 100 elementary reactions. The decomposition and isomerization reactions of the C6H,3 radicals (23-28) are of particular irnportan~e'~' 4H1. Calculations showed that the ratio n-hexyl/l-methylpentyl has a decisive influence on the ratio of the products ethylene/propene, and that in the temperature range studied the rate of the isomerization of n-hexyl ---t 1-methylpentyl is considerably higher than the rate of decomposition of the radicals. Therefore the equilibrium, which is shifted in the direction of I-methylpentyl is established rapidly, and the most important decomposition reaction is (24). Application of the p-rule to the decomposition reactions of the C6Hf 3 radicals indicates a ratio of reactions (23-26) of about 1 : 14: 1.5:9. This also confirms the importance of (24). 339 Fig. 7. Pressure rise in the pyrolysis of n-hexane (. experimental values; ---- calculated curve). This reaction is explosive in a particular p T range. The reaction model used consists of 24 elementary reactions, selected mainly from chemical considerations. Most of the kinetic constants were taken from the literature (Table 5). For regions showing a nonexplosive behavior, computations were camed out for a flow reactor assuming a steady state and taking into account heat flow across the reactor wall. The concentrations of the radicals, the reaction rates of the elementary reactions, and the rate of heat formation inside the reactor were calculated. The heat flow was determined as a function of temperature, pressure, and the geometry of the reactor. The explosion limit was defined as the temperature at which the heat formed in the reaction cannot be completely removed. The calculations indicate that the “critical” temperature is highly dependent on traces of water. An increase in water content from 1 to 10 to 100 ppm lowers the critical temperature from 1100 to lo00 to 900 K. Further the calculations showed that the influence of water occurs essentially via the reaction chains (4), ( 5 ) and (4), (12), (18) (Table 5). 4.6. Reactions in the Atmosphere 4.5. Oxidation of Carbon Monoxide Pilling and Noyes have done some work on the oxidation A number of publications deal with photochemically iniof~interest are ~~~~ ]. of carbon monoxide by oxygen in the presence of ~ a t e r [ ~ ~ , ’ ~ I . tiated reactions in the a t m o ~ p h e r e ~Topics Table 3. Reaction model of the pyrolysis of n-hexanc (computer diagram). a: Kinetic parameters evaluated by analogy. Reaction Elementary reaction No. C6H I4 C6H 14 C6H I4 C6HI4 + H C6H14 + H C6H14 + H C6H14 + CH3 C6Hl4 + CH3 C6H14 + CH3 C6H14 + C2H5 C6H14 + C2H5 C6H14 + C2H5 C6H14 + C3H7 C6H14 + C3H7 C6HI4 + C3H7 C2H5 C3H7 C3H7 C4H9 C4H9 C5Hl I C5HI 1 I-C6H13 2-C6H 13 3-C6H13 3-C6H I 3 I-C6H13 2-C6H I3 CH3 + CH3 CH3 + C2H5 C2H5 + C2H5 H + H C3H7 + C3H7 C4H9 + C4H9 C5Hll + CSHII I-C6H13 + I-C6H13 2-C6H13 + 2-C6H13 3-C6H13 + 3-C6H13 340 * + A * + -f -+ + * -+ -+ + + * + + -+ -+ -f -+ -+ -+ * + + + -+ CH3 + C5Hll C2H5 + C4H9 C3H7 + C3H7 (17.2/353.4) (16.6/343.3) (16.0/340.8) H2 + I-C6H13 H2 + 2-C6H13 H2 + 3-C6H13 ( 9.0/ 42.3) a (10.3/ 39.4) (10.W 39.4) a a CH4 + I-C6H13 CH4 + 2-C6H13 CH4 + 3-C6H13 ( 8.7/ 48.1) ( 8.9/ 40.2) ( 8.6/ 40.2) (13) C2H6 + I-C6H13 C2H6 + 2-CbH13 C2H6 + 3-C6H13 ( 7.8/ 52.8) ( 8.0/ 43.5) ( 7.7/ 43.5) a a C3H8 + I-C6H13 C3H8 + 2-C6H13 C3H8 + 3-C6H13 ( 7 . 6 f 52.81 ( 7.8/ 43.5) ( 7.5/ 43.5) a C2H4 + H C3H6 + H C2H4 + CH3 C3H6 + CH3 C2H4 + C2H5 C3H6 + C2H5 C2H4 + C3H7 C2H4 + C4H9 C3H6 + C3H7 C5H10 + CH3 C4H8 + C2H5 ( I 3.5/ 170.4) (16) (13.8/159.1) (16) (13.6/138.6) (l2.l/ll3.5) (13.6/121.4) (12.5/116.8) (12.6/120.2) (12.6/120.2) (13.5/108.9) (12.6/120.2) (14.3/113.5) 2-C6HI 3 I-C6H13 (ll.O/ 57.4) ( l l . l / 69.5) C2H6 C3H8 C4H10 H2 (10.6/ (10.6/ (10.6/ Tar Tar Tar Tar Tar Tar (ll.O/ 0.0) 0.0) 0.0) 0.0) (10.6/ 0.0) (10.6/ 0.0) (10.6/ 0.0) (10.6/ 0.0) (10.6/ 0.0) (10.61 0.0) a a (41) (13) (13) a a a (16) (16) (16) ( ( 5) 5) a (15) a a ( ( 5) 5) a a a a a a a a a a Angew. Chem. Ini. Ed. Engl. 19. 333-343 (1980) by “mathematical modeling” is an important tool in chemical kinetics. In mathematics this is called the “inverse problem”. A non-linear least-squares fit in dependence of the reaction parameters is performed using possibly weighted experimental data and values calculated from the model ( W,=weighting factors): Table 4. Comparison of the product distribution (16)of the pyrolysis of n-hexane: 700 K, 1 3 x lo4 Nm (0.13 bar), 1516 conversion. ’ 1 Product exp, - H2 CH 4 C2H4 C2H6 C3H6 C3H8 C4H8 C5H10 20.83 22.18 13.76 26.64 1.39 12.95 2.24 calc. 1.42 23.51 23.53 11.18 26.78 1.59 11.45 0.46 L (5) Function F should be minimized. I Table 5 . Reaction model of CO oxidation (computer diagram). a: Kinetic parameters evaluated by analogy. ~ ~~ ~ Reaction lOgloA/E, Ref. ( 8.35 /251.2 ( 8.00 / 0.0 (51) (52) (53) Elementary reaction No. * co2 + 0 C O + O + M H20 + 0 CO + OH 02 + H OH + 0 H20 + H H2 + 0 H2 + OH OH + OH O H + H + M -+ C02 2OH C02 02+H+M H 0 2 + OH H02 + H H02 + H H02 + H H02 + 0 CO + H02 H02 + H02 H202 + OH H202 + H H202 + H H202 + 0 H202 + M -+ co + 02 * * * * * * * -+ * * * -+ * -+ * -+ * * * * * + M +H m + o 02 + H H2 + OH OH + H H20 + H H20 + 0 H20 + M H02 + M H20 + 0 2 H20 + 0 H2 + 0 2 20H 0 2 + OH CO2 + OH H202 + 0 2 H20 + H02 H2 + H02 H20 + OH H02 + OH 20H + M (10.76 / ( 8.49 / (11.35 / (10.40/ (10.92 / (10.24 / (10.36 / ( 9.76 / (11.00 / ( 9.51 / (10.70 / (10.70/ (10.40 / (11.40 / (10.70/ 4.2 2.9 8.0 4.2 ( 8.00 / 41.9 (10.00 / 4.2 (10.00 / 7.5 (10.37 / 38.5 (11.50 / 37.7 (l0.70/ 5.0 (14.07 /190.5 the influence of particular substances emitted from the earth’s surface on the ozone layer in the upper atmosphere, and the formation of smog. One of the most interesting examples was the formation of the Los Angeles smog, which was studied using a model of 92 elementary In this work, the main problem is the formulation of suitable reaction models. The calculated time dependence of the reaction systems agreed relatively well with the experimental data obtained in a smog chamber. The computations are rather complex because of fluctuations in the radiation density and because of the need to include diffusion processes and convective transport phenomena in the calculations, which are difficult to describe exactly. Though experimental results can be only partly explained by model computations, important information has been attained on the participation of particular substances in smog formation, and there is no doubt that these methods will help to increase our understanding of chemical reactions in the atmosphere and their further consequences. 5. The “Inverse Problem” 5.1. Evaluation of Kinetic Constants Besides the “simulation” of chemical reactions, the evaluation of kinetic constants of particular elementary reactions Angew. Chem. Inr. Ed. Engl. 19, 333-343 (1980) ) ) 75.4 ) 2.5 ) 70.3 ) 0.0 ) 84.2 ) 39.6 ) 21.8 ) 3.8 ) 0.0 ) -4.2 ) 4.2 ) ) ) ) ) ) ) ) ) ) ) ) (55) (54) (55) (53) (53) (55) (53) (56) (54) (57) (57) (57) (57) (57) (57) (57) (54) (54) (54) a (54) The usual mathematical methods for the solution of such optimization problems can be divided into two groups: 1) Direct search methods: For each experimental value a corresponding value has to be calculated according to a model. The set of optimal parameters is obtained by a systematic search of certain parameter sets, frequently expressed as lattice points in higher dimensions. 2) Descent methods: These methods require in addition to the function F also the partial derivatives of F with respect to the parameters. These derivatives are used to find the steepest descent to the minimum of F. Bremermann and Milstein tested a number of different methods from each A direct search method with random vectors proved to be most effective and robust. The direction indicating the variation of the parameter vector is chosen by a random number generator. The “best” parameter set in this direction is determined by a special iteration method. This new parameter set is the starting point for the next random vector, and the calculation is continued until the new random vectors do not result in further improvements of the parameter sets. This method was applied to experimental data from the Calvin photosynthesis cyclef6’].The reaction model consisted of 22 reactions with 18 species. Optimal reaction constants were calculated for all reactions. The model calculation agreed very well with the experimental time behavior of the 18 species within the course of the reaction. 341 Another method, which is often more effective although more cumbersome, is the frequently used “trial-and-error’’ approach. The “improvement step”, i. e. the selection of the new parameter set, is executed “manually” taking into account all the previous results. In the numerical optimization procedure it is the “stiff solver” subroutine which demands the greatest proportion of computer time. Therefore, the greatest potential for improving the method lies in this area. The steady-state assumption decreases the computational effort considerably. Various parameters or relations between parameters can then often be evaluated and described by simple linear regression. Thus it is hardly surprising that nearly all parameter fittings in complex reaction systems assume QSSA conditions. It may happen that in such a parameter optimization procedure the minimum of the function F can be determined, but no satisfactory consistency between the experimental and the calculated values is obtained. This is a clear indication that the chemistry of the reaction model is incorrect. rectly. This objection is fundamentally wrong, as can be clearly demonstrated by an error analysis of the individual parameters. This is best done by evaluating a sensitivity matrix, which contains the variations of the concentrations of the different chemical species in dependence of the variations of the parameters of the individual elementary reactions. Such sensitivity matrices also depend on temperature, initial conditions, and reaction time. Table 6 shows the sensitivity matrix of the model for hexane pyrolysis (Table 3) at 718 K and initial concentrations of 0.0023 mol/l within the steady-state region. The columns contain the stable compounds, and the rows contain the elementary reactions. The elements of the matrix represents the percentage change in concentration of a particular species if the reaction constant of the corresponding elementary reaction is increased by a factor of 2. (In fact, however, the coefficients were calculated with the assumption of a differential variation of the system.) Negative values indicate that the concentrations decrease. The matrix of Table 6 was calculated assuming a steady state, which remarkably reduced the computational effort. A comparison shows that in this case the agreement of the QSSA values with the “exact” values is better than 0.1% under the conditions chosen. The matrix demonstrates that there are only a few reaction constants which significantly influence the concentrations of most of the products. Most rate constants are “insensitive” and their absolute values have little influence on the product 5.2. Sensitivities of Kinetic Parameters If for a reaction mechanism one of the optimization methods yields a set of parameter values which fits the experimental results well, an often heard objection is that almost any conceivable set of data can be simulated by so many parameters (even the shape of an eIephant!), and that ail these sets of parameters will describe the measurements more or less cor- Table 6. Sensitivity matrix of the reaction model of the pyrolysis of n-hexane (see Table 3) (computer diagram) Product 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 342 H2 CH4 C2H4 C2H6 C3H6 C3H8 18 24 9 18 23 9 18 23 9 I8 24 9 17 24 9 18 23 9 18 24 9 18 23 9 35 47 18 18 23 9 35 47 18 0 0 0 0 -0 0 0 -0 -0 0 0 3 -3 6 37 -32 9 65 9 3 -3 6 36 -31 6 46 4 -4 9 53 -51 10 73 -0 4 -4 9 52 - 50 10 71 4 -4 9 53 -51 10 73 4 -4 9 53 -51 10 73 -1 -1 -0 -0 0 -0 0 2 -2 4 0 1 0 -0 -0 0 -1 0 0 0 1 -I - -2 -6 40 0 -1 -0 -0 -0 12 7 50 31 -0 -0 -0 -0 -2 34 66 -61 -4 4 -4 10 -4 6 -4 I0 1 -1 1 0 -0 -0 0 -0 0 -0 -0 0 -1 2 -1 0 1 -1 -2 -6 -40 0 -1 -0 0 0 -0 -0 -0 0 0 1 -0 -0 -0 -0 -1 0 -0 -0 0 -1 1 0 0 -0 2 -2 0 0 1 0 0 -0 0 -0 0 0 -0 -0 0 -1 3 -2 0 0 1 -1 -2 -6 -40 0 -0 -0 0 -2 -6 -40 0 -1 -0 -2 -6 -40 0 -1 -0 -0 0 0 -0 -0 -0 -0 0 -1 0 0 -0 -0 0 -1 3 -2 0 1 -1 -2 -6 - 40 0 0 -0 -0 0 2 -3 4 5 -3 8 61 27 4 -4 -89 0 -1 -0 0 -0 -0 0 C4H8 C4H10 C5H10 Tar C6H14 1 -0 -1 12 6 47 29 -0 -0 1 4 1 1 1 -0 -1 23 -3 -22 -0 0 0 1 -1 1 2 0 -0 -0 -0 0 -0 0 -1 0 0 -0 0 0 0 0 1 1 -0 0 0 0 4 -4 9 53 -51 10 73 -0 0 3 -2 4 0 1 0 -0 -0 0 8 -8 18 107 - 102 20 145 0 -0 -0 -8 8 -9 -179 0 -0 -0 0 -1 -97 98 -2 5 -4 0 -0 2 -2 -5 -12 -! -1 -2 -6 -40 0 -2 -6 -40 0 -1 -5 -12 19 1 -1 -2 -6 - 40 0 0 -2 -1 98 -0 -0 0 -0 0 -0 0 -0 0 0 0 -1 -0 0 -0 -0 0 0 -0 -0 0 -0 -0 0 -80 0 0 0 Angew. Chem. Inr. Ed. Engl. 19, 333-343 (1980) distribution. From different sensitivity matrices evaluated on the basis of different experimental parameters a new matrix can be calculated containing the variances (standard deviations) and covariances of the reaction pararneterd6O1.These variances depend again on the experimental conditions (number of measurements, initial conditions, temperature, etc.) and on the accuracy of the experimental data. Thus, if a given kinetic parameter is to be evaluated with greater accuracy from experimental data, the experimental conditions must be such that the variance becomes as small as possible. This can also be done by selecting special initial conditions, e.g. by the addition of reactants which increase the sensitivity of the elementary reaction in question. Another possibility is the selection of another reaction system in which this particular elementary reaction has u priori a high sensitivity. 6. Outlook Work on the kinetics of complicated chemical reactions with the aim of solving all the details of a given mechanism and its dynamics were considered to be more or less hopeless until recently, because no methods were available for evaluating more general quantitative conclusions from the experimental materials. Today the very opposite seems to apply. Computational methods for kinetic treatment of even complicated chemical systems have been developed to a high standard, and they will surely be further improved. As computer capacity becomes much more readily available in the future, these methods will become more generally applicable. The quality of the reaction model will then be the limiting factor. The quality of models, however, largely depends on the accuracy with which the kinetic constants of the particular elementary reactions are known. Therefore methods of experimental determination and evaluation of kinetic constants will become more important. The application of such methods will initially be focused upon more special problems such as the effect of additional substances on reaction rates and product distributions. Within this kind of work the “inverse problem” and sensitivity analysis play an important part in obtaining quantitative information. The scope of these methods will be enhanced if they can be extended to the computation of heterogeneous reactions too. For this type of reaction, absorption, desorption, and diffusion processes have to be included in the reaction scheme. This does not mean that the reaction models become more complicated. On the contrary, in catalytic processes the number of “important” elementary reactions is comparatively small. Although the methods outlined above may be viewed with skepsis, perhaps due to the complicated mathematics involved, they nevertheless supply detailed information about the mechanism and time dependence of complicated reactions and permit a fast and reliable simulation of the influences of different parameters on chemical reactions. This is a significant progress in chemical kinetics. Received: October 31. 1979 [A 315 IE] German version: Angew. Chem. 92, 331 (1980) [ l ] M. Bodensiein. Z . Phys. Chem. 85, 329 (1913). 121 E. Slieffel: Einfuhrung in die numerische Mathematik. Teuhner, Stuttgart 1970. [3] S. W. Benson: The Foundations of Chemical Kinetics. McGraw-Hill, New York 1960 Angew. Chem. Ini. Ed. Engl. 19,333-343 (1980) 141 Y N. Kondraiieu: Chemical Kinetics of Gas Reactions. Pergamon Press, Oxford 1964. [5] D. L. Allara. D. Edelson. Int. J. Chem. Kinet. 7, 479 (1975). 161 J. Wei, C. W. Kuo, Ind. Eng. Chem. Fundam. 8, 114, 124 (1969). 171 Y. Ozawa, Ind. Eng. Chem. Fundam. 12, 191 (1973). [8] S. Y Golikeri, D. Luss, Chem. Eng. Sci. 29. 845 (1974). [9] S. M. Jacob. B. Gross, S. E. Volfz. V. W. Weekman. Jr.. AIChE J. 22. 701 (1976). [lo] G. L. Pratt: Gas Kinetics. Wiley, London 1969. 1111 M. R. Mulcahy: Gas Kinetics. Nelson, London 1973. 1121 L. H. Bamford, C. F. H. Tipper: Comprehensive Chemical Kinetics. Vol. I . Elsevier, Amsterdam 1969 I131 1.A. Kerr, M. J. Parsonage: Evaluated Kinetic Data on Gas Phase Hydrogen Transfer Reactions of Methyl Radicals. Butterworth. London 1976. [t4] J. A. Kerr. M. J. Parsonage: Reactions of Atoms and Radicals with Alkenes. Alkynes, and Aromatic Compounds. Butterworth. London 1972. 1151 Y N. Kondrafieut Rate Constants of Gas Phase Reactions. NSRDS-NBS. Washington 1972. 1161 S. W. Benson, H. E. O’Neal: Kinetic Data on Gas Phase Monomolecular Reactions. NSRDS-NBS 21, Washington 1970 1171 S. W. Benson: Thermochemical Kinetics. Wiley, New York 1976. [18] S. W. Benson, Int. J. Chem. Kinet. Symp. 1, 359 (1975). (191 Y N. I. Chan. I. Birnbaum. L. Lupidus. Ind. Eng. Chem. Fundam. 17, 133 (1978). 1201 W. C. Gear: Numerical Initial Value Problems in Ordinary Differential Equations. Prentice-Hall, New York 1972. [21] B. A. Goiiwald: Preprints Discussion Meeting: Kinetics of Physicochemical Oscillations. Aachen 1979. 1221 P. Deuflhard, personal communication. [23] R. E O’Malley: Introduction to Singular Perturbation Theory. Academic Press, New York 1974. 1241 A . H. Nayfeh: Perturbation Methods. Wiley, New York 1973. (251 C. dt Prima Richard: Modern Modeling of Continuum Phenomena. Vol 16. American Mathematical Society, Providence, Rhode Island 1977. 1261 W. Schneider: Mathematische Methoden der Storungsmechanik. Vieweg, Braunschweig 1978. [27] C. E. Pearson: Handbook of Applied Mathematics. Van Nostrand Reinhold. New York 1974, Chap. 14. [28] K. H. Eberf, H. J. Ederer, G. Isbarn, Preprint No. 31. SFB 123, Universitat Heidelherg 1979. [29] E. Hessruedf, 6. Hou, I . S. A . fsaksen, Int. J . Chem. Kinet. fa.971 (1978). [3Ol F C. Hoppenstedt, personal communication. 1311 D. L. Allara, D. Edelson. K. C. Irwin, Int. J . Chem. Kinet. 4, 345 (1972). 1321 K. M. Sundaram, G. F. Fromeni, Ind. Eng. Chem Fundam. 17, 174 (1978). 1331 K. H . Eberi, H. J. Ederer, P. S. Schmrd/. ACS Symp Ser. 65, 313 (1978). [34] H. J. Ederer, Dissertation. Universitat Heidelberg 1975. I351 K. Huyermann, A . W. Preuss, H. Cg. Wagner, Ber. Bunsenges. Phys. Chem. 79, 156 (1975). 1361 F E. Frey. H . J. Hepp, Ind. Eng. Chem. 25. 441 (1933) 1371 J. Chrysochoos, W. A. Bryce, Can. J. Chem. 43, 2092 (1965). 1381 L. Szepesy. V. Illes, K. Welrher. J. Simon, Acta Chim. Acad. Sci. Hung. 78, 341 (1973). 1391 V Illes, K. Weirher, I. Pleszkals, Acta Chim. Acad. Sci. Hung 78, 357 (1973). [40] Y Illes, I . Plesrkais, L. Srepesy, Acta Chim. Acad. Sci. Hung. 7Y. 259 (1973). 141) W. Tsang, J . Phys. Chem. 76, 143 (1972). 1421 A. Kossrakofl. F. 0. Rice, J. Am. Chem. SOC.65, 590 (1943). 1431 A. S Gordon, J. R. McNesby, J. Chem. Phys. 3 f . 853 (1959). 1441 M. J. Pearson, B. S. Rabinouiich, J. Chem. Phys. 42, 1624 (1965). 1451 F. Doue, G Guiochon, J. Phys. Chem. 73. 2804 (1969). 1461 D. C. Tardy, B. S. Rabinouiich, C. W. Larson, J. Chem. Phys. 45. 1163 ( 1966). I471 V. M. Rybin, Yu. P. Yampolskii. Pet. Chem. USSR 16, 167 (1976). [481 D.A . Leathard, J. H. Purnell, Annu. Rev. Phys. Chem. 21, 197 (1970). [491 M. J. Pdling, R. M. Noyes, Int. J. Chem Kinet. I / , 821 (1979). [SO] D.Edelson, Int. J. Chem. Kinet. If, 687 (1979). I511 C. H. Yang, A . L. Berlad, J . Chem. SOC.Faraday Trans. I 70, 1661 (1974). 1521 D.L. Baulch, D.D. Drysdale, A. C. Lloyd. High Temperature Reaction Rate Data, No. 1. Leeds University 1968. [531 D.L. Baulch, D.D. Drysdale, A. C. Lloyd. High Temperature Reaction Rate Data, No. 2. Leeds University 1968. 1541 D. L. Baulch, D.D.Drysdale, A. C. Lloyd: High Temperature Reaction Rate Data, No. 3. Leeds University 1969. 1551 W. E. Wilson, J. Phys. Chem. Ref. Data 1, 535 (1972). [561 W A. Trush, Prog. React. Kinet. 3, 63 (1972). 1571 A. C. Lloyd, Int. J. Chem. Kinet. 6 , 169 (1974). I581 S. W. Benson, Int. J. Chem. Kinet. Symp. l(1975). 1591 J. Milsfein, H. J. Bremermann, J . Math. Biol. 7, 99 (1979). [60] J. Milstein, J. Appl Math 35, 479 (1978) 1611 J. A. Bassham, M. Calvin: The Path of Carbon in Photosynthesis. Prentxe Hall, Englewood Cliffs 1957. 343

1/--страниц