Development of a Reactive Force Field for ... - ACS Publications

Subscriber access provided by UNIVERSITY OF ADELAIDE LIBRARIES

Article

Development of a ReaxFF for Hydrocarbons and Application to Iso-octane Thermal Decomposition Endong Wang, Junxia Ding, Zongjin Qu, and Keli Han Energy Fuels, Just Accepted Manuscript • DOI: 10.1021/acs.energyfuels.7b03452 • Publication Date (Web): 14 Dec 2017 Downloaded from http://pubs.acs.org on December 14, 2017

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 free 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 accessible to all readers and 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.

Energy & Fuels 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 22 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

Energy & Fuels

Development of a ReaxFF for Hydrocarbons and Application to Iso-octane Thermal Decomposition Endong Wang,†‡ Junxia Ding,† Zongjin Qu, † Keli Han†* †

State Key Laboratory of Molecular Reaction Dynamics, Dalian Institute of Chemical Physics, Chinese Academy of Sciences, 457 Zhongshan Road, Dalian 116023 (China) ‡

University of Chinese Academy of Sciences, Beijing 100049 (China)

Abstract Reactive force field (Reaxff) is a powerful method, which employs bond order/bond length formulism in order to describe the bond breaking and bond re-formation. In this work, a modification to the bond order formula was made and Reaxff parameters were re-optimized. The underlying idea of these modifications is to improve the energy gradient. Better agreements of the bond dissociation potential curves with the QM curves were obtained based on the aforementioned changes. Reaxff simulation was carried out to gain the understandings of the isooctane pyrolysis. The new parameters based apparent rate constants of the iso-octane decomposition fit well with the experimental results. The simulation results are in agreement with the existing experimental results. A maximum of C2-hydrocarbons were found to have the largest percentage. And distribution of the iso-octane decomposition pathway was illustrated.

ACS Paragon Plus Environment

1

Energy & Fuels 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 22

1 Introduction The real jet fuels, diesel or gasoline often contain up to hundreds of distinct chemical components. A study on the combustion or pyrolysis of these real fuels is not feasible and may be out of reach for the present laboratorial devices. To tackle with this, researchers often use a surrogate fuel to control the system within a manageable size.1 A common approach is to categorize all the surrogate fuels into as many types as needed based on their structural characters. The commonly used classes to distinguish the surrogates are chain-alkanes, branchedalkanes, aromatics, polycyclic-alkanes, olefins, naphthenes, and oxygenated hydrocarbons.2 Many experimental and computational works regarding mechanisms of surrogate fuel combustion and pyrolysis have been reported.2 Among all the surrogate fuels, iso-octane is a valuable branched-hydrocarbon because it is one of the primary reference fuel components used in the knock rating of gasoline and it is also viewed as a single component surrogate for gasoline in homogenous charge ignition engines engine modeling. Many experimental works have already focused on the iso-octane pyrolysis and combustion. Davis and Law studied the flame speeds of iso-octane.3 Davidson et al. also have done a lot of works on iso-octane, including obtained the shock-tube-derived ignition delay time data,4, 5 measured species time-histories for OH during iso-octane oxidation,6 and detected the methyl radical concentration time-histories during the oxidation and pyrolysis using UV laser absorption7. Besides, many kinetic works in iso-octane have also been reported. Slavinskaya et al. presented a simplified chemical kinetic work to describe the combustion of n-heptane, isooctane, and their mixtures in air.8 Curran et al. developed a mechanism of iso-octane oxidation in a jet-stirred reactor, flow reactors, shock tubes and in a motored engine.9 Some other chemical kinetic models were also developed for surrogate fuels of Jet A with iso-octane as a


2

Page 3 of 22 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

Energy & Fuels

component.10-12 Dooley et al. presented an implicit methodology, usually called first generation surrogate, to formulate jet aviation fuel surrogate.1 Then, a modification to the 1st generation surrogate was made and was found better agreement in predicting the formation of CO and CO2. Also improvements were seen in predicting species histories from flow reactor oxidation experiments.13 However, the modified 1st generation surrogate is still unable to predict concentration profiles of some species.14 And it was argued that the failures in predicting those species profiles are due to the empirical rules based estimated rate constants in the kinetic models. Thus, we try to study the iso-octane pyrolysis using a different method in this work. In recent years, Reaxff15, 16 has received more and more attentions. This method has been successfully used in many research areas. A thorough review article has been published recently.17 Reaxff was firstly proposed by Goddard and van Duin et al.16 and later improved by Chenoweth et al.15 In this method, bond order is used to bridge the gap between the quantum mechanics and classical mechanics, thus allowing the system to evolve in a reactive manner without the expensive quantum mechanical (QM) costs. Notably, this method is initially developed for hydrocarbons17 and has been widely used in hydrocarbon combustion and pyrolysis.15, 18-22 In addition, this method has also been used in the modelling of coal pyrolysis,2327

carbonaceous materials formation,28, 29 cellulose pyrolysis,30 crystal simulation31-35 and other

condensed state systems. Besides, Reaxff has also been used to simulate catalysis system.36 This method has been widely used in various fields. However, people also pointed out that there still exist some deficiencies in the widely used Chenoweth C/H/O parameter set37 and the results could be improved by tuning Reaxff parameters.38 Thus, in this work, efforts were devoted to explore ways to improve the accuracy of Reaxff. In out attempts, the bond length range with respect to bond order was enlarged to better fit the nature of chemical bonding and


3

Energy & Fuels 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 22

the previous research. In addition, the bond order to bond distance formula was modified to smooth the double bond and triple bond dissociation energy curve. Based on these changes related to bond order, parameters used in simulations were re-optimized using QCISD level of theory of results as QM training set. Detailed parameters were supplied in the electronic supplementary information (ESI). The new parameters were used to simulate the thermal decomposition of iso-octane. The rest of the manuscript was organized as follows. Firstly, theoretical method and computational detail is given. Secondly, results and discussions are supplied in the results and discussion section. And then this manuscript is summarized. 2 Theoretical method and computational details In Reaxff, the total energy can be divided into several terms as listed in the following formula. Esystem = Ebond + Eval + Etors + Econj + EvdWaals + ECoulomb + EH −bond + Elp + Eover + Eunder + E pen + Ecoa + EC 2

In the above formula, those terms represent bond energy, angle energy, torsion angle energy, four body conjugation energy, van der Waals energy, coulomb energy, hydrogen bond energy, lone pair energy, over-coordination penalty energy, under-coordination penalty energy, three body conjugation term and energy correction for C2-molecules respectively. Reaxff is essentially an empirical force field method employing energy-bond order-coordinates relationship which allows energy to change smoothly during bond breaking and bond formation. The heart of the method, to the best of our knowledge, lies in bond order formulism. In nonreactive molecular dynamics, two atoms bonded together can only oscillate back and forth within a range of bond length making the molecular topology remain unchanged. Instead, in Reaxff, specific energy terms were added to cope with the bond-breaking and bond-formation. The add of these energy terms is accomplished through the introduction of bond order as the basis.


4

Page 5 of 22 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

Energy & Fuels

Additionally, all the energy terms except coulomb energy and van der Waals energy lean on the bond order explicitly making these terms change along bond length implicitly. This is another distinct difference compared with the non-reactive force field. More technical details about this method are well summarized elsewhere.15-17 We modified the bond order formula and re-optimized the C/H parameters to better reproduce potential energy curves from QM calculations. The details will be discussed below. Several typical potential energy curves using the newly optimized parameters are shown below and the rest of the potential energy curves can be found in the supporing information. Parameters involving bond dissociation energy curves, bond bend energy curves, bond torsion energy curve, charge and some important reaction paths, like β-scission reactions and H abstraction reactions were re-optimized. All the QM calculations are performed by Gaussian09-d01 using quadratic configuration interaction method (QCISD) and the Pople 6-311+G(d,p) basis set. Reaxff molecular dynamics simulations were done using velocity Verlet integrator with a time step of 0.1 fs. Increasing density could save amounts of computational costs, thus 20 isooctane molecules were put in a 30Å × 30Å × 30Å box with a density of 0.14 g/cm3. The pressure ranges from around 3.063MPa to 30.63 MPa with temperature increases from 300 K to 3000 K. The pressure may be a bit high for typical combustion or pyrolysis experiments, but can easily be achieved in many reactive processes, including rocket engines and shock-induced reactions.18 Canonical ensemble was used in all simulations. Temperature was controlled by Berendsen thermostat using a temperature damping constant of 10 ps. Initially, 10 ps simulation at 300 K was conducted to reach pre-equilibrium. And then temperature was increased from 300 K to 3000 K in another 10 ps. Finally, simulations were run for 2 ns to collect data. All the simulations were performed using reax/c package in freely available code “Large-scale


5

Energy & Fuels 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

Atomic/Molecular

Massively Parallel

Simulator”

Page 6 of 22

(LAMMPS)

from

Sandia

National

Laboratories.39 3 Results and discussion 3.1 Bond order formula modification and force field parameters optimization The aim of Reaxff is to gain understanding on the reactive processes of the complex system, thus the potential energy curves depicted by Reaxff are innegligible. According to the inter-atomic distance, all the potential energy curves in Reaxff could be divided into two types: the first includes those describing bond breaking/forming and the second are those describing the equilibrium structure. Because the advantage of Reaxff is to allow the bond dissociation, more attentions should be paid on the first type on the basis of maintaining a reasonable equilibrium structure. Potential energy curves of bond dissociation/forming can be evaluated from two aspects: the first is the bond dissociation energy and the other is how closely the Reaxff energy curve matches with QM curve during the changing of bond length. Because the bond dissociation energy is important to the final results and the fitness of the Reaxff curve can influence the kinetics of the simulation. Based on some trials, the parameter set proposed in 200815 could be further improved in the above two aspects, as illustrated by Figure 1. It could be easily seen in Figure 1 that deviations exist along with the increasing of the bond length, as well as the bond dissociation energies. Thus, attempts were made to improve the Reaxff potential energy curves from above two respects.


6

Page 7 of 22 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

Energy & Fuels

(a)

(b)

(c)

Figure 1 Reaxff C-C bond dissociation energy curve of (a) isobutene, (b) ethylene, (c) acetylene using the parameters in reference paper.15 Unless otherwise specified, the old ff reported in this paper was obtained from reference paper,15 the new ff below referred to the newly optimized parameters. As stated previously, the bond order concept is the core stone of Reaxff. It establishes the relations between changes of energy and changes of bonds. Thus, given the importance of bond order and the comparison of the bond order curves over atom distances in the two formularelated papers,15, 16 we find that the derivative of bond order over atomic distance and the range of bond order would have a significant impact on the tendency potential energy curves. All chemical bonds of hydrocarbons can fall into single bond including carbon-carbon single bond (C-C), carbon-hydrogen single bond (C-H)and hydrogen-hydrogen single bond (H-H); carboncarbon double bond (C=C) and carbon-carbon triple bond (C≡C). Each bond type has its own term in the bond order to bond length formula as in equation (1). The three terms of equation (1) represent single bond, double bond and triple bond respectively.  r  BOijn = BOijσ + BOijπ + BOijππ = exp  pbo1 ⋅  ijσ    ro  

pbo 2   r  + exp  p ⋅  ij   bo3  roπ  

  

pbo 4   r  + exp  p ⋅  ij   bo5  roππ  

  

pbo 6    

(1)

In our several failed attempts, no ideal results of better double bond breaking energy curves and triple bond breaking energy curves were observed. Thus, to get better fitness of double bond and triple bond Reaxff curves, some changes were made to the latter two terms of the formula of bond order. The basic idea behind these changes is to tune the energy gradient in order to get better potential energy curves because the gradient is calculated as follows.

dE dE dBO = • dq dBO dq


7

Energy & Fuels 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 8 of 22

where, E represents energy, q represents atom coordinates, and BO means the bond order of this bond. In non-reactive force field, the gradient is just derivative with respect to the atom coordinates. However, as mentioned above, bond order is the bridge between energy and atom coordinates. All the energy terms except coulomb energy and van der Waals energy lean on the atom coordinates implicitly. Thus, one simple way to improve the results is to pay attention to the bond order, including the bond order formula and the corresponding force field parameters. The modified bond order formula is shown in equation (2) below. And this modification makes the carbon-carbon double bond and carbon-carbon triple bond smoother.   rij BOijn = BOijσ + BOijπ + BOijππ = exp  pbo1 ⋅  σ  ro  

  

pbo 2  1 1 + + pbo 3 p  r r     bo 5 ij  1+ 1 +  ij     pbo 4   pbo 6 

(2)

Besides, this manuscript also represents an attempt of re-optimizing the parameters used in the Reaxff molecular dynamics simulations. GARFfield, a genetic algorithm-based reactive force field optimizer,40 was used to optimize the parameters. To re-optimize the parameters, a training set to train the parameters is needed. Here, the training set contains results of energies, geometries, charges from QM calculations. We begin the optimization procedure with the previously published parameters by Chenoweth et al.15. The strategy used for the parameters optimization also involves adjustment of bond length range with respect to bond order. Pauling41 proposed that a bond of any order n could be calculated through a simple relationship with its bond length. Based on that, Johnston et al.42 introduced an empirical correlation of dissociation energy versus bond length for all the possibilities of carbon-carbon bonds. It seems that the carbon-carbon bond dissociates at around 3.0 Å given the bond energy decrease to 0 kcal/mol at around 3.0 Å and this is consistent with the QM results. Thus, the upper limit of bond length with respect to bond order in Reaxff is increased to around 3.0 Å from previously used 2.5 Å


8

Page 9 of 22 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

Energy & Fuels

which turns out to be able to give better dissociation energy and to provide the dissociation energy curves more closely to the QM results for C-C single bond and H-H single bond. While every effort is made, there still exists a deficiency in this work, i.e. a better C-H bond dissociation potential energy curve is not achieved due to the complex coupling among the parameters. And this is part of our ongoing work and we anticipate that it could be improved through the continuing improvement of the method. The changes of C-C and H-H bond order along with bond length is depicted in Figure 2, the C-H bond order along with bond length is kept unchanged. Representative figures of comparisons of each kind of potential energy curves and charges are listed in Figure 3 and the rest are given in Figure S1.

(a)

(b)

Figure 2 Atom distance dependency of (a) C-C bond order, (b) H-H bond order. In this figure, old ff_single means the first part of equation (1) described using parameters from the reference paper.15 Likewise, new ff_single means the first part of equation (2) described using the newly optimized parameters. The rest are similar.


9

Energy & Fuels 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 10 of 22

(a)

(b)

(c)

(d)

Figure 3 Reaxff C-C bond dissociation energy curves of (a) isobutene, (b) ethylene, (c) acetylene and (d) hydrogen using both old ff and new ff. In addition to potential energy curves, apparent rate constants of three individual runs of iso-octane decomposition were also provided to demonstrate the validity of the newly optimized force field. All the three runs (see the Convergence Test section) were performed using the identical parameters but the random input structures. The apparent rate constants were calculated using data collected from 2200K-3000K. For a clear view, the apparent rate constants obtained from the three runs together with the experimental rate constants are shown in Figure 4.


10

Page 11 of 22 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

Energy & Fuels

Reasonable agreements can be found in Figure 4. Results of fitted Arrhenius Parameters using newly optimized parameters along with the experimental result are displayed in Table S1.

Figure 4 Apparent rate constants obtained from the three test runs and experimental work. Reference a: Davidson et al.7. 3.2 Convergence Test Molecular dynamics calculation is essentially a statistical tool that coupled statistical mechanics together with classical physics. Thus, to get a more unbiased result and meaningful statistical sampling, it is necessary to check the convergence of some important properties obtained from the simulation results. In this work, the multi-molecular Reaxff molecular dynamics calculations were performed to verify the convergence. 20 iso-octane molecules were placed in a cubic box. Three runs were performed individually to test the convergence. All the three runs shared the identical input parameters, including the random seeds for the initial speeds of the atoms, except the starting atoms configurations. The number of molecules over time and the number of species produced each step were drawn in Figure 5. From Figure 5(a) and Figure 5(b), consistent results could be found among the three runs. Thus, the three trajectories could give converged results for the following analysis and all the three runs were used in the analysis for multi-molecular simulations.


11

Energy & Fuels 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

(a)

Page 12 of 22

(b)

Figure 5 Selected properties for the three test runs. (a) number of molecules over time, (b) number of species produced each step. 3.3 Overview of initial C8H18 decomposition To get an overview of the pyrolysis of iso-octane, the averaged number of main products obtained from three simulations at 3000K over time is shown in Figure 6. As shown in Figure 6(a), the main species with large concentrations reached equilibrium after around 300ps and fluctuated within a limited range thereafter. To investigate how the main species change in the early time of decomposition, the same figure was plotted in Figure 6(b) between 0-300ps. According to Figure 6(b), the decomposition of iso-octane initiates nearly at 20 ps. It takes around another 20ps to consume all the iso-octane molecules. Since the initiation of iso-octane decomposition, the number of all the others species begin to increase. At approximately 30 ps, most of the iso-octane are about to decompose while the counts of CH3, C3H6 and C4H8 reach their maximum which implies these species might be the most abundantly initial products, that is to say, generations of those three species might be directly correlated to the consumption of the reactants. This is consistent with others.7 D.F. Davidson et al.7 stated that the rapid conversion to CH3 and other small radicals and stable intermediates is responsible for the initial fuel consumption. After around 30ps, the number of CH3 began to decrease. Meanwhile, the amounts of CH4 began to increase. At around 50ps, most of the CH3 depleted accompanied by the number


12

Page 13 of 22 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

Energy & Fuels

of CH4 reached to its maximum. Thereafter, CH4 began to deplete. From Figure 6(a), we could conclude that nearly no CH4 existed at 2ns which is consistent with the recent experimental results.13 After 30 ps, the amount of CH3, C3H6 and C4H8 begin to decrease with the increasing number of C2H2, C2H3 and C2H4. After around 300ps, the tendency of main products became stable and began to fluctuate with time within a limited range, as shown in Fig. 5(b). At 2ns, the most abundant species are H2, C2H4 and C2H3, C2H2 come the next, and other species only take up a little percentage. This reflects the nature of dissociating into small intermediates for the thermal decomposition reactions. According to Figure 6(c), the number of new types of species increases constantly from 0 ns to 2 ns which mean that new species are still being generated via reactions till 2ns. This means that new reactions constantly happen though most of the major products equilibrate at around 300 ps. To the best of our knowledge, this reflects the complex nature of combustion.

(a)

(b)

(c)

Figure 6 (a) Counts of main species from 0ns to 2ns, (b) counts of main species from 0 ns to 200ps. Changes of numbers of major compounds over time were obtained from Reaxff molecular dynamics during iso-octane decomposition at 3000K. Figures were drawn using averaged values from the three runs. (c) Counts of new species along with time. Although several experimental works have been devoted to gain understandings of characteristics during pyrolysis/oxidation processes, information like the percentage of the


13

Energy & Fuels 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 14 of 22

reactants decomposition pathway is still rare. In this subsection, we focus the C8H18 decomposition process in a statistical fashion through tracing the chemical process of each C8H18 molecule. To achieve this, a scheme to discriminate the molecular topology was proposed as illustrated in Scheme 1.

Scheme 1 Scheme to identify C8H18 molecules. The number side the atom is the atom id of the corresponding atom; the two lines of numbers are “identification” proposed of this chain molecule. Other detailed explains about this scheme could be found in the text. This molecule is divided into two parts: one is the backbone part; the other is the branch part. Two atom types are involved in this molecule. The first type is carbon atom and is specified as 1. The second type is hydrogen atom and is specified as 2. In this scheme, the first number, 26, beneath the C8H18 molecule of Scheme 1 stands for the total atoms number of this molecule. The next 8 and 18 mean that there are 8 carbon atoms and 18 hydrogen atoms in this molecule. 3 mean that there are three atoms on the backbone of this molecule. The next underlined three numbers are the atom id of the three backbone atoms and their absolute values are pointless here. The next ‘1 1 1’ are to indentify the atom type of the atom 5, atom 6 and atom 9, i.e. they are carbon atoms. The next ‘4 4 4’ are to specify that the overall coordination number of each of the three backbone atoms is four. The next three numbers ‘4 2 3’ mean that there are four carbon


14

Page 15 of 22 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

Energy & Fuels

atoms, two carbon atoms and three carbon atoms coordinate with atom 5, atom 6 and atom 9, respectively. The next three numbers ‘0 2 1’ mean that there are zero hydrogen atoms, two hydrogen atoms and one hydrogen atoms coordinate with atom 5, atom 6 and atom 9. The next five underlined numbers are the atom id of the branched atoms and again their absolute values are pointless. The ‘1 1 1 1 1’ are to indicate the atom type of these branched atoms and atom type of 1 means the five branched atoms are carbon atom. The next five numbers, i.e. ‘4 4 4 4 4’, indicate the overall coordination number of these five branched atoms. The next ‘1 1 1 1 1’ and ‘3 3 3 3 3’ indicate that each of these five branched atoms coordinates with one carbon atom and three hydrogen atom. The next five numbers ‘1 1 1 3 3’ represent the five branched atoms which are directed coordinated to the backbone of this molecule. The rest hydrogen atoms are omitted in this methodology and could be added based on the above information. With this methodology and a script based on it, dissociation pathway of each C8H18 is given and is shown in Figure 7.

Figure 7 Percentage of each C8H18 dissociation pathway. The atom labels listed along with the atoms are shown without labeling other equivalent atoms (e.g. the three -CH3 coordinating with atom 4 are equivalent) for simplicity. Figure 7 demonstrates the distribution of each dissociation point. The carbon-carbon bond fission clearly dominates the iso-octane decomposition and few C8H18 molecules dissociate


15

Energy & Fuels 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 16 of 22

following the C-H bond fission. This might be caused the larger C-H bond dissociation energy and less polar of this molecule. The percentage of dissociation of bond 1-2, bond 1-3, bond 3-4 and bond 4-5 increase which is in accordance with the decreasing tendency of the bond dissociation energies in the same sequence at the CASPT2-(2e,2o)/6-31+g(d,p)//CAS(2e,2o)/631+G(d,p) level of theory which again demonstrates the validity of the newly optimized parameters14. As seen in Figure 7, the primary decomposition product of C8H18 is CH3 which fits well the results in Figure 6(b) in which CH3 reaches its peak first. 3.3 Products at 2ns The main products at 2ns are listed in Figure 8. Only those having a large number percentage (larger than 0.01) are presented. As seen in Figure 8, hydrogen has the biggest percentage followed by C2H4 and C2H3. This figure implies that fuel decomposition at high temperature is a process causing successive hydrogen abstraction and increasing the number of small radicals and intermediates. This is consistent with a similar research43 in which the author found increasing temperature would contribute to the formation of hydrocarbons with 6 or more carbon atoms during thermal decomposition of methane at 0.1 g/cm3. The author also claimed that intermediates with 6 or more carbon atoms are only found at temperatures of 3250K which is also agree with our results. In Figure 9, counts of intermediates with the same carbon atom number changing with time are listed. From this figure, hydrogen and H radical have the largest number. The second is the C2 intermediates. The next is C3 intermediates. This kind of intermediates gradually decreases from their peak and then gets stable. The peak might be related to the decomposition of C8H18 as stated above. In a similar study of methane pyrolysis,44 the authors also found a maximum in the concentration of the C2-hydrocarbons during the pyrolysis.


16

Page 17 of 22 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

Energy & Fuels

Figure 8 Main Products at 2ns. Percentage was calculated and averaged from all the three runs.

Figure 9 The number of intermediates with the same carbon atom number as a function of time. 4 Summary Reaxff method is an important tool used in combustion field and many other fields. A modification to the bond order formula was proposed and the force field parameters were reoptimized. All these changes rely on the significant role of bond order in this method and are centered on energy gradients in Reaxff. Better fitness with the QM potential energy curves were found based on these modifications. The newly optimized parameters based reaxff were used to extend the knowledge of thermal decomposition of iso-octane. The apparent rate constants calculated using the new parameters is consistent to the experimental results. Under the


17

Energy & Fuels 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 22

simulated conditions, all the iso-octane decomposed around 20 ps later. While the number of species with larger concentrations reached stable within 300 ps, new types of chemical reactions are still happening until a longer time later. At 2ns, a maximum of C2-hydrocarbons during the decomposition was found. At last, the distribution of iso-octane decomposition pathway was given and the carbon-carbon bond cleavage was found to dominate the iso-octane thermal decomposition. ASSOCIATED CONTENT Supporting information This supporting information is available free of charge on the ACS Publications website. Contents of the supporting information: Force field parameters used in the simulation Figure S1 Additional potential energy curves and charge changes described by parameters of Ref. 15 and the newly optimized parameters. Table S1. Fitted Arrhenius Parameters of iso-octane Pyrolysis AUTHOR INFORMATION Corresponding Author E-mail: [email protected] ORCID Keli Han: 0000-0001-9239-1827 Notes The authors declare no competing financial interest. ACKNOWLEDGMENT


18

Page 19 of 22 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

Energy & Fuels

This work was supported by the National Natural Science Foundation of China (No. 21403221 and 91441106). REFERENCES 1.

Dooley, S.; Won, S. H.; Chaos, M.; Heyne, J.; Ju, Y.; Dryer, F. L.; Kumar, K.; Sung, C.-J.;

Wang, H.; Oehlschlaeger, M. A.; Santoro, R. J.; Litzinger, T. A. Combust. Flame 2010, 157, 2333-2339. 2.

Westbrook, C. K.; Pitz, W. J.; Herbinet, O.; Curran, H. J.; Silke, E. J. Combust. Flame

2009, 156, 181-199. 3.

Davis, S. G.; Law, C. K. Combust. Sci. Technol. 1998, 140, 427-449.

4.

Davidson, D. F.; Gauthier, B. M.; Hanson, R. K. Proc. Combust. Inst. 2005, 30, 1175-

1182. 5.

Davidson, D. F.; Oehlschlaeger, M. A.; Herbon, J. T.; Hanson, R. K. Proc. Combust. Inst.

2002, 29, 1295-1301. 6.

Oehlschlaeger, M. A.; Davidson, D. F.; Herbon, J. T.; Hanson, R. K. Int. J. Chem. Kinet.

2004, 36, 67-78. 7.

Davidson, D. F.; Oehlschlaeger, M. A.; Hanson, R. K. Proc. Combust. Inst. 2007, 31,

321-328. 8.

Slavinskaya, N. A.; Haidn, O. J. J. Propul. Power 2003, 19, 1200-1216.

9.

Curran, H. J.; Gaffuri, P.; Pitz, W. J.; Westbrook, C. K. Combust. Flame 2002, 129, 253-

280. 10.

Chaos, M.; Kazakov, A.; Zhao, Z.; Dryer, F. L. Int. J. Chem. Kinet. 2007, 39, 399-414.

11.

Eddings, E. G.; Yan, S.; Ciro, W.; Sarofim, A. F. Combust. Sci. Technol. 2005, 177, 715-

739.


19

Energy & Fuels 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 20 of 22

12.

Zhang, H. R.; Eddings, E. G.; Sarofim, A. F. Proc. Combust. Inst. 2007, 31, 401-409.

13.

Malewicki, T.; Comandini, A.; Brezinsky, K. Proc. Combust. Inst. 2013, 34, 353-360.

14.

Ning, H.; Gong, C.; Li, Z.; Li, X. J. Phys. Chem. A 2015, 119, 4093-4107.

15.

Chenoweth, K.; van Duin, A. C. T.; Goddard, W. A. J. Phys. Chem. A 2008, 112, 1040-

1053. 16.

van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard, W. A. J. Phys. Chem. A 2001, 105,

9396-9409. 17.

Senftle, T. P.; Hong, S.; Islam, M. M.; Kylasa, S. B.; Zheng, Y.; Shin, Y. K.; Junkermeier,

C.; Engel-Herbert, R.; Janik, M. J.; Aktulga, H. M.; Verstraelen, T.; Grama, A.; van Duin, A. C. T. Npj Computational Materials 2016, 2, 15011. 18.

Chenoweth, K.; van Duin, A. C. T.; Dasgupta, S.; Goddard Iii, W. A. J. Phys. Chem. A

2009, 113, 1740-1746. 19.

Ding, J.; Zhang, L.; Zhang, Y.; Han, K.-L. J. Phys. Chem. A 2013, 117, 3266-3278.

20.

Liu, L.; Bai, C.; Sun, H.; Goddard, W. A. J. Phys. Chem. A 2011, 115, 4941-4950.

21.

Wang, Q.-D.; Wang, J.-B.; Li, J.-Q.; Tan, N.-X.; Li, X.-Y. Combust. Flame 2011, 158,

217-226. 22.

Zheng, M.; Li, X.; Liu, J.; Guo, L. Energy & Fuels 2013, 27, 2942-2951.

23.

Li, G. Y.; Ding, J. X.; Zhang, H.; Hou, C. X.; Wang, F.; Li, Y. Y.; Liang, Y. H. Fuel 2015,

154, 243-251. 24.

Mo, Z. Mol. Sim. 2015, 41, 13-27.

25.

Zheng, M.; Li, X.; Liu, J.; Wang, Z.; Gong, X.; Guo, L.; Song, W. Energy & Fuels 2015,

28, 522–534. 26.

Zhang, T.; Li, X.; Qiao, X.; Zheng, M.; Guo, L.; Song, W.; Lin, W. Energy & Fuels 2016,


20

Page 21 of 22 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

Energy & Fuels

30, 3140-3150. 27.

Wang, J.-P.; Li, G.-Y.; Guo, R.; Li, A.-Q.; Liang, Y.-H. Energy & Fuels 2017, 31, 124-

132. 28.

Zhang, C.; Zhang, C.; Ma, Y.; Xue, X. Phys. Chem. Chem. Phys. 2015, 17, 11469-11480.

29.

Xue, X.; Meng, L.; Ma, Y.; Zhang, C. J. Phys. Chem. C 2017, 121, 7502-7513.

30.

Zheng, M.; Wang, Z.; Li, X.; Qiao, X.; Song, W.; Guo, L. Fuel 2016, 177, 130-141.

31.

Liu, Y.; Zybin, S. V.; Guo, J.; van Duin, A. C.; Goddard III, W. A. J. Phys. Chem. B 2012,

116, 14136-14145. 32.

Wen, Y.; Xue, X.; Zhou, X.; Guo, F.; Long, X.; Zhou, Y.; Li, H.; Zhang, C. J. Phys. Chem.

C 2013, 117, 24368-24374. 33.

Zhou, T.; Song, H.; Liu, Y.; Huang, F. Phys. Chem. Chem. Phys. 2014, 16, 13914-13931.

34.

Wen, Y.; Zhang, C.; Xue, X.; Long, X. Phys. Chem. Chem. Phys. 2015, 17, 12013-12022.

35.

Xue, X.; Wen, Y.; Long, X.; Li, J.; Zhang, C. J. Phys. Chem. C 2015, 119, 13735-13742.

36.

Zou, C.; Duin, A. C. T. V.; Dan, C. S. Top. Catal. 2012, 55, 391-401.

37.

Jensen, B. D.; Wise, K. E.; Odegard, G. M. J. Comput. Chem. 2015, 36, 1587-1596.

38.

Christopher Browne, N. O., Ale Strachan The Summer Undergraduate Research

Fellowship (SURF) Symposium, Purdue University, West Lafayette, Indiana, USA 2014. 39.

Steve, P.; Paul, C.; Aidan, T. LAMMPS (Sandia National Labs, Albuquerque, NM, 2012)

2009. 40.

Jaramillo-Botero, A.; Naserifar, S.; Goddard, W. A. J. Chem. Theory Comput. 2014, 10,

1426-1439. 41.

Pauling, L. J. Am. Chem. Soc. 1947, 69, 542-553.

42.

Johnston, H. S.; Parr, C. J. Am. Chem. Soc. 1963, 85, 2544-2551.


21

Energy & Fuels 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 22 of 22

43.

Lümmen, N. Phys. Chem. Chem. Phys. 2010, 12, 7883-7893.

44.

Arutyunov, V.; Vedeneev, V.; Moshkina, R.; Ushakov, V. Kinet. Catal. 1991, 32, 234-240.


22

Development of a Reactive Force Field for ... - ACS Publications

Recommend Documents