Large-Scale Reactive Molecular Dynamics Simulation and Kinetic

Article pubs.acs.org/JPCB

Large-Scale Reactive Molecular Dynamics Simulation and Kinetic Modeling of High-Temperature Pyrolysis of the Gloeocapsomorphaprisca Microfossils Chenyu Zou†,# and Sumathy Raman*,† †

ExxonMobil Research and Engineering Co. (EMRE), 1545 Route 22 East, Annandale, New Jersey 08801-0998, United States

Adri C.T. van Duin‡ ‡

Department of Mechanical Engineering, The Pennsylvania State University, University Park, State College, Pennsylvania 16802, United States ABSTRACT: The ability to predict accurately the thermal conversion of complex carbonaceous materials is of value in both petroleum exploration and refining operations. Modeling the thermal cracking of kerogen under basinal heating conditions improves the predrill prediction of oil and gas yields and quality, thereby ultimately lowering the exploration risk. Modeling the chemical structure and reactivity of asphaltene from petroleum vacuum residues enables prediction of coke formation and properties in refinery processes, thereby lowering operating cost. The chemical structure-chemical yield modeling (CS-CYM) developed by Freund et al.1 is more rigorous, time-consuming, and requires a great deal of chemical insight into reaction network and reaction kinetics. The present work explores the applicability of a more fundamental atomistic simulation using the quantum mechanically based reactive force field to predict the product yield and overall kinetics of decomposition of two biopolymers, namely, the Kukersite and Gutternberg. Reactive molecular dynamics (RMD) simulations were performed on systems consisting of 104 to 105 atoms at different densities and temperatures to derive the overall kinetic parameters and a lumped kinetic model for pyrolysis. The kinetic parameters derived from the simulated pyrolysis of an individual component and the mixture of all four components in Guttenberg reveal the role of cross-talk between the fragments and enhanced reactivity of component A by radicals from other components. The Arrhenius extrapolation of the model yields reasonable prediction for the overall barrier for cracking. Because simulations were run at very high temperature (T > 1500 K) to study cracking within the simulation time of up to 1 ns, it, however, led to the entropically favored ethylene formation as a dominant decomposition route. Future work will focus on evaluating the applicability of accelerated reactive MD approaches to study cracking.

1. INTRODUCTION Natural organic matter (NOM) and anthropogenical organic compounds (AOC) comprise the bulk of organic carbon present in terrestrial and aquatic biospheres. NOM represents portions of degraded and partially degraded biopolymers (e.g., proteins and polysaccharides) subjected to a variety of physical, chemical, and biochemical transformations. Kerogen and black carbon constitute a class of NOM (referenced as carbonaceousNOM or C-NOM) characterized by enriched carbon and aromatic contents. Kerogen, the portion of sedimentary organic matter not soluble in organic solvents, is a complex mixture of biological lipids and polymers, which survive relatively intact during early diagenesis, and of condensed fragments from less stable biochemicals that have been altered by chemical and microbial processes. It is often called a geopolymer, but it has no defined structure composed of specific monomers; rather, it is a heterogeneous macromolecule composed of semirandomly cross-linked core structures that may be aliphatic, naphthenic, © 2014 American Chemical Society

aromatic, and/or heteroatomic in nature. As such, it can be characterized only by chemical or physical properties of the bulk structure or by molecular distributions that are representative of the core and cross-linking moieties and of the chemical nature of the linkages. Behar and Vandenbroucke2 developed 3D representations of Type I, II, and III immature kerogens and Type II kerogens at different stages of maturity. These structural models were synthesized from various analytical data, including elemental analysis, 13C NMR spectroscopy, thermogravimetry, field ionization mass spectrometry, functional group analysis, and pyrolysis. Faulon3 developed computer programs to construct 3D structural models of asphaltene and chemical descriptions of the entire heavy oil fraction. These structural models help to Received: February 24, 2014 Revised: May 12, 2014 Published: May 12, 2014 6302

dx.doi.org/10.1021/jp501925a | J. Phys. Chem. B 2014, 118, 6302−6315

The Journal of Physical Chemistry B

Article

method in reproducing the relevant energetic and hence its applicability in studying the initial steps of a complex pyrolysis phenomenon. However, to the best of our knowledge, the system sizes of the existing simulated pyrolysis studies are limited to several thousands of atoms due to limitations of computational resources. The duplication of the molecular model of the materials is limited to up to 20. However, Ding et al.25 have demonstrated the necessity of sufficient duplications of reactant molecular model in obtaining statistically meaningful results. In a relevant study of the oxidation of a coal model, Castro-Marcano33 performed simulations with a system containing ∼30 000 atoms. Thus, one can be confident in using this method to perform large-scale MD simulations, where the accuracy of the underlying chemical process is validated not by repeating the simulations at different initial conditions but by the statistics ensured by the large-scale systems containing ample duplication of molecules. In this work, we report the MD simulations of thermal decomposition process of Gloeocapsomorphaprisca microfossils using relatively large systems containing 104 ∼ 105atoms. Systems at this scale allow for a meaningful analysis of pyrolysis initiation steps, pyrolysate composition, and kinetic modeling with reasonable statistical accuracy. In addition, systems with a density of nearly 1.0 g/cm3 are studied to address the highpressure effect on the pyrolysis as commonly encountered during geologic burial. Furthermore, a lumped kinetic model is built on the basis of the kinetic information extracted from the simulation. Extrapolation of this model toward a real geological application shows a satisfactory agreement with available data. The next section highlights the molecular model of kerogen together with the methodology of the atomistic simulation employed in this study together with a brief highlight on the features of ReaxFF. In the subsequent section, analysis of the simulation results is presented to provide the elementary steps involved in the initiation stage along with key reaction steps leading to H2, CO, H2O, and small oxygenates in the pyrolysate composition. This is then followed by a lumped composition analysis of products into C1, C2−C5, C6−C14, and C15+ groups and deduction of first-order rate constants from the temporal profiles of the reactant molecules and the activation barrier from the temperature dependence of the first-order rate constant.

calculate thermodynamic parameters that are difficult to measure directly or determine experimentally. Understanding the thermal decomposition process of organic fossil matters is of great interest, because such a process initiates the generation of valuable crude oil and natural gas. Interests in such processes include not only the detailed pyrolysis pathways and product but also the underlying chemical kinetics information. Such information is difficult to obtain accurately due to the geological conditions of low temperature, high pressure, and geological time scales. Typical experimental studies in this area employ mainly the pyrolytic methods4 or the controlled oxidative degradation approaches5 to infer about the detailed molecular structure of the starting material and its degradation process under thermal stress by analyzing the product composition. Though the lowtemperature, high-pressure geo-conditions is not necessarily reproduced in a practical laboratory environment, results from these studies are valuable as they can be fed into other modeling approaches and provide validation for model predictions. One example of such an approach can be found in ref 6, where the proposed mechanistic model for predicting the pyrolysis of coal was based on an experimentally constrained structural model of coal. More recent studies employ a stochastic modeling approach to account for complex structures and properties of fossil matters.1,7 This approach treats the complex fossil material as an ensemble of molecules, of which the elementary reactions and the chemical kinetic parameters are predetermined by experimental observations. Results from this modeling approach support the extrapolation of high-temperature pyrolysis data. Recent decades have seen great improvements in the computational modeling approaches rooted in quantum chemistry and their application to chemical processes. Such first-principle-based modeling approaches provide fundamental information concerning the thermodynamic properties of reacting species, intermediates, products of chemical reactions, and kinetic information on chemical processes. However, practical application of such models to thermal processes of fossil matters is limited by the system size that it can handle at a reasonable computational cost. Accurate atomistic descriptions of the complex structure of such materials may require systems containing thousands of atoms, whereas quantum chemical approaches are limited to a couple of hundreds of atoms. For this reason, we employ the ReaxFF reactive force field developed by van Duin et al.8 to manage the complexity of material structures and the diversity of reaction pathways. The parameters of ReaxFF force field are trained against ab initio data at the density functional theory (DFT) level. Hence, the potential can nearly reproduce the accuracy of the ab initio calculation while maintaining a much lower computational cost. The potential was developed to reproduce fundamental information on materials such as accurate molecular structure and heats of formation of species as well as the reaction energies and barriers of associated chemical processes. This method has been applied to investigate the reaction dynamics information on various chemical systems involving transition metal oxide9 and hydride,10 alloy systems,11,12 as well as combustions13 and catalytic processes.14,15 Prior to this study, the method has also been applied to study the thermal decomposition of various materials including polymers,16−19 fuel type hydrocarbons,20−27 energetic,28 and sedimentary materials.29−32 These simulated pyrolysis studies using the ReaxFF potential have demonstrated the accuracy of the

2. MODELING METHODOLOGY 2.1. Molecular Model. The atomistic model of the fossil materials studied in this paper is proposed by Blokker et al.5 The chemical structures of Gloeocapsomorphaprisca microfossils were characterized using methods including RuO4 chemical degradation, Fourier transform/infrared spectroscopy (FTIR), and flash pyrolysis-GC/MS. The proposed 2D structure of Gutternberg and Kukersite were demonstrated in Figure 10 in the original paper. The Kukersite model (referred to hereafter as Kukersite) contains one large piece of n-alkyl resorcinol based monomer; the Gutternberg model (referred to hereafter as Gutternberg) consists of four different monomers of similar type. The zigzag next to the oxygen in the resorcinol represents the site of polymerization. In the simulation, these sites were terminated by a hydrogen atom, as demonstrated in Figure 1. Polymerization sites in Kukersite were treated in a similar fashion. Besides the resorcinol groups substituted with C8, C10, C15 and C17 aliphatic chains, these models also contain oxygenate functional groups like ethers and ketones. 6303



Article

while the system evolves from reactant to intermediates and then to products. The connectivity-dependent interactions (chemical bond, valence angle, and torsion angle) are bondorder dependent, and their contribution vanishes upon bond dissociation. The potential also includes the nonbonded interactions (van der Waals and Coulomb) calculated between each pair of atoms. The charge-transfer calculation employs the electronegativity equalization method (EEM),36 which is very similar to the charge-transfer ionic potential proposed by Rappe and Goddard,37 as well as another formulizm by Streitz and Mintmire.38 The force field parameters were derived from a large data set of quantum mechanical predictions of heats of formation, reaction energies, and transition state calculations. This parametrization allows nearly quantum accuracy of this empirical method while maintaining orders of magnitude lower computational costs. This feature makes it feasible to simulate decomposition of systems containing up to 106 atoms over a wide range of temperatures, densities, and compositions. This study used the C/H/O ReaxFF parameters reported in ref 12. 2.3. Reactive Molecular Dynamics Simulation. The general MD system preparation procedures follow the description in Figure 2. The 2D representation of the model was rebuilt into 3D using MOLDEN and was optimized at the B3LYP/6-31G level using Gaussian09. Then we built a periodic simulation box containing large numbers of duplications at a relatively low density using PACKMOL. 39 A detailed description of all the systems studied in this work is tabulated in Table 1. All the systems were prepared with the following procedure: First the system underwent the NVT energy minimization at 2K for 1.5 ps to optimize the intermolecular interactions. Thereafter the temperature was increased from 2 K to 300 K in 25 ps and equilibrated at the final temperature for 1.5 ps; after the room temperature equilibration, the simulation box was compressed by rescaling the box dimension to a density of 1.38 g/cm3, corresponding to the density at geological conditions under high pressure. Then the systems were equilibrated at a pressure of 1 bar and 300 K using the isothermal−isobaric (NPT) ensemble. This procedure causes a system expansion to

Figure 1. Atomistic model in 2D representing the n-alkyl resorcinolbased polymer composing the G. prisca microfossils from I (Gutternberg) and II (Kukersite).

The pyrolysis of individual fragments of Gutternberg, namely, A, B, and D, as well as its mixture and Kukersite were studied using reactive molecular dynamics simulations as a function of temperatures and densities. It is worth mentioning that the difference between fragment A and fragment C is mainly the length of the aliphatic side chains. This difference is expected not to affect the initiation process and hence the kinetic properties of the two structures but indeed will generate products with different carbon number distribution and different structures. 2.2. ReaxFF Reactive Force Field. The ReaxFF8 reactive force field potential employs the bond length/bond order and bond order/bond energy relationships.34,35 This allows smooth transitions between bonded systems to nonbonded systems

Figure 2. Protocol of simulated pyrolysis: MD, reactive molecular dynamics; NVT, canonical ensemble, number of atoms (N), volume (V), and temperature (T) are fixed during the simulation. 6304



Article

Table 1. System Details for the Varies NVT −MD Simulations for Different Models models fragment A fragment B fragment D Gutternberg Kukersite

no. of atoms per fragment 95 127 597 938

no. of duplications

densities (g/cm3)

temperature (K)

100 100 150 100 of each 100

0.6/0.99/1.38 0.6/0.95/1.38 0.6/1.0/1.38 0.6/1.0/1.38 1.0

1800−2500 1800−2500 1500−2200 1500−2200 1500−2200

procedures. Lastly, the system underwent NVT simulation at different elevated temperature (1800 K ∼ 2500 K, 100 K as the interval) for fragment A and B; Fragment D was heated isothermally as well at 1500 K ∼ 2200 K. We also studied the simulated pyrolysis of Gutternberg as a mixture of the aforementioned three fragments and Kukersite. Many previous studies, as listed in the introduction, have implemented a pyrolysis process of similar materials and have explored the reaction pathways at a temperature around 2000 K due to a nanosecond time-scale limit. In our study as well, a relatively high temperature range was chosen for the same reason. The duration of all the simulations (250 ps) is orders of magnitude lower than that in a practical experimental condition (several hours to tens of days) of similar investigation or that in

a density of about 1 g/cm3. For the purpose of investigating the density dependence of the activation energy, systems with a density of 0.6 g/cm3 were also prepared following similar

Figure 3. Carbon number distribution of the radicals in three systems at 25 ps (left figure of each panel) and 125 ps (right figure of each panel); panel A, fragment C35H56O4 (T = 2200 K, density = 0.99 g/cm3); panel B, fragment C44H80O3 (T = 2200 K, density = 0.95 g/cm3); panel C, fragment C230H363O22 (T = 1500 K, density = 1.0 g/cm3). 6305



Article

Figure 4. Initiation sites of the three fragments (indicated by arrows): monomer A, benzylic bond in C35H56O4; monomer B, phenoxylether bond in C44H80O3; monomer D, phenoxylether bond in C230H364O23.

Figure 5. Carbon number distribution analysis for alkanes from fragment C35H56O4 at 25 ps (the left figure of each panel) and 125 ps (the right figure of each panel): panel A, density = 1.0 g/cm3; panel B, density = 1.38 g/cm3. For both cases, T = 2200 K.

on the number of carbon atoms in the molecules. These analyses reveal key information on the initiation step of the material under thermal stress as well as on the product composition. The results of the simulation are analyzed below to infer on the chemistry in different stages of pyrolysis. 3.1. Pyrolysis Initiation. For the analysis of the Pyrolysis process of fragments A, B, and D in Gutternberg, we traced the population of radicals being formed at early stages of the process. Figure 3 demonstrates the radical species’ distribution in terms of the carbon number in the structure at early (25 ps) and later (125 ps) stages of decomposition. The radical is determined from the chemical formula of the species: if a species has a general molecular formula of CXHYOZ, then it is

a real geological conditions (millions of years). Nevertheless, it has been demonstrated that qualitative agreement can be achieved between simulations and experiment on the initiation reaction pathways of these fossil materials.29,30 Moreover, the results from the simulated pyrolysis at different densities contribute to our understanding of the thermal degradation of these materials under geological high-pressure conditions.

3. RESULTS AND DISCUSSION Due to system sizes, we focus our analysis of the process in a statistical fashion. Instead of tracing the chemical process of each single molecule, we analyzed the system composition primarily based on the general formula of a molecular class and 6306



Article

Figure 6. Carbon number distribution analysis for product alkenes from fragment C35H56O4 at 25 ps (the left figure of each panel) and 125 ps (the right figure of each panel): panel A, density = 1.0 g/cm3; panel B, density = 1.38 g/cm3. For both, T = 2200 K.

Figure 7. Mechanism of methyl (CH3) radical formation: intramolecular hydrogen transfer and β cleavage (indicated by the black spline; intramolecular hydrogen transfer indicated by the arrow).

initiate pyrolysis in other fragments via bimolecular hydrogen abstraction. The extent of conversion will also depend on density as it influences the diffusion/mobility of the radicals to sustain a long-lived radical chemistry. A further consequence of this change in the initiation step is the possibility to explore the dynamics at a relatively low T, 1500 K as opposed to >1800 K, in Gutternberg and Kukersite, as discussed in the sections below. 3.2. Pyrolysate Composition. The pyrolysis following the initiation step is mainly the formation of benzylic radicals and long aliphatic chain radicals. The thermal depolymerization of the aliphatic chain leads to the formation of a variety of saturated and unsaturated hydrocarbons via the free radical reaction mechanism. This process has been thoroughly studied for different fuel type of hydrocarbons20,25 using the same methodology, and the observed chemical processes agree with the high temperature pyrolysis mechanism. Figures 5 and

recognized as a free radical if it has a fractional value for the “degree of unsaturation”, which is estimated as (2x + 2 − y)/2. Molecular formulas corresponding to the species with relatively high molar fractions in each figure are provided. Although a variety of radicals were formed during the initiation stage, the radicals with high molar fractions suggest that the initiation of three fragments mainly occurs at the benzylic position in fragment A and at the phenoxyether bond in fragments B and D, as indicated by the arrows in Figure 4. A similar initiation mechanism is observed for the simulated pyrolysis of Kukersite due to a comparable molecular structure. These observations agree with the relative strength of the following types of chemical bonds: benzylic C−C bond (C6H5CH2−CH3) > C− C bond in an aliphatic chain (−CH2−CH2−) > phenoxyether bond (C6H5O−CH3). This result further suggests that pyrolysis in Gutternberg could be initiated in fragment B/D by bond scission, and the radicals generated from this step could further 6307



Article

Figure 8. Conversion of total carbon elements (in C35H56O4) after 250 ps NVT simulation into C2H4 (panel A) and CH4 (panel B) at three different densities.

ref 30. In contrast, the conversion of the total carbon element into methane remains relatively low (the highest conversion observed in the simulation is 1.3%) in the density and temperature range investigated. Pyrolysis of sedimentary materials in nature will result in an increase in the degree of unsaturation and the carbon to hydrogen ratio due to successive hydrogen abstraction. In our simulations, hydrogen molecules were formed due to the same mechanism. Figure 9 displays the hydrogen production from

Figure 6 provide the mole fraction distribution of alkanes and alkenes, respectively, in terms of carbon number of the pyrolysates from fragment A. (Species with a number in the bracket is the actual molar number of that species without scaling to the y axis.) A significant amount of ethylene is formed through the β cleavage mechanism, as expected at high temperature. Besides this reaction channel, the following reactions have been observed resulting in the formation of methane and higher alkanes and alkenes: Figure 7 depicted the formation of a methyl (CH3−) radical. First, a C14H29− radical detached from the reactant molecule as a result of the cleavage of a benzylic C−C bond. The radical releases a C2H4 molecule via β cleavage. The resulting C12H25− radical isomerizes twice via intramolecular hydrogen transfer to form a more stable secondary radical. The following β cleavage reaction led to the formation of a CH3− radical and a C11H22 molecule, and the former eventually forms methane via hydrogen abstraction from a reactant (C35H56O4) molecule. The decomposition process of the resulting secondary radical C35H55O4 involves a series of β-scission process starting by releasing a C10H21− radical. After releasing two ethylene molecules, the propagation process was terminated by forming a C6H14 molecule via hydrogen abstraction. Moreover, the formation of light compounds seems to be inhibited at higher density as expected. However, the formation of ethylene from the β cleavage of the radicals seems to dominate at the high temperatures employed in these simulations as it involves an increase in entropy and hence a decrease in the free energy changes of the reaction. Under the simulation conditions, the competing reaction of a radical addition to a double bond leading to molecular growth is entropically disfavored in spite of being an almost barrierless reaction. A rise in reaction temperature not only increases the reaction rate constants but also changes the ratio of low barrier to high barrier reaction rate constants and hence has difficulty in predicting the product selectivity of a low temperature process. Figure 8 provides the conversion analysis of the total carbon elements into ethylene and methane. The results reveal a clear temperature and pressure dependence of the conversion into ethylene molecules. Increasing pressure inhibits the conversion to ethylene at all temperatures ranging from 1800 K to 2500 K. A sudden increase in the conversion into ethylene after 2000 K is observed. Similar dependence was proposed in

Figure 9. Hydrogen formation from fragment A at a density of ∼1.0 g/ cm3 at different temperatures.

fragment A at a density of around 1 g/cm3 as a function of temperature. A dramatic increase in the hydrogen formation is observed when the temperature increased to 2300 K. Therefore, we closely examined the species formed at this temperature in terms of the degree of unsaturation at different simulation time for fragment A and fragment B. The Y axis in Figure 10 indicates the percentage of a particular group of species (identical degree of unsaturation) in all the heavy (skeleton carbon number greater than five) organic species. Notice we exclude the gaseous species in this analysis because they are predominant, as previously demonstrated in Figure 6. Numbers in the brackets indicate the actual population of that particular molecule. At the initial stage of the simulation, the predominant species is the reactant indeed, with a degree of unsaturation of 8 and 5, respectively. As the simulation proceeds to around 60 ps, resorcinol moieties starts to form due the cleavage of C−C bond at the benzylic position in 6308



Article

Figure 10. Heavy (C5+) unsaturated species generated during the simulated pyrolysis at 2300 K of fragment A and fragment B at a density around 1.0 g/cm3.

fragment A, as discussed in the previous section. At the later stage of the pyrolysis, olefins begin to form due to the degradation of the aliphatic chains. Nevertheless the formation of species with high degree of unsaturation is indeed observed in our simulation. A close examination of these species indicates that their formation is mainly due to the addition reactions between two radicals carrying aromatic rings. We do not observe the formation of high unsaturated species from the dehydrogenation of aliphatic chains. The latter process is believed to be among the predominant low-temperature pyrolysis. However, using this methodology, Lummen23 studied the noncatalytic thermal decomposition of methane and observed the formation of a six-member ring. The ring structure survived for a short period of time (14.5 ps) after its formation. The simulation was carried out at a temperature of 3500 K for about 788.5 ps until the observation of this ring structure. This study suggests that the formation of species with a high degree of unsaturation from alkanes can indeed be simulated by using ReaxFF. A similar event has been observed in the atomistic simulation of phenolic resin to carbon process,16 where the formation of fused ring structure was reported at a temperature of 3500 K. However, these species will tend not to survive the high temperature condition necessary for the simulation. Furthermore, the chemical process that leads to the formation of such species is not predominant compared to other high-temperature process such as β-scission. In addition, water formation was observed via an interchain dehydration mechanism. This mechanism has been observed in a previous pyrolysis study of phenolic resin.16 We did not observe a significant amount of H radical or OH radical generation at the temperature range investigated. Thus, water formation through these channels is rarely observed in the simulations of this work. In the case of fragment B, carbon monoxide formation as a function of temperature is displayed in Figure 11, especially at high temperature. This observation is against what has been studied of the carbonyl functional group under thermal stress.29 The author did not observe the thermal degradation of a ketone group leading to the formation of carbon monoxide. The reason for this discrepancy may be the

Figure 11. CO formation at the density of ∼1.0 g/cm3 as a function of temperature for fragment B.

position of the functional group as well as the duration and temperature of the simulation. In that study, the ketone group is substituted to a benzyl moiety, and thus its thermal stability is enhanced. In contrast, the ketone group in our structure is in an aliphatic chain and thus is more prone to thermal degradation. Furthermore, the simulation in that study was carried out for 50 ps, and the temperature is up to 2200 K. The longer simulation time and the higher temperature in our simulation thus validate the observation of CO formation. The thermal degradation of the aliphatic chain substituted by the carbonyl functional group results in the formation of small oxygenate species, as expected. 3.3. Organic Species Mass Distribution. Figure 12 demonstrates the evolution of different products from fragment A at two different temperatures. The products were lumped on the basis of the number of carbon atoms in the structure. This approach is frequently used in experimental studies on kerogen cracking. The pie diagram shows the mass fraction at the end of the simulation. At a relatively low temperature (1900 K) after 250 ps, the mass fraction of condensed matters (C15+) is higher than that of gas-phase compounds, whereas at higher temperature (2500 K), the heavy compounds decompose into light substances at the later stage of the simulation. Similar 6309



Article

Figure 12. Mass fraction evolution of four different organic compounds at 1900 K (panel A) and 2500 K (panel B) of fragment C35H56O4 at 0.99 g/ cm3. The mass fraction at the end of the simulation is shown in the pie diagram in each panel.

Figure 13. Comparison of mass fraction distribution between the simulation and a selected experimental study. Panels I-A and I-B: mass fraction distribution of fragment C35H56O4 at 1900 K after 250 ps at different densities. Panels II-A and II-B: mass fraction distribution of a type I kerogen in an open system and a closed system, respectively.40

trends were observed for fragment B and fragment D. A more detailed classification of the C15+ fraction reveals the generation of the C35+ compound, indicating the ability of the reactive force field to bring out the addition reactions of radicals to unsaturated hydrocarbons at the early stage of the simulation. The latter is thought to be the dominant reaction in lowtemperature cracking of dense carbonaceous systems. However, the extent of contribution of this reaction cannot be estimated

from a reactive molecular dynamics simulation performed at a high enough temperature to accelerate the reaction time to within a nanosecond. In addition, we ran a few modified simulations, including the following: (a) a simulated annealing simulation wherein the system is periodically ramped up to a high T followed by rapid quenching and simulation run at a lower temperature to see the progress of radical chemistry; (b) a low T simulation of fragment A with explicit radicals formed 6310



Article

in the initiation stage of the high-temperature run to improve the prediction on low T product selectivity. Because the system is very dense, neither of these attempts led to any significant conversion of reactants. Figure 13 shows a comparison of the mass fraction distribution from the simulation and one selected experimental study on a different type I kerogen. Panels I-A and I-B of Figure 13 combined with the pie diagram in panel A of Figure 12 reveal the pressure effect on the simulated pyrolysis. The highdensity high-pressure condition will inhibit the degradation process of the heavy compound, resulting in a decrease in the formation of the light compound. It can be seen from panel I-B in Figure 13 that there is no methane formation at the end of the simulation at 1900 K, 1.38 g/cm3, which is the highest density addressed in this study. This observation is consistent with the pressure effect on the formation of ethylene, as discussed in the previous session. The experimental data are obtained from the pyrolysis study at several hundred degrees for tens of days. Despite the distinction of time and temperature scales, qualitative agreement has been achieved upon the relative ratio of different mass compounds, as demonstrated in Figure 13. A detailed classification of the pyrolysates of Gutternberg and Kukersite at a density of 1 g/cm3 based on carbon number is provided in Figure 14. We observe fast degradation of the heavy compound (C30+) into the light compound at a temperature above 1800 K in both cases.

Figure 15. Temperature dependence of the time evolution of the relative concentration of reactant C35H56O4 at a density of ∼1.0 g/cm3.

K for all the species to degrade within 250 ps.The temporal profile of reactants at a given temperature enabled the determination of the rate constant with the assumption of a first order decay process. The temperature dependency of the rate constant further allows for the calculation of apparent activation energy for thermal cracking. Similar analysis has been done on all fragments of Gutternberg and Kukersite. The results are summarized in Table 2. Table 2. Arrhenius Parameters Predicted from the Simulation and Comparison with Selected Experimental Data system monomer

dens. (g/cm3)

fragment A

0.60 0.99 1.38

fragment B

0.60 0.95 1.38

fragment D

1.0

Kukersite

1.0

A (s−1) 5.7 9.2 2.1 1.3 2.1 1.2 1.3 5.1 4.3 2.1 1.8 5.3 4.7 5.6 5.6

× × × × × × × × × × × × × × ×

1015a 1013b 1016a 1014b 1016a 1014b 1015a 1014b 1016a 1015b 1016a 1015b 1015a 1014b 1014

Ea (kcal/mol) 57 37 62 39 64 41 45 41 50 46 57 56 44 41 37

± ± ± ± ± ± ± ± ± ± ± ± ± ± ±

2a 2b 3a 2b 3a 2b 2a 3b 2a 1b 3a 2b 3a 2b 2

a

: indicating that the data was estimated from the simulations of a system containing single species. b: indicating that the data was estimated from the simulations of a system containing multiple species.

Simulations of fragments A, B, and D of Gutternberg were studied as an isolated fragment and as a mixture of all three components at three different densities, namely, 0.6, 1.0, and 1.38 g/cm3. It is worthwhile to note from Table 1 that the size of the system varies almost an order of magnitude between simulations of an isolated fragment and a mixture. The calculated activation energy of the fragment A from the simulation of Gutternberg is significant lower compared to that from simulations of the isolated fragment itself. This drop in the activation energy indicates a different initiation mechanism of such species when present together with other molecules. In Gutternberg, fragment A is initiated by the bimolecular

Figure 14. Pyrolysates mass distribution of different compound of Gutternberg (panel A) and Kukersite (panel B) at a density of 1 g/cm3 after 250 ps.

3.4. Kinetic Modeling. The temperature dependence of the thermal decomposition of the reactant is demonstrated in Figure 15. This example shows the time evolution of fragment A in Gutternberg under thermal stress at a density of ∼0.95 g/ cm3. It is observed that at the lowest temperature (1800 K) investigated in this work, only 20% of the reactant had reacted within 250 ps. The temperature needed to be raised above 2200 6311



Article

magnitude difference in rate constants, the estimated values are within this uncertainty in Ea. Actual rate constants used in the model are provided in Table 3.

hydrogen abstraction reaction of its own molecule with the radicals formed from the thermal initiation of other molecules (fragment B and D in this case) with lower activation energy, whereas in the system with isolated fragment A, bond scission due to thermal stress initiates the radical formation. Wang et. al27 demonstrated a similar effect of fuel additives such as a diethyl ether molecular on the pyrolysis of n-decane. Furthermore, in both the simulations, the rate constant decreases with an increase in density, as is expected from reduced mobility of the dense system. From our experience, we found Arrhenius pre-exponential factors generally differ by an order of magnitude. These results are likely because NVT simulations have a constant increase in pressure as simulation continues, whereas experimental results are often the result of a constant pressure setup. Alternatively, the noise in the initial temporal profile where the system decays rapidly could introduce significant offset to the frequency factor. One of our goals is to evaluate the applicability of atomistic simulations in developing a lumped kinetic model that can describe the time evolution of changing composition under the pyrolysis condition. Such an ability provides unique strength to this tool as it avoids the laborious and time-consuming effort of constructing the detailed reaction mechanism involving thousands of reactions among hundreds of species and, in addition, estimates the rate parameters for all the reactions in CS-CYM approach. The model is based on the kinetic data predicted for fragment A as a representative model structure. Five lumped compounds are considered in the model: C1, C2− C5, C6−C14, C15−C34, C35, and C35+. The five reactions considered in the model are demonstrated in Figure 16.

Table 3. Kinetic Parameters for the Five-Step Lumped Kinetic Modela

a

k1 (s−1)

A1

A2

A3

A4

2.1 × 1016 × exp(62/RT) k2 0.4 × k1 k3 0.2 × k2 k4 0.4 × k3 k5 0.2 × k4

0.1463 A5 0.164 A8 0.1 A11 0.1

0.2432 A6 0.2726 A9 0.3 A12 0.9

0.5025 A7 0.5633 A10 0.6

0.108

R: the universal gas constant; T: temperature.

This assumption could be validated and corrected through further investigation of the simulation. However, comparisons between the model and the simulations (Figure 17) indicate that errors due to these assumptions and approximations are within acceptable range. Furthermore, assuming a geological thermal burial history of approximately 3 °C/million years (which is in reasonable range of the proposed hypothetical investigations41 and the references therein), our model (Figure 18) predicted that the model kerogen molecule can sustain a thermal stress up to ∼150 °C and the oil yield window is between 170 to 210 °C, which is within reasonable qualitative agreement with the predictions from experiment42 (the temperature of the onset of hydrocarbon generation is lower compared to this study, corresponding to the lower activation energy of the relevant type of kerogen). Furthermore, the temperature corresponding to the maximum oil generation (C6−C14) is determined to be 196.3 °C and will shift to a higher value at a faster heating rate within the framework of our model (198.2 °C at 4 °C/million year and 199.8 °C at 5 °C/ million year). This trend agrees with the general observations in the experimental study.43

Figure 16. Simplest five-step first-order model constructed to predict the thermal decomposition of fragment C35H56O4.

4. CONCLUSIONS The ReaxFF MD simulations were performed to investigate the high-temperature thermal decomposition process of Gloeocapsomorphaprisca microfossils. The method can accurately predict the initiation stage of such species under thermal stress. The main chemical process observed in these systems during the simulated decomposition includes: • Scission of thephenoxy-ether group and the cleavage of the benzylic carbon−carbon bond yielding aliphatic chain radicals. • Defunctionalization of the hydroxyl group in the resorcinol groups leading to the generation of water molecules via interchain dehydration. • Depolymerization of aliphatic chain yielding various gas phase hydrocarbons and condensed compound. • Production of ethylene via successive β cleavage is observed to be predominant at temperature range investigated. • Besides, intramolecular hydrogen transfer leading to the formation of methane and higher alkenes are observed.

The mass stoichiometry coefficient for the first step is estimated using the data from the initial 20 ps, where the pyrolysis is approximately dominated by this step. The coefficients of the second step are estimated in such a way that the ratios among the coefficients are identical to those of the first step; the coefficient for the remaining steps are purely estimated from approximation, as all the step progress simultaneously at the later stage of the pyrolysis, and hence, an accurate estimation from simulation becomes difficult. The rate constant k1 of the first step is estimated from the Arrhenius parametrization; however, k2-k5 is estimated on the basis of the assumption that these rate constants should be in a descending order. As a rough check, we made NVT simulations on hexadecane, one of the saturated components of C15−C34 fraction, and we arrived at a value almost close to that of k1. Because the simulation of low T (∼373 to 423 K) cracking chemistry is more sensitive to activation barrier and a difference of 2 kcal/mol in Ea at 423 K could lead to an order of 6312



Article

Figure 17. Comparison between the results from the simulation (Panel A) and the prediction of the model (Panel B).

complex heterogeneous structure of these geological fossil materials makes the ab initial modeling approach unfeasible. The ReaxFF reactive force field method provides an alternate approach in studying these systems. The method roots in the fundamental chemistry (empirical parameters were trained against information concerning the potential energy surface of a specific system) and provide chemical dynamics information at relatively high temperature. Comparison between these RD simulations and experimental results contributes to the understanding of the thermal process of these materials in nature. Recently, there are several studies concerning the feasibility of accelerated MD simulations: parallel replica dynamics44 and temperature accelerated dynamics.45 These studies aim at increasing the time-scale limits of a classical MD simulation and hence make plausible the simulation of a low-temperature process within practical simulation time. A recent study from our group46 introduced algorithms to couple parallel replica dynamics with the ReaxFF potential and implemented it in LAMMPS. This has enabled simulation of the pyrolysis of 1heptene at a temperature as low as 1350 K and brought the simulation time scale to microseconds. Progress in this area is believed to enhance the accuracy and reliability of RMD simulations addressing the problems of simulating the chemically relevant low-temperature process.

Figure 18. Model prediction of the temperature dependence of different compounds (assuming a thermal history of 3 °C/my).

• Meanwhile, higher alkanes could also form via hydrogen abstraction of an alkyl radical. • Species with a high degree of unsaturation formed via successive hydrogen elimination is not favored at the temperature range investigated in our study. These observations agree with previous studies of similar fossil materials. Furthermore, kinetic modeling based on the Arrhenius relation reveals a positively correlated relation between the activation energy and the system density. A similar trend is reported in the study of the pyrolysis of ndodecane24 at a similar density range. Finally, a simplified kinetic model based on the results of RMD simulation was built to predict the time rate of change of different compounds under thermal stress. Despite the simplification and approximation, the model has good agreement with the results from simulations. Furthermore, the extrapolation of the model under geological conditions yields reasonable predictions of the thermal stability of kerogen and the onset generation of crude oil. The prediction of the maturation process of kerogen is difficult because the process occurs in a time and temperature scale beyond the capability of practical experimental study; furthermore, inorganic materials such as minerals and water could have a more complicated effect on the thermal decomposition process. However, in the meantime, the

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Present Address

# Department of Mechanical Engineering, The Pennsylvania State University, University Park, State College, Pennsylvania 16802, United States.

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS

■

REFERENCES

We thank Howard Freund and Clifford Walters for their insightful comments and constructive suggestions.

(1) Freund, H.; Walters, C. C.; Kelemen, S. R.; Siskin, M.; Gorbaty, M. L.; Curry, D. J.; Bence, A. E. Predicting Oil and Gas Compositional Yields via Chemical Structure-Chemical Yield Modeling (CS-CYM):

6313



Article

Part 1 − Concepts and Implementation. Org. Geochem. 2007, 38, 288− 305. (2) Behar, F.; Vandenbroucke, M. Chemical Modeling of Kerogens. Org. Geochem. 1987, 11, 15−24. (3) Faulon, J. L.; Vandenbroucke, M.; Drappier, J. M.; Behar, F.; Romero, M. 3D Chemical-Model for Geological Macromolecules. Org. Geochem. 1990, 16, 981−993. (4) Jackson, K. J.; Burnham, A. K.; Braun, R. L.; Knauss, K. G. Temperature and Pressure Dependence of N-Hexadecane Cracking. Org. Geochem. 1995, 23, 941−953. (5) Blokker, P.; van Bergen, P.; Pancost, R.; Collinson, M. E.; de Leeuw, J. W.; Damste, J. S. S. The Chemical Structure of Gloeocapsomorpha prisca Microfossils: Implications for Their Origin. Geochim. Cosmochim. Acta 2001, 65, 885−900. (6) Solomon, P. R.; Hamblen, D. G.; Carangelo, R. M.; Serio, M. A.; Deshpande, G. V. General-Model of Coal Devolatilization. Energy Fuels 1988, 2, 405−422. (7) Kelemen, S. R.; Freund, H.; Siskin, M.; Carry, D. J.; Xiao, Y.; Olmstead, W. N.; Gorbaty, M. L.; Bence, A. E. Chemical Structual and Composition Yields Model for Predicting Hydrocarbon Thermolysis Products. U.S. Patent No. US 7344889 B2, Mar 18, 2008. (8) van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard, W. A., III ReaxFF: A Reactive Force Field for Hydrocarbons. J. Phys. Chem. A 2001, 105, 9396−9409. (9) Aryanpour, M.; van Duin, A. C. T.; Kubicki, J. D. Development of a Reactive Force Field for Iron-Oxyhydroxide Systems. J. Phys. Chem. A 2010, 114, 6298−6307. (10) Zou, C. Y.; van Duin, A. C. T.; Sorescu, D. C. Theoretical Investigation of Hydrogen Adsorption and Dissociation on Iron and Iron Carbide Surfaces Using the ReaxFF Reactive Force Field Method. Top. Catal. 2012, 55, 391−401. (11) Vasenkov, A.; Newsome, D.; Verners, O.; Russo, M. F.; Zaharieva, R.; van Duin, A. C. T. Reactive Molecular Dynamics Study of Mo-Based Alloys under High-Pressure, High-Temperature Conditions. J. Appl. Phys. 2012, 112, 013511−013523. (12) Shin, Y. K.; Kwak, H.; Zou, C. Y.; Vasenkov, A. V.; van Duin, A. C. T. Development and Validation of a ReaxFF Reactive Force Field for Fe/Al/Ni Alloys: Molecular Dynamics Study of Elastic Constants, Diffusion, and Segregation. J. Phys. Chem. A 2012, 116, 12163−12174. (13) Chenoweth, K.; van Duin, A. C. T.; Goddard, W. A., III ReaxFF Reactive Force Field for Molecular Dynamics Simulations of Hydrocarbon Oxidation. J. Phys. Chem. A 2008, 112, 1040−1053. (14) Mueller, J. E.; van Duin, A. C. T.; Goddard, W. A., III Development and Validation of ReaxFF Reactive Force Field for Hydrocarbon Chemistry Catalyzed by Nickel. J. Phys. Chem. C 2010, 114, 4939−4949. (15) Zou, C. Y.; Van Duin, A. Investigation of Complex Iron Surface Catalytic Chemistry Using the ReaxFF Reactive Force Field Method. JOM 2012, 64, 1426−1437. (16) Jiang, D.-e.; van Duin, A. C. T.; Goddard, W. A., III; Dai, S. Simulating the Initial Stage of Phenolic Resin Carbonization via the ReaxFF Reactive Force Field. J. Phys. Chem. A 2009, 113, 6891−6894. (17) Chenoweth, K.; Cheung, S.; van Duin, A. C. T.; Goddard, W. A.; Kober, E. M. Simulations on the Thermal Decomposition of a Poly(Dimethylsiloxane) Polymer Using the ReaxFF Reactive Force Field. J. Am. Chem. Soc. 2005, 127, 7192−7202. (18) Desai, T. G.; Lawson, J. W.; Keblinski, P. Modeling Initial Stage of Phenolic Pyrolysis: Graphitic Precursor Formation and Interfacial Effects. Polymer 2011, 52, 577−585. (19) Saha, B.; Schatz, G. C. Carbonization in Polyacrylonitrile (Pan) Based Carbon Fibers Studied by ReaxFF Molecular Dynamics Simulations. J. Phys. Chem. B 2012, 116, 4684−4692. (20) Castro-Marcano, F.; van Duin, A. C. T. Comparison of Thermal and Catalytic Cracking of 1-Heptene from ReaxFF Reactive Molecular Dynamics Simulations. Combust. Flame 2013, 160, 766−775. (21) Liu, L.; Bai, C.; Sun, H.; Goddard, W. A., III Mechanism and Kinetics for the Initial Steps of Pyrolysis and Combustion of 1,6Dicyclopropane-2,4-Hexyne from ReaxFF Reactive Dynamics. J. Phys. Chem. A 2011, 115, 4941−4950.

(22) Chenoweth, K.; van Duin, A. C. T.; Dasgupta, S.; Goddard, W. A., III Initiation Mechanisms and Kinetics of Pyrolysis and Combustion of Jp-10 Hydrocarbon Jet Fuel. J. Phys. Chem. A 2009, 113, 1740−1746. (23) Lummen, N. ReaxFF-Molecular Dynamics Simulations of NonOxidative and Non-Catalyzed Thermal Decomposition of Methane at High Temperatures. Phys. Chem. Chem. Phys. 2010, 12, 7883−7893. (24) Wang, Q. D.; Wang, J. B.; Li, J. Q.; Tan, N. X.; Li, X. Y. Reactive Molecular Dynamics Simulation and Chemical Kinetic Modeling of Pyrolysis and Combustion of N-Dodecane. Combust. Flame 2011, 158, 217−226. (25) Ding, J.; Zhang, L.; Zhang, Y.; Han, K.-L. A Reactive Molecular Dynamics Study of N-Heptane Pyrolysis at High Temperature. J. Phys. Chem. A 2013, 117, 3266−3278. (26) Leininger, J. P.; Minot, C.; Lorant, F. Two Theoretical Simulations of Hydrocarbons Thermal Cracking: Reactive Force Field and Density Functional Calculations. TheochemJ. Mol. Struct. 2008, 852, 62−70. (27) Wang, Q.-D.; Hua, X.-X.; Cheng, X.-M.; Li, J.-Q.; Li, X.-Y. Effects of Fuel Additives on the Thermal Cracking of N-Decane from Reactive Molecular Dynamics. J. Phys. Chem. A 2012, 116, 3794−3801. (28) Zhou, T. T.; Huang, F. L. Effects of Defects on Thermal Decomposition of Hmx via ReaxFF Molecular Dynamics Simulations. J. Phys. Chem. B 2011, 115, 278−287. (29) Salmon, E.; van Duin, A. C. T.; Lorant, F.; Marquaire, P. M.; Goddard, W. A. Early Maturation Processes in Coal. Part 2: Reactive Dynamics Simulations Using the ReaxFF Reactive Force Field on Morwell Brown Coal Structures. Org. Geochem. 2009, 40, 1195−1209. (30) Salmon, E.; van Duin, A. C. T.; Lorant, F.; Marquaire, P. M.; Goddard, W. A. Thermal Decomposition Process in Algaenan of Botryococcus Braunii Race L. Part 2: Molecular Dynamics Simulations Using the ReaxFF Reactive Force Field. Org. Geochem. 2009, 40, 416− 427. (31) Yan, G.; Zhang, Z.; Yan, K. Reactive Molecular Dynamics Simulations of the Initial Stage of Brown Coal Oxidation at High Temperatures. Mol. Phys. 2013, 111, 147−156. (32) Zheng, M.; Li, X.; Liu, J.; Guo, L. Initial Chemical Reaction Simulation of Coal Pyrolysis via ReaxFF Molecular Dynamics. Energy Fuels 2013, 27, 2942−2951. (33) Castro-Marcano, F.; Kamat, A. M.; Russo, M. F., Jr.; van Duin, A. C. T.; Mathews, J. P. Combustion of an Illinois No. 6 Coal Char Simulated Using an Atomistic Char Representation and the ReaxFF Reactive Force Field. Combust. Flame 2012, 159, 1272−1285. (34) Tersoff, J. New Empirical-Approach for the Structure and Energy of Covalent Systems. Phys. Rev. B 1988, 37, 6991−7000. (35) Brenner, D. W. Empirical Potential for Hydrocarbons for Use in Simulating the Chemical Vapor Deposition of Diamond Films. Phys. Rev. B 1990, 42, 9458−9471. (36) Mortier, W. J.; Ghosh, S. K.; Shankar, S. Electronegativity Equalization Method for the Calculation of Atomic Charges in Molecules. J. Am. Chem. Soc. 1986, 108, 4315−4320. (37) Rappe, A. K.; Goddard, W. A. Charge Equilibration for Molecular-Dynamics Simulations. J. Phys. Chem. 1991, 95, 3358−3363. (38) Streitz, F. H.; Mintmire, J. W. Electrostatic Potentials for MetalOxide Surfaces and Interfaces. Phys. Rev. B 1994, 50, 11996−12003. (39) Martinez, L.; Andrade, R.; Birgin, E. G.; Martinez, J. M. Packmol: A Package for Building Initial Configurations for Molecular Dynamics Simulations. J. Comput. Chem. 2009, 30, 2157−2164. (40) Behar, F.; Vandenbroucke, M.; Tang, Y.; Marquis, F.; Espitalie, J. Thermal Cracking of Kerogen in Open and Closed Systems: Determination of Kinetic Parameters and Stoichiometric Coefficients for Oil and Gas Generation. Org. Geochem. 1997, 26, 321−339. (41) Lewan, M. D.; Ruble, T. E. Comparison of Petroleum Generation Kinetics by Isothermal Hydrous and Nonisothermal Open-System Pyrolysis. Org. Geochem. 2002, 33, 1457−1475. (42) Behar, F.; Roy, S.; Jarvie, D. Artificial Maturation of a Type I Kerogen in Closed System: Mass Balance and Kinetic Modelling. Org. Geochem. 2010, 41, 1235−1247. 6314



Article

(43) Guo, L.; Xiao, X.; Tian, H.; Song, Z. Distinguishing Gases Derived from Oil Cracking and Kerogen Maturation: Insights from Laboratory Pyrolysis Experiments. Org. Geochem. 2009, 40, 1074− 1084. (44) Voter, A. F. Parallel Replica Method for Dynamics of Infrequent Events. Phys. Rev. B 1998, 57, 13985−13988. (45) Sorensen, M. R.; Voter, A. F. Temperature-Accelerated Dynamics for Simulation of Infrequent Events. J. Chem. Phys. 2000, 112, 9599−9606. (46) Joshi, K. L.; Raman, S.; van Duin, A. C. T. Connectivity-Based Parallel Replica Dynamics for Chemically Reactive Systems: From Femtoseconds to Microseconds. J. Phys. Chem. Lett. 2013, 4, 3792− 3797.

6315


Large-Scale Reactive Molecular Dynamics Simulation and Kinetic

Recommend Documents