Article pubs.acs.org/IECR
Generalized Model of Hydrocarbons Pyrolysis Using Automated Reactions Network Generation Adam Karaba,*,†,‡ Petr Zamostny,*,† Jaromir Lederer,‡ and Zdenek Belohlav† †
Faculty of Chemical Technology, Department of Organic Technology, Institute of Chemical Technology Prague, Technicka 5, 116 28 Prague, Czech Republic ‡ Research Institute of Inorganic Chemistry, Revolucni 1521/84, 400 01 Usti nad Labem, Czech Republic ABSTRACT: A computer program for the automated generation of steam-cracking reactions network (RNG) was developed. The RNG algorithms are based on graph theory tools. Kinetic parameters of generated reaction steps are derived using the information about the structure in the neighborhood of the reaction center by the developed group contribution method. Total radical activity in a reaction mixture is assumed to be stationary along the reactor. A method for estimating this activity in dependence on reactant structure and reaction conditions was also developed. Group contribution parameters were optimized using the results of laboratory experiments involving several pure hydrocarbons (n-alkanes, i-alkanes, cycloalkanes). The laboratory pyrolysis reactor was close to isobaric with plug-flow and specified temperature profile. Results of optimized model were in good agreement with experimental data obtained on a laboratory pyrolysis reactor.
■
INTRODUCTION Steam cracking models are important tools for the industrial development, as well as for the theoretical aspects, like reaction mechanism validation, simplifying assumptions testing, and experimental data verification. Empirical models usually use formal kinetics while mechanistic models are based on the exact description of a reaction mechanism. Semiempirical models are empirical in their nature but use relationships similar to mechanistic models for describing some properties of the described object. Semimechanistic models use theoretical knowledge for modeling chemical reactions like mechanistic models, but they rely partially on formal chemical reactions to substitute certain parts of the reaction mechanism. Comparisons of different types of pyrolysis models and evaluating of their advantages has been recently published by Sedighi et al.1 Many authors developed mechanistic models of hydrocarbon pyrolysis. For example Van Damme et al.2 published a model devoted to description of C3 mixtures pyrolysis, Allara and Shaw3 recommended a set of radical reactions for description of C4 hydrocarbons cracking including kinetic parameters collected from various literature sources, Sundaram and Froment4−6 developed a model of pyrolysis of light hydrocarbons. The methods of pyrolysis kinetics research were established in this period as well as simplifying assumptions introduced to developed models such as the effect of radial temperature profile,7,8 steadystate of radicals concentration,9 and its influence on model accuracy.10 Recently, some new models of light hydrocarbons and their mixtures were developed.11−15 The present higher variability of steam cracking feedstocks forces models to be much more flexible. Design of decomposition mechanisms and manual creating of mathematical description is an extremely tedious process, which is also prone to many human errors. It is necessary to accept some simplifications in the case of higher-boiling components (e.g., neglect some parallel reactions16). The model should also roughly maintain an equal level of detail for any two © 2013 American Chemical Society
compounds that are approximately equally abundant in the reaction mixture. The only way to comply with said requirements is to generate a reaction mechanism automatically. The requirement to add the reaction mechanism of a new component can be quickly satisfied in this way. A new approach in model developmentthe automated reactions network generation (RNG)is very promising. Van Geem et al.10 described experimental data obtained during hexane pyrolysis by the RNG model. Matheu et al.17 provided a very detailed model of methane pyrolysis. Broadbelt et al. utilized this modeling approach for ethane and cyclohexane steam cracking.18 The model of pyrolysis of more complex feedstock molecules by this approach was also developed.19 Billaud et al. developed a software package for n-hexane steam cracking modeling using completed numerical engines.20 Marquaire et al. modeled hydrocarbons pyrolysis at low temperature (400 °C) and developed tools for primary mechanism automated generation.21 This group utilized CHEMKIN software as their computation engine, which allows general utilization of generated mechanisms. Kinetic parameters were determined using various correlations also grouping chemical reactions to classes by similarities. Authors exhibited simulation results on examples of n-hexane and n-tetradecane pyrolysis with a good match to experimental data. Application on thermal decomposition of ndodecane at high-temperature was also demonstrated including comparison to experimental data.15 A system utilizing automated generation of a reactions network was also successfully applied in models of other processes such as oxidative coupling of methane,22 partial oxidation of propane,23 biomass to biofuel conversion,24 or esters oxidation.25 Other papers were published where authors Special Issue: NASCRE 3 Received: Revised: Accepted: Published: 15407
March July 5, July 9, July 9,
1, 2013 2013 2013 2013
dx.doi.org/10.1021/ie4006657 | Ind. Eng. Chem. Res. 2013, 52, 15407−15416
Industrial & Engineering Chemistry Research
Article
Figure 1. Diagram of pyrolysis reactions network generating algorithm.
an inspiring source of initial knowledge for us. However, our work is aimed at using the reaction network in our standalone pyrolysis simulation software eventually, so that we used own implementation independent of any commercial software. Computer Representation of Reaction Species. In order to support automated generation of reactions networks, a suitable representation of hydrocarbon molecules structure in the computer is necessary. In this work, the molecules and radicals are represented by graphs, exploiting the possibility of using graph theory tools, methods, and libraries for most of the routine tasks performed during the RNG process. The atoms and bonds in the hydrocarbon molecule are represented by graph nodes and edges respectively. The radical center is represented by a loop (the edge connecting one and the same node). The information about the bond types is included in the adjacency matrix of the graph. The information on the chemical nature of all participating atoms is provided in a vector of symbols representing chemical elements in the same order as in the adjacency matrix. However, the “matrix-oriented” approach is used only for data storage purposes. All operations carried out with the graph represented reaction species within the reactions network, such as listing, sorting, comparing, editing, merging etc., including reaction network generation process itself, is done by an object oriented programming approach. Thus, the graphs are represented by class instances in the computer memory connected by pointers. This dynamic structure can be easily created from/transformed to a matrix representation. This object-oriented approach may cost more memory space; but, program source code is very clear for possible future modifications, and graph theory algorithms can be easily and directly implemented. Network Generation Algorithm. Each elementary reaction type, such as hydrogen abstraction, β-scission, isomerizations, recombinations, etc., is provided by two rules:
demonstrated using an RNG-based kinetic model for detailed mechanism validation,26 kinetic parameters estimation methods verification,27 or simple mixtures steam cracking.5,11 Typically, the models are focused to describe only a single component in the feedstock, usually having the kinetic parameters tuned up to describe its behavior in great detail. On the other hand, our goal was to develop a tool for fully automated generation of the reactions network for an arbitrary hydrocarbon molecule, thus developing a generalized model, capable of simulating general hydrocarbon pyrolysis in a given pyrolysis reactor. The goal was to make the model as predictive as possible, maintain a high level of generalization, and prevent any overfitting of the experimental data. Obviously, high generalization is attainable only if the number of kinetic parameters is kept as low as possible. Steam cracking of hydrocarbons is a process based on radical pyrolysis reactions, and as such, it is largely nonselective, so that pyrolysis of even a single hydrocarbon can produce a variety of products by a multitude of reactions, which would normally increase the number of parameters to the extreme. However, hydrocarbon structure variability is quite limited (branching, double bonds, rings, aromatics). Thus, a whole set of many different reactions concerning specific molecules or radicals that share the reaction mechanism can be generalized by a “template”. A complete reactions set then can be generated automatically for each component of any feedstock by using a single set of “rules”. Moreover, kinetic parameters of those reactions may be derived from the basic template by making corrections for the different structure around the reaction center.
■
AUTOMATED REACTIONS NETWORK GENERATION Many authors developed their own algorithms for automated reactions network generation,10,17−19,28,29 and their papers were 15408
dx.doi.org/10.1021/ie4006657 | Ind. Eng. Chem. Res. 2013, 52, 15407−15416
Industrial & Engineering Chemistry Research
■
Article
MATHEMATICAL MODEL The mathematical model of the laboratory pyrolysis reactor (microreactor Shimadzu PYR4A) was developed for using the reaction network above. There is an extensive set of experimental data available using the reactor.31 The laboratory reactor was modeled as an isobaric straight tube reactor in steady-state with specified temperature profile along the length coordinate (z). Carrier gas flow rate is measured at ambient conditions, and the flow rate at reaction conditions is calculated using the ideal gas equation of state. The mole balance and the temperature and flow rate profiles along the reactor are represented by the following equations:
specifying which structures can enter the reaction and changes in the molecular structure of reactants due to the reaction. Those rules above were formalized by graph algorithms. Molecules can participate in hydrogen-abstraction reactions, while radicals can undergo β-scission, radical isomerizations, or recombination. In the beginning of the process, only graphs of the reactant molecules are present in the list of reaction species. The entire procedure of reaction network generation is then started by executing a whole set of the rules for the first reactant (Figure 1). Executing the rule means that the program looks for the parts of the reaction species structure which match the rule eligibility for a reaction. Where eligible, the second part of the rule is executed, making a change to the reaction species structure, generating a new molecule or radical, which is added at the end of the reaction species queue. In order to decide, whether any generated product represents a new component of the reaction system, the resulting graphs are tested for isomorphism with all existing components. It is needed because newly produced graphs of a particular compound formed by two different reactions may not be identical due to the different atom (node) numbering produced by the RNG algorithms. Only those identified as nonisomorphic with any compound already present in the list are appended to the list of molecules and radicals. Simultaneously, the reaction object representing the rule execution is constructed and fitted into the reaction network. After all rules were executed for the first reaction species, the generation moves to the next reaction species in the queue. The generation continues until the processing queue of reaction species is empty. All relevant possibilities of decomposition for all components in the list are known after the algorithm ends, and the reaction pathways used are stored in the resulting data structure. After processing each reactant, there may be duplicate reactions in the generated network due to many intermediate products being formed by multiple reaction pathways. Multiple instances of the same reaction are lumped together, which is achieved by simply summing them up by stoichiometry. Streamlining the Network. It was necessary to develop an approach, how to evaluate significance of different elementary steps in order to simplify and streamline the overextensive reaction network. This step requires kinetic parameters, so as to calculate approximate rate of reactions. Reactions having extremely high activating energy, e.g., β-scission of double CC bond, are rejected from system by an absolute filter. Then, for each reaction species participating in more than one reaction, partial selectivity of reactions is evaluated. If the selectivity of a reaction is lower than selected threshold, such as 5%, that reaction is removed from the network. The temperature-dependent selectivities are calculated for 1050 K, which is considered being characteristic temperature of the steam-cracking process. Products formed exclusively by removed reactions are removed as well from the database. At the present state, the model does not include a mechanistic description of radical addition reactions. Instead, the generated network is supplemented by set of formal molecular reactions of C1−C4 molecules, substituting radical addition, Diels−Alder reactions, and hydrogenation reactions. The set of formal reactions including their parameters was adopted from our well-established steam-cracking model30 developed earlier. It was necessary to enable the simulation of the high-conversion experiments where secondary pyrolysis reactions are more significant.
dJi dz
=
∑ νi ,jrj j
dT = f (z ) dz
⎛ d J ⎞⎤ dV ̇ RS ⎡⎢⎜⎛ dT ⎟⎞ = Ji + T ∑ ⎜ i ⎟⎥ ∑ dz p ⎢⎣⎝ dz ⎠ i ⎝ dz ⎠⎥⎦ i
where Ji is the intensity of component molar flow per reactor cross-section, rj the rate of reaction j, and νi,j is the stoichiometric coefficient of component i in reaction j, T is the reaction temperature, p is the reaction pressure, and R is the gas constant. Reaction rate is calculated assuming all reactions being of the first-order to each reactant. Rate constants of chemical reactions are evaluated by Arrhenius equation using the actual temperature at the particular point of the reactor and using the reaction activating energy and frequency factor as parameters. Hydrogen abstraction is modeled as a pseudomonomolecular chemical reaction, assuming the concentration of active radicals in the reaction mixture is stationary and related to reaction conditions and the structure of the starting compound. Combining the generated reactions network and model of chemical reactor, equations of the mathematical model are generated. A resulting set of stiff ODEs is solved using Rosenbrock method.32,33
■
KINETIC PARAMETERS ESTIMATION Knowledge of kinetic parameters of each chemical reaction is critical for the accurate modeling of the process. Generated systems usually contain hundreds or thousands of chemical reactions, and therefore, their parameters must be estimated “on-the-fly”. A group contribution method is used for estimating kinetic parameters of a reaction, resulting from the type of reaction and structure modifiers of the reaction components. Categories and subcategories of chemical reactions were created for this purpose, including the following: • hydrogen abstraction • β-scission recognizing ◦ β-scission of C−C bonds in aliphatic chain ◦ β-scission of C−C bonds in naphthenic ring ◦ β-scission of C−H bonds • isomerizations ◦ isomerization of unsaturated (allyl-type) radicals ◦ isomerization on carbon chain between the 1 and 5 positions 15409
dx.doi.org/10.1021/ie4006657 | Ind. Eng. Chem. Res. 2013, 52, 15407−15416
Industrial & Engineering Chemistry Research
Article
Figure 2. Mechanism of cycloalkanes initiation.
◦ isomerization on carbon chain between 1 and 4 positions • recombination of nonfissile radicals (e.g., hydrogen or methyl radicals) by the hydrogen radical One general frequency factor value is assigned to each one of the (sub)categories. Estimating the activation energies for Habstraction, as well as for both types of β-scission, exploits the premise that the activation energy of bond scission must be correlated with the energy of the bond. Therefore, activation energies are estimated by the group contribution method. It reflects the bond type (C−C or C−H), the number of alkyl substituents in the bond neighborhood, and the existence and distance of multiple bond(s) in the neighborhood: • base activation energy for β-scission of C−C bond (substituted by 3 alkyls) ◦ increments for different counts of substituents1, 2, 4 ◦ increment for multiple bonds in the β position ◦ increment for multiple bonds in the α position • base activation energy for β-scission of C−H bond (on secondary carbon atom) ◦ increments for primary and tertiary carbon atoms ◦ increments for multiple bonds in the α and β position (shared with C−C bond scission) • base activation energy of H-abstraction (on secondary carbon atom) ◦ increments for primary and tertiary carbon atom ◦ increment for multiple bonds in the β position • base activation energy of isomerization reactions (allyl-type, 1−4, 1−5) • base activation energy of recombination of nonfissible radicals by hydrogen radical
Despite being controversial, in order to achieve generalization for many different hydrocarbon molecules and describing more complex mixtures in the future, we decided to apply this simplification. If all active radicals were considered as potential attackers in the H-abstraction reaction in combination to any C−H bond any molecule, it would lead to a very complex system of chemical reactions. Solving of the resulting equations system can be problematic, solving the model of complex substances or industrial mixtures would be more time-consuming, and the kinetic parameter optimization would be difficult. Using this assumption, the sum of radical activity presents a “constant” multiplier in the rate equations of H-abstraction reactions. However, this “constant” is a function of state and the reaction mixture composition. Using this approach, it is not needed to include abstraction reactions by all individual radicals, but some relationship must be established between the structure of reacting molecules and the total activity of radicals. A novel approach was developed for estimating the radical activity in this study. The sum of activity of radicals can be determined at given temperature and pressure from the equation aR • = aRREF FS FC •
where the activity of radicals during pyrolysis of referential material under given conditions (T, p) is defined as a referential activity. FS is the structural factor to characterize the structural differences between referential and given compound, and FC is the factor taking into account the influence of cyclic hydrocarbons. Activity of radicals during reference substance pyrolysis can be determined from equilibrium condition of homolytic scission of reference substance under given conditions. The equilibrium constant of this reaction can be found at some referential conditions and easily recalculated to specified conditions. The structural factor of a component in the mixture reflects the structural difference of the component from reference substance. This factor represents a ratio of probability of homolytic scission in the given component to that probability in the reference substance at some reference temperature. The factor for a molecule is calculated by dividing the sum of all reaction rates of homolytic C−C scission reactions at a reference temperature of 1050 K possible in the particular molecule by the sum of all reaction rates of homolytic C−C scission reactions in the reference molecule. The structural factor in the mixture is calculated as a weighted sum of structural factors of all components present in the mixture. This approach replaces rate equations of all initiation reactions in the model by equilibrium calculation. In cyclic hydrocarbons, the homolytic scission in the hydrocarbon cycle leads to a biradical that is favorably decomposed without forming an active radical (see Figure 2). The FC factor reflects the effect of cycle presence on the active radical concentration.
■
MODELING THE RADICALS ACTIVITY As stated above, hydrogen abstraction is considered as a pseudomonomolecular reaction. Specifically, hydrogen abstraction reactions of a hydrogen atom by all possible radicals are lumped into a single reaction expressing the abstraction by “active radical”. The active radical activity is calculated as the sum of radical concentrations weighted by the relative rate of hydrogen abstraction by the respective radicals. Radical concentrations are obtained using a pseudostationary principle, taking into account the rates of initiation and termination reactions as functions of hydrocarbon mixture composition. Therefore, the activity of radicals entering the reaction is modeled as a single variable being a function of the reaction conditions, feedstock, and equipment. This assumption greatly simplifies solving the model equations. An exact description would extremely increase the number of hydrogen abstraction reactions and made the model much more difficult to solve.3,16,19 15410
dx.doi.org/10.1021/ie4006657 | Ind. Eng. Chem. Res. 2013, 52, 15407−15416
Industrial & Engineering Chemistry Research
Article
Only one set of generalized kinetic parameters was shared for modeling the reactions of all compounds in general. The set included 47 kinetic parameters in the following structure: • frequency factors (8) ◦ hydrogen abstraction, β-scission (C−C, C−H bonds, C−C bonds in ring) ◦ recombination by hydrogen radical ◦ isomerizations (allylic, 1−4, 1−5) • activation energies (18) ◦ E* of hydrogen abstraction on secondary carbon - ΔE of hydrogen abstraction on primary and tertiary carbon - ΔE of hydrogen abstraction if multiple bond is in the α and β position ◦ E* of β-scission of C−H bond on secondary carbon - ΔE of C−H scission on primary and tertiary carbon ◦ E* of β-scission of C−C bond substituted by 3 alkyls - ΔE of C−C bond scission substituted 1−2 and 3−4 alkyles - ΔE of β-scission if multiple bond is in the α and β position ◦ E* of isomerizations (allylic, 1−4, 1−5) • parameters of radical activity estimation method (2) ◦ Reference activity of reference substance (nhexane) under reference conditions (810 °C, 4 bar) ◦ parameter in the cyclic factor FC expression • parameters of formal reactions (19) replacing radical mechanism of light components decomposition Adaptive random search method (including parameters relaxation on the level of 30%) using the LCG generator was utilized for optimization. If the optimization process had stopped in a local minimum of objective function, statistical and residual analysis were compiled for appropriate actions from operator side and the optimization process could continue. More than 2 million iterations (about 3 months of clear computation time) of optimization were done before the optimization process was stopped. The full record of obtained kinetic parameters was published recently.35
It represents the fact that the active radicals level in the reaction mixture is lower in comparison to the analogous alkane. This influences the hydrogen abstraction rate and finally the conversion under otherwise identical conditions. This factor therefore reflects the fact that only in the less preferred case when in the biradical state C−H bond fission happens, the reaction chain can propagate, which is a function of cycle size, given by the following equation. ⎛ NCC ⎞ NCC− 4 FC = ⎜b ⎟ ⎝ NC ⎠
where b is adjustable constant, NC is the count of carbon atom, and NCC is the count of carbon atoms in cycle(s).
■
PARAMETER OPTIMIZATION The proposed model contains several kinetic parameters in general form. These parameter values were originally obtained from statistical analysis of literature data3,6,19,34 describing the kinetics of involved reactions. The conversions of feedstock molecules were predicted quite well by the model, but important differences between simulated and experimentally obtained pyrolysis products selectivity were observed. Therefore, the kinetic parameters of the model were optimized. Substantial differences in published data were also identified.31 Moreover, the effect of some structure factors on kinetics parameters could not be found reliably in the literature, e.g. change of activating energy caused by double bond presence in the α or β position to the reaction center. Parameters of the model were optimized to fit the experimental data obtained in the laboratory reactor,32 including the results for compounds of substantial structure variability. Selected experimental runs used for the model optimization are summarized in Table 1. Table 1. Experimental Conditions31 Chosen for Model Optimization (Various Peak Temperature and Gas Flow in N milliliters per minute at Constant Pressure of 4 bar) substance
T = 700 °C
n-heptane 2-methylheptane 2,4-dimethylpentane 2,2-dimethylpentane cyclopentane cyclohexane methylcyclohexane 1,2-dimethylcyclohexane 1,3-dimethylcyclohexane
150 150 200 200
T = 750 °C
T = 815 °C 100 100 100
■
IMPLEMENTATION EXAMPLES Model implementation can be better illustrated using attached examples as follows. Computer Representation of Reaction Species. Structures of chemical compounds are represented by graphs.
50 50 50 100 50
The weighted sum of squared differences between simulated and experimental mass-fractions (in weight percent) of main pyrolysis products or products groups was used as an objective function for the optimization. The products and product groups were as follows: Individual light componentshydrogen, methane, ethane, ethylene, acethylene, propane, propylene, propadiene + propyne, n-butane, i-butane, 1-butene, 2-butenes, 1,3-butadiene Individual aromatesbenzene, toluene, ethylbenzene, xylenes, styrene AggregateC5−C6, C7−C12, pyrolysis oil (>C12)
Figure 3. Chemical structure representation of 2-methyl-propenyl radical. 15411
dx.doi.org/10.1021/ie4006657 | Ind. Eng. Chem. Res. 2013, 52, 15407−15416
Industrial & Engineering Chemistry Research
Article
Figure 4. Scission of 3-methylpentane-2-yl radical.
Table 2. Demonstration of Kinetic Parameters Estimation Based on the Contributory Method for an Example Reaction (See Figure 4) reaction 1 2 3 4
frequency factor
activating energy increments
A(C−C β-scission) E*(C−C β-scission base, C−C β-scission 1-subst inc) A(C−C β-scission) E*(C−C β-scission base) A(C−H β-scission) E*(C−H β-scission base, C−H β-scission primary carbon incr.) A(C−H β-scission) E*(C−H β-scission base, C−H β-scission tertiary carbon incr.)
Table 3. Values of Selected Kinetic Parametersa parameter
value
A(C−C β-scission) A(C−H β-scission) E*(C−C β-scission base3 substituents) E*(C−H β-scission basesecondary carbon) ΔE(C−C β-scission 1× substituent increment) ΔE(C−H β-scission primary carbon increment) ΔE(C−H β-scission tertiary carbon increment)
1.38 × 1019 1.16 × 1010 190.86 118.45 72.07 49.87 −36.21
Figure 5. Illustration of generated networks for different reaction speciescount of molecules □ and radicals ■ contained in networks (H = n-heptane, 2MH = 2-methylheptane, 24DMP = 2,4dimethylpentane, CH = cyclohexane, MCH = methylcyclohexane, 13DMCH = 1,3-dimethylcyclohexane).
a
The entire list is published in ref 35. Frequency factors (A) are in inverse seconds, and activation energies (E*, ΔE) are in kilojoules per mole.
For example the i-butenyl radical is represented according to Figure 3. Radical structure can be represented by incidence matrix shown in the figure. Atoms connected by a chemical bond are represented by the multiplicity of bond. A radical center is represented by a loopthe bond connecting the atom with itself, i.e. the matrix diagonal. Elementary composition is characterized by a vector C, organized in the same order as atoms are numbered in the adjacency matrix. On the basis of a convention for generation purposes, a hydrogen atom is represented by the value of 1, carbon by value of 2. However, the computer implementation is based on the dynamic structure of objects linked by pointers in the computer memory. Using this virtual representation of the graph, it is possible to quite readily formalize rules used for generation process, but it also limits possibilities of transferring generated networks to another tools such as CHEMKIN. Reaction Rules Application. As an example of a scission generation algorithm, 3-methylpentene-2-yl radical decomposition is provided. The radical center is present on carbon number 2. The algorithm progresses in the way that, first, all bonds at distance 2 from the reaction center are found with limiting condition, there are no multiple bonds on the route from the reaction center to the bond, and only carbon atoms are on the road. Then, potential possibilities of bond scission are generated, and kinetic parameters are calculated by the group contribution method, as indicated in Table 2. The full list of all kinetic parameters was published by Karaba.35 The excerpt of the data necessary for illustration is provided in Table 3.
Figure 6. Illustration of generated networks for different reaction speciescount of reaction types involved: (white) H-abstraction, (blue striped) beta-scissions, (red striped) radical isomerizations, (purple checked) recombinations, (black) others (starting substance see Figure 5).
For example, reaction no. 1 is the scission of C−C bond substituted by two alkyl-substituents. Therefore, A(C−C β-scission) is used as frequency factor, and activation energy is calculated by combination of base activation energy of C−C scission and increment for two substituents around the bond: E*(C−C β-scission base) + ΔE(C−C β-scission 1 substituent increment) as 15412
dx.doi.org/10.1021/ie4006657 | Ind. Eng. Chem. Res. 2013, 52, 15407−15416
Industrial & Engineering Chemistry Research
■
Article
RESULTS AND DISCUSSION
The model was validated by comparison of experimental and simulated data for different feedstocks’ pyrolysis at different conditions. Effect of Structure on Conversion. As stated above, the dependence of conversion on reactant structure at constant conditions is quite complex. Therefore, the first validation step is the comparison of simulated and experimental conversion achieved for different structures. The comparison for the homologous n-alkane series is shown in Figure 7. According to the figure, the influence of carbon chain length is very well predicted by the model. It is important to note that the data are consistent with expected growth of conversion by increasing molecule length, i.e., increasing the count of C−C bonds that can react by the initiation and the count of C−H bonds that can undergo hydrogen abstraction. The effect of reactant structure on conversion and product selectivity was simulated quite well qualitatively, but some quantitative deviations were observed. After optimization of kinetic parameters, it was possible to compare conversions of very different feedstocks (Figure 8). The diagram shows a very close match between the experimental and simulated conversions for all acyclic hydrocarbons, despite significant structure variability. However, cycloalkane behavior is only partially reflected by the Fc factor. In order to correct this behavior, an optimization of kinetic parameters had to be used. Figures 9 and 10 show simulated weight fractions of selected pyrolysis products on the outlet of reactor compared to the experimental values. There is very good agreement between the experiment and the simulation, but there are some outliers. Those belong to the substituted cycloalkanes pyrolyzed at high temperature. More detailed comparison of experimental and simulated data is provided for the example of a compound that is not described quite well by the model (Figure 11). Methane yield is
Figure 7. Dependence of n-alkane conversion on the count of carbon atoms under constant conditions: T = 815 °C, p = 400 kPa, F = 100 N mL/min; (●) experiment,31 () simulation.
it is written in Table 3. Another example is reaction no. 4, the rate of which is described by the A(C−H β-scission) frequency factor and E*(C−H β-scission base). Examples of Generated Networks. The generated networks for the reaction species are exemplified in Figures 5 and 6 to demonstrate the count of molecules, radicals, and individual reactions species contained. It is obvious that generated reaction networks cannot be clearly shown on available space in this paper except cases of very simple molecules such as ethane or propane. Obviously, the size of the generated network increases with the size of the molecule as well as the molecule complexity.
Figure 8. Comparison of experimental31 and simulated conversions of different feedstocks: n-heptane (H), 2-methylheptane (2MH), 2,2dimethylheptane (22DMH), 2,4-dimethylheptane (24DMH), cyclopentane (CP), cyclohexane (CH), methylcyclohexane (MCH), 1,2dimethylcyclohexane (12DMCH). Conditions: p = 400 kPa, (●) 700 °C, 150 or 200 N mL/min, (■) 750 °C, 50 or 75 N mL/min, (▲) 810 °C, 100 N mL/min (flow rates selected to achieve approximately equal conversion for all hydrocarbons). 15413
dx.doi.org/10.1021/ie4006657 | Ind. Eng. Chem. Res. 2013, 52, 15407−15416
Industrial & Engineering Chemistry Research
Article
Figure 9. Comparison experimental31 and calculated fractions of selected pyrolysis products: (●) ethylene, (▲) propylene. For conditions, see Figure 8.
Figure 10. Comparison experimental31 and calculated fractions of selected pyrolysis products: (■) methane, (⧫) sum of C4, (○) butadiene. For conditions, see Figure 8.
discussion as it illustrate the benefits of the generalized approach. If the kinetic parameters within the present kinetic scheme were optimized, using only 1,3-dimethylcyclohexane data, the model would fit the data well and simulated butadiene yields would be comparable to experimental values. However, the parameters would not simulate the other hydrocarbons well. Therefore, the generalized model reveals that there is a systematic deficiency of reaction mechanisms represented in the model that prevents correct simulation of substituted
strongly underrated in the case of 1,3-dimethylcyclohexane. Unexpectedly strong demethylation of cycloalkanes has just been noted in the literature,36 but no significant reaction pathways or mechanisms are proposed that would be specific for substituted cyclic hydrocarbons. There is actually not a problem in the model itself, but the model rather reveals deficiency in the theoretical insight to radical reactions in substituted cyclic hydrocarbons. Unexpectedly high proportion of 1,3-butadiene in products of 1,3-dimethylcyclohexane pyrolysis is another example worth 15414
dx.doi.org/10.1021/ie4006657 | Ind. Eng. Chem. Res. 2013, 52, 15407−15416
Industrial & Engineering Chemistry Research
Article
Figure 11. Comparison of experimental31 and simulated weight distribution of selected pyrolysis products of 1,3-dimethylcyclohexane: (□) experiment, (■) simulation. Conditions: p = 400 kPa, T = 700 °C, F = 50 N mL/min.
cycloalkenes cracking. It can be only revealed by the generalized approach as the model is actually able to describe the single hydrocarbon data quite well. It was proved experimentally that the yields of C4 hydrocarbons cannot be attributed to secondary reactions of formed isoprene, by pure isoprene cracking experiments. Those experiments yielded almost no C4 products. A possible explanation involves probably reverse Diels−Alder reactions of some primary pyrolysis products (nonsaturated naphthenes).
Notes
CONCLUSION The developed model based on the automated reactions network generation can be implemented for compounds or mixtures that are practically important but still keeping the size of the reactions network within manageable limits. The established relationship between structure and reactivity describes behavior of alkanes and isoalkanes very well, but incorporation of further reaction mechanisms still holds potential for improving the cyclic hydrocarbon cracking simulation. Simulated feedstock conversion and distribution of main pyrolysis products are well-predicted by the developed model for very different feedstocks at widely varying conditions. It can be concluded, that the proposed modeling technique allows for a very good generalization, and it is promising to employ it for predictive calculations. One of the most important benefits of the described approach employing a universal set of generalized kinetic parameters that are common for arbitrary hydrocarbon structure is to be seen in its ability to evaluate the experimental results and compare them quantitatively with theory. Proposed cracking mechanisms can be assessed if they are viable or not, not only by comparison of simulated and experimental data, but also by comparison with results obtained for different species containing identical structure elements.
(1) Sedighi, M.; Keyvanloo, K.; Towfighi, J. Modeling of Thermal Cracking of Heavy Liquid Hydrocarbon: Application of Kinetic Modeling, Artificial Neural Network, and Neuro-Fuzzy Models. Ind. Eng. Chem. Res. 2011, 50 (3), 1536−1547. (2) Van Damme, P. S.; Narayanan, S.; Froment, G. F. Thermal cracking of propane and propane-propylene mixtures. Pilot plant versus industrial data. AIChE J. 1975, 21 (6), 1065−1073. (3) Allara, D. L.; Shaw, R. A compilation of kinetic parameters for the thermal degradation of n-alkane molecules. J. Phys. Chem. Ref. Data 1980, 9 (3), 523−559. (4) Sundaram, K. M.; Froment, G. F. Modeling of thermal cracking kinetics II. Cracking of iso-butane, of n-butane and of mixtures ethanepropane-n-butane. Chem. Eng. Sci. 1977, 32 (6), 609−617. (5) Sundaram, K. M.; Froment, G. F. Modeling of thermal cracking kinetics. I. Thermal cracking of ethane, propane and their mixtures. Chem. Eng. Sci. 1977, 32 (6), 601−608. (6) Sundaram, K. M.; Froment, G. F. Modeling of thermal cracking kinetics. 3. Radical mechanisms for the pyrolysis of simple paraffins, olefins, and their mixtures. Ind. Eng. Chem. Fundam. 1978, 17 (3), 174−182. (7) Sundaram, K. M.; Froment, G. F. Two dimensional model for the simulation of tubular reactors for thermal cracking. Chem. Eng. Sci. 1980, 35 (1), 364−371. (8) Van Geem, K. M.; Heynderickx, G. J.; Marin, G. B. Effect of radial temperature profiles on yields in steam cracking. AIChE J. 2004, 50 (1), 173−183. (9) Sundaram, K. M.; Froment, G. F. The accuracy of the pseudosteady-state approximation for radicals in thermal cracking. Int. J. Chem. Kinet. 1978, 10 (11), 1189−1193. (10) Van Geem, K. M.; Reyniers, M.-F.; Marin, G. B.; Song, J.; Green, W. H.; Matheu, D. M. Automatic reaction network generation
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS The publication was supported by the UniCRE project, funded by the EU Structural Funds and the state budget of the Czech Republic.
■
■
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected];
[email protected]. 15415
REFERENCES
dx.doi.org/10.1021/ie4006657 | Ind. Eng. Chem. Res. 2013, 52, 15407−15416
Industrial & Engineering Chemistry Research
Article
using RMG for steam cracking of n-hexane. AIChE J. 2006, 52 (2), 718−730. (11) Sabbe, M. K.; Van Geem, K. M.; Reyniers, M. F.; Marin, G. B. First principle-based simulation of ethane steam cracking. AIChE J. 2011, 57 (2), 482−496. (12) Chakraborty, J. P.; Kunzru, D. High pressure pyrolysis of nheptane. J. Anal. Appl. Pyrol. 2009, 86 (1), 44−52. (13) Dente, M.; Faravelli, T.; Pierucci, F.; Ranzi, E.; Barendregt, S.; Valkenburg, P. Primary distribution products from the pyrolysis of cyclo-alkanes. AIChE Spring National Meeting, Atlanta, April 10−15, 2005. (14) Xu, C.; Al, S. A. S.; Wang, C.; Carstensen, H.-H.; Dean, A. M. Kinetic Modeling of Ethane Pyrolysis at High Conversion. J. Phys. Chem. A 2011, 115 (38), 10470−10490. (15) Herbinet, O.; Marquaire, P. M.; Battin-Leclerc, F.; Fournet, R. Thermal decomposition of n-dodecane: Experiments and kinetic modeling. J. Anal. Appl. Pyrol. 2007, 78 (2), 419−429. (16) Karaba, A.; Zamostny, P.; Belohlav, Z. Modelling steam-cracking kinetics by the computer-generated network of reactions. 19th International Congress of Chemical and Process Engineering (CHISA), Prague, Aug 28−Sep 1, 2010. (17) Matheu, D. M.; Dean, A. M.; Grenda, J. M.; Green, W. H., Jr. Mechanism Generation with Integrated Pressure Dependence: A New Model for Methane Pyrolysis. J. Phys. Chem. A 2003, 107 (41), 8552− 8565. (18) Broadbelt, L. J.; Stark, S. M.; Klein, M. T. Computer Generated Pyrolysis Modeling: On-the-Fly Generation of Species, Reactions, and Rates. Ind. Eng. Chem. Res. 1994, 33 (4), 790−799. (19) Depeyre, D.; Flicoteaux, C. Modeling of thermal steam cracking of n-hexadecane. Ind. Eng. Chem. Res. 1991, 30 (6), 1116−1130. (20) Billaud, F.; Elyahyaoui, K.; Baronnet, F.; Kressmann, S. Chemical kinetic modeling of n-hexane pyrolysis of ACUCHEM, CHEMKIN and MORSE software packages. Chem. Eng. Sci. 1991, 46 (11), 2941−2946. (21) Bounaceur, R.; Warth, V.; Marquaire, P.-M.; Scacchi, G.; Dominé, F.; Dessort, D.; Pradier, B.; Brevart, O. Modeling of hydrocarbons pyrolysis at low temperature. Automatic generation of free radicals mechanisms. J. Anal. Appl. Pyrol. 2002, 64 (1), 103−122. (22) Simon, Y.; Baronnet, F.; Marquaire, P.-M. Kinetic Modeling of the Oxidative Coupling of Methane. Ind. Eng. Chem. Res. 2007, 46 (7), 1914−1922. (23) Hognon, C.; Simon, Y.; Marquaire, P.-M. Hydrogen Production by Homogeneous Partial Oxidation of Propane. Energy Fuels 2012, 26 (3), 1496−1508. (24) Rangarajan, S.; Bhan, A.; Daoutidis, P. Rule-Based Generation of Thermochemical Routes to Biomass Conversion. Ind. Eng. Chem. Res. 2010, 49 (21), 10459−10470. (25) Grana, R.; Frassoldati, A.; Cuoci, A.; Faravelli, T.; Ranzi, E. A wide range kinetic modeling study of pyrolysis and oxidation of methyl butanoate and methyl decanoate. Note I: Lumped kinetic model of methyl butanoate and small methyl esters. Energy 2012, 43 (1), 124− 139. (26) Cavallotti, C.; Polino, D.; Frassoldati, A.; Ranzi, E. Analysis of Some Reaction Pathways Active during Cyclopentadiene Pyrolysis. J. Phys. Chem. A 2012, 116 (13), 3313−3324. (27) Raman, S.; Carstensen, H.-H. Tree structure for intermolecular hydrogen abstraction from hydrocarbons (C/H) and generic rate constant rules for abstraction by vinyl radical. Int. J. Chem. Kinet. 2012, 44 (5), 327−349. (28) Clymans, P. J.; Froment, G. F. Computer-generation of reaction paths and rate equations in the thermal cracking of normal and branched paraffins. Comput. Chem. Eng. 1984, 8 (2), 137−142. (29) Di, M. F. P.; Lignola, P. G. KING, a Kinetic Network Generator. Chem. Eng. Sci. 1992, 47 (9), 2713−2718. (30) Belohlav, Z.; Zamostny, P.; Herink, T. The kinetic model of thermal cracking for olefins production. Chem. Eng. Process. 2003, 42 (6), 461−473.
(31) Zamostny, P.; Belohlav, Z.; Starkbaumova, L. A multipurpose micro-pulse reactor for studying gas-phase reactions. Chem. Biochem. Eng. Q. 2007, 21 (2), 105−113. (32) Kaps, P.; Rentrop, P. Generalized Runge-Kutta methods of order four with stepsize control for stiff ordinary differential equations. Numer. Math. 1979, 33 (1), 55−68. (33) Press, W. H.; Teukolsky, S. A.; Vetterling, W. T.; Flannery, B. P. Numerical recipes in C (2nd ed.): the art of scientific computing; Cambridge University Press: Cambridge, U.K., 1992. (34) Tsang, W. Chemical kinetic data base for hydrocarbon pyrolysis. Ind. Eng. Chem. Res. 1992, 31 (1), 3−8. (35) Karaba, A.; Zamostny, P. Simplifying Complex ComputerGenerated Reactions Network to Suppress Its Stiffness. Proc. Eng. 2012, 42, 1773−1782. (36) Bajus, M.; Vesely, V. Hydrocarbon pyrolysis. VI. C8 cycloalkane fractions. Ropa Uhlie 1976, 18 (3), 245−254.
15416
dx.doi.org/10.1021/ie4006657 | Ind. Eng. Chem. Res. 2013, 52, 15407−15416
