Ind. Eng. Chem. Res. 1995,34, 3749-3760
3749
Reduced Mechanisms for Propane Pyrolysis Alison S. Tomlin,**tMichael J. Pillingt John H. Merkin? and John Brindleyo Department of Fuel and Energy, School of Chemistry, and Department of Applied Mathematics, University of Leeds, Leeds, LS2 9JT, U.K.
Neil1 Burgess Research and Technology, ICI Chemicals and Polymers Limited, P.O. Box 90, Wilton, Middlesborough, Cleveland T S 6 8JE, U.K.
Arthur Gough" Department of Chemical and Process Engineering, University of Newcastle upon Tyne, Newcastle upon Tyne N E 1 7RU, U.K.
Reduced chemical mechanisms describing the pyrolysis of pure propane in a cracking tube with a uniform firebox temperature are presented. The full chemistry is taken from the model of Dente and Ranzi and consists of 422 reactions involving 48 species. By use of techniques for the identification of redundant species and the principal component analysis of the local rate sensitivity matrix, the mechanism is reduced to 122 reactions. The concentration profiles along the cracking tube calculated using this mechanism are within 1-5% of the profiles calculated using the full chemistry. The use of overall sensitivities and approximations involving fast reversible reactions enables the selection of a further 72 redundant reactions in the mechanism. The resulting 50-step scheme reproduces the concentrations of the major products and the main reactor features, although some minor products are poorly represented. The number of species in the model is further reduced by the application of the quasi-steadystate approximation (QSSA). The calculation of the instantaneous error induced by the QSSA and the species lifetime reveals a set of 12 QSSA species, and explicit analytical expressions for a subset of these species are presented. The substitution of these expressions results in a model containing only 12 coupled differential equations. The errors induced in the concentrations of non-steady-state species by the application of the QSSA are shown to be minimal. 1. Introduction
Models describing the simulation of pyrolysis reactors need t o combine complex kinetics with important physical features such as heat and mass transfer, firebox profiles and fluid dynamic characteristics. The central part of the model is the kinetic mechanism and the related rate constants. Typically, the chemistry component can consist of several hundred reactions involving large numbers of species, and this together with the coupling of the kinetic rate equations to the physical features can lead to computationally expensive calculations. The use of full kinetic mechanisms for on-line simulations such as plant optimization is therefore rarely possible, and empirical models are most often used. However, recent developments in the field of mechanism reduction have led t o the simplification and partial automation of reduction techniques. It is now possible t o select only the most important species and reactions from any full model by using sensitivity analysis (Turiinyi, 1989, 1990; Tomlin et al., 1992) and to achieve partial lumping of the equations using techniques such as the quasi-steady-state approximation (Turanyi et al., 1993). The dimension of the kinetic equations is then greatly reduced and extensive savings in computer simulation time can be achieved when these
* Author to whom correspondence should be addressed. E-mail:
[email protected]. ' Department of Fuel and Energy. School of Chemistry. 8 Department of Applied Mathematics. I ' E-mail:
[email protected].
*
0888-5885/95/2634-3749$09.0~l0
low order equations are incorporated into the full physical model. Such savings, when combined with faster computer hardware, may lead to the possibility of using detailed kinetic mechanisms for on-line optimization. Experiments used for the validation of empirical models would not then be necessary. Another important use for reliable cracking models is in the design and development of new reactors, where general predictions of product and by-product distribution are required, and where the optimum running conditions need to be established. For reactor design it may often be possible to use considerably reduced mechanisms describing the main concentration distribution features and temperature conditions. Such mechanisms will greatly reduce the time spent carrying out the repeated calculations often used in design. A strategy for reducing mechanisms is discussed below, and a low-dimensional set of equations is found which can be used to describe the high-temperature pyrolysis of propane in a cracking tube. The strategy contains three essential steps. Firstly sensitivity analysis and secondly assumptions based on fast-reversible reactions are used to remove redundant species and reactions from the full mechanism. Finally the number of species is reduced further by using the quasi-steadystate approximation (QSSA),thus reducing the number of differential equations needed t o describe the kinetics. Techniques for the selection of QSSA species are illustrated. Although the results presented apply t o a particular physical situation and direct migration of the reduced model t o different systems is not t o be encouraged, the techniques described can be used to similar effect for other types of problems.
0 1995 American Chemical Society
3760 Ind. Eng. Chem. Res., Vol. 34, No. 11, 1995 Table 1. List of All Species; Important and Necessary Species Are Marked by a Tick and Redundant Species by a Cross species no. species species status 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 39 40 41 42 43 44 45 46 47 48
steam CH4
3 3J X
3J J J X X X
X X X X
X
3 X
X X X X
3J J J J J J J X
3
J X X
X X X X
J X X
application of such reduction techniques to a model widely applied in the petrochemical industry. No attempt has been made to modify the rate parameters in line with current kinetic databases (Tsang and Hampson, 1986; Tsang, 1987, 1988, 1991). The model is applied to cracking propane at commercial conversions in a cracking tube of internal diameter 0.05 m, with a heated length of 32.9 m in the firebox, and a final adiabatic section of 3.5 m. The cracking process was modeled as a plug flow reactor (Reynolds number > 100 000) with steam as a diluent. We assume heat transfer to the tube according to Stefan’s law, and for the tube to the bulk process fluid heat transfer, a film heat transfer coefficient was calculated from the Nusselt number to give the process temperature. All reaction was assumed to occur in the bulk gas. The concentration profiles, gas temperature, gas velocity, hydrocarbon partial pressure, and total heat input must all be reproduced satisfactorily by the model based on the reduced chemistry. The inlet conditions were a total pressure of 3 bar with a hydrocarbon partial pressure of 1.233 bar and a temperature of 873 K. Preliminary model runs showed that nonuniform firebox temperatures have little effect on predicted yields a t a given conversion. Therefore, a range of uniform firebox temperatures between 1400 and 1500 K has been covered, leading t o gas temperatures between 900 and 1140 K. Pressure drop down the tube was approximately 1bar. A simulation including the full chemistry is first carried out in order to produce a set of concentrations, temperatures, and pressures at various points along the tube. The reaction mechanism can then be investigated using this information in order to find a reduced scheme. The reduction process is begun by employing local rate sensitivity analysis at each of the selected points. The calculations were carried out a t 12 points along the cracking tube as indicated in Table 2, so that a wide range of concentration and temperature conditions were covered. In using this technique we are essentially assuming constant temperature and good mixing at selected points along the tube. The mass balance of each species is therefore described by the following set of ordinary differential equations:
2. Full Model and Reaction Features
We take our example of a “full” kinetic mechanism describing the pyrolysis of pure propane from Dente and Ranzi (1983). It consists of 422 reactions involving the 48 species listed in Table 1. They claim that the model has been tested extensively against both experimental results and results from other modeling studies, showing good agreement with both. We have chosen this scheme over other mechanisms (Allara and Edelson, 1975; Sunderam et al., 1978; Edelson and Allara, 1980) because of its acceptance by industry for modeling full scale processes to high conversions. The mechanism is reproduced in the Appendix, where the rate parameters have been taken from the Dente and Ranzi scheme. Our objective is to demonstrate the objective reduction of a detailed mechanism to generate a reduced mechanism that reproduces, within acceptable limits, the yields of major product components and temperature profiles predicted by the full model. No attempt was made t o compare the predictions with data from a real furnace, but rather to compare the resulting reduced mechanisms with the full one. Our aim is to demonstrate the
where ci is the concentration of species i, kj the rate constant for reaction j , vu the stoichiometry of species i in reaction j , and Rj the rate of reaction j , where j signifies the reaction number as shown in the Appendix. 3. Reduction Procedure 3.1. Identification of Redundant Species. The first stage of the reduction process is to select the important species and features of the reaction. The reduced model must generate values for the important species and features which are in agreement, to within a specified accuracy, with those generated by the full model. In general the important species will include the initial reactants and main products of the reaction. For non-isothermal reactions, temperature may also be an important feature. The eight important species chosen in this case were C3H8 and an arbitrary set of the more commercially significant products of the py-
Ind. Eng. Chem. Res., Vol. 34, No. 11, 1995 3751 Table 2. Overall Sensitivity Orders for Various Extents of Conversion and Gas Temperatures; A Constant Firebox Temperature of 1450 K Is Assumed data point distance down tubelm extent of convlwt % gas temp/K overall sensitivity order 1 2 3 4 5 6 7 8 9 10
1 0.0 0 875
2 5.0 2.7 974
3 10.0 14.8 1023
4 15.0 30.3 1054
5 20.0 48.8 1072
60 106 107 55 61 67 27 221 128 152
152 125 67 27 102 153 59 129 241 38
67 60 27 106 61 59 30 13 152 129
13 55 67 61 27 130 175 74 59 107
67 27 38 59 125 129 152 60 30 55
rolysis, c&,cz&,CzH6, C3H6, C4H6, C6H6, and Hz. This choice of products will of course affect the species which will eventually be regarded as redundant. The next stage is the identification of redundant species, which may be excluded from the model, and necessary species, which must be included in order to produce accurate results for the important species. The necessary species will be coupled t o the important ones through significant reactions (Turanyi, 1990). Redundant species will be those which are not coupled in any significant way to the set of important and necessary species and so can be excluded from the model with only negligible direct and indirect effects on the important concentrations. A species may be considered redundant if a change in its concentration does not greatly affect f n , the rate of production of important and necessary species. A n element of the normalized Jacobian is a measure of this effect. It represents the fractional change in the rate of production of a species n caused by a fractional change in concentration ci. For an N-membered group of important and necessary species we take the sum of squares of the normalized Jacobian elements
as a measure of the overall effect on all important and necessary species by changes in concentration of species i. On the first iteration the larger the value of Bi the greater the direct effect of a nonimportant species on the set of important ones. At this stage a threshold must be defined. This threshold will be different for each model studied, but in most cases a natural limit is found with differences of several orders of magnitude between the Bi values for one species and the next; species with Bi above this limit now become necessary species. Since important species may be coupled through reactions involving several necessary species, the necessary species identified in this way have to be admitted to the N-membered group. The procedure must therefore be iterative in order to account for such indirect effects: on each iteration new necessary species can be included in the sum. Usually convergence is achieved after not more than four or five iterations and species still remaining below the threshold are considered redundant. The iteration must be carried out at various time points in order t o identify species which remain redundant throughout the reaction period. Table 1 shows the total of 19 necessary and important species which were selected a t the 12 points tested. The selected species are indicated by a tick and the redun-
6 25.0 66.4 1097
7 29.40 78.7 1120
reaction number 59 59 129 60 60 129 152 152 74 16 106 74 16 106 30 30 62 62 198 198
8 30.1 80.2 1113
9 30.8 81.3 1107
10 31.5 82.2 1103
11 32.2 83 1099
12 32.9 83 1095
67 27 55 30 13 38 62 59 61 63
67 27 55 38 13 30 62 102 63 101
30 62 59 63 55 60 13 129 152 76
59 60 129 152 16 74 106 30 198 55
30 62 55 13 63 67 27 102 76 61
dant species by a cross. C4H6 is not redundant only because it was selected as an important product. Its concentration does not affect the concentration of the other important and necessary species. The concentrations of the 29 redundant species do not have a significant effect on the important species either directly or indirectly via necessary species. These species are therefore not considered in the next stage of the reduction where reaction importance is investigated. Steam has been identified as a redundant species and therefore has little effect on the other species through its reaction in the gas phase. It is therefore possible to remove the third-body reactions (reactions 12 and 54) containing steam from the mechanism. It is of course essential to take account of the steam in pressure drop and heat transfer calculations. In real cracking tubes steam also plays an important part in controlling coking. 3.2. Identification of Redundant Reactions. 3.2.1. Background. The second stage of the reduction is the identification of redundant reactions giving the mechanism comprising only important and necessary species. We consider the simplest and most effective technique for achieving this t o be local rate sensitivity analysis (Turanyi, 1989, 19901, which measures the effect of a perturbation in rate parameter kj on the rate of change of species i. This method has an advantage over concentration sensitivity methods in that the lognormalized rate sensitivity matrix pij depends algebraically on reaction rates and is easily computed, viz.,
(3) where fi is the rate of production of species i which can be calculated from the concentration sets according to eq 1. The computing time required to calculate this matrix is therefore less than for the concentration sensitivity matrix which requires the solution of a set of coupled differential equations, and for large mechanisms this can greatly increase the computational speed of the reduction process. In the modeling of industrial processes there are further advantages of using such a local technique. It permits the provision in a large scale simulation including physical effects such as heat and mass transfer of sets of concentrations at various times and temperatures. The fi values can then be calculated for these concentration sets, and local rate sensitivity analysis can be used as a postprocessing technique requiring no further integrations. In the present study there are 48 species and 422 rate parameters, and so a large amount of information is contained in the rate sensitivity matrix which has to
3762 Ind. Eng. Chem. Res., Vol. 34, No. 11, 1995 be interpreted. The application of principal component analysis of a submatrix of Fv will quantify the importance of reaction groups for the set of 19 important and necessary species (Vajda, 1985; Turanyi, 1989, 1990) described in the previou: section. A principal component analysis of matrix F is based on the eigenvalueeigenvector decomposition of the cross-product matrix FTF.The eigenvalues measure the significance of their respective eigenvectorsin the overall mechanism. Each eigenvector then represents a set of coupled reactions whose relative contributions are recognized by the relative size of the eigenvector components. The eigenvector elements also provide information about the connections between reactions. Defining two thresholds, for significant magnitudes of eigenvalues and eigenvector components, therefore provides an automatic way of deciding which parameter groups can be eliminated. The thresholds which have been used in the present study are typically 1.00 x for the eigenvalues and 0.02 for the eigenvector components, although a more natural threshold may exist in specific cases. The overall rate sensitivity coefficients (Turanyi, 1989)
(4)
a
100
80
.,’ .-
40
20
0 0
5
10 15 20 Length along tube/m
25
30
3.01 2.5
‘6 1.5 1.0
provide a way of ordering the overall importance of the reactions selected by the principal component analysis a t different reaction times. The optimum way to select the set of important reactions is to take the set of necessary and important species a t each reaction time and to compute the principal components when only these species are included in the objective function. The sum of the reactions at each time must then be taken which provides a reduced mechanism applicable over the whole reaction period. 3.2.2. Results from Sensitivity Analysis. The 122 reactions which were selected by the principal component analysis are indicated by an asterisk (*) in the Appendix (Table 6). The mechanism was analyzed at all 12 data points for each firebox temperature using eigenvalue and eigenvector thresholds of 1.00 x and 0.02, respectively. If a reaction appeared above these thresholds at any reaction point, it was included in the reduced model. The overall sensitivity order changes with extent of propane conversion and temperature along the tube as shown in Table 2,-althoughthere are certain reactions which have high sensitivities for most of the conditions covered. The rate data for many reactions in the full mechanism are “averaged” in relation t o single reaction classes. For the most sensitive reactions, however, small errors in an averaged rate constant can lead to large errors in the rate of change of species concentration and therefore the final species concentration profile. Although a real assessment of parameter uncertainty should be made through the calculation of concentration sensitivities, information about the overall rate sensitivity order is very useful in forming priorities about which reaction rate data should be updated. Figures 1 and 2 compare species concentrations and reactor features along the cracking tube for the full and reduced chemistry. The 122-step reaction scheme is barely distinguishable, in most cases, from the full 422step reaction scheme. Table 3 compares the final
0.5
0.0
0
5
10 15 20 Length along tube/m
25
30
Figure 1. Comparison of concentration profiles for full and reduced schemes along the cracking tube with a uniform firebox temperature of 1450 K. The full scheme is represented by a solid line, the 122-step scheme by a long dashed line, and the 50-step scheme by a short dashed line. The long dashed line is completely covered by the solid line. (a) Propane and high-yield products; (b) low-yield products.
product yields for the full and reduced 122-step scheme, showing how they remain within 1-5% of the weight percent yields for the important products, even though less than a third of the original reactions remain. We have selected a firebox temperature of 1450 K for which to present the results, although the scheme was found to apply equally well over the range of firebox temperatures tested. 3.3. Overall Sensitivities and Fast Reversible Reactions. Any further reduction in the number of reactions was carried out in a less automatic way. The strict accuracy conditions required of the 122-step reaction scheme may be relaxed. For example if we wish to use the model for reactor control or design, it is not often necessary to model the final yields to within 5%for all product species: it may be sflicient to match the accuracy of the model to experimental results for major reaction features. Firstly, it is possible to use the overall sensitivities as calculated by eq 4 to find those reactions which were selected by the principal components, but which remain low in the overall sensitivity order. Such reactions can often be removed with very little loss of information. Again it is necessary t o investigate conditions all along the tube, since a reaction may have a low overall sensitivity for one set of temperature and concentration conditions, but not for another. Secondly, information
Ind. Eng. Chem. Res., Vol. 34,No. 11, 1995 3753
b
1.251
0
5
Length along tube/m
10
15
20
25
30
25
30
Length along tube/m
-I
1100
0
5
10
15
20
25
30
0
5
10
15
20
Length along tube/m
Length along tube/m
Figure 2. Comparison of reactor conditions for full and reduced schemes along the cracking tube with a constant firebox temperature of 1450 K. The full scheme is represented by a solid line, the 122-step scheme by a long dashed line, and the 50-step scheme by a short dashed line. (a) residence time; (b)hydrocarbon partial pressure; (c) gas temperature; (d) heat input. Table 3. Comparison of Final Product Yields for Full and Reduced Schemes; The Firebox Temperature Is as in Table 2 yield, wt % species
422 scheme
122 scheme
18.42 34.80 2.10 22.42 16.13 0.83 1.73 1.59
18.19 34.31 2.09 22.51 17.00 0.80 1.66 1.57
50 scheme 19.19 34.87 3.10 22.62 14.90 0.94 2.14 1.49
from the principal components and from a rate of production analysis can help to identify fast reversible reactions, which can be eliminated in pairs from the mechanism if their overall rate is much slower than the net rate of production for the species taking part in the reaction. Such reactions can be identified in pairs from the eigenvector elements of the principal components. Using these techniques, it is possible to reduce the number of reactions in the mechanism t o 50 without losing the general concentration features. The 50 remaining reactions invol-ving18 species are indicated by a bullet in the Appendix (Table 6). The concentration profiles and main reactor features for the 50-step model are again illustrated in Figures 1 and 2, and the final product yield is represented in Table 3. The weight percent of the major products such as C2H4, C3H6, and CHI are reproduced quite well by this skeleton scheme.
Some of the more minor products such as C2H6 and C6H6 agree less closely with the full model. The final yield for C2H6 differs by 30% from the full scheme. The reason for this is partly that some of the less important species have been removed completely from the scheme and therefore percentages have altered, although the individual concentrations may not differ by as much. The residence time and heat input to the reactor almost exactly match those from the full scheme, and the gas temperature differs a t most by 10 K. The hydrocarbon partial pressure shows the most significant errors although they are still small. This stage of the reduction process is less easily automated than the sensitivity analysis stage, but it is necessary as an aid to the solution of the steady-state equations in the latter stages of the reduction. 3.4. The Quasi-Steadystate Approximation or QSSA. 3.4.1. Comments on the Applicability of the QSSA. The above methods have enabled the removal of all redundant species and reactions from the original mechanism. Any further progress requires the possible lumping of reactions and species in order t o reduce the number of differential equations needed t o model the system. One method of combining reactions is via the quasi-steady-state approximation or QSSA. The concentrations of QSSA species cqs can be expressed algebraically in terms of the other species cnqs since it is assumed that their rates of change can be decoupled from the differential equations and the right-hand sides
3764 Ind. Eng. Chem. Res., Vol. 34, No. 11, 1995
set to zero, i.e., eq 1 can be replaced with the equation set dcinqs - f (')(ci,kj) dt
Table 4. Lifetimes and Estimated Single Species Errors for Each Species at a Sample Gas Temperature of 1045 K Showing 12 QSSA Species; The Firebox Temperature Is as in Table 2 ~
--
0 = f '2'(ci,kj)
order
(6)
The number of differential equations is therefore reduced, since some have been replaced by algebraic expressions. It is then, in principle possible t o substitute for cq8 in the remaining differential equations by analytic solution of the algebraic equations, and the solution of the remaining smaller set should require less computing time. The solution of such coupled sets of algebraic equations is however made difficult by the nonlinearities present, and therefore numerical solution of the coupled sets of algebraic-differential equations is often necessary. Such steady-state calculations are an integral part of the Dente and Ranzi SPYRO model although they provide no explanation of the reasons behind its successful application. An important aspect of applying the QSSA is the choice of QSSA species. The ultimate criterion must the accuracy of the resulting model, but pursuit of this accuracy is aided by employing objective procedures for identifying QSSA species. The difference in the QSSA species concentrations a t the starting time of the application of the QSSA, calculated from eq 1 and eqs 6 and 5, is the instantaneous error of the quasi-steadystate approximation. In particular the instantaneous error induced by the application of the QSSA to a single species ACg (Turanyi et al., 1993) indicates the possible steady-state species and is calculated by the following expression:
(7)
est error for single species, Ac:/c,
species
1 2
-2.970 x -1.437 x -7.831 -4.037 -1.714 1.000 4.736 x -3.155 x -1.130 x -1.234 x -1.096 x -8.148 x 1.578 x -1.376 x 7.845 x -7.109 x -6.502 x -2.601 x
I 11 12 13 14 15 16 17 18
lo1 10'
lo-' lo-' lo-'
lo-'' lo-''
lifetime = -l/Jii
(8)
it is clear that small errors require small lifetimes or slow rates of change for a species. It is useful therefore to examine species lifetimes since this can indicate possible QSSA species. In reality the QSSA is rarely applied to a single species in isolation, and error calculations must take account of the interactions between QSSA species. It is important therefore t o use eq 7 only as a first guide to the selection of QSSA species and then to calculate A e , the error resulting from the application of the QSSA to a group of species with small single errors, using the following equation:
where i and k run over all QSSA species in the group. The calculation of the group errors therefore requires the solution of a coupled set of linear algebraic equations. If the group errors remain small, then the chosen group of species can be considered as QSSA species. In our work we have calculated the fractional instan-
x 10'
x 10'
x lo-'
10-3 x x x x x x x x x
x lo-''
x
Table 5. Estimated Group Errors for the 12 Chosen QSSA SDecies ~
~~~
order
~
est group error, AcB/c,
species
1
4.572 -4.003 -1.117 -3.623 -3.135 -3.108 -3.100 -2.911 -1.785 -1.233 -8.030 -1.285
2 3 4 5 6 7 8 9 10 11 12
x
x x x x x x x x x
lo-'
lo-' lo-'
lo-' lo-' lo-' lo-' lo-'
lo-' lo-' lo-''
x x lo-''
taneous errors given by Aci/ci
where J i k = [afi(c,k)/&kl. Since the lifetime of a species is related to the reciprocal of the diagonal Jacobian element for that species, viz,
lifetime 0.1536 0.1008 0.4312 0.2202 0.1011 0.9765 0.1112 0.1794 0.1979 0.2057 0.6769 0.5254 0.4589 0.1356 0.3339 0.1431 0.9976 0.2506
(10)
so that dimensionless quantities can be compared for each species. 3.4.2. Results from Error Calculation. The lifetimes and fractional instantaneous errors for single species Acg/ci were calculated at each point along the tube and for all species with consuming reactions. An example calculation is shown in Table 4. The results are for a gas temperature of 1045 K, but are representative of the results obtained at other temperature and concentration conditions. There are 12 species which have lifetimes which are orders of magnitude smaller than for the species above them: they are potential QSSA species since the same group has the lowest fractional instantaneous errors Ace/+ When the QSSA is applied to the group of 12 species a-C4H7, CzHs, C4H91, diallyl, u - C ~ HCzH3, ~, C4H7-4, C3H7-2, H, C3H7-1, a-C3Hb and CH3, and the fractional error calculated via eq 9, we find that the group errors remain small for most species, as shown in Table 5. The species a-CdH7 and C2H5 have the largest fractional errors, but their lifetime is short implying that the decay of instantaneous error will be fast (Turanyi et al., 1993) and the overall error small. The influence of the errors on nonQSSA species cannot be assessed using this method but is likely to be minimal and so the QSSA should be a valid approximation for this group. All but one of the chosen species (diallyl) are radical species. It would have been wrong, however, to assume a priori that all
the radicals would naturally be QSSA species and that nonradicals should not be considered. The second point t o note concerning the QSSA is that steady-state behavior of the QSSA species, Le., slow rates of change, is not necessary for the application of the QSSA, since the error will be negligible if the lifetime of the species is small in comparison to the net rate of change of species concentration. For example, the net production rate of a-C3H5 at the time point tested in Table 4 is 2.80 x approximately 5 times smaller than the production rate for C2H6, although its fractional error is 13 orders of magnitude smaller than that for C2H6. 3.4.3. Application of Q S S k The algorithm presented in the previous sections enables the selection of QSSA species. The simplest and most common application of the QSSA then solves numerically the coupled set of differential and algebraic equations which results from assigning the right-hand sides of some of the original equations to zero. This is certainly a good way to test the selection of the group of QSSA species, and could result in computational savings because of the reduction in stiffness of the problem (Turdnyi et, al., 1993). To reduce the order of the problem in a significant way, however, the steady-state equations must be solved algebraically giving explicit expressions describing the concentrations of the QSSA species in terms of the non-QSSA species, which can be substituted into the remaining differential equations. The nonlinearities present in many of the QSSA equations make this a difficult task, and often the solution has to be restricted to a subset of the QSSA species. If there are many QSSA species, the explicit solution is not found easily and can rarely be achieved by hand. Algebraic manipulation packages such as REDUCE (Hearn, 1987) and MATHEMATICA (Wolfram, 1988) can be used t o aid the solution, and both are able t o produce the resulting expressions in the form of a FORTRAN code which can be used directly with a differential equation solver. The present problem contains coupling between many of the steady-state species, so that, although in theory it is possible to solve the coupled sets of algebraic linear equations, in practice the algebra becomes very complex. There are also nonlinear equations describing the concentration of species a-C3H5 and CH3. However, there are six species whose concentrations can be described simply in terms of the other species, and so a reduction of the order of the equations by the removal of a further six variables can easily be achieved. More complex iterative techniques may be needed to solve the remaining steady-state equations. The explicit expressions for these six species are shown below:
[C,H7-1] =
A+B '60
+ '59
making the algebra complex. Although in theory it is possible to solve these coupled equations, in practice the resulting expressions are so complicated as to outweigh the computational benefits gained from removing the species from the differential equations. Truncation of
3756 Ind. Eng. Chem. Res., Vol. 34, No. 11, 1995
the steady-state expressions would be necessary to simplify the calculations, but as yet there are no formal rules governing truncation procedures (Peters and Rogg, 1992). Including only the first 6 species and substituting the above expressions results in a set of 12 coupled differential equations. Such a model is approaching the kind of order required by an on-line optimization code. For the purposes of testing the application of the QSSA we have found numerical solutions for the coupled differentiayalgebraic set of equations, including all 12 species in the approximation. The concentrations and percentage yields calculated using the QSSA agree with those calculated from the set of differential equations describing the 50-step scheme to at least 3 figures for all QSSA species. The results have not therefore been reproduced since they are indistinguishable from those already discussed in section 3.3. 4. Conclusions
The techniques of sensitivity analysis and the calculation of redundant species coefficients and principal components have enabled the original 422-step scheme to be reduced to 122 reactions. The weight percentages of the main products along the cracking tube from this reduced chemistry are within 1%of the weight percentages from the model incorporating the full mechanism for some species, and within 5% for less important products. Overall sensitivities and assumptions based on fast reversible reactions led to a further reduction in the number of reactions and t o a model consisting of only 50 elementary steps. This short model reproduces the main concentration profiles and the reactor features and could therefore be useful in reactor design where repeated calculations are often necessary. Further reduction in the number of variables has been achieved through the application of the quasi-steady-state approximation. Error calculations revealed a set of 12 QSSA species, although explicit expressions describing the concentrations of QSSA species in terms of the remaining species in the model could only be found for half of the group. Numerical tests showed that the concentrations of the non-QSSA species were not affected up to three decimal places by the application of the QSSA to these 12 species. The above approach can be used not only to reduce existing mechanisms but also in the development of new mechanisms. Dente and Ranzi point out that simplification procedures are needed in order to “avoid an initial enormous number of species”. They add that “in establishing the reactions network, account is taken of the most relevant species and reactions” though they introduce no method of classifying the species. The redundant species calculations used in the present work provide an automatic method for the classification of species into important, necessary, and redundant categories. The size of the initial model becomes unimportant since redundant terms can easily and automatically be removed and the temptation to make unjustified simplifications avoided. The simulation time of the initial mechanism will be increased by including all possible reactions, but systematic reduction methods will lead to a more widely applicable reduced model. The initial development of the model requires not only the construction of the kinetic equations but also the establishment of a set of kinetic rate parameters. Often average parameters are used which relate to single reaction classes without any emphasis being placed on whether the model is sensitive t o errors in such estima-
tions. Sensitivity analysis not only provides a reduction in the number of reactions, but also indicates which of the remaining reactions are the most sensitive to changes in the rate parameters, hence showing which parameters need to be most accurately measured. It is important to produce a reduced scheme where the original rate parameters are still embedded in the equations (rather than global steps with empirical reaction rates) t o avoid repeating the reduction process each time a kinetic parameter is updated. As much flexibility as possible must be retained in a reduced mechanism. It is important to reemphasize our aim of demonstrating the application of reduction procedures using a well-established hydrocarbon-cracking mechanism. The version of the Dente and Ranzi model used here successfully models propane pyrolysis under industrial conditions but cannot be described as a fundamental chemical model. The same rate parameters have been used for different reactions of the same type; e.g., the activation energies for abstraction by CH3 have been equated. The rate constants differ substantially from those quoted by Tsang. The success of this early implementation of the model derives not from its adherence to the best available rate parameters for its component elementary reactions, but because it has been tested, and so some extent tuned, against operational crackers. The detailed mechanistic conclusions of our study and, in particular, the identity of the 50 key reactions and their overall sensitivities should not be over-interpreted, because of the approximations inherent in the original model. The procedure demonstrated retains both the accuracies and inaccuracies of the original model, because the criterion employed to test the validity of the reduced model is agreement with the full model for important species. In terms of the overall strategy, it is important to make comparisons of the reduced scheme with the full scheme rather than with experimental results. Any improvement in agreement with experiment on reduction of the full scheme might only be the result of a balancing out of inaccuracies. The resulting reduced scheme will not, therefore, lead to a better understanding of the processes involved. The application of the QSSA, sometimes known as the continuously-varying steady-state (CVSS (Dente and Ranzi, 1983)) approach, is often employed as a numerical simplification in pyrolysis simulation. Dente and Ranzi use such a technique in their SPYRO model for radical species. Their justification of the method appears to be based on experience and testing, however, rather than on well-founded a priori principles. The criteria described above provide a method for the selection of QSSA species and gives estimations of the error induced from using such an approach. It also explains why steady-state assumptions can be used in describing species with typically nonstationary profiles (Tomlin et al., 1992). Our experience has shown that the algebraic solution of the steady-state equations using manipulation packages such as REDUCE and MATHEMATICA provides greater savings in simulation time than the numerical solution of coupled sets of Merentid-algebraic equations. Clearly problems arise for nonlinear relationships, and the truncation or approximation of the steady-state expressions may be necessary. Further developments in model reduction must include lumping, where new species are defined which are functions of the original ones. So far the only lumping
Ind. Eng. Chem. Res., Vol. 34, No. 11, 1995 3767 Table 6. The Dente and Rami Scheme for Propane Pyrolysis (SeeAppendix) log(A) E 17.0 89 000 ., *
. * * , * *
* *
* ,. * ,. *
* * *
., .,
0 ,*
* *
* ., * ., * * * *
0.*
* * * * * *
* * *
., ., ., *
* * ., * ., * ,. * 0 .* * * * * * *
* * ., *
., ., *
*
*
., *
16.9 15.8 15.4 16.2 14.8 16.0 15.5 11.7 11.7 16.7 16.0 14.0 13.5 13.5 13.6 13.6 14.1 13.6 13.4 13.9 14.0 14.0 13.7 13.9 14.5 11.5 11.6 11.3 11.6 10.9 10.9 11.0 11.5 11.4 11.7 11.6 11.3 11.3 10.8 10.1 10.9 10.8 10.6 10.9 9.6 9.6 9.0 9.4 11.7 11.8 10.1 11.5 14.9 13.9 13.0 12.2 12.2 13.7 14.3 14.3 14.3 13.0 14.0 14.1 14.1 14.2 13.7 14.0 14.0 14.0 14.3 12.4 13.0 13.4 12.7 13.3 13.0 12.9 11.7
86000 73300 73 500 72 000 59000 72000 73500 30000 30000 83000 10 000 2 000 2 000 2 000 2 000 2 000 1900 2 000 2 000 2000 2500 2000 5000 5000 12000 7600 7600 7600 6 000 6 000 6 000 5 000 5 000 6 000 5 000 5 000 7 600 13 000 21 500 12 500 21 500 21 500 21 500 21 500 12500 12500 5000 5000
* *
0,*
* *
*
* .,*
0,* 0,*
* .,*
.,* .,* .,* * *
.,*
0 ,*
*
0,*
.,* .,* .,*
14500 25000 * 12000 * 13000 * 29 700 41 000 59 000 37 500 37 500 34 000 39 000 43 000 39 500 34 000 49 000 55 000 39 000 31000 * 41000 0 , * 41 000 0 , * 34000 * 38000 * 34000 .,* 34000 .,* 38000 .,* 49000 * 35000 * 0 0 0 0
log(A)
E
12.3 13.9 11.9 13.5 13.2 13.7 11.9 12.4 12.3 13.3 16.8 13.7 13.9 11.0 11.2 10.9 11.0 10.5 15.0 14.0 14.7 14.7 14.8 14.5 14.4 14.8 14.4 14.4 14.5 14.8 14.1 14.5 14.7 14.4 14.7 14.4 14.8 14.5 14.5 14.8 14.8 14.7 14.7 12.2 12.8 13.0 12.7 12.5 13.0 12.5 12.5 12.7 13.0 12.2 12.7 12.8 12.5 12.8 12.5 12.7 12.7 12.7 13.0 13.0 12.8 12.8 12.1 12.7 12.9 12.6 12.4 12.9 12.4 12.4 12.6 12.9 12.1 12.6 12.7 12.4
0 71 000 38 000 63 000 61 500 70 000 38 000 60 000 62 000 58 000 59 000 64 000 64 000 21 500 22 500 22 500 21 500 26 500 66 700 66 700 8 000 8 000 8 000 8 000 8 000 8 000 8 000 8 000 8 000 8 000 8 000 8 000 8 000 8 000 8 000 8 000 8 000 8 000 8 000 8 000 8 000 8 000 8 000 11 900 11 900 11900 11900 11900 11900 11900 11900 11900 11900 11900 11900 11900 11900 11900 11900 11900 11900 11 900 11 900 11 900 11900 11900 7 800 7 800 7 800 7 800 7 800 7 800 7 800 7 800 7 800 7 800 7 800 7 800 7 800 7 800
3768 Ind. Eng. Chem. Res., Vol. 34,No. 11, 1995 Table 6 (Continued) log(A)
240 C3H7-2 + CH4
+
C3H8 + CH3
E
12.9 7 800 12.7 7 800 12.7 7 800 12.1 14000 12.7 14000 12.7 14000 12.6 14000 12.4 14000 12.9 14000 12.4 14000 12.4 14000 12.6 14000 12.9 14000 12.1 14000 12.6 14000 12.7 14000 12.4 14000 12.7 14000 12.4 14000 12.6 14000 12.6 14000 12.6 14000 12.9 14000 12.9 14000 12.7 14000 12.7 14000 12.6 23500 13.2 23500 13.2 23500 13.4 23500 12.9 23500 13.4 23500 12.9 23500 12.9 23500 13.1 23500 13.4 23500 12.6 23500 13.1 23500 13.2 23500 12.9 23500 13.2 23500 12.9 23500 13.1 23500 13.1 23500 13.1 23500 13.4 23500 13.4 23500 13.2 23500 13.2 23500 12.1 7800 12.7 7800 12.7 7 800 12.9 7 800 12.6 7 800 12.9 7 800 12.4 7 800 12.4 7 800 12.6 7 800 12.9 7 800 12.1 7 800 12.6 7 800 12.7 7 800 12.4 7 800 12.7 7 800 12.4 7 800 12.6 7800 12.6 7 8 0 0 * 12.6 7 8 0 0 * 12.9 7 800 12.9 7800 12.7 7800 12.7 7 800 11.9 15500 12.5 15 500
log(A)
320 C4H7-4
+ C&Hs
- C4H8-1+ CyCsH5
. .
E
12.5 15500 12.7 15500 12.4 15500 12.2 15500 12.7 15500 12.2 15500 12.4 15500 12.7 15500 11.9 15500 12.4 15500 12.5 15500 12.2 15500 12.5 15500 12.2 15500 12.4 15500 12.4 15500 12.4 15500 12.7 15500 12.7 15500 12.5 15500 12.5 15500 11.9 15500 12.5 15500 12.5 15500 12.7 15500 12.4 15500 12.2 15500 12.7 15500 12.2 15500 12.4 15500 12.7 15500 11.9 15500 12.4 15500 12.5 15500 12.2 15500 12.5 15500 12.2 15500 12.4 15500 12.4 15500 12.4 15500 12.7 15500 12.7 15500 12.5 15500 12.6 24000 13.2 24000 13.2 24000 13.4 24000 13.1 24000 12.9 24000 13.4 24000 12.9 24000 12.9 24000 13.1 24000 12.6 24000 13.1 24000 13.2 24000 12.9 24000 13.2 24000 12.9 24000 13.1 24000 13.1 24000 13.1 24000 13.4 24000 13.4 24000 13.2 24000 13.2 24000 11.9 12200 12.5 12200 12.5 12200 12.7 12200 12.4 12200 12.2 12200 12.7 12200 12.2 12200 12.2 12200 12.7 12200 11.9 12200 12.4 12200 12.5 12200 12.2 12200
Ind. Eng. Chem. Res., Vol. 34, No. 11, 1995 3759
12.5 12.7 12.4 12.2 12.7 12.2 12.2 12.4 12.7 11.9 12.4 12.5 12.2 12.5 12.2 12.4 12.4
14000 14000 14000 14000 14000 14000 14000 14000 14000 14000 14000 14000 14000 14000 14000 14000 14000
techniques which have been applied to large systems such as pyrolysis have involved lumped species which are simple summations over groups of compound classes (Liguras and Allen, 1992). This form of lumping does not however allow for differences in rate constants within these classes, and the reactivity of the lump, which changes with time, cannot be expressed by a simple rate expression. Often the same reactivity is assumed for each species in the lump and a significant amount of information is therefore lost. The conditions for exact and approximate lumping, and general methods for its application, have been investigated (Li and Rabitz, 19891, although they have not yet been applied to large reaction systems. The development of such systematic lumping techniques will be an important factor in the future of kinetic modeling of pyrolysis systems. 5. Computations
Redundant species identifications, rate sensitivities, principal components, and species lifetimes were calculated by KINAL (Turbnyi, 19901, and REDUCE was
log(A)
E
12.4 12.7 12.5 12.5 11.9 12.5 12.5 12.7 12.4 12.2 12.7 12.2 12.2 12.4 12.7 11.9 12.4 12.5 12.2 12.5 12.2 12.4 12.4 12.4 12.7 12.7 12.5 12.5 12.6 13.2 13.2 13.4 13.1 12.9 13.4 12.9 12.9 13.1 13.4 12.6 13.1 13.2 13.2 12.9 13.1 13.1
14000 14000 14000 14000 14000 14000 14000 14000 14000 14000 14000 14000 14000 14000 14000 14000 14000 14000 14000 14000 14000 14000 14000 14000 14000 14000 14000 14000 32500 32500 32500 32500 32500 32500 32500 32500 32500 32500 32500 32500 32500 32500 32500 32500 32500 32500 32500 32500 32500 32500 32500
13.1
13.4 13.4 13.2 13.2
used for the explicit solution of the algebraic QSSA equations.
Acknowledgment We thank IC1 plc. for research support. We are grateful to T. Turanyi, M., Braithwaite, and I. Parker for helpful discussions.
Nomenclature ki = mth order rate constant, (molecules ~ m - ~s-l) ~ - ~ A = preexponential factor, (molecules E = activation energy, J mol-I R = gas constant, J K-l mol ci = species concentration, molecules cm-3 uy = stoichiometric number Rj = rate of reactionj, molecules cm-3 s-l t = time, s Appendix The Dente and Ranzi scheme for propane pyrolysis consisting of 422 reactions and 48 species is shown in
3760 Ind. Eng. Chem. Res., Vol. 34,No. 11, 1995
Table 6. Reactions contained in the 122-step scheme are indicated by asterisks (*) and those in the 50-step scheme by bullets (0). The rate constants for each reaction is in the form ki = AT" exp[-E/RTl with units (molecules ~ m - ~ ) s-l - - ~for an mth order reaction. In all reactions n = 0. The temperature dependencies and thermochemical data are also taken from the Dente and Ranzi scheme.
Literature Cited Allara, D. L.; Edelson, D. Znt. J . Chem. Kinet. 1975, 7, 479. Dente, M. E.; Ranzi, E. M. I n Pyrolysis: Theory and Zndustrial Practice. Albright, L. F., Crynes, B. L., Corcoran, W. H., Eds.; Academic Press: New York, 1983; p 133. Edelson, D.; Allara, D. L. Znt. J . Chem. Kinet. 1980,12,605-621. Hearn, A. C. RAND Publ. 1987, CP78,7. Li, G.; Rabitz, H. Chem. Eng. Sci. 1989,44, 1413. Li, G.; Rabitz, H. Chem. Eng. Sci. 1990,45, 977. Liguras, K. L.; Allen, D. T. Znd. Eng. Chem. Res. 1992,31, 45. Peters, N., Rogg, B., Eds. Reduced kinetic mechanisms for applications in combination systems; Springer-Verlag: Berlin, 1993. Sunderam, K. M.; Froment, G. F. Znd. Eng. Chem. Fundam. 1978, 17. 174.
Tomlin, A. S.; Pilling, M. J.; Turzinyi, T.; Merkin, J. H.; Brindley, J. Combust. Flame 1992, 91, 107. Tsang, W. J . Phys. Chem. Ref. Data 1987,16, 471. Tsang, W. J. Phys. Chem. Ref. Data 1988,27,887. Tsang, W. J . Phys. Chem. Ref. Data 1991,20,221. Tsang, W.; Hampson, R. F. J . Phys. Chem. Ref. Data 1986, 15, 1087. T u r h y i , T. J . Math. Chem. 1990, 5, 203. Turanyi, T. New J . Chem. 1990,14, 795. Turanyi, T. Comput. Chem. 1990,14, 253. Turanyi, T.; Berces, T.; Vajda, S. Znt. J . Chem. Kinet. 1989,21, 83. TurBnyi, T.; Tomlin, A. S.; Pilling, M. J. J . Phys. Chem. 1993,97, 163. Vajda, S.; Valko, P.; Turanyi, T. Znt. J . Chem. Kinet. 1985,17,55. Wolfram, S . Mathematica: a system for doing maths by computers; Addison-Wesley: Reading, MA, 1988.
Received for review December 13, 1994 Accepted May 22, 1995@ IE940735A
Abstract published in Advance A C S Abstracts, August 15, 1995. @
