Article pubs.acs.org/JPCA
Ab Initio Reaction Kinetics of Hydrogen Abstraction from Methyl Formate by Hydrogen, Methyl, Oxygen, Hydroxyl, and Hydroperoxy Radicals Ting Tan,† Michele Pavone,‡,⊥ David B. Krisiloff,† and Emily A. Carter*,‡,§,∥ †
Department of Chemistry, ‡Department of Mechanical and Aerospace Engineering, §Program in Applied and Computational Mathematics, and ∥Andlinger Center for Energy and the Environment, Princeton University, Princeton, New Jersey 08544-5263, Unites States S Supporting Information *
ABSTRACT: Combustion of renewable biofuels, including energydense biodiesel, is expected to contribute significantly toward meeting future energy demands in the transportation sector. Elucidating detailed reaction mechanisms will be crucial to understanding biodiesel combustion, and hydrogen abstraction reactions are expected to dominate biodiesel combustion during ignition. In this work, we investigate hydrogen abstraction by the radicals H·, CH3·, O·, HO2·, and OH· from methyl formate, the simplest surrogate for complex biodiesels. We evaluate the H abstraction barrier heights and reaction enthalpies, using multireference correlated wave function methods including size-extensivity corrections and extrapolation to the complete basis set limit. The barrier heights predicted for abstraction by H·, CH3·, and O· are in excellent agreement with derived experimental values, with errors ≤1 kcal/mol. We also predict the reaction energetics for forming reactant complexes, transition states, and product complexes for reactions involving HO2· and OH·. High-pressure-limit rate constants are computed using transition state theory within the separable-hindered-rotor approximation for torsions and the harmonic oscillator approximation for other vibrational modes. The predicted rate constants differ significantly from those appearing in the latest combustion kinetics models of these reactions.
1. INTRODUCTION Fossil fuels are limited in supply and produce greenhouse gases upon combustion. The growing demand for energy motivates great interest in alternative liquid fuels that can be used efficiently, especially for transportation. Key alternatives include biofuels, such as biodiesel, which is an energy-rich fuel comprised of a blend of monoalkyl esters of long chain fatty acids1 that can be derived from renewable sources such as vegetable oils, animal fats, or microalgae.2 Besides reducing our dependence on petroleum, biodiesel combustion also produces less net CO2, particulate matter, and CO.3 Given these advantages, the U.S. domestic biodiesel consumption increased from 10 to 220 million gallons from 2001 to 2010, and its use is expected to accelerate in the future.4 The rapid expansion of the biodiesel industry, as well as the need for greater efficiency and reduced pollution, are motivating intense research efforts to achieve a comprehensive understanding of biodiesel combustion chemistry. While macroscopic qualities of biodiesel fuels, such as emission performance, are known,5−7 their combustion mechanisms remain highly uncertain. Biodiesel’s large molecular size results in a multitude of combustion pathways that hinder experimental and theoretical studies. Thus, combustion investigations have mostly addressed biodiesel surrogates © 2012 American Chemical Society
comprised of smaller methyl esters. For example, Fisher et al.8 developed a detailed combustion mechanism for methyl formate (MF) and methyl butanoate, as model compounds for biodiesel fuel. Subsequent studies9,10 have improved these mechanisms. Methyl decanoate, a larger methyl ester closer in size to biodiesel fuel, has also been studied.11,12 Investigation of these small-to-medium size methyl ester surrogates has improved our understanding of the actual biodiesel combustion chemistry. MF is the smallest methyl ester surrogate of biodiesel. Its simple structure makes it possible to isolate the combustion behavior of the ester functional group. MF also plays a significant role in the atmospheric oxidation of dimethyl ether and dimethoxy methane.13 Recently, Dooley et al.14 constructed and validated a mechanism for MF combustion with a rich set of experimental data, later augmented with low-pressure flame experiments and kinetics modeling.15 However, the kinetics parameters in these mechanisms are often just approximated from analogous simpler reactions. Received: May 17, 2012 Revised: July 20, 2012 Published: July 25, 2012 8431
dx.doi.org/10.1021/jp304811z | J. Phys. Chem. A 2012, 116, 8431−8443
The Journal of Physical Chemistry A
Article
verification of experimental data. G2 also uses electronic energies from QCISD(T),34 which does not properly describe chemical processes with significant multiconfigurational character, i.e., when two or more electron configurations contribute significantly to the electronic ground-state wave function. Bond dissociation, transition states (TSs), and diradicals all exhibit multiconfigurational character, all three of which figure prominently in combustion chemistry. To obtain barrier heights, we use a recently developed scheme that accounts for multiconfigurational behavior and achieves chemical accuracy without resorting to empirical corrections: the complete basis set extrapolated multireference singles and doubles configuration interaction (CBS-MRSDCI)28 approach. Meanwhile, lower levels of theory provide sufficiently accurate Arrhenius pre-exponential factors, probably due to error cancellation in the partition functions.35 The partition functions used to compute the Arrhenius pre-exponential factors here are obtained at the hybrid density functional theory (DFT) level. In the present work, we report rate constants of MF hydrogen abstraction reactions by five radical species (H, O, OH, CH3, and HO2), obtained by TST within the separablehindered-rotor approximation. The reaction enthalpies and barrier heights of these reactions are predicted by the CBSMRSDCI composite scheme,28 including size-extensivity corrections. We compare against barrier heights estimated from data in a recently reported MF combustion mechanism14 as well as other theoretical predictions. These results also provide further validation of the size-extensivity-corrected CBSMRSDCI model chemistry approach to studying elementary chemical reactions related to biodiesel combustion chemistry.
Biodiesel combustion mechanisms contain a daunting number of possible reactions. Reaction pathway16 and sensitivity17 analyses can help to distinguish the most relevant pathways and the most important elementary reactions. At high temperatures, these analyses8,11,14 reveal that the dominant methyl ester oxidation pathway starts with hydrogen abstraction followed by β-scission. In this way, a large methyl ester breaks into small molecules and radicals.11,18 For example, kinetics simulations indicate that H abstraction reactions by H, OH, and O radicals are responsible for approximately 55.5% of the consumption of MF in a fuel-to-air ratio Φ = 1.0 flame in the region of 0−7 mm from the burner surface. 15 Consequently, H abstraction reactions from the methyl ester largely determine the ignition characteristics of the fuel as well as the distribution of products. Very accurate thermodynamic and kinetic data of the H-abstraction reaction are required for a definitive assessment of the MF combustion mechanism, but the elusiveness of the intermediate species has hindered thorough experimental investigations. Because of the significant role MF plays in atmospheric chemistry, a large number of studies have investigated H abstraction from MF by hydroxyl radical, whereas other abstracting species have attracted less attention. Le Calvé et al.19 studied H abstraction by OH using pulsed laser photolysis laser-induced fluorescence, reporting the rate constant over the temperature range 233−372 K. Willington et al.20 used flash photolysis resonance fluorescence, and Szilágyi et al.21 employed isothermal fast flow coupled with resonance fluorescence detection to measure the reaction rate of H abstraction by OH at room temperature. Donovan et al. measured the rate constant of H abstraction from MF by methyl radical using mass spectroscopy, deriving an activation energy of 10.3 ± 0.2 kcal mol−1 in the temperature range 385 to 499 K.22 Mori measured reaction rates for H abstraction from MF by ground state 3P oxygen atom at room temperature.23 To the best of our knowledge, no direct experimental observations have been reported for the specific reactions of MF with H and HO2. Considering the molecular structure of MF, there are two different positions where H abstraction can occur. Few experiments have measured the two reaction rates separately.24 Consequently, quantum chemistry calculations provide the most effective way to investigate these reaction pathways in detail. A rigorous theoretical method is needed to study combustion mechanisms. Transition state theory (TST)25 has been widely used to predict kinetics parameters, both pre-exponential Arrhenius factors and barrier heights.26,27 A high level of theory is required to compute reliable barrier heights; an error of only a few kcal mol−1 in energy barriers can lead to qualitatively incorrect predictions, especially at low temperatures. Because current model chemistry approaches are often able to predict thermochemical and kinetics parameters of molecular species within chemical accuracy (i.e., within an error of ∼1 kcal mol−1 with respect to experiments),28−31 they can be used to study the energetics of these combustion reactions. Only a few theoretical studies of MF H abstraction have been reported. Good and Francisco used the G2 methodology to predict reaction enthalpies and barrier heights for MF reacting with an OH radical32 and rate constants of H abstraction reactions of MF with H, CH3, and OH radicals in a subsequent study with QCISD(T).33 The G2 methodology relies on empirical correction terms, which limits its independent
2. METHODS AND COMPUTATIONAL DETAILS We predict both reaction energetics (barrier heights and reaction enthalpies) and rate constants for MF H abstraction in this work. The activation energy (barrier height) is defined as the energy of the TS minus the energy of reactants, including zero point energies (ZPEs). Reaction enthalpies are obtained as the energy of the products minus the energy of reactants, including ZPEs and enthalpic thermal corrections (TCs) calculated from molecular partition functions (using ideal gas, rigid rotor, and harmonic oscillator (IGRRHO) approximations). The detailed features of our approach were described recently.28 Briefly, we use hybrid DFT for computing the optimized structures and corresponding vibrational frequencies of reactants, products, and TSs. MRSDCI provides the electronic energy for each of the optimized structures. The MRSDCI energies are corrected for size-extensivity errors (vide infra) and then extrapolated to the CBS limit. This composite approach has been shown to achieve chemical accuracy when predicting bond dissociation energies of hydrocarbons,28 and here, it is tested for H abstraction barrier heights. The highpressure-limit rate coefficients are predicted by TST with the torsional motions treated as hindered internal rotors, while all other contributions to the partition functions are treated as before (IGRRHO approximations). All geometry optimizations were performed using quasiNewton searches and the quadratic approximation to determine step size and direction within the GAMESS-US code.36 The structures of each stationary point on the reaction potential energy surface were optimized with DFT using the B3LYP37 exchange-correlation (XC) functional, which has been shown to predict accurate geometries.38 The large cc-pVTZ basis set39 used ensures negligible basis set superposition error (BSSE), 8432
dx.doi.org/10.1021/jp304811z | J. Phys. Chem. A 2012, 116, 8431−8443
The Journal of Physical Chemistry A
Article
electrons to not be important for the ∼1 kcal mol−1 accuracy threshold, so no excitations of core electrons were included. In MRSDCI, the truncation of the CI expansion to single and double excitations leads to a size-extensivity error. Although the size extensivity error grows with molecular size, it causes nontrivial errors even in the small molecules considered here. Two kinds of size-extensivity corrections are widely used, the a posteriori Davidson-type corrections47,48 and the a priori selfconsistent multireference averaged coupled pair functional (MRACPF) method.49,50 In order to test and validate a correction scheme for this error, both the Davidson−Silver (DS) correction47 and the MRACPF method were used. The DS correction is evaluated a posteriori, i.e., after the MRSDCI calculation, and can be evaluated at essentially no additional cost. MRACPF, unlike the DS correction, is a self-consistent method to correct the MRSDCI energy:49 it compensates for missing excitations that would have made the wave function size-extensive by modifying the renormalization of the wave function.50 However, it tends to slightly overestimate correlation contributions. The MRACPF-2 approach uses a different renormalization that improves both the accuracy and stability of MRACPF.50 Therefore, MRACPF-2 (hereafter referred to as MRACPF for simplicity) was used in our work, and the calculations were performed with our TIGER-CI code.51−54 The MRSDCI+DS correction and MRACPF calculations were performed with the cc-pVDZ and cc-pVTZ basis sets; the resulting energies were then extrapolated to the CBS limit, denoted as, e.g., MRSDCI+DS/cc-pV∞Z. As in ref 28, we employ the two-point extrapolation scheme proposed by Truhlar.55 This scheme requires both a reference energy and a correlation energy to extrapolate. We used the CASSCF total energy and the energy difference between MRSDCI+DS (or MRACPF) and CASSCF total energies as reference and correlation energies, respectively, with extrapolation parameters α = 3.4 and β = 2.4. On the basis of barrier heights computed as described above, we predict the high-pressure-limit rate constants according to TST using the CANTHERM code.56 The IGRRHO approximations are applied for translational, rotational, and vibrational partition functions, respectively. However, low frequency modes play an important role in the vibrational partition function. Some low frequency torsional modes act as hindered internal rotations at high temperature,57 and we treat these modes, torsional rotations along a single bond, within the onedimensional, separable, hindered-rotor approximation. We obtain these rotational barriers by scanning the PES at fixed dihedral angles in 10° increments using DFT-B3LYP/cc-pVTZ. For equilibrium geometries, all degrees of freedom are optimized except the rotational dihedral angle being scanned. For TSs, in addition to the rotational dihedral angle being scanned, the lengths of the breaking and forming bonds are fixed at the TS values during optimization. The scanned potential is fitted to a fifth-order Fourier series in the torsional angle ϕ; V(ϕ) = ∑5m = 0Am cos(mϕ) + Bm sin(mϕ). Using the reduced moment of inertia58,59 of the optimized geometry, the one-dimensional rotational Schrödinger equation with the fitted potential is solved numerically. The canonical partition function thus obtained then replaces what the HO approximation would compute for these torsional modes. The partition function combines the coupled external and internal rotations and is calculated using the semiclassical Pitzer−Gwinn approximation.58 The final TST rate coefficients include the asymmetric
which is important in cases where bimolecular complexes form (vide infra). Restricted DFT was used for closed-shell molecules, whereas open-shell molecules were described with unrestricted (spin-polarized) DFT. Harmonic vibrational frequencies computed via diagonalization of numerical Hessians constructed from finite differences of analytic gradients were used to determine the nature of each stationary point, as well as to compute ZPEs and TC terms (to 298.15 K). The scale factor of 0.9691 recommended for DFT-B3LYP/cc-pVTZ40 was applied to the normal-mode frequencies in order to account for anharmonicity prior to computation of ZPEs and TCs. Hydrogen bonding in reactions involving HO2 and OH radicals results in reactant and product complexes with multiple conformers (vide infra). Intrinsic reaction coordinate calculations at the DFT-B3LYP/cc-pVDZ level of theory were performed to determine the reactant and product complex geometries. These structures were then further refined at the DFT-B3LYP/cc-pVTZ level. However, DFT-B3LYP predicted a barrierless process for OH abstracting the formyl H from MF. Approximate XC functionals, such as B3LYP, are well-known to underestimate reaction barrier heights41 due to self-interaction and delocalization errors.42 However, the M06-2X functional43 has been shown to predict TS structures with small mean errors.44 Consequently, we used M06-2X to obtain structures and frequencies of all species in the MF plus OH radical reaction. We continued using the B3LYP frequency scale factor (0.9691) for consistency; the recommended scale factor for M06-2X, 0.9735,45 is only slightly different and produces negligible energy changes (∼0.005 kcal mol−1). Apart from these exceptions, DFT-B3LYP predicts nearly identical TS structures as DFT-M06-2X (
