ReaxFF Molecular Dynamics Simulations of Oxidation of Toluene at

Sep 21, 2012 - The shielding functions have been used to avoid the excessive repulsion of atoms at close distances, and the Electronegativity Equaliza...
16 downloads 0 Views 1MB Size
Download PDF
Article pubs.acs.org/JPCA

ReaxFF Molecular Dynamics Simulations of Oxidation of Toluene at High Temperatures Xue-Min Cheng,† Quan-De Wang,‡ Juan-Qin Li,† Jing-Bo Wang,† and Xiang-Yuan Li*,† †

College of Chemical Engineering, Sichuan University, Chengdu 610064, People's Republic of China College of Chemistry, Sichuan University, Chengdu 610065, People's Republic of China

‡

S Supporting Information *

ABSTRACT: Aromatic hydrocarbon fuels, such as toluene, are important components in real jet fuels. In this work, reactive molecular dynamics (MD) simulations employing the ReaxFF reactive force field have been performed to study the hightemperature oxidation mechanisms of toluene at different temperatures and densities with equivalence ratios ranging from 0.5 to 2.0. From the ReaxFF MD simulations, we have found that the initiation consumption of toluene is mainly through three ways, (1) the hydrogen abstraction reactions by oxygen molecules or other small radicals to form the benzyl radical, (2) the cleavage of the C−H bond to form benzyl and hydrogen radicals, and (3) the cleavage of the C−C bond to form phenyl and methyl radicals. These basic reaction mechanisms are in good agreement with available chemical kinetic models. The temperatures and densities have composite effects on toluene oxidation; concerning the effect of the equivalence ratio, the oxidation reaction rate is found to decrease with the increasing of equivalence ratio. The analysis of the initiation reaction of toluene shows that the hydrogen abstraction reaction dominates the initial reaction stage at low equivalence ratio (0.5−1.0), while the contribution from the pyrolysis reaction increases significantly as the equivalence ratio increases to 2.0. The apparent activation energies, Ea, for combustion of toluene extracted from ReaxFF MD simulations are consistent with experimental results. intermediate and final product distributions obtained by experiment alone. Therefore, kinetic modeling and computational approaches can help us widen our sight and provide an efficient way to understand the true nature of the complex combustion reactions. Generally, a detailed reaction mechanism consists of numerous species and elementary reactions, such as the blended fuel model (toluene−butane) KBG composed of 97 species and 529 reactions,5 the detailed reaction mechanism for toluene oxidation developed by Dagaut composed of 120 species and 920 reactions,6 and the recent kinetic model employed by Qi et al. with 176 species and 804 reactions.7 However, it is difficult to apply these detailed reaction mechanisms to computational fluid dynamics (CFD) simulation. To obtain a reduced reaction mechanics for CFD simulation, an accurate, computationally feasible, atomistic method to describe the combustion of hydrocarbons is urgently needed. In theory, quantum mechanical (QM) calculations can provide accurate reaction rate constants and thermochemical parameters used for chemical kinetic modeling. However, it should be noted that QM calculations now can only deal with simple reactive systems for their computational expense limits

I. INTRODUCTION Aromatic hydrocarbon fuels are an important class of components in real jet fuel. A detailed investigation on combustion and pyrolysis reaction mechanisms of these species can provide important information for the development of combustion mechanisms of real fuels. Further more, aromatic hydrocarbon attracts particular attention for its particular molecular structure, high energy density, and antiknock rating.1 Additionally, it also can be formed during the oxidation of aliphatic hydrocarbon in the fuels, so a deeper understanding of their oxidation mechanisms is needed. As the simplest monosubstituted aromatic hydrocarbon, toluene is commonly used as a precursor fuel in gasoline reference fuels to study the combustion behaviors of aromatic components in practical experiments and theories, and it is a key alkylated aromatic intermediate leading to subsequent formation of polycyclic aromatic hydrocarbons (PAHs) and eventually to soot when incomplete hydrocarbon oxidation occurs.2 To reduce PAH and soot formation, a detailed understanding of reaction pathways of toluene under practical combustion conditions is necessary. Additionally, toluene has been used as an aromatic representative in surrogate mixtures for jet fuels.3,4 Because of the desired properties compared with other aromatics, toluene is one of the most experimentally investigated aromatic fuels. However, it is a great challenge to disentangle the various fundamental combustion steps from the © 2012 American Chemical Society

Received: April 26, 2012 Revised: September 20, 2012 Published: September 21, 2012 9811

dx.doi.org/10.1021/jp304040q | J. Phys. Chem. A 2012, 116, 9811−9818

The Journal of Physical Chemistry A

Article

employed to account for polarization effects. ReaxFF force field exhibits good performances on alkanes;9,16−18 however, it has been seldom utilized to the studies of aromatic hydrocarbon combustions. To establish whether ReaxFF is able to describe aromatic hydrocarbon combustions, a series of MD simulations have been carried out to study the oxidation mechanisms of toluene. In this work, a wide range of input conditions for the toluene/O2 mixtures have been studied. For the ReaxFF simulations, all structures are minimized at 5 K with a time step of 0.1 fs with a constant number of atoms (N) in a constant volume (V) and temperature (T), denoted as NVT-MD. Then these systems are equilibrated for 50 ps (ps) using NVT-MD simulation at 5 K with a MD time step of 0.1 fs. This time step was recommended for hydrogen oxidation at high temperature (2500 K) by Chenoweth et al.,11 and it gives reasonable descriptions for the oxidation reaction of hydrocarbons. C−O and H−O bond parameters are switched off during the equilibration simulations to prevent reactions occurrence. The equilibrated systems are then used in the NVT-MD simulations where the bond-order cutoff is 0.3 and the MD time step is 0.1 fs. During all the NVT-MD simulations, the system temperatures are controlled via a Berendsen thermostat23 with a 0.5 ps damping constant. The equivalence ratio (Φ) is determined by the ratio required for complete combustion of the toluene to carbon dioxide and water, and the system densities are set by the edge lengths of the cubic boxes. A summary of all the input conditions can be seen in Table 1. The system pressures

and practical applications. On the other hand, a molecular dynamics (MD) simulation is an important tool to investigate molecular properties that are difficult to measure or expensive to obtain by experimental methods directly. Traditional MD simulations can deal with large systems (thousands of atoms), but they generally fail to describe the processes of bond breaking and bond formation in the chemical reactions. Previously, Brenner et al. developed a force field for hydrocarbons that can describe chemical bonds broking, which is inspired by the introduction of bond order/bond distance relationship.8 In their force field, the connectivity of every atom in the system is determined at each iteration and is a function of bond order. However, the van der Waals and Coulombic interactions are not included in this force field, so dissociative PESs predicted by the Brenner force field can not be quantified. Consequently, these methods can not give a full description of breaking and forming bonds during complex chemical reactions. The most recent success attempt is that of van Duin and co-workers who developed a force field that includes Coulombic and van der Waals interactions explicitly at each iteration during the whole MD simulation.9 The parameters in ReaxFF are derived based on QM calculations and can reproduce the PES for specific systems. Therefore, the ReaxFF force field has good transferability and can be used in a wide range of systems. Recently, van Duin et al. have extended their force field training set to include transition-state and chemical reactivity data for hydrocarbon oxidation reactions.10,11 Successful applications using ReaxFF simulations to model the oxidation of hydrocarbons have been reported9−15 and help us to gain a deep understanding of the reaction mechanisms in the complex combustion phenomena. In the present work, we present systematic MD investigations on the high-temperature oxidation of toluene. The ReaxFF force field for C/H/O systems developed by van Duin et al.9 is used to simulate the oxidation of toluene. The hightemperature MD simulations are performed at equivalence ratios of 0.5, 1.0, and 2.0 under different temperatures and densities. The paper is organized as follows: first, the full descriptions of ReaxFF MD simulation details are presented; second, the detailed reaction mechanisms extracted from the computed trajectories of the simulation systems, together with the effects of temperatures, densities, and equivalence ratios on the reaction mechanisms of toluene oxidation are discussed; finally, the main conclusions about oxidation of toluene is presented. It is expected that this work will help us to gain a better understanding of the oxidation of toluene at high temperatures.

Table 1. Temperature, Density, and Equivalence Ratio of the Systems Studied equivalence ratio (Φ)

toluene/ oxygen

temp (K)

density (kg/ dm3)

ensemble

0.09 0.5 1 2

1/100 10/180 10/90 10/45

2500−3000 2500−3500 2500−3500 2500−3500

0.35 0.05−0.35 0.05−0.35 0.05−0.35

NVT NVT NVT NVT

corresponding to a series of chosen initial systems have been estimated by the van der Waals equation of state for the real gas,24 and they range from 8.65 MPa at the lowest temperature and density (T = 2500 K, ρ = 0.05 g/cm3, Φ = 2.0) up to 409.3 MPa at the highest temperature and density (T = 3500 K, ρ = 0.35 g/cm3, Φ = 0.5). The high temperatures and pressures used in MD simulations would not affect the reaction mechanisms but only influence the reaction rates, and these simulation conditions are not fully equivalent with that of reaction conditions. 2.2. Quantum Mechanical Calculations. To validate the reactive force field and its transferability, QM calculations are performed on the key intermediates and reactions observed during the ReaxFF MD simulations. Geometry optimization and the heats of formation for key reactions at 298 K are calculated from Gaussian-3 (G3),25 a high-level composite method, by using the Gaussian 03 program.26

2. COMPUTATIONAL DETAILES 2.1. ReaxFF MD Simulation Details. ReaxFF, a reactive force field based on the bond-order principle, considers the bond energy as a function of the bond length, valence angle, and torsion angle. It can give an accurate description for the oxidation and pyrolysis reactions of hydrocarbon as well as for other systems in recent years.13,16−21 The bond orders are updated at each simulation step, leading to the variation of whole system. Since the connectivity-dependent interactions such as valence angles and torsion angles are all bonddependent, their energy contribution changes when bonds break. The van der Waals and Coulomb interactions are also taken into account in ReaxFF. The shielding functions have been used to avoid the excessive repulsion of atoms at close distances, and the Electronegativity Equalization Method22 is

3. RESULTS AND DISCUSSION 3.1. Initial Reaction Mechanisms for Oxidation of Toluene. To investigate the initial reaction mechanisms of toluene oxidation, a series of NVT-MD simulations are carried out on unimolecular combustion models (containing 100 oxygen molecules and 1 toluene molecule) and multimolecular 9812

dx.doi.org/10.1021/jp304040q | J. Phys. Chem. A 2012, 116, 9811−9818

The Journal of Physical Chemistry A

Article

combustion models (containing 90 oxygen molecules and 10 toluene molecules). These simulations are performed at 0.35 g/ cm3 with temperatures ranging from 2500 to 3000 K in 100 K intervals. These independent simulations show benzyl is a dominant intermediate at different temperatures. This is reasonable and can be convinced by the facts that the methyl C−H bond is the weakest (89.2 kcal/mol)27 and the cleavage of this bond has the lowest enthalpy among all the reactions to form C7H7 radicals.28 Detailed analysis of toluene oxidation at each step observed in ReaxFF simulations at 2700 K have been listed in Scheme 1, where the dash cycles indicate the final

carbon dioxide. The extended simulation time of 1 ns (ns) for high temperature simulation (3000 K) does not result in further oxidation of carbon monoxide to carbon dioxide. This possibility is attributed to the longer simulation time demands for the further oxidation of carbon monoxide.29 In comparison with the unimolecular oxidation of toluene, we see that the initiation reactions of toluene oxidation occurring in the multimolecular oxidation reactions are very different. In addition to the H-abstraction reactions by oxygen molecules, many other reactions in the ReaxFF MD simulations have been observed and listed in Table 2. These reaction

Scheme 1. Toluene Partial Oxidation Pathway Observed in the ReaxFF NVT-MD Simulation at T = 2700 K

Table 2. Comparison of ReaxFF, G3, and Literature Values of Reaction Enthalpies (kcal/mol) for Reactions Related for Toluene Oxidation and Pyrolysisa reactions (A1) C6H5CH3+O2→ C6H5CH2+HO2 (A2) C6H5CH3+HO2→ C6H5CH2+H2O2 (A3) C6H5CH3+O→ C6H5CH2+OH (A4) C6H5CH3+OH→ C6H5CH2+H2O (A5) C6H5CH3+OH→C6H5OH+ CH3 (A6) C6H5CH3+H→ C6H5CH2+H2 (A7) C6H5CH3+CH2→ C6H5CH2+CH3 (B1) C6H5CH3→C6H5CH2+H (B2) C6H5CH3→C6H5+CH3

ΔHro0K (ReaxFF)

ΔHro298K (G3)

ΔHro (ref)b

20.57

11.06

41.12

−7.27

3.08

1.56

−28.32

−13.00

−13.04

−26.20

−27.47

−29.05

−6.32

−8.79

−9.00

−16.84

−14.15

−14.47

−27.94

−20.00

−20.69

92.73 95.95

90.51 104.44

89.73 103.72

a Molecular structures of the molecules and radicals (C6H5CH3, C6H5CH2, C6H5OH, and C6H5) are given as Supporting Information. b Literature values are taken from ref 35.

mechanisms can be divided into two classes: the H-abstraction reactions and the pyrolysis reactions. The H-abstraction reactions by oxygen molecules or other radicals are still the main initial reactions at lower equivalence ratio (Φ = 0.5). By contrast, pyrolysis reactions have major contributions to toluene consumption at high equivalence ratio conditions (Φ = 2.0). The ring-opening reaction of toluene have been observed in the simulations that can hardly happen in practical conditions because of high CC doublet bond energy. However, the reaction occurs in our MD simulations, which may be induced by the high temperatures used in MD simulations. These initial steps are consistent with the kinetic models reported by previous reports.6,30−34 This suggests that ReaxFF simulations can provide rational reaction mechanisms for the combustion of aromatic substances. 3.2. Force Field Validation. The force field used in the ReaxFF MD simulation has been reported by Chenoweth et al. for aromatic hydrocarbon;11 now we will focus here on the goodness of the force field specific to the toluene system. Table 3 compared the standard heats of formation of ReaxFF and literature data30,35,36 for the key intermediates of toluene oxidation. The geometries used in the ReaxFF MD simulation are obtained by the optimized geometry of G3. We observed that the heats of formation of ReaxFF match well with the values reported in literatures with a maximum deviation of ∼8.7 kcal for the alkoxy benzyl radical. The initiation reactions of toluene oxidation occurring in the multimolecular oxidation

products observed in the whole simulation. From this observation, the hydrogen abstraction by oxygen is the first step of toluene oxidation in most cases, while only one step was initiated by pyrolysis of toluene to form hydrogen and benzyl radicals (IM1). This is consistent with experimental results, namely, side chains on the ring are facilitated to react first.27 Next benzyl and hydroperoxide radicals combine with each other to form C7H8O2 molecule (IM2). The C7H8O2 molecule (IM2) exists for only a short time at lower temperature, and it can not be observed at higher temperatures (≥2800 K), implying that the weak O−O bond hardly exists at high temperatures. The high-energy O−O bond is cleaved to form hydroxyl and benzyl alcohol radical (IM3) species. The benzyl alcohol radical (IM3) undergoes hemolytic C−C bond cleavage to form phenyl radical (IM4) and formaldehyde. The oxygen molecule then reacts with the phenyl radical to form C6H5O2 radical (IM5). Next, one of oxygen radicals migrates to the adjacent carbon and thus leads to C6H5O2 radical (IM6). Within a small simulation time, the oxygen radical would react with the adjacent carbon lead to the C−C bond opens, and then forms carbon monoxide and C5H5O radical (IM7). The ring radical, C5H5O (IM7) can be observed in most simulations, and it is formed via dissociation of CO from C6H5O2 radical (IM6). At last the ring opens and is further oxidized to the final products. The final products observed from the NVT-MD simulations include water, carbon monoxide, and 9813

dx.doi.org/10.1021/jp304040q | J. Phys. Chem. A 2012, 116, 9811−9818

The Journal of Physical Chemistry A

Article

indicating that this force field should be suitable for describing toluene system. 3.3. Temperature Effect on the Oxidation of Toluene. To analyze the temperature effect on the oxidation of toluene, a series of NVT-MD simulations are performed at temperatures from 2500 to 3500 K with an interval of 200 K for a total simulation time of 1 ns (N = 240 (Φ = 2.0)). The time evolution of the number of total species, toluene molecules, and the main products coupled with the potential energies at density of 0.35 g/cm3 are shown in Figure 1. Figure 1a describes the potential energy changes against time. It is found that the potential energy decreases rapidly with the increasing of temperature. The rate of potential energy decrement depends on temperature, that is, the rate reduces more quickly with increasing of temperature, indicating that temperature accelerate toluene oxidation. The reason for this phenomenon can be explained by the heat release of system, which depends on the temperature. For the higher temperature (T = 3500 K), the potential energy increases within the first 15 ps but decreases as the simulation moves on. However, the potential

Table 3. Comparison of Standard Enthalpy of Formation at 0 K of Key Species of Hydrogen Combustion in ReaxFF and Literature Values name toluene (C6H5CH3) benzyl (C6H5CH2) benzyl hydroperoxide (C6H5CH2OOH) alkoxy benzyl (C6H5CH2O) phenyl (C6H5)

ΔfHmΘ0K (ReaxFF)

ΔfHmΘ298K (ref)

15.14 43.02 −12.36

17.5a 49.7b −8.2c

18.56 74.56

27.3c 80.78c

a Literature value is taken from ref 35. bLiterature value is taken from ref 36. cLiterature data are taken from ref 30.

reactions have been employed in this section. Table 2 shows the comparison of the heats of reaction between ReaxFF and as obtained from G3 and Baulch et al.35 A notable discrepancy of reaction (A1) in Table 2 shows a larger disagreement with QM and literature value. However, most of reactions predicted by ReaxFF are agreement with QM calculations and literature data,

Figure 1. Time evolution of the potential energies, total number of species, toluene, and main products at the temperatures 2500−3500 K in the interval of 200 K (NVT-MD simulation, ρ= 0.35 g/cm3, Φ = 2.0). 9814

dx.doi.org/10.1021/jp304040q | J. Phys. Chem. A 2012, 116, 9811−9818

The Journal of Physical Chemistry A

Article

Figure 2. Time evolution of the potential energies, total number of species, toluene molecules, and main products at density ρ = 0.05, 0.15, and 0.35 g/cm3 (NVT-MD simulation, T = 3500 K, Φ = 2.0).

carbon monoxide at equivalence ratios of 0.5 and 1.0, no tendency is observed at equivalence ratio of 2.0, probably because of the low concentration of oxygen molecules comparing with the simulations at the equivalence ratio of 0.5 and 1.0. Similar global trends are also obtained for simulation systems with density of 0.15 g/cm3 and 0.05 g/cm3, and these results are not shown. 3.4. Pressure/Density Effect on the Oxidation of Toluene. The pressure effect on the combustion of toluene is investigated by analyzing the NVT-MD simulations at the system densities ρ = 0.05, 0.15, and 0.35 g/cm3 (N = 240 (Φ = 2.0), T = 3500 K). Figure 2 shows the time evolutions of the potential energy, total number of species, toluene molecules, and representative products observed during the ReaxFF MD simulations. From the changes of the potential energies, we can see that the system at 0.35 g/cm3 is the most exothermic. The tipping point tc in these simulations are all obvious at 0.05, 0.15, and 0.35 g/cm3, and the reactions at 0.15 and 0.35 g/cm3 have reached thermodynamic equilibrium within 1 ns. The tendencies of the potential energies and the toluene consumption indicate that increasing the density lead to an earlier initiation of reaction and fast consumption of fuel

energy for the system at 2500 K decreases dramatically until 440 ps. Together with the total number species and the number of toluene molecules, we can be certain that the whole system begins with an endothermic reaction. We define the tipping point tc as the time of reaching the maximum of the potential energy. It can be seen that tipping point tc increases with the temperature decreasing, indicating that endothermic reactions need much more time along with decreasing of temperature. The total number of fragments indicates that the reaction proceeds faster as the simulation temperature rises, and the conversions of toluene molecule also increase rapidly. From the simulation results, the primary initial reaction is the pyrolysis of toluene and hydrogen-atom abstraction by oxygen. This is consistent with the fact that the initial reaction is endothermic. The initial reaction products and byproducts then react with oxygen to release tremendous amounts of energy, and finally lead to the dramatic decrease of potential energy. The final products observed from the MD simulations include CO2, CO, and H2O. From the map of product distributions, H2O and CO begin to increase as the temperature goes up and it reaches a maximum number at 3500 K. Although we have observed that CO2 has the same tendency with the distribution of water and 9815

dx.doi.org/10.1021/jp304040q | J. Phys. Chem. A 2012, 116, 9811−9818

The Journal of Physical Chemistry A

Article

Figure 3. Temperature evolution of the initiation time as well as the half-time period extracted from ReaxFF MD simulations (NVT-MD simulation).

molecules. The density effect also can be reflected by the final product distributions. The products CO and H2O increase with the system densities go up and they reach a maximum number at 0.35 g/cm3. Figure 2f shows that CO2 is first appeared at 0.35 g/cm3, and its maximum concentration is lower than that at 0.15 g/cm3. The product distributions of CO2 at various densities are always lower than that of CO distributions which is consistent with the flame experiments reported by Qi et al.7 These are primarily due to the incomplete combustion. 3.5. Equivalence Ratio Effect on Kinetic Properties. The equivalence ratio effects on toluene/O2 systems are analyzed by performing a series of NVT-MD simulations covering a range of equivalence ratios from 0.5 to 2.0 (N = 510 (Φ = 0.5), N = 330 (Φ = 1.0), and N = 240 (Φ = 2.0)). The NVT-MD simulations of toluene/O2 are run at two different densities (0.15 and 0.35 g/cm3) for temperatures ranging from 2500 to 3500 K (in 200 K interval). The time required for initiation as well as half-consumption of fuel molecules are determined for each simulation, which we call initiation time and half-time period in the paper. The initiation time and halftime period determined for NVT-MD simulations are shown in Figure 3. It can be seen that longer initiation time and half-time period are required at lower temperature, lower density, and higher equivalence ratio. This suggests that toluene is less reactive at low temperature, low density, and high equivalence ratio. The temperature and density effects have discussed above, so we will further focus on the effect of equivalence ratio on toluene consumption. As shown in Figure 3a, it appears that the initiation time is positively correlated with equivalence ratio at the same temperature and density, and a similar trend also exhibits for the half-time period. The presence of some irregular date is likely due to the reactive event associated with the initiation determined for each simulation. The global trend of toluene consumption indicates that the ignition delay time of toluene would become longer as the equivalence ratio increases. This is consistent with the experimental results,30,31 and the kinetic modeling studies also reveals this phenomenon.6 Our ReaxFF MD simulations capture the unusual combustion properties of toluene, and this confirms the validation of our present simulation results. The hydrogen abstraction reaction is the dominant initiation reactions of toluene oxidation under fuel-lean conditions, while the pyrolysis reactions of (B1) and (B2) contribute more along with the increase of the equivalence ratio. The reflected shock-tube technique has investigated the dynamic behavior of reactions (B1), (B2), and

(A1) and obtained the activation energy of these reactions is 79.8,32 99.5,32 and 17.333 kcal/mol, respectively. This explains why low equivalence ratio can accelerate the consumption of toluene. 3.6. Overall Activation Energy of the Toluene Oxidation Reactions. The overall oxidation of toluene can be represented as follows C7H8 + 9O2 → 7CO2 + 4H 2O

NVT-MD simulations of toluene/O2 systems were run at two different densities (0.15 and 0.35 g/cm3) for temperature ranging from 2500 to 3500 K (in 200 K interval) and with equivalence ratio ranging from 0.5 to 2 (N = 240, 330, and 510). The overall rate constant k of reaction (R1) at each temperature can be simplified by using eq 1, i.e. k=−

d[C7H8] 1 dt [C7H8][O2 ]

(1)

where [C7H8] and [O2] are the average concentrations of toluene and oxygen, respectively. In the present work, these values are obtained from five simulations at each simulation temperature. Each five simulations were run for 0.8 ns with different initial configurations. The overall activation energies Ea are then calculated by using the Arrhenius rate equation (eq 2), i.e. kov = A exp( −Ea /RT )

(2)

We obtain the overall activity energy values ranging from 47.62 to 86.76 kcal/mol as shown in Table 4. Vasudevan et al. previously report an overall activation energy Ea (55.2 kcal/ Table 4. Fitted Arrhenius parameters ReaxFF

density(g/cm3)

equivalence ratio

Ea (kcal/mol)

R2

0.15

0.5 1.0 2.0 0.5 1.0 2.0

86.76 47.62 64.63 70.69 58.70 71.16 61.1a 55.2b

0.9494 0.9470 0.9110 0.9176 0.9198 0.9279

0.35

exp

a

Experimental result is taken from ref 31. bExperimental result is taken from ref 34. 9816

dx.doi.org/10.1021/jp304040q | J. Phys. Chem. A 2012, 116, 9811−9818

The Journal of Physical Chemistry A

Article

mol) by using the shock tube technique,34 while Mittal and Sung give Ea (61.1 kcal/mol) by experiment conducting in the rapid compression machine.31 There are no existing experimental results at our input conditions; however, the apparent activation energies extracted in this work are in accord with these reported experimental results.31,34

(7) Li, Y. Y.; Cai, J. H.; Zhang, L. D.; Yuan, T.; Zhang, K. W.; Qi, F. Proc. Combust. Inst. 2011, 33, 593−600. (8) Brenner, D. W. Phys. Rev. B 1990, 42, 9458−9471. (9) van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard, W. A. J. Phys. Chem. A 2001, 105, 9396−9409. (10) Agrawalla, S.; van Duin, A. C. T. J. Phys. Chem. A 2011, 115, 960−972. (11) Chenoweth, K.; van Duin, A. C. T.; Goddard, W. A. J. Phys. Chem. A 2008, 112, 1040−1053. (12) Chenoweth, K.; van Duin, A. C. T.; Dasgupta, S.; Goddard, W. A. J. Phys. Chem. A 2009, 113, 1740−1746. (13) Wang, Q. D.; Hua, X. X.; Cheng, X. M.; Li, J. Q.; Li, X. Y. J. Phys. Chem. A 2012, 116, 3794−3801. (14) Li, X. Y.; Wang, Q. D.; Wang, J. B.; Li, J. Q.; Tan, N. X. Combust. Flame 2011, 158, 217−226. (15) Liu, L. C.; Bai, C.; Sun, H.; Goddard, W. A. J. Phys. Chem. A 2011, 115, 4941−4950. (16) Jeon, B.; Sankaranarayanan, S. K. R. S.; van Duin, A. C. T.; Ramanathan, S. J. Chem. Phys. 2011, 134, 234706. (17) Ojwang, J. G. O.; van Santen, R. A.; Kramer, G. J.; van Duin, A. C. T.; Goddard, W. A. J. Chem. Phys. 2009, 131, 044501. (18) Valentini, P.; Schwartzentruber, T. E.; Cozmuta, I. J. Chem. Phys. 2010, 133, 084703. (19) van Duin, A. C. T.; Merinov, B. V.; Han, S. S.; Dorso, C. O.; Goddard, W. A. J. Phys. Chem. A 2008, 112, 11414−11422. (20) Zhang, L. Z.; Zybin, S. V.; van Duin, A. C. T.; Dasgupta, S.; Kober, E. M.; Goddard, W. A. J. Phys. Chem. A 2009, 113, 10619− 10640. (21) Weismiller, M. R.; van Duin, A. C. T.; Lee, J.; Yetter, R. A. J. Phys. Chem. A 2010, 114, 5485−5492. (22) Mortier, W. J.; Van Genechten, K.; Gasteiger, J. J. Am. Chem. Soc. 1985, 107, 829−835. (23) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684−3690. (24) Mortimer, R. G. Physical Chemistry; Academic Press: London, 2008; pp 2−15. (25) Henry, D. J.; Radom, L. In Quantum-Mechanical Prediction of Thermochemical Data; Cioslowski, J., Ed.; Kluwer Academic Publishers: Norwell, MA, 2001; Vol. 22 of Understanding Chemical ReactiVity, pp 161−197. (26) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, revision B.05; Gaussian, Inc.: Wallingford, CT, 2004. (27) Emdee, J. L.; Brezinsky, K.; Glassman, I. Symp. Combust. 1991, 23, 77−84. (28) Luo, Y. R. Comprehensive Handbook of Chemical Bond Energies; CRC Press: Boca Raton, FL, 2007. (29) Gardiner, W. C. Gas-Phase Combustion Chemistry; SpringerVerlag Press: New York, 2000. (30) Bounaceur, R.; Da Costa, I.; Fournet, R.; Billaud, F.; BattinLeclerc, F. Int. J. Chem. Kinet. 2005, 37, 25−49. (31) Mittal, G.; Sung, C. J. Combust. Flame 2007, 150, 355−368. (32) Sivaramakrishnan, R.; Michael, J. V. Proc. Combust. Inst. 2011, 33, 225−232.

4. CONCLUSIONS MD simulations employing the ReaxFF reactive force field have been performed to investigate the high-temperature oxidation of toluene at different temperatures, densities, and equivalence ratios. For the unimolecular oxidation process, the most important reaction at the initial stage is hydrogen abstraction reaction by the oxygen molecule. The pyrolysis of toluene is very rare and only occurs once in our all unimolecular oxidation simulations. However, in the multimolecular model, the pyrolysis of toluene becomes more competitive when the equivalence ratio increases. The effects of temperature, density, and equivalence ratio on high-temperature oxidation of toluene are further investigated. We have found that increasing the density and simulation temperatures of the systems can greatly accelerate the oxidation rate of toluene. Unlike the effect of temperature and pressure, our ReaxFF MD simulations demonstrate that the reaction rate of toluene consumption increases along with the decrease of equivalence ratio, which coincides with experimental conclusions and kinetic modeling results. Finally, the apparent kinetics is used to describe the oxidation of toluene, and the fitted apparent activation energies agree well with the experimental results.

■

ASSOCIATED CONTENT

S Supporting Information *

The molecular structures of the molecules and radicals in Table 3, the detailed kinetic analysis of the ReaxFF MD simulation results, and the ReaxFF force field parameters used in the work. This information is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank Dr. Adri van Duin at Pennsylvania State University for his help on the use of the ReaxFF code. This work is supported by the National Natural Science Foundation of China (Grant Nos. 91016002 and 21103117).

■

REFERENCES

(1) Lovell, W. G.; Campbell, J. M.; Signaigo, F. K.; Boyd, T. A. Ind. Eng. Chem. 1934, 26, 475−479. (2) Agafonov, G. L.; Smirnov, V. N.; Vlasov, P. A. Proc. Combust. Inst. 2011, 33, 625−632. (3) Bieleveld, T.; Frassoldati, A.; Cuoci, A.; Faravelli, T.; Ranzi, E.; Niemann, U.; Seshadri, K. Proc. Combust. Inst. 2009, 32, 493−500. (4) Cancino, L. R.; Fikri, M.; Oliveira, A. A. M.; Schulz, C. Proc. Combust. Inst. 2009, 32, 501−508. (5) Klotz, S. D.; Brezinsky, K.; Glassman, I. Symp. Combust. 1998, 27, 337−344. (6) Dagaut, P.; Pengloan, G.; Ristori, A. Phys. Chem. Chem. Phys. 2002, 4, 1846−1854. 9817

dx.doi.org/10.1021/jp304040q | J. Phys. Chem. A 2012, 116, 9811−9818

The Journal of Physical Chemistry A

Article

(33) Oehlschlaeger, M. A.; Davidson, D. F.; Hanson, R. K. Combust. Flame 2006, 147, 195−208. (34) Vasudevan, V.; Davidson, D. F.; Hanson, R. K. Proc. Combust. Inst. 2005, 30, 1155−1163. (35) Baulch, D. L.; Bowman, C. T.; Cobos, C. J.; Cox, R. A.; Just, T.; Kerr, J. A.; Pilling, M. J.; Stocker, D.; Troe, J.; Tsang, W.; Walker, R. W.; Warnatz, J. J. Phys. Chem. Ref. Data 2005, 34, 757−1397. (36) Douglas, J. D.; Robert, T. M., Jr.; Warren, J. H. J. Am. Chem. Soc. 1980, 102, 3334−3338.

9818

dx.doi.org/10.1021/jp304040q | J. Phys. Chem. A 2012, 116, 9811−9818