Catalytic Conversion of Alkenes on Acidic Zeolites: Automated Generation of Reaction Mechanisms and Lumping Technique.

Acid-catalyzed hydrocarbon transformations are essential for industrial processes, including oligomerization, cracking, alkylation, and aromatization. However, these chemistries are extremely complex, and computational (automatic) reaction network generation is required to capture these intricacies. The approach relies on the concept that underlying mechanisms for the transformations can be described by a limited number of reaction families applied to various species, with both gaseous and protonated intermediate species tracked. Detailed reaction networks can then be tailored to each industrially relevant process for better understanding or for application in kinetic modeling, which is demonstrated here. However, we show that these networks can grow very large (thousands of species) when they are bound by typical carbon number and rank criteria, and lumping strategies are required to decrease computational expense. For acid-catalyzed hydrocarbon transformations, we propose lumping isomers based on carbon number, branch number, and ion position to reach high carbon limits while maintaining the high resolution of species. Two case studies on propene oligomerization verified the lumping technique in matching a fully detailed model as well as experimental data.

Introduction

Catalytic conversion of alkenes over acidic zeolites finds application in a wide variety of refining and petrochemical processes including oligomerization, alkylation, and aromatization.[1] Acidic zeolites are porous, three-dimensional, aluminosilicate frameworks that contain catalytically active Brønsted acid sites. These active sites catalyze the transformation of hydrocarbons in a confined environment (zeolite pores), which allows for shape-selective effects.[2,3] The regular geometry and the microporous structure of acidic zeolites lead to the high selectivity of specific reactants and products as well as impact the stabilization of the reaction intermediates. Furthermore, acidic zeolites are nontoxic and environmentally friendly catalysts compared to homogeneous analogues, with exceptional physical properties such as excellent mechanical and thermal stability.[4,5]

Acidic zeolite-catalyzed light alkene transformations proceed through extensive and highly interconnected reaction networks. Generally, these transformations can be broadly categorized as:

Oligomerization. Oligomerization of light alkenes is an economically beneficial process to produce liquid linear and branched higher alkenes, which provide key intermediates for manufacturing high-octane gasoline and products such as detergents, oil additives, and petrochemicals.[6,7] The increased availability in recent years of shale gas resources that contain significant percentages of ethane and propane has directed considerable emphasis on oligomerization over acidic zeolites.[8]

Alkylation. The acid-catalyzed alkylation of alkanes with alkenes is crucial in gasoline reformulation processes where light hydrocarbons obtained from catalytic cracking are converted to a complex mixture of branched alkanes, called alkylate, which is used as a blending component to increase gasoline octane number.[9,10]

Aromatization. The production of aromatics (benzene, toluene, xylenes) from light hydrocarbons is pivotal for the manufacture of fine chemicals and plastics.[11] The conversion of alkanes to aromatics over acidic zeolites proceeds through the formation of intermediate alkenes. The addition of a metal promoter to the zeolite favors the dehydrogenation of alkanes into alkenes.[12−15]

The processes described above typically involve complex feedstocks that consist of several hundreds or thousands of molecular species. Furthermore, the acid-catalyzed chemistry results in the formation of protonated surface species as intermediates that are interconnected by thousands of elementary reactions steps. A mechanistic approach to model the kinetics of these complex reacting systems consists of the automated development of reaction networks and microkinetic models that provide a detailed description of the reaction pathways, including each possible elementary step and reaction intermediate.[16] As an example, in an oligomerization network, Table quantifies the potential number of species that are involved in the reaction networks for alkene conversion on acidic zeolites. Networks composed of this level of complexity and specificity of species defy manual construction; rather automated network generation algorithms are necessary computational tools to construct detailed networks and elucidate the complexity of chemical transformations.

Table 1

Number of Gaseous and Protonated Structural Alkene Isomers for a Given Total Carbon Length

alkene carbon length	number of gaseous structural isomers	number of protonated structural isomers (no primary)
3	1	1
6	13	9
9	153	122
12	2281	1819
13	5690	4852
14	14 397	11 602
15	36 562	29 623
16	93 646	76 037

However, kinetic models must balance the specificity of the species that are included with computational expense; if too large a span of carbon number is covered, the model may be too unwieldy to be solved or used. Returning to Table , the possible number of alkene isomers increases exponentially as the carbon number increases. Once a carbon length of C15 is reached, a fully detailed microkinetic model already includes a total of about 60 000 unique gaseous species and about 50 000 unique protonated intermediates, for a total of about 110 000 unique species. Similarly, the hundreds of thousands of reactions between each unique isomer create a system of equations too large and stiff to be solved molecularly that likely includes numerous kinetically insignificant pathways. Rather, lumping similar isomer groups together and assuming each lump is at equilibrium can allow for higher carbon lengths to be reached.

Lumping techniques have been used in the past, and lumps are often a compromise between the capabilities of the experimental analytics to characterize and the computational complexity to replicate.[17] For example, one of the first lumped models was proposed by Nace and Weekman in 1971 to model fluid catalytic cracking (FCC) that contained only three lumps corresponding to the three distillation cuts: unreacted gas oil, gasoline, and gas + coke.[18] Since then, experimental analytics and computational power have improved, leading to lumping approaches based on carbon number and number of branches or carbon number and number of rings that have been applied by a variety of kinetic modeling endeavors to model catalytic hydrocracking[19−21] and catalytic reforming,[22] to name a few.[17] However, these previous lumping studies mainly focus on lumping based on carbon number alone and are assumed because large groups of intermediates are considered in equilibrium already (such as during hydrocracking). In this work, more detailed kinetic lumping is required for processes like oligomerization, where lumping should depend on not only carbon number but also the degree of branching and the protonated surface intermediates covering the catalyst within the reactor.

In this work, a library of detailed reaction networks was constructed using an automated network generator to describe several industrially relevant alkene transformations over acidic zeolites, namely, light alkene oligomerization, alkylation, and aromatization. The underlying mechanisms for these transformations can be described based on a limited number of reaction families, and the complexity of the reaction mechanisms can be tuned by selecting appropriate termination criteria for the network generation algorithm. We also proposed a lumping strategy based on carbon number, degree of branching, and ion type for gaseous species and protonated intermediates to create kinetic models that can reach high carbon numbers while keeping a high level of mechanistic detail. The proposed reaction networks are sufficiently general to be tailored to the different acidic zeolite frameworks of interest via combination with an appropriate set of kinetic descriptors. The aim of the work is to provide ready-to-use tools for the mechanistic description of alkene transformations on acidic zeolites and the construction of microkinetic models with different degrees of complexities.

The computational details for the automated generation of reaction networks and lumping strategy are presented in Section . Section describes the generated reaction networks. Finally, we discuss the proposed lumping strategy based on two case studies applied to propene oligomerization (Section ).

Computational Details

Automated Network Generation and Visualization

Automated network generators are used to develop networks of elementary reactions with different degrees of complexities depending on the reaction types, or reaction families, and the reaction rules that are defined by the user.[16,23]

In the present implementation, chemical species are represented using graphs and converted into mathematical expressions using the concept of the bond and electron (BE) matrix.[24] The diagonal element [BE] of a BE matrix represents the number of free valence electrons of atom i. The off-diagonal element [BE] of the BE matrix represents the order of the bond between atoms i and j with [BE] = 0 if i and j are not adjacent. Types of chemical reactions are often referred to as reaction families and represent what bonds break and form when a species undergoes a chemical reaction. A reaction family operator can be given as an input to generate a product from a reactant or it can be defined as the difference of the BE matrices that represent the product(s) and the reactant(s) of the chemical transformation. This concept is illustrated in Figure that shows an exemplary oligomerization step between ethene and ethoxy (that is represented for simplicity as a carbenium ion).

Figure 1

Construction of the oligomerization reaction operator. The oligomerization of ethene with ethoxy is chosen as an exemplary reaction step. First, the reactants and products are expressed as BE matrices based on the atom connectivity. Second, the reactant BE matrices are combined, and reduced BE matrices that include only the atoms that are affected by the chemical reaction are identified for both reactant and product. Third, the reaction operator is obtained as a difference of product and reactant BE reduced matrices.

The BE matrices for the reactants ethene and ethoxy are combined into one matrix and permuted to move the atoms that take part in the chemical reaction to the top. This step helps identifying reduced BE matrices that include only atoms for which either bonds break or form. Following this procedure, a reaction operator can be defined for an exemplary elementary reaction of each reaction family.

The same reaction operators can be used to describe multiple unique reactions that belong to the same general family by applying the same reaction operator to a variety of reactant species BE matrices. Moreover, complex networks can be created from a small set of starting species and reaction families by allowing products in one run to be reactants in the next.

The reaction family operators are implemented in NetGen, a software developed by Broadbelt and co-workers.[25,26] NetGen automatically generates the reaction products by adding the reaction operators to the reduced BE matrices that represent the reactant species. Graphical networks in this work have been created via Cytoscape 3.9.0 software[27] to visualize the reaction mechanisms and their growth as a function of the termination criteria.

Kinetic Parameters

The kinetic parameters of these models are typically estimated based on a solid theoretical basis, and no a priori assumptions about the rate-determining step(s) are needed.[17] As a result of their robust structure, microkinetic models are powerful tools to design industrial processes, optimize the operating conditions of chemical reactors, and inspire novel catalyst development.

The Arrhenius equation was used to calculate rate constants k, where A is the pre-exponential factor, Ea is the activation energy, R is the universal gas constant, and T is the temperature (eq ).Further, Evans–Polanyi relationships[28] were used to calculate the activation energies (eq )where E0 is the intrinsic energy barrier and ΔHR is the enthalpy of the reaction. For all cases, the α transfer coefficient (0 ≤ α ≤ 1) was set to 0.1 and 0.3 for oligomerization and protonation, respectively, and 0.5 for all isomerization reactions.[29,30]

The enthalpy of reaction is dependent upon ΔHf, which denotes gaseous enthalpy of formation, ΔHphy, which is the physisorption energy for a physisorbed species within a specific framework, Δq, which is the stabilization energy for a protonated intermediate within a specific framework, and υ, which is the species’ stoichiometric coefficient (eq ).The kinetic parameters of kinetic models for acid catalysis are a combination of “kinetic/thermodynamic” and “catalyst” descriptors.[31] Kinetic/thermodynamic parameters include gas-phase enthalpies of formation, which are independent of the catalyst. The gas-phase formation enthalpies of the species considered here were calculated using group additivity values based on the formulation by Benson.[32] A value of 365.7 kcal/mol was considered for the enthalpy of formation of the proton in the gas phase.[33] Catalyst descriptors depend on the specific framework of the zeolite, which account for the impact of the zeolite topology on the kinetics (physisorption enthalpy and stabilization enthalpy). Physisorption enthalpy trends were obtained from hydrocarbon-focused density functional theory (DFT) studies published by De Moor et al.,[34−36] as well as Kostetskyy et al.[37] In those works, physisorption enthalpies are reported as linear scaling relationships dependent on carbon number. Stabilization energy trends were also taken from DFT studies.[38−41] The combination of the contributions from the gas-phase reaction enthalpy with the catalyst descriptors allows the estimation of the reaction enthalpy on the zeolite surface for each of the elementary steps included in the network.

Carbon and Rank Termination Criteria

The automated generation of a reaction network consists of repeatedly applying the reaction operators to the reactants and their progeny. However, this algorithm results in the possibility of generating an infinite number of reactions and species. To produce a reasonably sized reaction network, a termination criterion can be applied to the generation algorithm. In this work, a “carbon and rank” termination criterion[25,42] has been applied, which consists of limiting the maximum carbon number of the species that are allowed to react as well as the rank. The rank of a species reflects its associated order of appearance in the reaction network. Reactants have rank 0. The rank of the products increases by 1 depending on their order of appearance such that primary products are associated with rank 1, secondary products with rank 2, and so forth.[42] The increase of the species rank is only associated with the formation of molecules and not with the formation of protonated species. In this work, the criteria used to terminate the network generation algorithm have been indicated using the expression CR, where C represents the maximum carbon number and R represents the maximum rank of the species allowed to react.

Oligomerization Networks

An acid-catalyzed alkene-oligomerization network converts lower carbon number alkenes to both linear and branched oligomers. The first step of the oligomerization pathway is the activation of the reactant alkene. This consists of the physisorption of the alkene from the gas phase to the pores of the zeolite, followed by protonation to a chemisorbed protonated intermediate.[43] The protonated intermediate has been proposed to take the form of either a covalent alkoxide that is connected to an oxygen atom within the zeolite framework or as an ion pair (carbenium ion) with the negatively charged zeolite framework.[44] All protonated intermediates are represented in this work as carbenium ions for simplicity.

Deprotonation is the reverse step that regenerates a physisorbed alkene and returns a proton to the zeolite. The number of physisorption and desorption reactions are both equal to the total number of alkenes that are supplied and generated in the network. The protonated species can increase their hydrocarbon chain length through an oligomerization step, representing the addition of a physisorbed molecule to an alkoxide or carbenium ion. The reverse step is β-scission that returns a smaller alkene and a chemisorbed species.

The skeletal rearrangement of the chemisorbed species proceeds through four types of isomerization steps, namely, hydride shift, methyl shift, as well as α- and β-protonated-cyclopropane-mediated (PCP) branching.

The list of reaction families that are involved in the overall oligomerization network and the exemplary elementary steps are listed in Table with the associated reaction operators. To note, this network assumes only alkenes and no alkanes are present. This alkene-oligomerization network can be assumed when only small alkenes are present with few abstractable hydrogens, limiting the hydride-transfer reaction from producing alkane products as well as alkenyl and allylic intermediates. By assuming the network only depends on alkene chemistry, the complexity and size of the network are greatly reduced, which lends to a much more manageable kinetic model. For this reason, it is important to consider if the alkene feedstock and products remain small (

Table 2

List of Reaction Families with Exemplary Elementary Steps and Associated Reaction Operators that Are Included in the Alkene-Oligomerization Networka

The suffixes (g) and (p) represent the gaseous phase and pores of the zeolite, respectively.

Alkylation Networks

Alkenes and alkanes can kinetically interact by introducing a hydride-transfer step in the list of reaction families. The mechanism of hydride transfer involves the attack of a physisorbed molecular species by a protonated intermediate resulting in the transfer of a hydrogen atom.

The transfer of a hydride from the alkyl group of an alkene to a protonated intermediate forms an alkenyl intermediate, which deprotonates and desorbs as a diene. Similarly, the transfer of a hydride from an alkane to a protonated intermediate that is otherwise saturated results in the formation of a different alkane.

Allylic isomerization was also included in the network to represent the different resonance forms of an allylic ion and the corresponding reaction products, given that charge is associated with a single atom for bookkeeping purposes. The reaction families that complete the list presented in Table to generate the alkylation network are depicted in Table .

Table 3

List of Reaction Families with Exemplary Elementary Steps and Associated Reaction Operators to be Added to the List in Table to Generate the Alkylation Network

Cyclization and Aromatization Networks

Cyclization is the transformation of a linear protonated intermediate into a cyclic protonated intermediate. These can either undergo a deprotonation step resulting in the formation of a cycloalkene or a cyclic hydride-transfer step to form a cycloalkane. Four cyclization families that lead to the formation of 1,5 endo, 1,5 exo, 1,6 endo, and 1,6 exo cyclization protonated intermediates were included in this network. Additionally, the reverse reactions of the cyclization steps were included in the network as a form of β-scission. Aromatics can be formed based on a series of deprotonation and hydride-transfer steps applied to a cyclohexene alkoxide or carbenium ion. The reaction families that complete the list presented in Tables and 3 to generate the cyclization and aromatization networks are depicted in Table .

Table 4

List of Reaction Families with Exemplary Elementary Steps and Associated Reaction Operators to be Added to the List in Tables and 3 to Generate Cyclization and Aromatization Networks

Lumping Strategy

Lumping involves grouping isomers, which may be done through various means.[48] In our case, lumps are based on unique combinations of carbon number, degree of branching, and ion type (primary, secondary, or tertiary) (Figure ). Within the group of isomers or lump, species are considered to be in thermodynamic equilibrium. In other words, internal isomer reactions that do not change carbon number, degree of branching, nor ion position are assumed to be much faster than reactions that affect the carbon backbone or ion type. This resolution of different ion types advances our lumping technique beyond previously used approaches by retaining the differences in ion stabilities on reaction rates, which can be accounted for with Evans–Polanyi relationships. Pathways between lumps with distinct kinetic parameters are still retained in the model. Each lump is represented by a reference species, which is used to calculate the equilibrium concentrations of all other molecules within the lump. The equilibrium term added to the rate expression is often referred to as a “lumping coefficient” and is determined by thermodynamic properties.

Figure 2

Examples of hexene lumped species. Each lumped species has a unique combination of carbon number, branching degree, and ion type. Ion types include primary alkoxide (1), secondary alkoxide (2), tertiary ion (3), or gaseous species (0). Highlighted species represent example reference species for the entire lump.

Lumping is done after network generation to ensure all pathways are still considered and to ensure kinetic rate constants remain specific to the unique reactant. Most often, generating the large network through matrix addition is far easier than modeling kinetics by solving (stiff) differential equations. Moreover, kinetic models must be solved multiple times for optimizations, while the underlying network is only generated once. In other words, network complexity generally impedes developing the kinetic modeling rather than automatic network generation. However, there are works that approach lumped network generation as well when species get very large.[20]

By applying lumping to an initially fully defined system, the reaction pathways as well as rate constants still retained all specificity of the initial system, prior to lumping, in which all isomers are considered. Only in calculating the concentration of the species, to be multiplied by the rate constant for the final rate calculation, was lumping influential (Figure ). Concentrations are rewritten assuming equilibrium with respect to a reference species. Species in thermodynamic equilibrium with each other can be simplified into algebraic expressions, effectively reducing the stiffness of the differential equations in the model. Additionally, internal isomerization reactions in which the reactant and product are in the same lump can be removed, further reducing the number of rates.

Figure 3

Schematic exemplifying the reduction realized by lumping of rates and differential equations for flow of gaseous species on the protonation/deprotonation of linear 2- and 3-hexene. In this example, h+ represents a linear C6 protonated intermediate with h+,2 and h+,3 being a secondary or tertiary ion, respectively, which can produce either linear 2-hexene or 3-hexene. H+ represents a proton from the acidic support. Bolded species, h+,2 and 2-hexene, are the reference species for each lump. When lumping, equilibrium (Keq) is assumed within a lump and the concentrations in the rate expression are rewritten based on the lump reference species (C2-hexene or Ch). The rate constants, such as pre-exponential factors and activation energies, are not changed and are still calculated specifically for each unique isomer. Isomerization rates within a lump can be removed as the relationship is maintained with an equilibrium expression. Overall, only the reference lump differential flow rates need to be solved, reducing the model stiffness and complexity while maintaining high mechanistic detail.

This type of lumping technique is notably applicable in modeling catalytic reactions in which species can grow in carbon number, and the network could theoretically grow to infinite size. Limiting the model to too few carbon numbers can lead to inaccurate accumulation of species at the carbon limit, and such termination effects have been reported already in previous propene oligomerization models.[29] Thus, being able to extend an oligomerization model to lengths past what is experimentally reported is likely necessary to accurately capture kinetics capable of oligomerization.

Results and Discussion

This section describes the reaction networks obtained using the methodology presented above. The library of reaction networks described in this section is available on the GitHub repository that is linked to this paper (https://github.com/JosephColin1/supporting-documents). The networks were catalogued based on the carbon and rank termination criteria (indicated as CR) that were used for the automated generation process. Each elementary step is associated with a reaction family and assigned a gas-phase reaction enthalpy.

This section presents the networks generated using the reaction families listed in Table by applying the termination criteria CR0 with i ranging from 2 to 14. The networks were generated by considering ethene and the zeolite proton as the only rank 0 species. However, these can be used to simulate the oligomerization chemistry of any other alkene included in the network. An exponential growth trend was observed for the number of generated alkenes and protonated intermediates (Figure a) as well as elementary reactions (Figure b) as a function of the carbon termination criterion. The exponential increase of the isomerization reaction steps included in the network corresponds to an increase in the number of possible isomers. This provides a complete description of the chemistry of the process, potentially allowing for a detailed prediction of process performance and distribution of isomers. However, this increased complexity of the reaction network is difficult to handle, and the solution of the associated models is computationally expensive or impossible without the implementation of appropriate reduction techniques, as discussed in Section .

Figure 4

Number of (a) molecules and protonated intermediates and (b) elementary reactions included in reaction networks for alkene oligomerization. The networks were generated using ethene as the initial reactant and applying the termination criteria CR0 with i ranging from 2 to 14.

Figures and 6 show exemplary acid-catalyzed oligomerization networks obtained using ethene as the initial reactant and different termination criteria. Red edges indicate protonation and deprotonation steps, green edges indicate oligomerization and β-scission steps, and blue edges indicate isomerization steps. Black nodes refer to molecular species, and white nodes refer to protonated species. The complete reaction networks including the list of elementary reactions and the Cytoscape files for visualization are available in the section “Oligomerization Networks” of the GitHub repository.

Figure 5

Projected reaction network to describe the acid-catalyzed alkene oligomerization. The network was obtained using ethene (“spe2” highlighted in red) as the initial reactant and the C6R0 termination criterion. Red edges indicate protonation and deprotonation steps, green edges indicate oligomerization and β-scission steps, and blue edges indicate isomerization steps. Black nodes indicate molecular species, and white nodes indicate protonated species. The graphs of molecular and protonated species are listed in Table S1.

Figure 6

Exemplary projected reaction networks that describe the acid-catalyzed alkene oligomerization. The networks were obtained using ethene as the initial reactant and by applying the following termination criteria: (a) C8R0, (b) C9R0, and (c) C10R0. Red edges indicate protonation and deprotonation steps, green edges indicate oligomerization and β-scission steps, and blue edges indicate isomerization steps.

In this section, we present a library of reaction networks that were generated using the reaction families presented in Tables and 3. The alkylation networks were obtained using the termination criteria CR with i ranging from 2 to 12 and j ranging from 0 to 3. The carbon termination criteria were limited to C8 for rank termination criteria R. Each network was generated by considering ethene and the zeolite proton as the only rank 0 species. The presence of alkanes in the product mixture is guaranteed by the occurrence of hydride-transfer steps, as previously discussed. Figure shows an exemplary reaction network describing the acid-catalyzed alkylation of isobutane with propene to 2,2-dimethylpentane. Additional alkylation products can be obtained through isomerization and skeletal rearrangements of the C7 protonated intermediate.

Figure 7

Exemplary core catalytic cycle that is part of the reaction network of the acid-catalyzed alkylation of isobutane with propene to 2,2-dimethylpentane.

Figure shows the exponential growth of the number of gas-phase molecules and protonated intermediates (a) and elementary steps (b) included in the networks as a function of the carbon termination criteria for a R0 rank termination criterion. As evident from a comparison between Figures and 8, the inclusion of hydride-transfer steps results in a considerable increase in the size of the reaction network. As an example, for a C12 carbon termination criterion, the inclusion of hydride-transfer steps results in an increase of a factor of 6 in the number of generated species, and of a factor of 8 in the number of generated elementary steps.

Figure 8

Number of molecules and protonated intermediates (a) and elementary reactions (b) included in the reaction networks for alkane/alkene alkylation networks. The networks were generated using ethene as the initial reactant and applying the termination criteria CR0 with i ranging from 2 to 12.

The generated molecular species include alkanes, alkenes, and dienes with the distribution depicted in Table for a carbon termination criterion ranging from C8 to C12. Figure a,b summarizes the results of the network generation procedure, in terms of the number of generated species and elementary steps, for different combinations of the carbon and rank termination criteria. The complete reaction networks described in this section including the list of elementary reactions and the Cytoscape files for visualization are available in the section “Alkylation Networks” of the GitHub repository. The same reaction networks presented in this section can be used to describe the acid-catalyzed cracking process of heavy alkene and alkane feeds.[49−51] In these processes, hydride-transfer steps need to be considered because of the presence in the feed of large alkenes with many abstractable hydrogens.

Table 5

Number of Gas-Phase Molecules for the Alkane/Alkene Alkylation Network Obtained by Imposing the Termination Criteria CR0 with i Ranging from 8 to 12

	formula	C8	C9	C10	C11	C12
alkanes	CnH2n+2	39	74	149	308	643
alkenes	CnH2n	116	269	646	1560	3841
dienes	CnH2n–2	128	390	1138	3232	9112

Figure 9

Total number of generated species (a) and elementary steps (b) for the alkane/alkene alkylation network as a function of the carbon (C) and rank (R) termination criteria. The networks were generated using ethene as the initial reactant.

This section presents the reaction networks that were generated based on the complete set of reaction families presented in Tables –4. Figure shows the number of generated species (a) and elementary steps (b) as a function of the carbon termination criterion for a R0 rank termination criterion. The species distribution of the networks presented in Figure is highlighted in Table , including alkanes, alkenes, dienes, cycloalkanes, and cycloalkenes.

Figure 10

Number of molecules and protonated intermediates (a) and elementary reactions (b) included in cyclization and aromatization networks. The networks were generated using ethene as the initial reactant and applying the termination criteria CR0 with i ranging from 2 to 12.

Table 6

Number of Gas-Phase Molecules for the Alkane/Alkene Cyclization/Aromatization Network Obtained by Imposing the Termination Criteria CR0 with i Ranging from 8 to 12

	formula	C8	C9	C10	C11	C12
alkanes	CnH2n+2	39	74	149	308	643
alkenes	CnH2n	116	269	646	1560	3841
dienes	CnH2n–2	128	390	1138	3232	9112
cycloalkanes	CnH2n	22	62	173	477	1315
cycloalkenes	CnH2n–2	92	309	992	3069	9286

The presence of aromatics in the reaction network can be obtained using the same reaction families listed in Tables –4 by varying the rank termination criterion. This is exemplified in Figure that shows a possible aromatization pathway of ethene to benzene. First, three molecules of ethene undergo oligomerization to 1,4-hexadiene, which is a rank 2 species. A hydride-transfer step and the following cyclization and deprotonation result in the formation of 1,3-cyclohexadiene as a rank 3 species. Benzene is finally formed from an additional hydride-transfer step followed by deprotonation. Although this is only one possible aromatization path of ethene to benzene, it is worth noting that a rank termination criterion R≥3 is required to generate benzene from ethene. Figure shows the simultaneous impact of the carbon and rank termination criteria on the total number of generated species (a) and elementary steps (b) using ethene as the initial reactant. The presence of branched aromatics heavier than benzene can be obtained by increasing the maximum carbon number of the species that are allowed to react. For example, the termination criterion C7R3 results in the formation of toluene, and the termination criterion C8R3 results in the formation of xylenes.

Figure 11

Aromatization pathway for the formation of benzene from ethene. R indicates the rank of the molecular species in the reaction network.

Figure 12

Total number of generated species (a) and elementary steps (b) for the cyclization and aromatization network as a function of the carbon (C) and rank (R) termination criteria. The networks were generated using ethene as the initial reactant.

The complete reaction networks described in this section including the list of elementary reactions and Cytoscape files for visualization are available in the section “Cyclization and Aromatization Networks” of the GitHub repository.

Lumping in Microkinetic Modeling

In the next sections, we use propene oligomerization to test the validity of the lumped modeling technique. The oligomerization network generated using the termination criterion C9R0 and described in Section was used to generate both the full model and the lumped model, ensuring both models have the same underlying reaction networks and calculated thermodynamics. The full network is available in the GitHub repository that is linked to the paper, and more details of the model can be found in the Supporting Information.

Case 1. Full Model vs Lumped Model at Equilibrium

When applying lumping techniques, we assumed that internal isomerization reactions among species with the same carbon number, degree of branching, and ion type occurred very fast and reached thermodynamic equilibrium. To ensure isomer congruity, we compare the lumped isomer (“lump”) and fully defined (“full”) models when isomerization is manually set to be very fast. In the case of fast isomerization reactions, assumed in the lumped model and applied in the full model, both results should converge to similar isomer fractions.

Isomerization reactions were set to be very fast by increasing the pre-exponential factor to 5 × 1012 s–1 and lowering the intrinsic activation energy to 5.0 kcal/mol (Table S3). Product mole fractions with respect to the entire product distribution are compared in Figure . A low conversion of about 1% was used in the comparison.

Figure 13

Product mole fractions for full and lump model simulations. Fractions were calculated with respect to the entire product distribution (C4–C9). Conversion for both simulations was about 1%.

As expected, both models converge to very similar results, with the minor differences that can be explained by the inherent differences in the model network (irreversible reactions now assumed reversible), tolerances imposed to solve the stiff set of model equations, or reaction path degeneracy. In both cases, the models predict that about 45 mol % of hexenes produced are linear, followed by about 35 mol % methylpentenes, and 2 mol % dimethylbutenes. Isomers of butene are produced in about even amounts (note that mole fractions are calculated across the entire product distribution and not within each carbon number). Overall, isomer fractions between the full and lumped model show similar mole fractions, suggesting that the models converge with fast isomerization or assuming isomerization equilibrium.

Case 2. Experiments vs Full Model vs Lumped Model

In this section, we compared the results of the fully defined and lumped isomer models against actual experimental data. The experiments were taken from the previously mentioned propene oligomerization study[29] and are summarized in Figure , and the experimental conversion and selectivities are reported in Table S4.[29] Initial time on stream data was used due to catalyst deactivation, and only up to C9 species were measured in the experiments. In this case, we want to show that the lumped model can still match experimental data, even when the ion distribution is not at equilibrium, as assessed by the full model.

Figure 14

Parity plots of conversion and selectivity comparing the full model and the lumped model to propene oligomerization experimental values. Kinetic parameters can be found in Tables S4–S6.

First, a C3–C9 model was developed, following the procedure described elsewhere.[29] The fully defined model encompassed 628 species and 2615 reactions. Compared to the experimental data, the full model results in an RMSE value of 24.8 (Figure A, kinetic parameters in Table S5). However, the high production and even overproduction of C9, the maximum carbon length, suggest that flux is artificially stuck at this limit as no further oligomerization pathways are present. As a result, C9s accumulate and are forced to crack into smaller products. This, in turn, can lead to incorrect kinetic parameters upon optimization such that the experimental data is captured, which results in a detailed kinetic model with limited utility. The lumped model does sufficiently well in describing the experiments, with an RMSE of 33.7 (Figure B, kinetic parameters in Table S6). A significant amount of this error comes from underproducing C5 and C7, which are difficult to produce with a carbon limit of C9, as either butenes (C5 + C4) or ethene (C2 + C7) must be produced, which is often energetically unfavorable.

Comparing kinetic parameters, both the full model and the lumped model share similar rates of protonation and deprotonation, as well as oligomerization and β-scission rate constants. Rate constants for isomerization, specifically for methyl and hydride shift, change more significantly. This difference is most likely due to the lumped model assuming most of these isomerization reactions are in equilibrium as well as the model being less sensitive to changes in isomerization reactions, which are fast, compared to oligomerization, which are slower and often rate-determining. To note, the lumped model is significantly smaller than the full model, containing only 79 reference species and 974 reactions. Although the full model has a lower RMSE compared to the lumped model, the solution of the lumped model is significantly computationally less expensive than the full model.

Next, a C3–C12 model including 8832 species and 42 097 reactions was developed to better describe the experimental data. Due to the extended carbon limit, the solution of this model failed past about 2.8% conversion (Figure C). On the other hand, the lumped model was able to be solved and optimized to fit the experimental data. The extended lumped model included 152 species and 4814 reactions and drops the RMSE to 31.7 (Figure D, kinetic parameters in Table S7). In this extended lumped model, C12s are being produced and more easily cracked to match the higher selectivities of C5 and C7. The absence of >C9 species in the experimental data does not exclude the formation and accumulation of these heavier species in the pores of the catalyst, causing the reported deactivation. Nevertheless, these heavier species being produced in either the experiments or the model remain relatively low for the mass balances to close ≥ 95%. Additionally, the rate of oligomerization in the C3–C12 model is lower than the C3–C9 model and is likely more accurate. In the previous C3–C9 case, an overestimation of the oligomerization rate would only provide more carbon into the C9 fraction, whereas in this larger C3–C12 model, an overestimation of oligomerization rates will direct carbon away from C9 and into C12.

Despite the RMSE being slightly higher than the full C9 model, the lumped C12 model likely provides a better description of the intrinsic kinetics of the system, without the false accumulation of the maximum carbon species and more viable paths to produce C5 and C7.

One of the benefits of using a mechanistic approach in kinetic modeling is being able to distinguish products of unique carbon number and branch degree. However, many times, experimental analysis only reports groups based on carbon number, such as in these propene experiments. Both the full and lumped models suggest slightly different branching isomer fractions; however, it is impossible to know which is closest to the experimental “truth” until such detailed experiments are carried out. Intriguingly, the availability of models that capture this level of detail as offered in the present work provides impetus for additional experimentation and even provides guidance about what species may dominate to aid the interpretation of complex spectroscopic analyses.

Conclusions

Computational automatic network generation is a powerful tool to create large networks encompassing many industrially relevant acid-catalyzed hydrocarbon chemistries including oligomerization, alkylation, and aromatization. These complex chemistries can be summarized in succinct sets of reaction families with chemical species written as a bond and electron matrix. Carbon and rank limits are required to provide bounds and focus on relevant carbon numbers and products. Examples of reaction networks and their reaction families, with Cytoscape visualization, show how these networks grow in complexity and can quickly become extremely extensive. To better utilize extensive kinetic models, lumping is required to reduce model size and computational load. Our lumping approach uniquely uses a combination of carbon number, branching degree, and ion position to retain high product specificity for both gaseous and protonated intermediates. The addition of ion position allows for a better description of the catalytic surface and preserves ion stability differences. This novel lumping technique was explored in three case studies, which confirms congruence with a fully defined kinetic model. Additionally, the application of lumping was confirmed to capture experimental oligomerization kinetics to a high degree of fidelity. Moving forward, this lumping technique could be applied to a wide variety of chemistries allowing for more robust and expansive models in the future.