488
Energy & Fuels 2006, 20, 488-497
Density Functional Theory Calculations of the Energetics and Kinetics of Jet Fuel Autoxidation Reactions§ Steven Zabarnick*,† and Donald K. Phelps‡ UniVersity of Dayton Research Institute & UniVersity of Dayton Mechanical & Aerospace Engineering Department, 300 College Park, Dayton, Ohio 45469, and Air Force Research Laboratory, Propulsion Directorate, 1790 Loop Rd. N., Wright-Patterson Air Force Base, Ohio 45433 ReceiVed October 21, 2005. ReVised Manuscript ReceiVed December 16, 2005
Density functional theory calculations of the energetics and kinetics of important reactions for jet fuel oxidation are reported. The B3LYP functional along with 6-31G(d) and larger basis sets are used for calculation of peroxy radical abstraction reactions from hydrocarbons and heteroatomic species, the reaction of sulfides, disulfides, and phosphines with hydroperoxides to produce nonradical products, and the metal catalysis of hydroperoxide decomposition. Reaction enthalpies and activation energies are determined via DFT calculations of the structures and energies of stable species and transition states. The peroxy radical abstraction study shows the high reactivity (Ea’s of 6-11 kcal/mol) of the H atoms which are weakly bonded to heteroatoms, including nitrogen, oxygen, and sulfur. These species, at part-per-million levels, are able to compete for peroxy radicals with the bulk fuel hydrocarbon species. Benzylic hydrogens on aromatic hydrocarbons are shown to be significantly more reactive (by 4 to 5 kcal/mol) than paraffinic hydrogens with the result that the aromatic portion of fuel sustains the bulk of the autoxidation process. Sulfides and disulfides are found to react readily with fuel hydroperoxides (Ea’s of 26-29 kcal/mol) to produce alcohols and the oxidized sulfur species. Triphenylphosphine reacts with hydroperoxides with a very low activation energy (12.9 kcal/mol). The metal catalysis of hydroperoxide decomposition is calculated to occur through the formation of a complex with subsequent decomposition to form radical species without regeneration of the metal ion. The reaction pathways found and activation energies calculated can be used to improve chemical kinetic models of fuel autoxidation and deposition.
Introduction Prior to being combusted for propulsion, jet fuel is heated during passage through aircraft fuel system components. This heating occurs incidentally while passing through fuel pumps, but it is encouraged via heat exchangers, particularly in advanced military aircraft, to remove excess heat from numerous aircraft subsystems. Systems which may require cooling include avionic, hydraulic, lubrication, and environmental control. The use of fuel to cool fuel system and combustion components is an enabling technology for advanced military aircraft because of the large quantity of excess heat produced. Unfortunately, the heat absorbed by the fuel is not always innocuous. When fuel temperatures approach ∼140 °C, the fuel begins to react via an autoxidation chain mechanism with the small amount of dissolved oxygen (65-75 ppm wt/wt)1 present from exposure to air. These autoxidation reactions ultimately result in formation of detrimental surface deposits and bulk insolubles.2 These deposits can plug narrow passageways in valves, filters, and * To whom correspondence should be addressed. § Disclaimer: The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright notation thereon. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of Air Force Research Laboratory or the U.S. Government. † University of Dayton. ‡ Air Force Research Laboratory. (1) Striebich, R. C.; Rubey, W. A. Prepr.-Am. Chem. Soc., DiV. Pet. Chem. 1994, 39, 47-50. (2) Hazlett, R. N. Thermal Oxidation Stability of AViation Turbine Fuels; ASTM: Philadelphia, 1991.
nozzles and can inhibit desired heat transfer in heat exchangers. Numerous techniques have been investigated to limit the formation of deposits, including: fuel system designs to minimize fuel temperatures, fuel additives to inhibit autoxidation and deposit formation, fuel deoxygenation, fuel system surface coatings, and inclusion of “sacrificial” or coke tolerant components.3 Each of these solutions has corresponding advantages and disadvantages which are dependent on the application, but no single method is able to eliminate the deposition problem under all current and proposed aircraft fuel system conditions. In recent years, chemical kinetic models have been developed which simulate the major autoxidation pathways that occur in jet fuels.4-6 Development of a widely applicable autoxidation mechanism, which enables the prediction of deposit formation, would greatly aid the fuel system design process and allow more efficient use of the fuel as a heat sink.7 As jet fuels consist of thousands of individual species, which vary in their identity and concentration in different fuel samples, it is impractical to build detailed chemical kinetic mechanisms. Grouped or lumped mechanisms, sometimes referred to as “pseudo-detailed” mech(3) Heneghan, S. P.; Zabarnick, S.; Ballal, D. R.; Harrison, W. E. J. Energy Res. Technol. 1996, 118, 170-179. (4) Zabarnick, S. Ind. Eng. Chem. Res. 1993, 32, 1012-1017. (5) Zabarnick, S. Energy Fuels 1998, 12, 547-553. (6) Kuprowicz, N. J.; Ervin, J. S.; Zabarnick, S. Fuel. 2004, 83, 17951801. (7) Balster, L. M.; Zabarnick, S.; Ervin, J. S.; Striebich, R.; DeWitt, M. J.; Doungthip, T.; Kuprowicz, N. J. Predicting the Thermal Stability of Jet Fuel: Analytical Techniques Toward Model Validation. Presented at the 8th International Conference on Stability and Handling of Liquid Fuels, Steamboat Springs, CO, 2004.
10.1021/ef050348l CCC: $33.50 © 2006 American Chemical Society Published on Web 01/20/2006
Jet Fuel Autoxidation Reactions
anisms, have been used to simulate the most important reactive pathways, including the effects of antioxidants and catalytic surfaces.4-6 In addition, these mechanisms have been combined with computational fluid dynamics techniques with the goal of simulating the complex time and temperature variation during fuel flow in aircraft fuel system components.8,9 Most recently, initial efforts aimed at including global deposit formation reactions in these mechanisms have been performed.7,10 Unfortunately, experimental kinetic parameters for autoxidation reactions are often not available for inclusion in these chemical kinetic mechanisms. It is desirable to have a consistent and reliable method to estimate such parameters. Recently, researchers have begun to apply semiempirical quantum mechanical11 and density functional theory12 methods to determine Arrhenius parameters for autoxidation reactions. These computational techniques promise to provide methods to readily estimate the rate parameters and evaluate the plausibility of candidate autoxidation and deposit-forming reactions. In the current work, we explore the use of density functional theory (DFT) quantum chemistry methods, specifically the B3LYP hybrid functional, for the calculation of the energetics and kinetics of autoxidation reactions. These DFT techniques promise to provide an accurate and consistent method for estimating the reaction energetics and Arrhenius parameters of fuel autoxidation reactions. Computational Methodology Calculations were performed using Gaussian 03 software (Gaussian, Inc.) on the USAF Aeronautical Systems Center Major Shared Resource Center computational facilities. All frequency and thermodynamic calculations were performed at 25 °C and one atmosphere of pressure. Transition state geometries were calculated using the synchronous transit-guided quasi-Newton (STQN) method or via saddle-point geometry optimizations starting with educated geometry estimates. All transition states reported were confirmed to contain a single imaginary frequency. Animations of the normal mode coordinate corresponding to the imaginary frequency of the transition state and intrinsic reaction coordinate (IRC) calculations were performed to confirm the connection between the reactants and products. Activation energies were calculated from the difference between the zero-point and thermal energy-corrected enthalpies of the transition state and the reactants. All calculations used unscaled frequencies, as the scaling factors for the functional/basis set combinations reported here result in negligible changes in the resulting thermodynamic and kinetic parameters. Unrestricted calculations were performed on open shell species (i.e., radicals and transition states), and restricted calculations were performed on closed shell (stable) species (i.e., the default behavior of Gaussian 03 was used). Basis set superposition error corrections were performed on the transition state complexes of the peroxy radical abstraction reactions by the counterpoise correction technique described below. Basis Set Superposition Error. It is well-known that errors occur in ab initio and DFT calculations of the energetics of intermolecular interactions in complexes because of the effect of the basis functions of one molecule augmenting the other molecule. The error occurs because the individual molecule energetics are usually calculated without the presence of the basis functions of the additional species. This so-called basis set superposition error (8) Ervin, J. S.; Zabarnick, S. Energy Fuels. 1998, 12, 344-352. (9) Doungthip, T.; Ervin, J. S.; Zabarnick, S.; Williams, T. F. Energy Fuels 2004, 18, 425-437. (10) Doungthip, T. Ph.D. Thesis, University of Dayton, Dayton, OH, 2004. (11) Nikolaeva, E. V.; Shamov, A. G.; Khrapkovskii, G. M.; Kharlampidi, K. E. Russ. J. Gen. Chem. 2002, 72, 748-759. (12) Denisova, T. G.; Emel’yanova, N. S. Kinet. Catal.. 2003, 44, 441449.
Energy & Fuels, Vol. 20, No. 2, 2006 489 Table 1. Calculated Bond Dissociation Enthalpies for the O-O Bond in n-Butylhydroperoxide level of theory/basis set
bond dissociation enthalpy (kcal/mol)
HF/6-31G(d) B3LYP/6-31G(d) B3LYP/6-31G(d,p) B3LYP/6-311G(d,p) MP2/6-31G(d) CBS-4M G3B3 exptl (based on analogous reactions)18,19
0.3 40.0 40.0 37.9 48.8 45.0 43.5 44-45
(BSSE) results in an incorrect lowering of the overall energy of the complex. Techniques of correcting for BSSE’s, such as the counterpoise correction of Boys and Bernardi,13 have been used extensively for the calculation of intermolecular complexes but have not generally been employed for transition states despite the fact that calculations of transition state interactions should lead to similar errors. Recently, counterpoise corrections of BSSE’s have been used for a limited number of transition state calculations with good success.14,15 Here, for the peroxy radical abstraction reactions studied, we use the counterpoise correction method contained within Gaussian 03, which utilizes the techniques of Boys and Bernardi13 and Simon et al.16
Results and Discussions Justification for Level of Theory and Basis Set Used and Technique Validation. The selection of the quantum chemical technique and basis set used for a given calculation represents a compromise between computational cost and accuracy. We would like to select the method that yields the most accurate results in a reasonable period of time. In addition, techniques which include electron correlation are required for the bond breaking and forming processes studied here. Unfortunately, the relatively large molecules which constitute jet fuel preclude the use of high accuracy compound methods, such as G3 theory and complete basis set methods, for calculation of transition states. Density functional theory, and in particular the B3LYP hybrid functional, is increasingly being proven to yield very good geometries and energies with relatively modest computational costs.17 Thus, it is particularly attractive for the relatively large molecular systems of interest in fuel autoxidation. Initial validation runs for selecting the DFT functional/basis set combination were conducted on the calculation of the bond dissociation enthalpy of the O-O bond in n-butylhydroperoxide. Hydroperoxides are known to be important initiators of autoxidation and intermediates in the overall chain process. nButylhydroperoxide was chosen as a small model molecule of the larger hydroperoxides that are produced in fuel. Previous work indicates that the bond dissociation enthalpies of alkyl hydroperoxides, essentially independent of the identity of the alkyl group, is 44-45 kcal/mol.18,19 Calculated bond dissociation enthalpies at various levels of theory are shown in Table 1. The results indicate, as expected, that using methods which include the effects of electron correlation (all except the HF result) are (13) Boys, S. F.; Bernardi, F. Mol. Phys. 1970, 19, 553. (14) Kobko, N.; Dannenberg, J. J. J. Phys. Chem. A 2001, 105, 19441950. (15) Sordo, J. A. THEOCHEM 2001, 537, 245-251. (16) Simon, S.; Duran, M.; Dannenberg, J. J. J. Chem. Phys. 1996, 105, 11024-11031. (17) Barckholtz, C.; Barckholtz, T. A.; Hadad, C. M. J. Am. Chem. Soc. 1999, 121, 491-500. (18) Sebbar, N.; Bockhorn, H.; Bozzelli, J. W. Phys. Chem. Chem. Phys. 2002, 4, 3691-3703. (19) Reints, W.; Pratt, D. A.; Korth, H.-G.; Mulder, P. J. Phys. Chem. A 2000, 104, 10713-10720.
490 Energy & Fuels, Vol. 20, No. 2, 2006
Zabarnick and Phelps
essential to obtain accurate values. The B3LYP DFT functional is similar in computational cost to HF but yields more accurate values. Ab initio (MP2), complete basis set (CBS-4M), and compound techniques (G3B3) yield even better results at a significantly higher computational cost. It was concluded that the DFT method resulted in the best compromise between accuracy and computational cost for the relatively large reaction systems that need to be studied in the autoxidation chemistry of jet fuel. All subsequent calculations reported here employ the B3LYP level of theory. The transition state calculations performed here are particularly costly and thus a relatively modest basis set, 6-31G(d), was selected for most of the calculations. We have performed limited studies to assess the change in accuracy of these calculations with increasing basis set size (e.g., 6-311G(d,p) and 6-311G(3df, 2p)), as reported below. Peroxy Radical Abstraction Reactions. One of the most important reactions in autoxidation is the propagation reaction in which peroxy radicals abstract a hydrogen atom from a fuel species, producing a hydroperoxide and a hydrocarbon radical
RO2 + RH f ROOH + R
(reaction 1)
This reaction tends to be the rate-limiting step in autoxidation, as the addition of O2 to hydrocarbon radicals, R + O2 f RO2, is a fast reaction. In addition, many antioxidants function by competing for peroxy reactions via
RO2 + AH f ROOH + A
(reaction 2)
where AH is the antioxidant species and A is an antioxidant radical that does not propagate the autoxidation chain. As the strength of the O-H bond (85-86 kcal/mol) in hydroperoxides tends not to vary significantly with the identity of the R group,18 the reaction energetics are primarily determined by the strength of the R-H or A-H bond. It has been shown that antioxidants function by rapidly intercepting peroxy radicals, which occurs readily because of their relatively weak A-H bonds.20,21 In addition to synthetic antioxidants which are added after fuel processing, fuels can contain naturally occurring species which may have antioxidant properties. These species, such as phenols, thiols, amines, indoles, carbazoles, etc., though only present at low levels (ppm), determine the oxidation rate and the tendency of the fuel sample to produce deposits upon oxidation. Thus it is important to know the energetics of these reactions as well as the R-H and A-H bond strengths. Unfortunately, very few measurements or calculations have been performed on the energetics or activation energies of reactions 1 and 2, although there are many measurements and calculations of bond strengths of some relevant R-H and A-H species.17,22,23 Recently, Nikolaeva et al.11 used semiempirical quantum chemical methods (PM3) and DFT to estimate the energetics of some of these reactions. Denisova and Emel’yanov12 used the intersecting parabolas method and DFT methods to estimate transition state geometries and energetics of peroxy radical reactions with various species. In the present paper, we utilize DFT for the calculation of the activation energies of peroxy radical abstraction reactions from hydrocarbons and heteroatomic species. These species were (20) Heneghan, S. P.; Zabarnick, S. Fuel 1994, 73, 35-43. (21) Scott, G. Chem. Ind. 1963, 271-281. (22) Nishiyama, T.; Suzuki, T.; Hashiguchi, Y.; Shiotsu, S.; Fujioka, M. Polym. Degrad. Stab. 2002, 75, 549-554. (23) Bordwell, F. G.; Zhang, X.; Cheng, J.-P. J. Org. Chem. 1991, 56, 3216-3219.
chosen to represent classes of species commonly found in jet fuels and include representatives of the following classes: alkanes, alkyl-substituted aromatics, alkyl-substituted phenols, indoles, carbazoles, and thiols. We have chosen the 2-butylperoxy radical as the abstracting radical in this study. This peroxy radical was chosen as a compromise between the larger hydrocarbon peroxy radicals (e.g., those with 9-16 carbon atoms) that likely dominate jet fuel autoxidation and the smaller radicals, such as methylperoxy, which may be too small and polar to approximate the reactivity of larger peroxy radicals. Analysis of the results obtained below demonstrates that the peroxy radicals most likely encountered in jet fuels are obtained from the reaction of fuel aromatic species, but it is likely that the hydrocarbon group of the peroxy radical has little effect on its reactivity, as the bond strength of the O-H bond in various hydroperoxides is largely independent of the identity of the hydrocarbon group.18 We show below, for a single abstraction reaction, that the identity of the abstracting peroxy radical (2butylperoxy vs benzylperoxy) has no effect on the activation energy obtained. In addition, we have also used DFT for the calculation of the bond strengths of the weakest bonds to the hydrogen atom for relevant fuel species. We have calculated the reaction enthalpies and activation energies with the B3LYP functional using the 6-31G(d) basis set for the reaction of 2-butylperoxy with a number of model fuel species, as shown in Table 2. Initial results indicated a problem with the values of the activation energies obtained for many of these reactions. As shown in the table, for many of these reactions, which tend to be slightly endothermic, the calculated activation energies were slightly lower than the calculated reaction enthalpies using the same level of theory and basis set. Obviously, a reaction cannot have an activation energy that is lower than the reaction enthalpy. Here we find that basis set superposition error (BSSE) can account for the activation energies that are too low. BSSE is well-known in the calculation of energies of intermolecular complexes (e.g., hydrogen-bonded complexes) but has only recently been applied to transition states.14,15 These errors occur when the basis functions of one species in the complex or transition state augments the basis functions of the other species, resulting in a lower overall energy of the complex or transition state. BSSE results in activation energies that are too low relative to the values obtained in the absence of basis function augmentation. BSSE can be corrected by the application of the counterpoise correction, in which the energies of the individual reactants are calculated in both the absence and presence of the augmented basis set resulting from other species.13 The error caused by this augmentation can then be assessed and subtracted from the system. In the present work, we have chosen to use the counterpoise correction as implemented in Gaussian 03, with the use of three molecular fragments. The three fragments chosen for each of the peroxy abstraction reactions are (1) the hydrogen being abstracted, (2) the peroxy radical, and (3) the resulting radical being formed. For example, in the case of the reaction of 2-butylperoxy and 3-ethylphenol, the three fragments used in the counterpoise correction are the phenol hydrogen, the 2-butylperoxy radical, and the 3-ethylphenoxy radical. Table 2 lists the activation energies calculated with the inclusion of the counterpoise correction for the BSSE. The use of the counterpoise correction results in an increase in the calculated activation energies from 3.1 to 4.7 kcal/mol relative to values obtained without the use of the correction. The activation energies calculated with the counterpoise correction are now more reasonable, being greater than the calculated reaction
Jet Fuel Autoxidation Reactions
Energy & Fuels, Vol. 20, No. 2, 2006 491
Table 2. Calculated Reaction Enthalpies, Activation Energies, and Bond Strengths Using B3LYP/6-31G(d) for the Abstraction of Hydrogen from Model Fuel Species by 2-Butylperoxy Radicals reacting species (hydrogen atom being abstracted)
∆H (kcal/mol)
n-butane (secondary hydrogen) ethylcyclohexane (ethyl CH2 hydrogen) ethylcyclohexane (3-ring CH2 hydrogen) ethylbenzene (benzylic hydrogen) cumene (benzylic hydrogen)
20.9 20.5 21.2 9.2 7.5
hydrocarbon species 20.1 19.3 21.6 15.4 14.8
23.9 (20.9)a 23.6 25.0 18.9 (16.5)a 18.9
3-ethylphenol (phenolic hydrogen) 3-ethylphenol (phenolic hydrogen) reaction with methylbenzylperoxy 2,6-di-tert-butyl-4-methylphenol (phenolic hydrogen) 3-ethylbenzenamine (amine hydrogen) 2-methylindole (indole hydrogen) 3-methylindole (indole hydrogen) carbazole (carbazole hydrogen) 2-methylindoline (indole hydrogen) 2,5-dimethylpyrrole (pyrrole hydrogen) ethylphenylthiol (thiol hydrogen)
3.4 3.4
heteroatomic species 1.8 1.6
6.0 (5.9)a 6.0
78.9 78.9
6.9 (6.7)a
70.8 (79.7-82.6)b
Ea without BSSE correction (kcal/mol)
-4.6
2.2
9.9 10.4 8.7 9.4 3.6 7.6 -2.4
5.6 6.8 6.5 5.6 2.4 6.7 4.3
Ea with BSSE correction (kcal/mol)
9.6 11.0 10.1 9.8 6.6 11.2 7.4
X-H bond strength (kcal/mol) 96.4 (98.3)b 96.0 96.6 84.6 (85.4)b 83.0
85.4 85.8 84.2 84.9 (93.6)b 79.0 83.0 73.1
a Activation energy values in parentheses are single point energies calculated via B3LYP/6-311G(d,p) DFT theory using B3LYP/6-31G(d) geometries with the counterpoise correction for the transition state. b Bond strength values in parentheses are experimental values tabulated by Luo.25 All imaginary frequencies for these reactions were between -1626 and -1935 cm-1.
enthalpy for all but one reaction (3-ethylbenzenamine). The cause of the benzenamine discrepancy is not known at this time. Presumably, the counterpoise correction would be successful in correcting previous calculations in which activation energies obtained appear too low or even negative.24 The magnitude of the BSSE correction obtained here is relatively small compared to many reactions which exhibit moderate to large activation energies (>20 kcal/mol) and thus the correction may be unnecessary for such transition state calculations and has not been employed for the other reaction types studied in this paper (see sections below). In contrast, for the slightly endothermic reactions being studied here, the correction is quite significant relative to the reaction enthalpy. The first five species in Table 2 are hydrocarbons which represent the types of these species present in fuel. n-Butane represents the normal and branched alkane fuel components that along with cycloalkanes and aromatic species compose the bulk of jet fuel. The n-butane reaction studied is the abstraction of one of the secondary hydrogens by the peroxy radical. The cycloalkanes are represented by ethycyclohexane, and calculations were performed on the abstraction of both an ethyl CH2 hydrogen and a CH2 hydrogen on the 3-ring carbon. The fuel aromatics are represented by ethylbenzene and cumene (isopropylbenzene), which contain relatively reactive benzylic hydrogens. The enthalpy of the reaction calculations show that the enthalpies are similar for the n-butane and cyclohexane species but are substantially lower for ethylbenzene and cumene. The results reflect the well-known decrease in C-H bond strength through the series: secondary hydrogen, benzylic hydrogen, tertiary benzylic hydrogen. Our calculated bond strengths shown in the table also support this interpretation and agree well (within 2 kcal/mol) with the experimentally measured values for these C-H bonds.25 The results show that the activation energies for n-butane and ethylcyclohexane are similar, while carbon bonds to benzylic hydrogens exhibit significantly lower activation energies. The changes in activation energy through the series of hydrocarbon species is reflected
in the corresponding C-H bond strengths and enthalpies of reaction. Also shown in the table is the counterpoise-corrected activation energies for a few of the hydrocarbon reactions using the B3LYP/6-31G(d) geometries with a single-point energy calculation with a larger basis set, B3LYP/6-311G(d,p). The larger basis set lowers the activation energies by 2.4 and 3.0 kcal/mol for the two hydrocarbons shown. Calculations using even larger basis sets (not shown), such as 6-311G(3df,2p), yield only slightly lower (
