Molecular Dynamics Simulation on the Pyrolysis Process of PODE3-5

Abstract

This paper investigates the pyrolysis of PODEn (n = 3, 4, 5) using ReaxFF molecular dynamics simulation. A large-scale model, which contains 2000 PODEn molecules, is simulated at 3000 K. The higher frequencies of the initial PODEn decomposition reaction at α or β C-O bond show that the α or β C-O bond in PODEn is not easy to break, which is approximately half the number of the other type of C-O bond dissociation. Furthermore, the bond dissociation energies (BDEs) are calculated using the ReaxFF method. The BDE of α or β C-O bond is higher than that of the other C-O bond, ~3–11 kcal/mol, indicating that BDE is one of the factors causing the different proportions of bonds broken. The evolution of pyrolysis products is also investigated. The results reveal that long-chain pyrolysis products from the initial PODEn decomposition are prone to further reaction, while a large amount of CH3O and CH3 remains in the system. This helps explain the difficulty in α and β C-O bond dissociation reactions. The results of the pyrolysis products are consistent with the result in further chemical kinetic simulation. The C2 species in pyrolysis products is relatively low, especially for C2H4 and C2H3, which is around zero. This supports the ability of PODEn to reduce soot emission.

1. Introduction

Diesel engines have been widely used in vehicles and engineering machinery due to the advantages in thermal efficiency and output torque [1]. However, soot emission due to the incomplete combustion of conventional hydrocarbon fuels, including diesel and gasoline, is still a serious problem which degrades the air quality [2]. Such problems have attracted more and more attention from governments all around the world, who have introduced a series of strict environmental protection regulations [3]. Reducing engine soot emissions motivates the development of fuel additives or alternatives. Among the fuel additives, oxygenated fuels are proven to reduce soot emissions, due to their extra oxygen content [4,5,6].

As one type of oxygenated fuel, Poly(Oxymethylene) Dimethyl Ethers (CH₃O(CH₂O)_nCH₃, n ≥1), referred to as PODEn, are promising additives that could be used in diesel engines. PODEn have high oxygen content that is proven to reduce the emissions in the compression ignition (CI) engine [7,8]. In addition, PODEn have a high cetane number and similar physical-chemical properties to those of diesel [9], which enables PODEn to be a suitable additive. Moreover, the production of PODEn could be from coal or biomass, which benefits the commercial application of PODEn [10].

PODEn have been studied by many researchers in experimental and numerical studies. Lumpp et al. [11] reported that PODEn could dramatically reduce the mass concentration of the particulate matter as an additive. Liu et al. [12] conducted a series of diesel engine experiments with diesel-PODE3 blends (10, 20, 30 vol%) and proved PODEn to be able to alter the combustion process in the engine. Our previous study also showed that PODE3 can reduce both the number and mass concentration of the soot particles due to the extra oxygen content in PODEn [13]. Moreover, Sun et al. [14], and He et al. [15] have developed a detailed mechanism for PODE3 as an analogy to the DMM (PODE1) mechanism. Furthermore, Cai et al. [16] developed a detailed mechanism for PODE2-4 with a class-based automatic reaction alternator (CLARA) which is also based on the DMM mechanism. ReaxFF molecular dynamics simulation is also widely used to investigate the chemical reaction system. It has shown its capabilities to unravel the reaction pathways and networks of hydrocarbon fuels [17,18]. In addition, the computational cost is not expensive as compared to the quantum mechanics calculation [19,20]. Xin, Liu, Liu, Huo, Li, Wang and Cheng [18] studied the decomposition of a series of hydrocarbons by ReaxFF-based MD simulation. The decomposition of PODE3 and its effect as an additive in diesel engines has also recently been studied using ReaxFF MD simulation on PODE3 [21]. The decomposition study on PODE3 was simulated with a 20 PODE3 system.

According to Zheng, Tang, Wang, Liao and Wang [9], only PODE3, PODE4, and PODE5 are suitable fuel additives that can be used in diesel engines. Therefore, in the following section, PODEn only refer to PODE3, PODE4, and PODE5. In this study, the PODEn (n = 3, 4, 5) were studied based on ReaxFF MD simulation to understand the thermal decomposition process of PODEn. The system was a large-scale MD simulation, which contained 2000 PODEn molecules. To the knowledge of the author, it is the first large-scale MD simulation on PODE3-5. Specifically, we attached great importance to the initial C-O bond dissociations reactions in PODEn. The bond dissociation energies (BDE) of C-O bonds in PODEn were also calculated using the ReaxFF method. The trajectory of PODEn pyrolysis products was investigated to understand the decomposition mechanism of PODEn. Also, the effect of chain length of PODEn on the decomposition process was carefully studied by comparing the simulation of the PODEn. The simulation results can be a good reference for further developing the detailed chemical kinetic mechanism of PODE3-5 pyrolysis.

2. Computational Procedures

2.1. Density Functional Theory (DFT) Calculation

In this study, the initial geometry of PODE3-5, and the pyrolysis products, including CH3 and CH3O, radicals were optimized at the M06-2X/6-311++G(d,p) level of theory [22,23,24]. The internal rotations were scanned in 10-degree increments as a function of the dihedral angle of every C-O bond in the above-mentioned species to make sure that all the geometry-optimized species were in their global lowest point of energy. All of the quantum chemical calculations were conducted using Gaussian 16 [25].

2.2. ReaxFF Molecular Dynamic (MD) Simulation

The molecular dynamic simulation in this study used the ReaxFF force field to simulate the reaction process in the PODEn systems, which adopts bond order to correlate bond distance and bond energy. The ReaxFF force field is designed for a large-scale reactive chemical system in MD simulation which was first developed by Van Duin, Dasgupta, Lorant and Goddard [20]. For this force field, the non-bond interactions between atoms were identified by Coulomb and van der Walls potentials and the atomic charges were calculated by the electronegativity equilibrium method (EEM) [26,27]. The Equation (1) describes the system energy E_system:

$$E_{system}=E_{bond}+E_{over}+E_{under}+E_{lp}+E_{val}+E_{tor}+E_{vdW}+E_{Coulumb}$$

(1)

where E_bond, E_over, E_under, E_lp, E_val, E_tor, E_vdW, E_coulumb represents bond energy, energy penalty terms for over-coordination and under-coordination, lone pair electron energy, valence energy, tortional energy, van der Waals energy, and Coulumb energy, respectively. In this study, the ReaxFF parameters of CHO-2016 [28], as an extended version of CHO-2008 [19], were applied. Compared with CHO-2008, CHO-2016 is able to accurately describe the C1 chemistry and oxidation initiation reaction, as well as the pyrolysis and oxidation of the species, which contains C, H, and O atoms.

This study focused on the initial pyrolysis of PODE3-5. In a typical system, 2000 PODEn (n = 3, 4, 5) molecules are randomly distributed in a periodic cubic box with a side length of 500 Å. The monomer geometry was optimized using the DFT method. The randomization of the molecules was implemented by PACKMOL [29]. The initial systems were equilibrated at 300 K for 20 ps with the Nose-Hoover thermostat (NVT) to avoid atom overlapping and achieve energy optimization. After that, the configurations obtained were used to further study the kinetics of PODEn. In Figure 1, the optimized structures of PODE3-5 are presented as an example, and the initial configuration of the PODE3 reaction system is shown via OVITO [30].

With the equilibrated configurations, all the PODEn systems were heated from 300 K to 3000 K. The heating process lasted for 2 ps and was controlled by a NVT thermostat. After that, the reaction systems were controlled at 3000 K for 200 ps to evaluate the pyrolysis of PODEn. Considering the expensive computational cost of a large-scale MD simulation, the time scale was ps, which is much smaller than the real chemical reaction time scale, especially for the reaction under engine conditions. Therefore, the temperature of the reaction system was increased to 3000 K, which is much higher than the temperature range in CI engine, to accelerate the reaction rate of PODEn decomposition. Such a strategy has been validated by many researchers [18,31,32], and the results show that 3000 K is suitable for PODEn to completely decompose in 200 ps within an acceptable computational time.

3. Results and Discussions

3.1. Initial C-O Bond Dissociation of PODEn

Since there are no C-C bonds in PODEn species, the initial decomposition reactions of PODEn were accompanied by the bond dissociations of the C-O bonds. Also, as shown in Figure 1, PODE3-5 are all self-symmetric molecules, and the external symmetries are all 2. Therefore, for example, there are four different C-O bonds in PODE3, which are labeled as α, β, γ, and δ C-O bonds in sequence from terminal to the middle. For PODE4 and PODE5, there are five and six different C-O bonds, respectively, and these are labeled in a similar way.

The dissociations of the C-O bonds in PODEn were carefully studied. Because of the different position of C-O bonds, the dissociation of the bonds led to different decomposition products of PODEn. In the MD simulation, the different decomposition reactions could be post-processed using ReacNetGenerator [33], and thus it was possible to identify the reaction pathway of the initial decomposition reactions of PODEn to reveal the tendency of C-O bonds dissociation at different positions. Figure 2 shows the four possible C-O bond dissociations in PODE3. As shown, the proportion of α and β C-O bond dissociations are close, which are 17.36%, and 16.55%, respectively, while the proportion of γ, δ C-O bond dissociations are closely dominated in 34.94%, and 31.15%, respectively. This indicates that the C-O bond dissociations mainly occur at a position close to the central C atom.

The decomposition reactions of PODE4 and PODE5 are depicted in Figure 3 and Figure 4, respectively. Although the chain length of PODEn varied, the C-O bond dissociations in PODEn followed a similar pattern, that is, the proportions of α and β C-O bond dissociations were close, while the proportions of the other C-O bond dissociations were close. The proportion ratio of α or β C-O bond dissociation and other kind of C-O bond dissociation was also around 1:2, showing that the C-O bonds that were near the middle point of the PODEn molecule break more frequently in the MD simulation. Furthermore, as the chain length of PODEn increased, the types of C-O bonds increased, which is n + 1 (‘n’ refers to the n value in PODEn), and the proportion of α and β C-O bond dissociations among all possible C-O bond dissociations actually dropped following the dependence of ~1/n. Therefore, at the beginning of PODEn decomposition, most of the C-O bonds that were close to the central C atom are broken and it was expected that a larger amount of the corresponding pyrolysis products would be produced.

3.2. Bond Dissociation Energies of PODEn

The bond dissociation energies (BDEs) of the C-O bonds in PODEn were then calculated using the ReaxFF method to better understand the decomposition of PODEn. The results of BDEs of all aforementioned C-O bonds in PODEn are shown in

Table 1. Bond dissociation energies of C-O bonds in PODEn calculated using ReaxFF.

C-O Bond	ReaxFF Method (kcal/mol)
PODE3 α	71.71
PODE3 β	67.96
PODE3 γ	64.61
PODE3 δ	62.67
PODE4 α	71.79
PODE4 β	72.82
PODE4 γ	64.73
PODE4 δ	63.03
PODE4 ε	62.53
PODE5 α	71.94
PODE5 β	73.87
PODE5 γ	64.46
PODE5 δ	67.54
PODE5 ε	62.30
PODE5 ζ	62.54

. As shown, the BDEs of α and β C-O bonds were close, which were all around 70 kcal/mol, while the BDEs of the other kinds of C-O bonds were close and all less than 65 kcal/mol, which were smaller than that of α and β C-O bonds. This is consistent with the results in Section 3.1, indicating the intramolecular BDE is an important factor that affects bond dissociation.

The result of BDEs calculated using the ReaxFF method favored the result in Section 3.2, which is that α and β C-O bond dissociations are less frequently to happen than the C-O bond dissociations of other positions. Nevertheless, among all the BDEs of C-O bonds of PODEn calculated using the ReaxFF method, the difference is not significant. However, the result in Section 3.1 reveals that the proportion of α or β C-O bonds broken was only about half the proportion of other kinds of C-O bonds. Therefore, the difference in BDEs of the C-O bonds may not be the only factor affecting the C-O bond dissociations. To understand the bond dissociation process of PODEn, it is necessary to analyse the time evolution of PODEn and pyrolysis products from the initial decomposition of PODEn.

3.3. Pyrolysis of PODEn and Its Initial Decomposition Products

In this section, PODEn and their initial decomposition products, in the molecule form of CxHyOz, are traced with the evolution of time (0–202 ps, including the heating process). To avoid the ambiguity of the molecular formula, the identical SMILES forms of CxHyOz species are provided in Table S1 in Supplementary Material. Therefore, in this paper, the CxHyOz molecular formula only indicates the species in Table S1.

Figure 5 presents the time evolution of PODE3 and its pyrolysis products during the MD simulation. During the simulation period, PODE3 was completely decomposed (99.5%). It is shown that the main pyrolysis products were CH3, CH3O, and C2H5O. Among them, at the beginning of the simulation, the production rate was CH3O > C2H5O > CH3. However, as the simulation progressed, the amount of CH3O and C2H5O started to drop after they reach their peaks. This can be attributed to their further reaction to other species through β-scission or H-abstraction. At the same time, CH3 continued to be produced during the whole simulation time and became the species of the largest proportion at the end of the simulation. As for C2H5O2 and C3H7O2, their species profiles are similar, and both dropped after they reached their maximum amount. However, the maximum amounts were still much lower than that of the initial PODE3 amount, around 5% and 6%, respectively, and these amounts were around zero at the end. For C3H7O3, C4H9O3, and C4H9O4, their amounts remained at a very low level throughout the simulation, which was around zero. Even the maximum amount among them was just 0.11% of the initial PODE3 amount.

The time evolution of PODE4 and PODE5 and their pyrolysis products during the MD simulation are shown in Figure 6 and Figure 7. Generally, the species profiles of PODE4 and PODE5 have a similar pattern to PODE3, thus the discussions on PODE4 and PODE5 decomposition will be more concise to avoid repetition. The top three pyrolysis products were still CH3, CH3O, and CH2O during both PODE4 and PODE5 decomposition, and they have the same mode of evolution as compared to PODE3. For the production process of larger radicals, they follow the basic pattern that CxHyOz species first evolve to a peak point of the amount in the first quarter of the simulation then drop to around zero at the end of the simulation.

The results in this section are consistent with the results in Section 3.1. Among all the products of the initial decomposition of PODEn, only CH3 and CH3O still remained in the reaction system in large amounts at the end of the simulation, which prevented the further α and β C-O bonds from breaking, the products of which contain CH3 and CH3O, respectively. While for the other PODEn pyrolysis products, their concentrations were kept at a low level, or even close to zero, throughout the simulation process, which supports the accordingly bond dissociation reactions to happen. Therefore, the dominant factor in bond dissociation of PODEn at high temperature was no longer the bond dissociation energy, but the further reaction of the initial pyrolysis product.

3.4. C0-C2 Species in the Pyrolysis of PODEn

The evolution of C0-C2 species (CH2, CH3O, CH2O, CH4, C2H2, C2H3, C2H4, C2H5, C2H6, CO, CO2, OH, H, H2, HO2) during the PODEn pyrolysis were studied in this research. On the one hand, according to the study by Sun, Wang, Li, Zhang, Yang, Yang, Li, Westbrook and Law [14], the main pyrolysis products of PODE3 are CH3, CH3O, and CH2O through C-O bond dissociation. The number concentrations of these three products are well studied in this part. On the other hand, according to Tan et al. [34], the formation of the first benzene ring mainly follows the C3 pathway and C2 + C4 pathway and PODE3 rarely provides C3 or C4 species. Therefore C2 species, as well as CH3, have been well studied to understand the reason why PODEn could reduce soot formation in combustion. Moreover, to be comprehensive, other C0-C1 species were also included to evaluate the PODEn pyrolysis.

The time evolution of C0-C2 species during PODE3 pyrolysis is shown in Figure 8. The top three pyrolysis products at the end of the simulation were CH2O, CH3, and CH3O. The final amount of CH2O was slightly over 6000, which is about three times the initial reactant amount. The amount of CH3 grew monotonically throughout the MD simulation and reached about 87% of the initial PODE3 amount. The ratio of the amount of CH3O over the initial PODE amount, however, grew to over 70% at the beginning, then fell back to about 50% at the end of the simulation. The PODEn decomposition generally follows Equation (2) [35]. This formula indicates that after one PODEn decomposes, the final products should consist of one methoxy radical (CH3O), one methyl radical (CH3), and n formaldehyde (CH2O), which is related to the chain length of PODEn. The result of CH2O production in pyrolysis was consistent with this formula. However, the final amount of CH3 was slighter less than the expected amount, while the final amount of CH3O was significantly underestimated according to the formula. Results of MD simulation such as these indicate that CH2O is the most stable, and CH3 is the second with a small amount of consumption, while CH3O undergoes rapid consumption after being formed.

Apart from the products from C-O bond dissociation, other products are also presented in Figure 8. However, at the beginning of PODEn pyrolysis, the amount of these products was much lower than the amount of C1 species. Among them, H, H2, and OH were the top three products in sequence. The high content of H revealed that the pyrolysis of PODE3 was accompanied, though not dominated, by C-H dissociation. The existence of H2 indicated the H-abstraction in progress and OH was necessary for further ignitability. During the simulation, 3.75% CO and 2.4% CH4 (divided by the amount of initial PODE3) also formed, which could be converted by CH2O and CH3, respectively.

For C2 species, the amounts of C2H6 and C2H5 were close, and each one was less than 2% of that of CH3, which is the source of the formation of C2H6 and C2H5. Furthermore, the amount of C2H4 was even lower, while C2H3 and C2H2 are rarely formed in the simulation. This indicates that, during the pyrolysis of PODE3, a small amount of C2 species was produced, which makes it difficult to further form C2H2. That could explain why PODE3 could reduce soot formation in the combustion process.

$$\mathrm{PODEn}=\mathrm{C}\mathrm{H}3\mathrm{O}•+\mathrm{C}\mathrm{H}3•+\mathrm{n}\mathrm{C}\mathrm{H}2\mathrm{O}$$

(2)

Time evolution of C0-C2 species during PODE4 and PODE5 pyrolysis are shown in Figure 9 and Figure 10. When the chain length increased, the patterns of the number of C0-C2 species were still similar, however, the distributions of different products varied. The amount of CH2O still followed Equation (2), which is about four and five times the initial PODEn amount. The evolution of CH3O and CH3 during PODE4 and PODE5 pyrolysis had the same trend as that during PODE3 pyrolysis. Nonetheless, the production rates of these three pyrolysis products increased when PODEn had a longer chain length. Such results reveal that increasing the chain length of PODEn cannot increase the total amount of CH3 and CH3O production but increases the CH2O. However, it accelerates the pyrolysis process. The evidence of the speed-up of the PODEn decomposition is also shown in the distribution of the other pyrolysis products. The total amounts of H, H2, OH, and other products also increased as the chain length of PODEn grew, though, the amounts of C2 species in PODE4 and PODE5 pyrolysis products were still very low. The productions of C2H2 and C2H3 were still neglectable. This indicates that PODE4 and PODE5 are both expected to reduce soot formation in combustion as an analogy to PODE3.

3.5. Chemical Kinetic Study on PODE3 Pyrolysis

In this section, a chemical kinetic study was conducted on PODE3 in a 0-D reactor with CHEMKIN-PRO [36] as a comparison with MD simulation. The detailed mechanism of PODE3 developed by He, Wang, You, Liu, Wang, Li and He [15] was used. The reactor in this simulation was controlled at a constant temperature (3000 K) and volume to keep the same condition as that of MD simulation. The initial pressure of the reactor also remained the same as MD simulation, which is ~33 bar. Correspondingly, the presented time scale of this simulation was also 0–200 ps. The initial pyrolysis products of PODE3 are compared and shown in Figure 11. During the reaction time scale, four observable species were CH2O, CH3, CH3O and C2H5O, while the other pyrolysis products decomposed rapidly after their generation. Furthermore, a large amount of CH2O was produced, which was much higher than the amount of CH3, CH3O and C2H5O. Despite the production rate of CH3O and C2H5O being higher than that of CH3 at the beginning, both CH3O and C2H5O were consumed rapidly. It can be concluded that these consumption rate trends are consistent with the MD simulation. However, the absolute values of the consumption rates of these pyrolysis products were underestimated in the MD simulation. This indicates that although MD simulation can recognize the tendency of the decomposition reactions of PODE3, the reaction rates of the reactions still need further calculation, e.g., quantum-chemical calculation.

4. Conclusions

In this study, we investigated the pyrolysis processes of PODEn (n = 3, 4, 5) by using ReaxFF molecular dynamics (MD) simulation. The MD simulation model is on a large scale which contains 2000 PODEn molecules and remains at 3000 K.

The frequencies of the initial PODEn decomposition reaction were carefully recorded. The results reveal that the proportion of the α or β C-O bond dissociation was only about half the proportion of the other kind of C-O bond dissociation, which indicates that it was relatively difficult to break the α and β C-O bond. This phenomenon was attributed to the different bond dissociation energies (BDEs) of different C-O bonds in PODEn. The BDEs of all C-O bonds in PODEn were therefore calculated using the ReaxFF method. The relative difference in the BDEs of C-O bonds calculated by the ReaxFF method was not evident as compared to the difference in the frequencies of C-O bonds dissociations. Therefore, we propose that the difference in BDEs of C-O bonds is not the only factor affecting the bond dissociation reactions of PODEn at high temperature.

The time evolution of the pyrolysis products of PODEn was also presented in this study. The evolution of the initial pyrolysis products due to the C-O bond dissociation of PODEn shows that apart from CH3 and CH3O, other CxHyOz species quickly decompose further into other smaller species after their generation. The remaining CH3 and CH3O in the reaction system prevented the α and β C-O bond from breaking further, which leads to the low frequencies of α and β C-O bond dissociation reactions. The abundance of H and H2 also suggests that ongoing C-H breaking and H-abstraction occurs, though not dominated. In addition, the results show that PODEn decomposition produces a large amount of OH, which could promote ignition in combustion. As compared with CH3 and CH3O, the number concentration of C2 species is very low, especially for C2H4 and C2H3, which were also close to zero throughout the simulation. This supports the point that PODEn can reduce soot particle emission when it is blended with other conventional fuel during combustion.

The PODE3 MD simulation was also compared with a chemical kinetic simulation. The results show that MD simulation can capture the pyrolysis products of PODE3 correctly, though the reaction rates still need further calculation. This indicates that MD simulation could help identify the key reactions of the PODEn reaction system and help develop or optimize a detailed chemical kinetic mechanism.