A Fragment-Based Mechanistic Kinetic Modeling Framework for

Oct 1, 2018 - A Fragment-Based Mechanistic Kinetic Modeling Framework for Complex Systems. Kehang Han and William H. Green. Ind. Eng. Chem...
0 downloads 0 Views 2MB Size
Subscriber access provided by UNIV OF LOUISIANA

Kinetics, Catalysis, and Reaction Engineering

A Fragment-Based Mechanistic Kinetic Modeling Framework for Complex Systems Kehang Han, and William H. Green Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.8b02870 • Publication Date (Web): 01 Oct 2018 Downloaded from http://pubs.acs.org on October 4, 2018

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 27 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 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

A Fragment-based Mechanistic Kinetic Modeling Framework for Complex Systems Kehang Han and William H. Green∗ Department of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, United States E-mail: [email protected] Phone: +1 (617)-253-4580. Fax: +1 (617)-258-8992

Abstract There is great interest in developing chemical kinetic models with predictive power. A predictive kinetic model maintains high fidelity to the true chemistry usually via a full-detail molecule representation which distinguishes all the molecules and isomers. Unfortunately such full-detail representation quickly becomes intractable for large molecules, because the number of isomers grows very rapidly with size. Several less detailed representations have been developed, including the popular StructureOriented Lumping (SOL) method, which allow one to model larger molecules. In this paper we first examine scalability for two existing representations (full-detail and SOL) and propose a new modeling methodology tailored for heavy systems: fragment-based modeling, which promises the best scalability. We show via a case study of heavy molecule pyrolysis that the new approach creates a much smaller reaction network but with similar prediction accuracy on feedstock conversion and products’ molecular weight distribution compared to the corresponding model using full-detail representation. An open source modeling package AutoFragmentModeling using this methodol-

1

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 27

ogy has been developed. Some suggestions are presented for future work based on this new modeling approach.

I. Introduction Facing with changing feedstock compositions and operating conditions, the oil and gas industry has been looking for reliable methodologies to generate kinetic models for various reacting systems and hoping to extrapolate well to new situations. With improvements in computer power and the successes in method development of molecular dynamics, quantum mechanics, numerical simulations, the past decades have seen the emergence of predictive kinetic modeling. A predictive kinetic model usually represents species in a very detailed way which distinguishes all the molecules and isomers. That full-detail representation allows high fidelity to the true chemistry, and so the corresponding models have superior extrapolation potential. The full-detail approach has achieved many successes in modeling relatively small, simple reacting systems, e.g., natural gas combustion, 1 butanol pyrolysis and combustion, 2 steam cracking of n-hexane, 3 n-alkane combustion up to n-hexadecane. 4 However it quickly becomes intractable for large molecules in applications such as crude oil upgrading, vacuum residue cracking, coal pyrolysis, etc., which are of great interest to the industry. To illustrate the poor scalabilty, consider a simple feedstock composed of molecules with linear backbones of M carbon atoms: up to L functionalized carbon atoms (e.g., carbonyl, carboxyl, etc.) and all the rest unsubstituted CH2 or CH3 carbon atoms. Suppose there are S different carbon types for each functionalized carbon to choose from. The number of distinct molecules is

Nf ull−detail

 L  X M = SN N N =0

2

ACS Paragon Plus Environment

(1)

Page 3 of 27 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 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

The last term in the summation is

ln Nf ull−detail ≥ ln

M L



S L , which gives a lower bound for Nf ull−detail .

 M  L

  S L = ln

M !S L  (M − L)!L!

(2)

Usually vacuum gas oil can have 20-45 carbons, and vacuum residue has more than 45 carbons. If we have a large molecule with 30 carbon atoms (M = 30), and functionalize half of the carbons (L = 15) from 8 functional group types (S = 8), then ln Nf ull−detail would be larger than 50: it requires e50 ≈ 1022 distinct molecules to represent the feedstock. Even 1 byte per each molecule requires total memory of 1022 bytes = 1010 Tb, which is 30 times the entire Internet data of 2011. 5 Unfortunately, many applications that have significant impact on energy, environment and society are complex chemical processes which may have even more complex feedstocks with more than 30 carbon atoms and including branched and cyclic species, such as crude oil, coal, and materials derived from biomass and municipal waste. To model these feedstocks, less-detailed representation methods are needed. Over the past decades, people have tried to understand those systems and model the essential chemistry. Lumping strategies have been widely employed where in early practice molecules were lumped based on molecular weights (or boiling points). See an example of a 3-lump model for the catalytic cracking process in Figure 1. Despite simplicity, the lumping strategy introduces an intrinsic error since the lumps contain molecules with different reactivities; they don’t all react at the same rate as assumed in the lumped model. Consequently, the composition of a lump changes as the reaction proceeds (e.g. certain molecules get enriched) so experimentally the lump does not have fixed properties nor follow first order kinetics. The kinetic parameters are often fitted from a very limited set of pilot plant experiments. These models typically have poor accuracy in extrapolation. When the feed or the operating condition changes, new experiments are usually required and model parameters refitted (sometimes even model structure needs to be revised to account for new chemistry). To overcome this problem, Structure-Oriented Lumping (SOL) 7,8 was invented, which distinguishes molecules with different functional groups (i.e. different reactions). SOL pre3

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 27

Figure 1: A simple lumped model for catalytic cracking process by Weekman and Nace 6 defines a list of functional groups, and represents a molecule using a vector of which each entry corresponds to a group and records the count of that group existing in the molecule. It assumes the functional groups are chemically independent in a molecule and ignores the connectivity between them. Its scalability can be illustrated in the same framework set by the aforementioned example: we have S pre-defined functional groups and up to L carbon atoms are functionalized. The number of distinct molecules is essentially the number of distinct solutions (vi ) to the following inequality: S X

vi ≤ L

(3)

i=0

where vi is the ith entry of the vector (length S) that represents a molecule in SOL. If L is sufficiently large (e.g. L > S), the number of physically meaningful SOL vectors can be approximated by the volume of the region in the positive quadrant which satisfies Eq.(3). When the equality holds, Eq.(3) is the definition of the hyperplane cleaving a hypercube in two, so this is half the volume of a hypercube with length L:

NSOL ≈

LS 2

(4)

For the case considered above, M = 30, L = 15 and S = 8, NSOL ≈ 109 , which is much less than the full-detail representations but still quite large. As a consequence, SOL simulations are also usually restricted to simplified unimolecular chemistry. From Eq.(4) we also see the number of species in SOL grows rapidly with the number of functional groups 4

ACS Paragon Plus Environment

Page 5 of 27 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 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

allowed per molecule (L, which scales with the molecular weight) as well as the number of distinct types of functional groups (S). Aimed to model reacting systems with large and complicated feeds, we present a novel strategy with better scalability: fragment-based approach, which stems from the idea of pseudospecies (distinct functional groups) proposed by Mehl et al. 9 where they modeled combustion chemistries by creating a reaction network of functional groups present in fuels. The fragments in our approach, in contrast to the pseudospecies, may consist of multiple functional groups as needed and preseve their original connectivity. The larger chemical context allows more accurate estimation of kinetic parameters. To standardize the modeling workflow, we designed the fragment-based framework with four components: fragmentization of model compounds, fragment reaction generation, model parameter estimation and reattachment of fragments (Section II), the last component of which enables converting predictions in fragment space back to molecule space. We demonstrate the methodology using a case study of pyrolysis of a C18 hydrocarbon: phenyldodecane (PDD) in Section III. Comparing the results using the fragment method against those from a full-detail method, we find fairly good agreement in feedstock conversion and product molecular weight distribution (Section IV). For heavy feed reacting systems, generating the reaction networks containing hundreds to thousands of reactions by hand is often time-consuming and error-prone. The past decades have seen many efforts devoted to automatic mechanism generation for the full-detail approach, among which are MAMOX, 10 EXGAS, 11 NetGen, 12 RMG 13 and RMG-Py. 14 To help reaction network generation for fragments, we built the first version of an automatic mechanism generation package AutoFragmentModeling, where we adopted several core features from RMG-Py e.g., chemical graph representation, tree-based kinetic parameter estimation. The source code is made available at https://github.com/ReactionMechanismGenerator/ AutoFragmentModeling.

5

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 27

II. Fragment modeling framework Fragment In fragment modeling, we aim to build kinetic models that focus on the key reactive regions of molecules and treat them as independent entities (Figure 2), which are named fragments. We see a large molecule as a combination of many molecular fragments. Each fragment consists of multiple functional groups and preserves their original connectivity in the molecule, which makes fragment a middle layer entity between molecule and functional groups. Ideally, a fragment should be much smaller than the original molecule but still provides large enough environment that is crucial to kinetics. Unlike SOL assuming little spacial interaction between functional groups, the fragment framework accounts for the interaction by accommodating related functional groups within the fragment, which may help improve model fidelity.

Figure 2: In fragment modeling, a model deliberately focuses on key reactive regions and sees a molecule-based reaction as a reaction between related reactive regions The fragment concept makes a much more compact representation for feedstock than aforementioned methods. Using the previous example, suppose fragments have K carbon atoms on average, the total number of distinct species in the fragment representation is

Nf ragment = S K 6

ACS Paragon Plus Environment

(5)

Page 7 of 27

Since fragments are much smaller than molecules, we choose K = 5 (we are good at dealing with chemistry with five carbons) with same parameter setup as previous (S = 8), Nf ragment = 33, 000. This is small enough that kinetic simulations could be solved using modern stiff ODE solvers. That is, arbitrarily complicated reactions of and between these fragments can simulated. Eq. 5 also shows Nf ragment doesn’t increase when dealing with larger molecules as long as modellers keep fragment size and functional groups fixed. Scalability comparison is made in Figure 3. 5-carbon fragments provide a smaller number of species than SOL for C10 or larger feeds. For heavy feeds, e.g., >C20, even 8-carbon fragments provide a more compact representation than SOL. 120

Full-detailed SOL Fragment with size 5 Fragment with size 8

100

lnNdistinctmolecules

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 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

80 60 40 20 0

0

10

20 30 40 Carbon number in feedstock

50

60

Figure 3: Scaling behaviors of number of distinct species with feedstock size for three different methods, assuming 8 distinct carbon arom types, and that at least half of the carbon atoms in the large molecules are not functionalized (i.e., L = M/2). In AutoFragmentModeling, we represent fragments using chemical graphs, with vertices representing atoms and edges representing bonds. This allows graph isomorphism algorithms to perform structure comparison between fragments, facilitating duplicate identification and fragment property estimation (e.g., group additivity method for thermochemistry). We also designed two text-based fragment representations: adjacency list and SMILES (Simplified Molecular-Input Line-Entry System) chemical identifier for users to interact with the AutoFragmentModeling’s internal representation: fragment chemical graph. Figure 4 shows these three representations of an example fragment. 7

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 49 50 51 52 53 54 55 56 57 58 59 60

Figure 4: Three fragment representations for ethyl-R

Fragmentization The fragmentization is a modeling step of converting feedstocks to fragments. It has three working scenarios. • A) systems where individual molecule structures are currently not possible to characterize but molecular weight distribution can be obtained as well as concentrations of functional groups. • B) systems where individual molecule structures can be characterized • C) systems where modellers have selected model compounds (with known molecular structures), typically as a surrogate for the real fuel or feedstock In scenario A, one needs to predefine fragment types based on characterized functional groups and solve for initial fragment concentrations by minimizing the deviations between model and characterization data, which is very similar to the regular feedstock reconstruction procedure. 8,15,16 In scenarios B and C, with original molecule structures at hand, we fragmentize the feedstock molecules with two competing considerations: making fragments small to reduce computation cost, while enlarging fragments to keep important reactivities locally. This paper includes a case study of phenyldodecane pyrolysis, of which the fragmentization step 8

ACS Paragon Plus Environment

Page 8 of 27

Page 9 of 27 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 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

is shown in Figure 5. The R labels, named fragmentization labels, indicate the positions we chose for fragmentization.

Figure 5: In the case study, phenyldodecane is fragmentized into three types of fragments

Fragment reaction generation Similar to molecules, a fragment or a pair of fragments is allowed to react by a certain reaction family if it contains the corresponding reactive region (substructure). For instance, Figure 6 shows the reactant fragment (left) can undergo β-scission since it contains the whole reactive region made up of γ, δ,  carbons.

Figure 6: Example β-scission reaction for fragment In practice, we integrated AutoFragmentModeling with RMG-Py (an automatic full-detail mechanism generation package), 14 to access all the RMG-Py’s reaction families and generate all the unimolecular and bimolecular reactions corresponding to each pair of fragments. Similar to fragments, we also designed a text-based representation (Figure 7) for fragment reaction. If using customized fragment names, users should provide a dictionary mapping between fragment name and SMILES for each reactant and product in the text representation (see the example dictionary and reaction text below). # This is a fragment input snippet # which maps fragment name to fragment SMILES 9

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 49 50 51 52 53 54 55 56 57 58 59 60

Figure 7: Conversion between two reaction representations in AutoFragmentModeling # # 1st part is for stable fragments # the format below is

# fragment name: fragment SMILES ArC(CCCR)CCCR: c1ccccc1C(CCCR)CCCR ArCC: c1ccccc1CC ArCCC: c1ccccc1CCC H2: [H][H] RCC: RCC RCCC: RCCC RCCCCC=CC: RCCCCC=CC ...

# 2nd part is radical fragments # the format below is # fragment name: fragment SMILES

RCCC*: RCC[CH2] ArC*: c1ccccc1[CH2] 10

ACS Paragon Plus Environment

Page 10 of 27

Page 11 of 27 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 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

RCC*: RC[CH2] ArCC*: c1ccccc1C[CH2] RC*: R[CH2] ArCCC*: c1ccccc1CC[CH2] ... # This is a fragment reaction input snippet # In order to understand the reaction texts below, # AutoFragmentModeling needs users to provide # fragment input like the above snippet. # # Beta-scission reactions

ArC*CCCR ArC=C + RCC*

ArCC*CCR ArCC=C + RC*

# Disproportionation or Reverse Disproportionation reactions

ArC=C + C=CC=CC ArC*C + C=CC=CC*

ArC=C + ArCCCCR ArC*C + ArC*CCCR ... Thermochemistry and kinetics estimation It’s not straighforward to define a fragment’s thermodynamic properties such as heat of formation since a fragment is not a physically independent entity like a molecule. However, in kinetic modeling, we usually evaluate thermodynamic property change between two sides 11

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 49 50 51 52 53 54 55 56 57 58 59 60

of a reaction and rarely need the absolute values of those properties. That allows us to choose some practical approach to estimate thermochemistry as long as the thermodynamic property change for a reaction is physically consistent. In AutoFragmentModeling, we create a representative molecule for a fragment whose thermochemistry is to be estimated by consistently concretizing fragmenation labels R into certain groups. The molecule’s thermochemistry will be treated as the fragment’s. Figure 8 shows the two-step estimation process: 1) R label is replaced with an ethyl group to generate a representative molecule. 2). Benson-style group additivity method is used via relevant RMG-Py algorithms to estimate thermochemical parameters (∆Hfo , S o and Cp ) of the representative molecule, which eventually is assigned to the original fragment. More details on how the functional groups are organized and estimation logic carried out can be found in RMG-Py’s main paper. 14

Figure 8: Fragment thermochemistry estimation scheme in AutoFragmentModeling Similarly for kinetics estimation (Figure 9), AutoFragmentModeling extracts reaction templates (describing reactive sites) of any given fragment reactions and accesses rate rules (available via RMG-Py) which map reaction templates to corresponding kinetic parameters, of which original data are usually from experiments and ab inito calculations (e.g., CBS-QB3 level). The rate rules are organized in a hierarchical tree of pre-defined functional groups from generic groups on top to specific ones at bottom. More details on the kinetic tree and estimation logic can also be found in RMG-Py’s main paper. 14 To ensure thermodynamic consistency, we select one direction of a reversible reaction to estimate forward rate coefficient kf using rate rules. The rate coefficient of reverse reaction,

12

ACS Paragon Plus Environment

Page 12 of 27

Page 13 of 27 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 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

kr , is computed using the equation below:

kr = kf /Kc

where Kc is the equilibrium constant in concentration units (computable from thermochemistry). The direction we select to use rate rules needs to satisfy the following requirements: • the reaction template of the direction has rate rule data • if both directions have data, choose the one with better quality data (e.g., experimental data is regarded better than quantum calculations, high level of theory is regarded better than low level of theory, etc.) • if both directions have similar quality of data, choose the exothermic direction

Figure 9: Fragment kinetics estimation scheme in AutoFragmentModeling

13

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 49 50 51 52 53 54 55 56 57 58 59 60

III. Case study In this section we show a case study to demonstrate the fragment-based modeling workflow. In addition, by comparing prediction results of a fragment-based model with those of a detailed model by RMG, we were able to assess feasibility of the proposed framework. The target system was purposefully chosen to be large: low-temperature pyrolysis of the C18 H30 molecule phenyldodecane (thereafer PDD), which is one of the largest systems ever modeled using the full-detail representation by RMG-Py. It also has been studied experimentally by several groups. 17–19

fragmentization In a PDD molecule, the carbons on the long alkyl chain experience aromatic influence, which decreases with distance from the benzene ring. The first four carbons, α, β, γ and δ carbon, are most influenced, which leads us to create the first fragment ArCCCCR that contains these carbons shown in Figure 5. Second fragment RCCCCR represents the remaining alkyl carbons that behave as regular carbons on a long alkyl chain. The terminal carbon is separately represented by the third fragment RC.

fragment reaction generation In order to capture the key chemistry of PDD pyrolysis, we selected four elementary reaction types (Figure 10) for the case study. • bond-fission / radical-recombination • β-scission / addition to multiple-bond • hydrogen-abstraction • disproportionation / reverse disproportionation

14

ACS Paragon Plus Environment

Page 14 of 27

Page 15 of 27 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 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 10: Example reactions for four selected reaction types: bond-fission / radicalrecombination, β-scission / addition to multiple-bond, hydrogen-abstraction and disproportionation / reverse disproportionation AutoFragmentModeling provides three options for reaction generation with different modeller input. • option 1: modeller inputs a set of fragments and their reactions (see the input format below). This option gives modellers highest control over what specific reactions to create. AutoFragmentModeling is used to facilitate the model construction process by automatic bookkeeping, duplicate checking, thermochemistry and kinetics estimation. • option 2: modeller inputs a whole set of fragments and desired reaction families. This option leaves modellers with control on which fragments are of interest but AutoFragmentModeling generates all the reactions between each possible pair of fragments and estimates required thermochemistry and kinetics. • option 3: modeller inputs a small set of fragments (usually just those present in feedstock) as well as desired reaction families and lets AutoFragmentModeling expand the initial set to have all the important fragments for the applications. This option achieves automatic reaction network generation for fragments. AutoFragmentModeling relies on the flux-based algorithm 20 to select important fragments and generates cross-reactions as well as estimates required thermochemistry and kinetics.

15

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 49 50 51 52 53 54 55 56 57 58 59 60

In this case study, we used the second option to construct the fragment mechanism.

IV. Results and Discussions A final model was generated using AutoFragmentModeling v1.0.0. It has 41 fragment species and 157 reactions. Compared with the full-detail model generated by RMG, it reduces species and reactions by 4 and 42 times, respectively. We conduct kinetic simulations in Cantera’s 21 homogeneous batch reactor module at 673 K and 350 bar using both models: one is full-detail model by RMG-Py v2.0.0, the other is fragment-based model by AutoFragmentModeling v1.0.0.

feedstock conversion The feedstock composition for both simulations is 100% PDD. In the fragment-based model PDD is not represented as a whole molecule but instead by its fragments: 26.7 mol% ArCCCCR, 46.7 mol% RCCCCR, and 26.7 mol% RC. When calculating feedstock conversion, we choose ArCCCCR’s conversion as an approximation for PDD’s conversion since most ArCCCCR belongs to PDD, especially at low conversions.

Figure 11: Agreement between feedstock conversions predicted by full-detail model (RMG v2.0.0) and fragment-based model (this work) 16

ACS Paragon Plus Environment

Page 16 of 27

Page 17 of 27

As Figure 11 shows, feedstock conversions predicted by these two models agree well with each other when conversion < 0.6, after which the two predictions start to deviate. One major source causing the discrepancy is that ArCCCCR appears in many other products in late conversion of PDD, which makes it no longer an accurate proxy for PDD. In Figure 12, we see a better agreement between moles of ArCCCCR predicted by the fragment-based model and the summed moles of ArCCCCR-containing species (e.g., pentylbenzene, hexylbenzene, ..., PDD) in the full-detail model (Figure 12). Savage et al.’s pyrolysis experiment data 22 is also added, which shows the predictions of both models achieve similar accuracies. 1.0 ArCCCCR moles / initial PDD mole

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 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

0.8 0.6 0.4 0.2

Detailed model Fragment model Savage Experiment

0.0 10-1

100 Time / hr

101

Figure 12: Agreement between moles of ArCCCCR predicted by fragment-based model (this work) and moles of ArCCCCR-containing species (e.g., pentylbenzene, hexylbenzene, ..., PDD) predicted by the full-detail model (RMG v2.0.0)

selected products comparison Molar yields (defined as moles produced per initial PDD mole) of a few selected products are also compared between the two models (Figure 13). Toluene and ethylbenzene are the two major products from PDD pyrolysis, whose trends are both captured by the fragment-based model. It also predicts heavy products, which are typically generated from radical recombination (e.g., bottom right in Figure 13, ArC(CCCR)C(CCCR)Ar) or multiple bond addi17

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 27

tion followed by hydrogen abstraction (e.g., bottom left in Figure 13, ArC(CCCR)CCCR). ArC(CCCR)C(CCCR)Ar from the fragment-based model is compared with all the full-detail molecules such as molecule (a) in Figure 14 that contain ArC(CCCR)C(CCCR)Ar, and ArC(CCCR)CCCR from the fragment-based model is compared with those such as molecule (b) in Figure 14 that contain ArC(CCCR)CCCR. The reason we see the trend difference (fragment model prediction of moles doesn’t go down) for ArC(CCCR)C(CCCR)Ar is that current fragment model doesn’t included decomposition pathways for ArC(CCCR)C(CCCR)Ar.

R

R

R

R

Figure 13: Major light products (toluene and ethylbenzene) and heavy products are predicted by full-detail model via RMG v2.0.0 and fragment-based model via AutoFragmentModeling v1.0.0 with their molar yields compared.

molecular weight distribution Simulations in fragment space predict fragment distributions, which make only limited predictions on the unlumped molecule level, as shown in Figure 11 and 13, since the same 18

ACS Paragon Plus Environment

Page 19 of 27 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 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Figure 14: Two dominant heavy products predicted by the RMG full-detail model. Molecule (a) is the recombination product of two PDD α radicals, while molecule (b) is from PDD α radical addition to double-bond of undecene followed by hydrogen abstraction. fragment appears in multiple molecules. To link fragment distribution back to molecule distribution, we designed a reattachment algorithm that converts the lumped representation (fragments) into the unlumped full-detail representation molecules. However, the fragmentization step has lost connectivity information between fragments in the original feedstock molecule. If fragments are allowed to merge randomly, reattachment would frequently be nonphysical, e.g., reattachment between ArCCCCR and another ArCCCCR since the R-labels of ArCCCCR are originally connected with other types of fragments e.g., RCCCCR, RC. That greatly impairs molecular prediction accuracy. Thus, we modified the fragmentization step (Figure 5) to retain information of fragment connectivity, as shown in Figure 15 by annotating fragments with pairing fragmentization labels; if two fragments are originally connected, one of them has R-label and the other L-label. The new mechanism has slightly larger size (105 species and 1094 reactions), but is still much smaller than the full-detail RMG-Py model. Most of the nonphysical reattachments can be easily prevented with the the new fragmentization by requiring two merging fragments to have compatible label types; reattachment is only allowed when one fragment contains an R-label and the other an L-label (Figure 16). Since any creation of a R-label is always accompanied by creation of a L-label during fragmentization step, R and L-labels are in equal moles initially. The moles of R-labels and L-labels stay constant due to the fact that the labels don’t react. Thus, when reaching 19

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 49 50 51 52 53 54 55 56 57 58 59 60

Figure 15: Fragmentization of PDD with two fragmentization label types: R-label and Llabel

Figure 16: Compatible reattachment example (left panel) and incompatible reattachment example (right panel) reattachment step, we have same amount of compatible fragments in the system. Figure 17 illustrates a systems with three R-labeled fragments and five L-labeled fragments, both of which have concentration of 26 (arbitary unit). A random process makes sure each R-labeled fragment is matched to one or more L-labeled fragments and figures out the amount for each matching. For instance, the first R-labeled fragment has 10 units in total, 5 units of which match the first L-labeled fragment and the other 5 units go to the second L-labeled fragment. One reattachment run provides one realization which comprises of merging each fragment with only one another compatible fragment, generating one possible molecular product distribution. Statistically, the variance of the realizations (by running reattachment multiple times) can be reduced by introducing more pairs of specific fragmentization labels, as it helps retain more information of fragment connectivity in the original molecules. But it also creates more distinct fragments, increasing computation workload during model construction and simulation. To balance accuracy with computation cost, in this case 20

ACS Paragon Plus Environment

Page 20 of 27

Page 21 of 27

Figure 17: Illustration of matching up compatible fragments study we choose to have one pair of fragmentization labels. Its reconstruction effectiveness is shown in Figure 18; the fragment-based model accurately reproduces the products’ molecular weight distribution computed with the detailed model constructed by RMG. 1.0 Detailed Mechansim Fragment Mechanism 0.8 Mole Fraction

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 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

0.6 0.4 0.2 0.0

550

Figure 18: Molecular weight distributions of products of PDD pyrolysis agree reasonably well between full-detail model by RMG v2.0.0 and fragment-based model by AutoFragmentModeling v1.0.0.

V. Conclusion In this paper, we present a new kinetic modeling framework based on fragment, a middle layer entity between molecule and functional groups. We present several innovative fragment 21

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 49 50 51 52 53 54 55 56 57 58 59 60

modeling schemes, which make convenient interfacing with the existing cheminformatics functionalities. For instance, in order to use Benson-style group additivity methods available in RMG-Py, we represent fragments as chemical graphs and create representative molecules to compute thermochemical parameters. Similarly for kinetics parameter estimation, we extract reaction templates from fragment reactions to access rate rule trees in RMG-Py’s database. This framework has several important features for kinetic modeling of systems containing large molecules: 1) the fragment representation of feedstock molecules drastically reduces the number of species and reactions needed to describe the chemical systems compared with full-detail and SOL modeling, 2) the model parameters are estimated via RMG-Py’s existing functionalities, of which the data are either from quantum mechanics calculations or experiments, 3) the reattachment of fragments allows predictions of molecular weight distribution and, if desired, detailed molecular composition. The Fragment representation has superior scalability for modeling systems that contain large molecules over other popular methods such as SOL. A proof-of-concept case study of phenyldodecane pyrolysis has been carried out and reported here. Compared with RMG’s detailed model, the fragment-based model has much smaller size (42 times fewer reactions), but gives similar predictions of feedstock conversion, major product molar yields and molecular weight distribution. This demonstrates the promising potential of the fragment-based framework in modeling large systems which full-detail approaches such as RMG cannot handle well. While the current results are promising, further work is needed to test and develop various possible algorithms for automatically fragmentizing large molecules, periodically re-attaching fragments, and keeping track of changes in the molecular weight distribution, to minimize the loss of information associated with any lumping method. In addition, the feedstock size effect on model prediction is currently only expressed as initial concentration effect; a feedstock molecule with larger size is fragmentized into same set of fragments but with higher concentrations, compared to smaller feedstock molecules. It is worthwhile investigating if

22

ACS Paragon Plus Environment

Page 22 of 27

Page 23 of 27 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 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

such concentration effect is capable of capturing the major size effect. For systems whose feed compositions are too complicated to characterize, this new approach should be coupled with feedstock-reconstruction algorithms to generate representative feed molecules, which would go through fragmentization to provide the initial conditions of fragments for the kinetic simulations. Some (simpler) systems may be modeled more accurately and efficiently using other methods, such as SOL; further work is needed to clearly define when the Fragment approach is the best approach.

VI. Acknowledgments We gratefully acknowledge financial support from the MIT Energy Initiative for funding the research presented here. We also thank Mark J. Goldman for constructive criticism of the manuscript.

Supporting Information Available Here is a list of supporting filed provided with this article: • both the one-labeled and two-labeled fragment mechanisms for phenyl-dodecane pyrolysis are provided in chemkin format • concentration profile predictions of selected products in Figure 13 are provided in csv format • molecular weight distribution predictions in Figure 18 are provided in csv format

References (1) Frenklach, M.; Wang, H.; Yu, C.; Goldenberg, M.; Bowman, C.; Hanson, R.; Davidson, D.; Chang, E.; Smith, G.; Golden, D. Gri-mech-1.2, an optimized detailed chemical 23

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 27

reaction mechanism for methane combustion. Gas Research Institute 1995, (2) Harper, M. R.; Van Geem, K. M.; Pyl, S. P.; Marin, G. B.; Green, W. H. Comprehensive reaction mechanism for n-butanol pyrolysis and combustion. Combust. Flame 2011, 158, 16–41. (3) Van Geem Kevin M.,; Reyniers MarieFrancoise,; Marin Guy B.,; Song Jing,; Green William H.,; Matheu David M., Automatic reaction network generation using RMG for steam cracking of nhexane. AIChE J. 2005, 52, 718–730. (4) Westbrook, C. K.; Pitz, W. J.; Herbinet, O.; Curran, H. J.; Silke, E. J. A comprehensive detailed chemical kinetic reaction mechanism for combustion of n-alkane hydrocarbons from n-octane to n-hexadecane. Combust. Flame 2009, 156, 181–199. (5) How

much

data

is

on

the

Internet?

https://www.quora.com/

How-much-data-is-on-the-Internet, Accessed: 2018-04-16. (6) Weekman, V. W.; Nace, D. M. Kinetics of catalytic cracking selectivity in fixed, moving, and fluid bed reactors. AIChE J. 1970, 16, 397–404. (7) Quann, R. J.; Jaffe, S. B. Structure-oriented lumping: describing the chemistry of complex hydrocarbon mixtures. Ind. Eng. Chem. Res. 1992, 31, 2483–2497. (8) Jaffe, S. B.; Freund, H.; Olmstead, W. N. Extension of Structure-Oriented Lumping to Vacuum Residua. Ind. Eng. Chem. Res. 2005, 44, 9840–9852. (9) Mehl, M.; Pitz, W. J.; Sarathy, S. M.; Westbrook, C. K. Modeling the combustion of high molecular weight fuels by a functional group approach. Int. J. Chem. Kinet. 2012, 44, 257–276. (10) Ranzi, E.; Faravelli, T.; Gaffuri, P.; Sogaro, A. Low-temperature combustion: Automatic generation of primary oxidation reactions and lumping procedures. Combust. Flame 1995, 102, 179–192. 24

ACS Paragon Plus Environment

Page 25 of 27 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 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

(11) Battin-Leclerc, F. Development of kinetic models for the formation and degradation of unsaturated hydrocarbons at high temperature. Phys. Chem. Chem. Phys. 2002, 4, 2072–2078. (12) 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, 790–799. (13) Song, J. Building robust chemical reaction mechanisms : next generation of automatic model construction software. Thesis, Massachusetts Institute of Technology, 2004; Thesis (Ph. D.)–Massachusetts Institute of Technology, Dept. of Chemical Engineering, 2004. (14) Gao, C. W.; Allen, J. W.; Green, W. H.; West, R. H. Reaction Mechanism Generator: Automatic construction of chemical kinetic mechanisms. Comput. Phys. Commun. 2016, 203, 212–225. (15) Neurock, M.; Nigam, A.; Trauth, D.; Klein, M. T. Molecular representation of complex hydrocarbon feedstocks through efficient characterization and stochastic algorithms. Chem. Eng. Sci. 1994, 49, 4153–4177. (16) Ahmad, M. I.; Zhang, N.; Jobson, M. Molecular components-based representation of petroleum fractions. Chem. Eng. Res. Des. 2011, 89, 410–420. (17) Lewan, M. D. Sulphur-radical control on petroleum formation rates. Nature 1998, 391, 164–166. (18) Behar, F.; Lorant, F.; Budzinski, H.; Desavis, E. Thermal Stability of Alkylaromatics in Natural Systems: Kinetics of Thermal Decomposition of Dodecylbenzene. Energy Fuels 2002, 16, 831–841.

25

ACS Paragon Plus Environment

Industrial & Engineering Chemistry Research 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 49 50 51 52 53 54 55 56 57 58 59 60

(19) Savage, P. E.; Klein, M. T. Discrimination between molecular and free-radical models of 1-phenyldodecane pyrolysis. Ind. Eng. Chem. Res. 1987, 26, 374–376. (20) Susnow, R. G.; Dean, A. M.; Green, W. H.; Peczak, P.; Broadbelt, L. J. Rate-Based Construction of Kinetic Models for Complex Systems. J. Phys. Chem. A 1997, 101, 3731–3740. (21) Goodwin, D. G.; Moffat, H. K.; Speth, R. L. Cantera: An Object-oriented Software Toolkit for Chemical Kinetics, Thermodynamics, and Transport Processes. 2016; http: //www.cantera.org, Version 2.2.1. (22) Savage, P. E.; Klein, M. T. Discrimination between molecular and free-radical models of 1-phenyldodecane pyrolysis. Industrial & Engineering Chemistry Research 1987, 26, 374–376.

26

ACS Paragon Plus Environment

Page 26 of 27

Page 27 of 27 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 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

ACS Paragon Plus Environment