Theoretical Studies of the Hydrogen Abstraction from Poly(oxymethylene) Dimethyl Ethers by O2 in Relation with Cetane Number Data.

Poly(oxymethylene) dimethyl ethers (PODME n , n = 2-6) are novel oxygenated compounds that can be used as promising candidates for new-generation fuels because of their excellent combustion performance. The oxidation of PODME n (n = 2-6) is essential for the understanding of the combustion process. It is necessary to study the relationship between kinetic parameters and cetane number (CN) of PODME n (n = 2-6). In order to predict initiation step rate constants for the oxidation of PODME n (n = 2-6), quantum mechanical calculations are performed using M06-2X/6-311G(d,p) and B3LYP/6-311G(d,p) methods. Structural, energetic, thermodynamics, and kinetics of the automatic ignition process are assessed. The kinetic model based on the conventional transition state theory is used to calculate the initiation step reaction rate constants at 1000 K. In both M06-2X/6-311G(d,p) and B3LYP/6-311G(d,p) methods, the calculated rate constants of the hydrogen abstraction process have an excellent correlation with the experimental CN of PODME n (n = 2-6). Our methodology presented here can be used to simulate chemical kinetics for other fuel additives.

Introduction

Cetane number (CN) is an important index to evaluate the ignition performance of diesel, whose value has important influence on diesel engine cold start, emission performance, and combustion noise.[1−3] The ignition quality of the diesel fuel is closely related to the chemical composition.[4−6] Oxygenated compounds (molecules which contain oxygen in their structure) can be used as additives or substitutes to diesel fuels, which could greatly reduce nitrogen oxides (NO) and particulate pollutants in tail gas and improve the combustion performance of diesel fuels.[6−8] As a new type of diesel additive, poly(oxymethylene) dimethyl ethers (PODME) can significantly improve the combustion performance of diesel fuels, effectively improve the thermal efficiency, and greatly reduce NO and soot emissions.[9−12] In particular, the average CN of PODME (CH3O(CH2O)CH3 with 2 ≤ n ≤ 6) is more than 76, the oxygen content is up to 47–50%, the boiling point is more than 160 °C, and all the indices fall within the range of diesel indices.[9−12] Furthermore, PODME have a good miscibility with an ordinary diesel fuel and are an environmentally friendly diesel blend components with great application prospects.[11,13] The combustion and emission characteristics of the (diesel + PODME3–4) mixtures with the PODME3–4 volume blending ratios 10 and 20% have been investigated in a light-duty direct injection diesel engine.[14] The results show that PODME3–4 can significantly improve the engine efficiency and reduce soot emissions.[14] The effects of (diesel + PODME3–6) mixtures on the combustion and emission characteristics of a heavy-duty diesel engine have been studied with the PODME3–6 volume blending ratios 15 and 25%.[15] Compared with diesel without PODME3–6, the soot emissions could be reduced by 88 and 95% with the addition of 15 and 25% PODME3–6, respectively.[15] Therefore, PODME (n = 2–6) have potential to be an environmentally friendly additive that significantly reduces soot formation in diesel engines. The general trend of CN value increasing with the chain length is still valid for PODME (n = 2–6).[9−12]

To understand the effect of CN improver on diesel combustion, it is important to understand the initiation step that the CN improver creates free radicals qualitatively.[16,17] Therefore, a great deal of research has been done on this topic.[17−21] In general, for the hydrogen abstraction reactions between organic compounds and oxygen, shell opening oxygen induces the formation of reactive free radicals, followed by branched chain oxidation.[17−21] The reactive free radicals produced in the initiation step are generated by the reaction of fuel (RH) with oxygen.[4−6,17,21] The mechanism of hydrogen abstraction reactions is simple, which can be illustrated clearly by quantum chemistry theory.[20,21] Numerous related computational studies have been carried out in the literature ranging from the Hartree–Fock (HF) theory to hybrid density functional theory (DFT), such as the popular B3LYP, BH&HLYP, M06-2X, and M08-HX methods.[4−6,16,22] Kaliaguine et al. had studied the theoretical kinetic on a few simple compounds, such as n-paraffins (C2–C6), primary linear alcohols (C2–C5), dimethyl ether, diethylether, and methylal.[4−6] The results were consistent with the hypothesis that the rate constants in the initiation step of the autoignition process of compounds being studied were correlated with the experimental CN of these oxygenates.[4−6]

To our knowledge, theoretical kinetic data of the autoignition process of the PODME and its correlation with the corresponding experimental CN data have not been reported in the literature. Hence, the objective of this work is to study the theoretical kinetic data of the autoignition process of PODME (n = 2–6) and to study the relationship between quantum calculation results and the experimental CN of these compounds.

Results and Discussion

Structures and Energetics

In the combustion or general oxidation of RH, a bimolecular gas-phase reaction is considered as the initiation step, given as follows[4,17]

In this reaction, an oxygen molecule reacting by H-abstraction with the initial organic compounds leads to the formation of the radicals R• and HO2• via a molecular transition state (TS) [R···H···O2].

The optimal geometric structures of reactants, products, and TSs that connect reactants and products are searched. The geometrical optimizations of PODME (n = 2–6) molecules show the distorted helix structure and minimum energy, similar to the gg-helix of polyoxymethylene.[23,24] The geometry optimized at the B3LYP/6-311G(d,p) level of theory is very similar to the structure obtained with M06-2X/6-311G(d,p). Figure shows the structure of PODME5 using the M06-2X/6-311G(d,p) method as an illustration. The Mulliken charges of atoms of PODME5 are shown in Table . Although some hydrogen atoms bonded to the same carbon atom are not rigorously equivalent to each other, the energy differences between the TSs of the H-abstraction reactions are <0.2 kcal/mol, which could be ignored. Hence, it is reasonable to treat them as identical. For example, according to the positions of hydrogen atoms of PODME5, these hydrogen atoms can be divided into four groups. For comprehensive study of the reaction mechanism, all O2-attack possibilities have been taken into account in Table based on the different hydrogen sites. Table S1 (Supporting Information) presents the total electronic energies [including zero-point energy (ZPE)] of the various reactants, TSs, and products at the B3LYP and M06-2X levels of the theory used.

Figure 1

Configurations of PODME5 using the M06-2X/6-311G(d,p) method.

Table 1

Mulliken Charges of Atoms of PODME5 Using B3LYP and M06-2X Methods

		Mulliken charge
label	atom	B3LYP	M06-2X
1	C	–0.124	–0.197
2	H	0.116	0.140
3	H	0.109	0.137
4	H	0.091	0.119
5	O	–0.357	–0.368
6	C	0.144	0.084
7	H	0.117	0.150
8	H	0.095	0.130
9	O	–0.374	–0.380
10	C	0.137	0.073
11	H	0.116	0.151
12	H	0.116	0.152
13	O	–0.370	–0.376
14	C	0.136	0.072
15	H	0.115	0.150
16	H	0.115	0.150
17	O	–0.370	–0.376
18	C	0.137	0.073
19	H	0.116	0.152
20	H	0.116	0.151
21	O	–0.374	–0.380
22	C	0.144	0.084
23	H	0.095	0.130
24	H	0.117	0.150
25	O	–0.357	–0.368
26	C	–0.124	–0.197
27	H	0.116	0.140
28	H	0.109	0.137
29	H	0.091	0.119

Table 2

Possible Reaction Pathways of H-Abstraction of PODME (n = 2–6) by O2

systems	possible O2-attack pathways
PODME2	CH3O(CH2O)2CH3 + O2 → CH3O(CH2O)2CH2• + HO2•	R1a
	CH3O(CH2O)2CH3 + O2 → CH3OCH2OC•HOCH3 + HO2•	R1b
PODME3	CH3O(CH2O)3CH3 + O2 → CH3O(CH2O)3CH2• + HO2•	R2a
	CH3O(CH2O)3CH3 + O2 → CH3O(CH2O)2C•HOCH3 + HO2•	R2b
	CH3O(CH2O)3CH3 + O2 → CH3OCH2OC•HOCH2OCH3 + HO2•	R2c
PODME4	CH3O(CH2O)4CH3 + O2 → CH3O(CH2O)4CH2• + HO2•	R3a
	CH3O(CH2O)4CH3 + O2 → CH3O(CH2O)3C•HOCH3 + HO2•	R3b
	CH3O(CH2O)4CH3 + O2 → CH3O(CH2O)2C•HOCH2OCH3 + HO2•	R3c
PODME5	CH3O(CH2O)5CH3 + O2 → CH3O(CH2O)5CH2• + HO2•	R4a
	CH3O(CH2O)5CH3 + O2 → CH3O(CH2O)4C•HOCH3 + HO2•	R4b
	CH3O(CH2O)5CH3 + O2 → CH3O(CH2O)3C•HOCH2OCH3 + HO2•	R4c
	CH3O(CH2O)5CH3 + O2 → CH3O(CH2O)2C•HO(CH2O)2CH3 + HO2•	R4d
PODME6	CH3O(CH2O)6CH3 + O2 → CH3O(CH2O)6CH2• + HO2•	R5a
	CH3O(CH2O)6CH3 + O2 → CH3O(CH2O)5C•HOCH3 + HO2•	R5b
	CH3O(CH2O)6CH3 + O2 → CH3O(CH2O)4C•HOCH2OCH3 + HO2•	R5c
	CH3O(CH2O)6CH3 + O2 → CH3O(CH2O)3C•HO(CH2O)2CH3 + HO2•	R5d

Due to the TSs of the combustion reaction of PODME (n = 2–6) molecules show similar structural characteristics at the corresponding H-abstraction position, there is just showing the structure of TSs of the combustion reaction of PODME5 in Figures and 3 as an illustration. All the C–H bonds without breaking in the methyl group and methylene group of PODME5 have a length of ∼1.09 Å. In the reaction with the hydrogen atom, the breaking C–H bond is extended by 26–47%. The breaking C–H bond lengths (Å) decrease in the order R4a (1.602/1.591) > R4b (1.532/1.509) > R4c (1.449/1.376) ≈ R4d (1.450/1.384) for B3LYP/M06-2X, while the forming O–H bond lengths increase in the order R4a (1.053/1.050) < R4b (1.071/1.076) < R4c (1.098/1.148) ≈ R4d (1.098/1.144) for B3LYP/M06-2X. B3LYP prediction of the C–H bond is slightly longer than M06-2X but less than by 7%. In all the TS structures, the formed O–H bond is shorter than the broken C–H bond, which indicates that the TS structure is closer to the products than the reactants, which is consistent with the endothermic properties of these reactions. Moreover, O2 can form product complexes with PODME5 via hydrogen bonds between the O atom of O2 and H atom of methylene or methyl hydrogen in PODME5. All TS structures formed are also stabilized by hydrogen bonding, which is between 2.0 and 2.6 Å in this work. Additionally, with an increase in the steric hindrance with H-abstraction position moving to the middle of the molecular chain, the angles of C–H–O increase in the order R4a (126.5°/121.8°) < R4b (141.7°/130.2°) < R4c (148.9°/143.8°) ≈ R4d (148.8°/143.7°) for B3LYP/M06-2X.

Figure 2

Optimized TS structures in PODME5 + O2 at the B3LYP level of theory. (a) R4a; (b) R4b; (c) R4c; (d) R4d.

Figure 3

Optimized TS structures in PODME5 + O2 at the M06-2X level of theory. (a) R4a; (b) R4b; (c) R4c; (d) R4d.

The angles of H–O–O are around 109° for both B3LYP and M06-2X. These are the features of H-abstraction reactions from O2. Notably, with the exception of the H-abstraction from the methyl group and methylene group closed-end position, such as R2b, R3b, R4b, and R5b, H-abstraction from the methylene groups contained in the five molecules shows almost the same activation energy and TS structures.

As shown in Table S1, the reactants are almost stable as the products by about 14.9–17.6 and 20.8–25.3 kcal/mol at B3LYP and M06-2X levels, respectively. Overall, the difference between the total energy of the products and their corresponding TS is smaller at the B3LYP level. In contrast, the total energy of products is about 0.4–6.4 kcal/mol below their corresponding TS at the M06-2X level. The activation energies for H-abstraction reactions of oxygen with PODME (n = 2–6), denoted as R1–R5 (a–d), are reported in Table . The relative energies ΔEI and ΔEII represent the calculated value of activation barrier between the reactants and the TS at B3LYP and M06-2X levels, respectively. It is concluded that the B3LYP functional has lower activation barrier than M06-2X with the same H-abstraction processes. This may be the B3LYP underestimate the activation barrier on account of self-interaction and delocalization errors.[25] It can also be observed from Table that the activation energy of the hydrogen atoms at each end of the chain is higher than that of the other hydrogen atoms in the molecular chain. For almost all reaction pathways of PODME (n = 2–6), H-abstraction is more likely to occur on methylene hydrogen than on methyl hydrogen. The Mulliken population analysis is also used to confirm this observation. The net charge of an atom in a molecule gives a rough idea of the electrostatic force between hydrogen and carbon atoms. Table shows the Mulliken population analysis of atomic net charges of PODME5 using B3LYP and M06-2X methods. In PODME5, because the electrophilic oxygen atoms decrease the electron density on carbon atoms, these carbon atoms of the methyl group and methylene group are negatively charged and positively charged, respectively. As a result, the C–H bonds of the methyl group have stronger and more attractive global charge than that of the methylene group.

Table 3

TS Barriers (in kcal/mol) for Hydrogen Abstraction Reactions

reaction	transition state	ΔEI	ΔEII
R1a	[CH3O(CH2O)2CH2···(H···O2)]#	17.35	25.03
R1b	[CH3OCH2OCH···(H···O2)OCH3]#	15.44	21.57
R2a	[CH3O(CH2O)3CH2···(H···O2)]#	17.42	25.09
R2b	[CH3O(CH2O)2CH···(H···O2)OCH3]#	15.63	23.33
R2c	[CH3OCH2OCH···(H···O2)OCH2OCH3]#	14.93	20.88
R3a	[CH3O(CH2O)4CH2···(H···O2)]#	17.52	25.20
R3b	[CH3O(CH2O)3CH···(H···O2)OCH3]#	15.82	23.56
R3c	[CH3O(CH2O)2CH···(H···O2)OCH2OCH3]#	15.06	20.84
R4a	[CH3O(CH2O)5CH2···(H···O2)]#	17.54	25.18
R4b	[CH3O(CH2O)4CH···(H···O2)OCH3]#	15.83	23.64
R4c	[CH3O(CH2O)3CH···(H···O2)OCH2OCH3]#	15.07	20.83
R4d	[CH3O(CH2O)2CH···(H···O2)O(CH2O)2CH3]#	15.10	20.83
R5a	[CH3O(CH2O)6CH2···(H···O2)]#	17.58	25.24
R5b	[CH3O(CH2O)5CH···(H···O2)OCH3]#	15.88	23.59
R5c	[CH3O(CH2O)4CH···(H···O2)OCH2OCH3]#	15.05	20.90
R5d	[CH3O(CH2O)3CH···(H···O2)O(CH2O)2CH3]#	15.12	20.84

Quantitative Study: Kinetic Modeling and Predicting CNs

The rate constant calculation for the bimolecular reactions of the autoignition process is implemented within formalism of the TS theory (TST). All the energy parameters of the reactants and the TS are shown in Table . The rate constant is calculated at 1000 K, which is used to simulate the temperature in the combustion chamber of the diesel engine. Consequently, the rate constant for these initiation steps k is computed using conventional TST[4,6,26] in the followingwhere κ is the transmission coefficient of the corrected tunneling effect, kB is the Boltzmann constant, h is the Planck constants, R is the gas constant, and ΔG0 is the standard state Gibbs free energy of activation at temperature T and pressure P0 = 1 atm. In this work, ΔG0 has been replaced by ΔE in the calculation process of k. Tunneling effects may be described quantitatively as the Wigner correction factor[26,27]where v⧧ (cm–1) is the magnitude of the imaginary frequency of the TST. However, the tunneling effect has little effect at high temperature, such as 1000 K. Furthermore, the imaginary frequencies of all reactions studied here are exactly similar, and the correction of the tunneling effect has no significant effect on the observed relative reactivity. Hence, it makes sense to set κ ≈ 1. The results of H-abstraction reactions are given in Table . It clearly shows that the order of the rate constants of these reactions (R1–R5) is the same in both B3LYP and M06-2X methods. Hence, the rate constants of H-abstraction reactions with PODME (n = 2–6) are in accordance with the order PODME2 < PODME3 < PODME4 < PODME5 < PODME6, which is well correlated with the experimental CN of these pure oxygenated compounds.[9−12]

Table 4

Computed (m × k) Rate Constants in cm3/mol/s for the Bimolecular Reaction Pathways (R1–R5) at 1000 Ka

	m × ki,TST
reaction	B3LYP	M06-2X
R1a (m = 6)	1.79 × 105	4.23 × 10–1
R1b (m = 4)	2.99 × 106	0.98 × 102
R2a (m = 6)	1.59 × 105	3.84 × 10–1
R2b (m = 4)	2.18 × 106	4.95
R2c (m = 2)	3.42 × 106	1.57 × 102
R3a (m = 6)	1.34 × 105	3.20 × 10–1
R3b (m = 4)	1.58 × 106	3.37
R3c (m = 4)	5.71 × 106	3.35 × 102
R4a (m = 6)	1.30 × 105	3.29 × 10–1
R4b (m = 4)	1.56 × 106	2.96
R4c (m = 4)	5.55 × 106	3.38 × 102
R4d (m = 2)	2.65 × 106	1.69 × 102
R5a (m = 6)	1.22 × 105	2.99 × 10–1
R5b (m = 4)	1.43 × 106	3.24
R5c (m = 4)	5.74 × 106	3.02 × 102
R5d (m = 4)	5.14 × 106	3.34 × 102

m is the number of hydrogen atoms with the same chemical environment in PODME (n = 2–6).

Results of Table also indicate that the H-abstraction by O2 from the middle methylene of PODME (n = 2–6) are the dominant paths among all reaction pathways attributed to its lowest barrier height. The additive group contribution approach is used to calculate the global rate constant, which includes the contribution of all H-abstract reactions for PODME (n = 2–6).[28] The global rate constant of the additive rule is expressed aswhere m is the number of hydrogen atoms with the same chemical environment in PODME (n = 2–6) and k is the individual rate constant computed for each category of H-abstraction. The computed total rate constants kglobal and their logarithmic form log(kglobal) for PODME (n = 2–6) combustion reactions (R1–R5) are reported in Table . It is indicated that the kinetic parameters obtained at the M06-2X/6-311G(d,p) level are far less than those at B3LYP/6-311G(d,p).

Table 5

Global Rate Constants Calculated at 1000 K with B3LYP and M06-2X, with kglobal(T) = ∑mk

	log kglobal
systems	B3LYP	M06-2X	CNexp.
PODME2	6.50	1.99	63
PODME3	6.76	2.21	78
PODME4	6.87	2.53	90
PODME5	7.00	2.71	100
PODME6	7.09	2.81	104

There is a significant positive linear correlation between the log(kglobal) calculated at 1000 K for reactions (R1–R5) and the reported CN values for the corresponding oxygenated compounds.[9−12] The linear regression function based on the least-squares method is used and expressed at B3LYP and M06-2X by eqs and 6, respectively

In Figure , the calculated results are compared with the experimental CN data for each oxygenated compound in B3LYP and M06-2X methods. It is clear from Figure that the calculated results are in good agreement with the experimental results[9] for both examined cases. The standard deviation values in B3LYP and M06-2X methods (, where N is the number of experimental data) are 0.80 and 0.67, respectively. The average relative deviation values in B3LYP and M06-2X methods (, where N is the number of experimental data) are 1.91 and 1.58%, respectively. Both correlations exhibit a higher CN value with PODME6 than the reported CN value. It is possible that there is a negative deviation in the experimental value with PODME6.

Figure 4

CN data as a function of log(kglobal) for reactions (R1–R5): (a) at B3LYP, the solid line represents the linear regression function of eq and (b) at M06-2X, the solid line represents the linear regression function of eq .

Conclusions

The oxygenated PODME (n = 2–6) are attractive as clean diesel fuel or fuel additive. In this work, the theoretical kinetic data of the autoignition process of the PODME in their helix configuration are calculated using conventional TST with B3LYP/6-311G(d,p) and M06-2X/6-311G(d,p) methods at 1000 K. Our study shows that the TSs of the combustion reaction of PODME (n = 2–6) molecules show similar structural characteristics at the corresponding H-abstraction position. The B3LYP functional has lower activation barrier than M06-2X. For almost all H-abstraction pathways of PODME (n = 2–6), the activation energy decreases as the H-abstraction position moves toward the middle of the molecular chain. H-abstraction occurs preferentially on the hydrogen of the methylene closed center rather than on the closed-end position. Hydrogen bonds are formed between the O atom of O2 and H atom of methylene or methyl hydrogen. These hydrogen bonds could stabilize the TS structures. The linear relationship between the rate constant of the autoignition process of pure compound and the experimental CN data is established. Our calculated rate constants are in good agreement with experimental CN data. It is indicated that the calculation is reliable, and a relationship between kinetics and CN is established. The methodology presented here can be used to generate theoretical kinetic data for other fuel additives.

Computational Methods

Large varieties of computational methods indicate that the HF methods have been known to overestimate barrier heights, and pure DFT methods tend to underestimate the energy barrier.[22,26,29−31] Hybrid DFT methods perform better for the computation of activation barriers than the pure DFT because of variation of the exact HF exchange energy term.[22,32,33] In particular, B3LYP is one of the most popular density functionals for a reasonable computational cost,[34,35] and the M06-2X functional has been shown to predict TS structures with small mean errors.[22] In this work, kinetic parameters of H-abstraction reactions will be determined using the hybrid DFT methods B3LYP and M06-2X. The present study provides new insights into the autoignition process of PODME (n = 2–6).

All the electronic structure calculations were performed using the Gaussian 16 software package.[36] Calculation results of the geometrical optimizations and harmonic vibrational frequency for reactants, TSs, and products were performed with the hybrid DFT of B3LYP[37−39] and M06-2X[40] coupled with the 6-311G(d,p) basis set. All TSs were obtained by using the routine TS (Berny) method.[41,42] On the basis of the TST,[4,6,26] the bimolecular rate constants of all reactions were calculated by using parameters and energy barriers obtained by the B3LYP and M06-2X methods without being scaled by the factor anymore. All the TS structures were confirmed by only one imaginary frequency, and intrinsic reaction coordinate calculations were also carried out to verify that the TS connects the corresponding reactants and products. The molecular oxygen O2 participated in the reaction in the form of ground triplet state, and the HO2 radical and R radical were considered as doublet states. The ZPE and thermal correction values were used for the evaluation of the reaction activation energies at 1000 K, which represented the operating temperature of diesel engine during the combustion process.[4−6,16] The kinetic model for combustion reactions of PODME (n = 2–6) based on the conventional TST was used to calculate the rate constants at 1000 K.