Article pubs.acs.org/EF
Theoretical Modeling of Thermal Decomposition of Methylnaphthalene Derivatives: Influence of Substituents Robin Chaudret,* Andreas Bick, and Xenophon Krokidis Scienomics, 16 rue de l’Arcade, 75008 Paris, France ABSTRACT: The kinetics and thermodynamics of thermal decomposition of naphthalene and several mono- and dimethylnaphathalene derivatives have been established through reactive molecular dynamics simulations using ReaxFF. These results were compared to previous theoretical and experimental studies and also compared to density functional theory calculation results. This work demonstrates that the kinetics and thermodynamics of the initial activation reaction is directly impacted by the position and number of methyl substituents. Subsequently, the activated molecules react to form either small organic or large char molecules. The char formation mechanism is shown to occur in three steps: activation, dimerization/ trimerization, and condensation. Finally, temperature effects on the char formation reaction were also studied.
■
INTRODUCTION Polycyclic aromatic hydrocarbons (PAHs) are organic compounds formed of multiple aromatic rings bound or fused together, where electron delocalization is very high. Such compounds can be found in fossil fuels, such as crude oils, but can also be formed during high-temperature reactions. Understanding the thermal stability of PAH is a central question in several industrial processes, such as pyrolysis or high-temperature reactivity1−3 or in the reactivity of crude oils in natural reservoirs.4 More particularly, understanding the formation of char or soot compounds during these processes is of high importance and has been widely studied through pyrolysis experiments of PAH compounds.5−10 Naphthalene derivatives are small PAH compounds that are often used as models of more complex PAH molecules. It has been shown experimentally that the position and number of methyl substituents have an influence on the properties of naphthalene derivatives.5,6 It is therefore important to understand, from a molecular point of view, the effect of the substituents on the stability of the PAH molecule and how that influences its decomposition thermal dependencies. ReaxFF reactive force field approach, developed by van Duin et al.,11, enables the modeling of bond formation and breaking using bond-order-dependent parameters initially introduced by Tersoff.12 A longer presentation of the different force field terms can be found in ref 1. ReaxFF simulations have been applied to many types of systems and reactions, such as reactivity of high-energy materials13,14 or metal and metal oxide reactivity.15−20 One of the main application areas of ReaxFF remains however the thermal stability, reactivity, or combustion of fuels, PAH, or organic molecules.21−26 In a previous study from Leininger et al.,27 the influence of the temperature on the properties of 1-methylnaphthalene was inspected using ReaxFF. The goal of the present study is to understand the influence of an alkyl side chain on the thermal stability of PAH. To obtain a systematic understanding, it was decided to focus on the influence of the position and number of methyl substituents on the thermal stability of naphthalene derivatives, which are good models for more complex PAH molecules. The © XXXX American Chemical Society
naphthalene derivative activation and char formation mechanisms were more specifically investigated.
■
METHOD
Notations. Table 1 is listing the different notations used in the paper.
Table 1. List of Notations Used in the Paper notation Naph 1-MNaph 2-MNaph 1,2-MNaph 1,3-MNaph 1,4-MNaph MNaphD MNaphD-nH C20+ C@C20+ Cmax C@Cmax
meaning naphthalene 1-methylnaphthalene 2-methylnaphthalene 1,2-dimethylnaphthalene 1,3-dimethylnaphthalene 1,4-dimethylnaphthalene any of the methylnaphthalene derivatives any of the methylnaphthalene derivatives from which n hydrogen atoms have been abstracted (n = 1 or 2 in our study) molecules containing 20 carbon atoms or more number of C atoms belonging to C20+ molecules heaviest molecule in the simulation box at a given time number of carbon atoms belonging to the Cmax molecule
The numbering of the different naphthalene derivatives corresponds to the position of the methyl radicals, as shown in Figure 1. Figure 1 also shows that, with one substituent only, positions 1 and 4 as well as positions 2 and 3 are equivalent. Force Field Simulation Setup. All systems were set up following the procedure from Leininger et al.27 For every MNaphD molecule considered in this study, amorphous systems composed of 20 molecules at a density of 0.809 g cm−3 were built using the amorphous builder of Scienomics’ MAPS platform28 (see Figure 2). Note that this density follows ref 27 and also corresponds to the typical density of a crude oil in its reservoir. Additionally, it was shown in ref 27 that system density influences the overall reaction kinetics but not the kinetic constant or reaction order. Received: May 4, 2016 Revised: July 16, 2016
A
DOI: 10.1021/acs.energyfuels.6b01062 Energy Fuels XXXX, XXX, XXX−XXX
Article
Energy & Fuels
Geometry optimizations were performed using unrestricted B3LYP33,34 DFT together with Def2-SV(P)35 basis sets and the resolution of identity approximation using NWChem36 software. The following convergence parameters were achieved: (i) convergence of the maximum gradient to 4.5 × 10−4 Ha/Bohr, (ii) convergence of the root-mean-square (RMS) gradient to 3.0 × 10−4 Ha/Bohr, (iii) convergence of the maximum displacement to 1.8 × 10−3 Bohr, and (iv) convergence of the RMS displacement to 1.2 × 10−3 Bohr. These energies were corrected from the zero point energy. The dissociation energy, defined as EDiss = E(MNaphD) − E(MNaphD-H) − E(H), where E(MNaphD), E(MNaphD-H), and E(H) denote the energies of MNaphD, MNaphD-H, and the hydrogen atom, was computed.
■
Figure 1. Representation of the (a) Naph with the atom numbering of (b) 1-MNaph, (c) 2-MNaph, (d) 1,2-MNaph, (e) 1,3-MNaph, and (f) 1,4-MNaph molecules.
RESULTS AND DISCUSSION The initial step of the decomposition of naphthalene derivatives involves the breaking of the C−H bond into a hydrogen atom and the “activated naphthalene derivative” (MNaphD-H). The system then reacts further to form either larger char molecules or light hydrocarbons (mainly C2H2 and C2H4) and H2 (see Figure 3 for evolution of the different species from 1-MNaph at
Figure 2. Initial amorphous configuration of the 1-MNaph system. Periodic boundary conditions were used throughout this study. The molecular systems were initially optimized using 5000 steps of conjugated gradient geometry optimization after typing with TeamFF force field.29 Subsequently, ReaxFF force field parameters (CHO combustion force field)30 were assigned to the system for molecular dynamics simulation. The system was run for 100 ps with a time step of 0.1 fs in a NVT ensemble, and the particle mesh31 approach was used to treat electrostatic interactions. This simulation protocol was performed at four different temperatures: 3000, 3500, 3800, and 4000 K. All geometry optimizations and molecular dynamics simulations were performed using LAMMPS software.32 These very high temperatures were used to accelerate the processes that usually occur at geologic time at lower temperatures. Earlier studies27 showed that increasing the ReaxFF simulation temperature would mainly have an influence on the light compounds formed during the naphthalene derivative decomposition and more specifically on the amount of H2, while neither the decomposition of the derivative itself nor the heavy product formations were found to be impacted. MAPS ReaxFF analysis tool was used to analyze the output of these ReaxFF simulations and to derive the evolution of molecule species during the simulation. To analyze the amount of char formed in the system, the evolution of the number of carbon atoms belonging to C20+ molecules (C@C20+) was computed together with the number of carbon atoms in the heaviest molecule at a given time (C@Cmax). The comparison between C@Cmax and C@C20+ can give an indication of the degree of condensation of the MNaphD compound. Density Functional Theory (DFT) Calculation Setup. A DFTbased analysis of the hydrogen dissociation energies of MNaphD molecules was conducted. The hydrogen atom to remove was carefully selected. For naphthalene, hydrogen was removed from either position 1 or 2 of the naphthalene ring, while for 1-MNaph and 2-MNaph, hydrogen was removed from the methyl substituent. In the case of dimethylated naphthalene derivatives (1,2-MNaph, 1,3-MNaph, and 1,4-MNaph), hydrogen was removed from either of the methyl substituents.
Figure 3. Evolution of the number of different molecules and number of C atoms belonging to a C20+ molecule along the reactive MD simulation of 1-MNaph at 3500 K.
3500 K). These results are in good agreement with several experimental studies, all of which found this type of product as a result of either 1-MNaph or 2-MNaph pyrolysis.5−7,10 We will now focus on the initial step of the reaction and the char formation mechanism. Initiation Reaction. Several mechanisms for the initiation reaction exist; however, two main mechanisms are apparent. (1) The dissociation of hydrogen from the naphthalene derivative, where one C−H bond from one derivative breaks and induces formation of an activated radical compound and a hydrogen atom: NaphD = NaphD-H + H. (2) The reversed radical disproportionation mutation, where two derivatives react to form two radical species, one where a hydrogen has been abstracted (noted NaphD-H) and another one where a hydrogen has been added. The latter species is so reactive that it almost immediately reacts to form mainly NaphD or NaphDH compounds. Therefore, in the following, only NaphD-H intermediate species will be discussed. These two mechanisms have been experimentally shown to be predominant in 1-MNaph pyrolysis.7 On the Naph compound, hydrogen is abstracted from the naphthalene ring, whereas, when at least one substituent is present, hydrogen is abstracted from the methyl group and the formed radical is stabilized through mesomeric effects. This B
DOI: 10.1021/acs.energyfuels.6b01062 Energy Fuels XXXX, XXX, XXX−XXX
Article
Energy & Fuels
Figure 4. Evolution of the number of MNaphD: (a) sum of all different forms of MNaphD, (b) MNaphD-H, (c) MNaphD-2H, and (d) (MNaphD + MNaphD-H + MNaphD-2H) along the reactive MD simulation at T = 3000 K.
Table 2. Dissociation and Relative Dissociation Energies for the Different MNaphD Compounds in kJ/mola Naph EDiss ΔEDiss
558.7 570.0 0.0 −11.3
1-MNaph
2-MNaph
439.0
434.0
−131.0
−136.0
1,2-MNaph
1,3-MNaph
416.6 423.6 −153.4 −146.4
405.2 414.2 −164.8 −155.8
1,4-MNaph 431.1 −138.9
a
When several energies are proposed, they correspond to the abstraction of hydrogen in different positions: for Naph, the abstraction of hydrogen from the ring in positions 1 and 2, and for 1,2-MNaph (1,3-MNaph), the abstraction of hydrogen from the methyl substituent in positions 1 and 2 (1 and 3). Note that, for 1,4-MNaph, the two abstractions are equivalent by symmetry. ΔEDiss corresponds to the relative dissociation energy computed as follow: ΔEDiss(X) = EDiss(X) − EDiss (Naph in position 1).
Both the number and position of methyl substitutions were shown to have an impact on the reactivity of methylnaphthalene derivatives, because they impact the degradation kinetics and reaction path. To further understand the formation and stability of the activated radical species, it was decided to compare their energy to their related naphthalene derivatives. To this end, DFT calculations of the hydrogen dissociation energy (EDiss) (Table 2) were performed. As expected, the Naph compound has a much higher dissociation energy (between 130 and 170 kJ/mol) than any other MNaphD compound. Additionally, all disubstituted compounds have a lower dissociation energy than monosubstituted compounds. Finally, it is interesting to notice that the abstraction of hydrogen from molecules containing a methyl substituent in position 1 (1-MNaph and 1,4-MNaph) has a slightly higher dissociation energy compared to the case of (at least) one substituent in position 2 or 3 (2-MNaph, 1,2MNaph, and 1,3-MNaph). To summarize, the ab initio results corroborate very well with the findings from ReaxFF simulations. These results are in good agreement with the results from Leininger et al.,27 who also found that, in the temperature range of 3000−4000 K, the predominant initiation reaction is C−H bond breaking. They also corroborate previous experimental studies, which showed Naph to be much more stable than 1MNaph or 2-MNaph with respect to the temperature, whereas the latter two exposed a similar behavior.5,6 These experiments
induces strong differences between Naph and MNaphD compounds, as shown in Figure 4. Additionally, further differences can found between mono- and dimethylated naphthalene. Indeed, with one substituent (1-MNaph and 2MNaph), 50% or less MNaphD is converted, whereas with two substituents (1,2-MNaph, 1,3-MNaph, and 1,4-MNaph) at least 75% MNaphD is converted. Similarly, the substituent appears to have an effect on the kinetics of the degradation of MNaphD. Indeed, for molecules containing one substituent in the position 2 or 3 (2-MNaph, 1,2-MNaph, and 1,3MNaph), MNaphD started being decomposed after less than 5 ps, whereas for the other molecules, MNaphD is only becoming decomposed after more than 10 ps (11.2 ps for 1,4-MNaph and 39.4 ps for 1-MNaph). The sum of all of the different forms of MNaphD (MNaphD + MNaphD-H + MNaphD-2H) was computed to estimate the global amount of MNaphD species (activated or not) remaining in the system at a given time step and, thereby, indicate how fast MNaphD reacts once it is activated. Figure 4b shows that MNaphD reacts almost immediately if it is substituted in position 2 (2-MNaph, 1,2-MNaph, and 1,3MNaph equivalent to 2,4-MNaph); it is degraded after less than 5 ps. On the contrary, MNaphD and its activated species remain quite stable if there is no substitution in position 2 (Naph, 1-MNaph, and 1,4-MNaph). Indeed, in such cases, even after 35 ps, none of the MNaphD species has reacted. The substitution in position 2 seems therefore to increase the kinetics of the degradation reaction and disadvantage the back reactions regenerating MNaphD. C
DOI: 10.1021/acs.energyfuels.6b01062 Energy Fuels XXXX, XXX, XXX−XXX
Article
Energy & Fuels
Figure 5. Number of carbon atoms in C20+ molecules (C@C20+, blue lines) and number of carbon atoms in the heaviest molecule (C@Cmax, orange triangle) as a function of the simulation step (ps) for (a) Naph, (b) 1-MNaph, (c) 2-MNaph, (d) 1,2-MNaph, (e) 1,3-MNaph, and (f) 1,4-MNaph compounds at 3500 K.
(dimers or trimers). During that step, the number of C@C20+ increases, while the number of C@Cmax remains small (less than 50). Finally, in the aggregation step, the low-weight chars react with each other and create some high-molecular-weight molecules. During this step, C@Cmax increases rapidly and reaches almost the same value as C@C20+. If the activation step is really fast, the dimerization/trimerization step is slow and depends upon the substituent number and position (Figure 5). The influence of the temperature was also studied for the char formation. If, for T = 3000 K, the char compounds are almost not created after 100 ps, when the temperature increases, the formation kinetics increases accordingly. It is worth noting that the temperature does not greatly influence the time required for the activation step, which corresponds to about 5 ps for the different temperatures (Figure 6). On the contrary, the time for the dimerization/trimerization step is greatly reduced when the temperature increases. This finding means that the formation of heavy char molecules is faster at high temperatures.
also showed that the activation step of MNaphD was usually hydrogen abstraction from the methyl group. Char Formation Mechanism. Another interesting aspect of the naphthalene derivative decomposition reaction is the formation of heavy char molecules. This reaction is crucial in many industrial processes as the char formation is usually an undesired side effect. At 3500 K, the number of carbon atoms in a char compound (molecules containing 20 or more carbon atoms, C@C20+) initially increases and then converges (see Figure 5) toward a value corresponding to up to 70% of the total number of carbon atoms in the system for the different MNaphD. These results are in very good agreement with Leininger et al. simulations;27 i.e., that more than 75% of the carbon atoms were belonging to C20+ molecules for 1-MNaph at 3800 K. It is expected, from the simulations performed, that longer time simulations may increase the final char formation yield, leading to the presence of only heavy char molecules and small hydrocarbons. Several experimental results corroborate these findings because it was found that thermal cracking MNaphD molecules induces the formation of soot and heavy char molecules.5−7 The char formation mechanism has, however, never been clearly identified. In the following, we will use results from the ReaxFF simulation to analyze in more detail the heavy char formation mechanism. The comparison between the number of carbon atoms belonging to the heaviest molecule (C@Cmax) of the system and C@C20+ shows that the formation of char appears in three steps. In the initiation step, the only reaction is the abstraction of hydrogen from the MNaphD compound. During this step, no C20+ molecule is created. In the dimerization/trimerization step, the different activated MNaphD-H molecules are reacting together to form many low-weight char precursor molecules
■
CONCLUSION
Modeling the thermal decomposition of PAHs is crucial in many industrial fields because it has been shown that such compounds were leading to the formation of undesired char molecules. In this study, 100 ps of ReaxFF molecular dynamics simulations were used to study the thermal decomposition of several methyl-substituted naphthalene compounds, with the goal being to model and understand the effect of the position and number of the methyl substituent on the reactivity of the molecule. To achieve this, we focused our analysis on the initial activation and final char formation reactions. D
DOI: 10.1021/acs.energyfuels.6b01062 Energy Fuels XXXX, XXX, XXX−XXX
Article
Energy & Fuels
(6) Yang, J.; Lu, M. Environ. Sci. Technol. 2005, 39, 3077−3082. (7) Leininger, J. P.; Lorant, F.; Minot, C.; Behar, F. Energy Fuels 2006, 20 (6), 2518−2530. (8) Smith, C. M.; Savage, P. E. Energy Fuels 1992, 6, 195−202. (9) Behar, F.; Budzinski, H.; Vandenbroucke, M.; Tang, Y. Energy Fuels 1999, 13, 471−481. (10) Lorant, F.; Behar, F.; Vandenbroucke, M.; McKinney, D. E.; Tang, Y. Energy Fuels 2000, 14, 1143−1155. (11) van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard, W. A., III J. Phys. Chem. A 2001, 105 (31), 9396−9409. (12) Tersoff, J. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 6991. (13) Zhang, L.; Zybin, S. V.; van Duin, A. C. T.; Dasgupta, S.; Goddard, W. A., III; Kober, E. M. J. Phys. Chem. A 2009, 113 (40), 10619−10640. (14) Strachan, A.; van Duin, A. C. T.; Chakraborty, D.; Dasgupta, S.; Goddard, W. A., III Phys. Rev. Lett. 2003, 91 (9), 098301. (15) Hong, S.; van Duin, A. C. T. J. Phys. Chem. C 2015, 119 (31), 17876−17886. (16) Psofogiannakis, G. M.; McCleerey, J. F.; Jaramillo, E.; van Duin, A. C. T. J. Phys. Chem. C 2015, 119 (12), 6678−6686. (17) Verners, O.; van Duin, A. C. T. Surf. Sci. 2015, 633, 94−101. (18) Fortunelli, A.; Goddard, W. A., III; Sementa, L.; Barcaro, G.; Negreiros, F. R.; Jaramillo-Botero, A. Chem. Sci. 2015, 6 (7), 3915− 3925. (19) Huang, L. L.; Gubbins, K. E.; Li, L. C.; Lu, X. H. Langmuir 2014, 30 (49), 14832−14840. (20) Cheng, H. Y.; Zhu, Y. A.; Chen, D.; Astrand, P. O.; Li, P.; Qi, Z. W.; Zhou, X. G. J. Phys. Chem. C 2014, 118 (41), 23711−23722. (21) Li, W.; Zhu, Y. M.; Wang, G.; Wang, Y.; Liu, Y. J. Mol. Model. 2015, 21 (8), 188. (22) Zhang, Y.; Wang, X. L.; Li, Q. M.; Yang, R.; Li, C. R. Energy Fuels 2015, 29 (8), 5056−5068. (23) Liu, X. P.; Zhan, J. H.; Lai, D. G.; Liu, X. X.; Zhang, Z. J.; Xu, G. W. Energy Fuels 2015, 29 (5), 2987−2997. (24) Li, G. Y.; Xie, Q. A.; Zhang, H.; Guo, R.; Wang, F.; Liang, Y. H. Energy Fuels 2014, 28 (8), 5373−5381. (25) Wang, Q. D.; Wang, J. B.; Li, J. Q.; Tan, N. X.; Li, X. Y. Combust. Flame 2011, 158 (2), 217−226. (26) Chenoweth, K.; Cheung, S.; Van Duin, A. C. T.; Goddard, W. A., III; Kober, E. M. J. Am. Chem. Soc. 2005, 127 (19), 7192−7202. (27) Leininger, J.-P.; Minot, C.; Lorant, F. J. Mol. Struct.: THEOCHEM 2008, 852, 62−70. (28) Scienomics. MAPS Platform, Version 3.4.2; Scienomics: Paris, France, 2015. (29) Sun, H. Fluid Phase Equilib. 2004, 217, 59−76. (30) Chenoweth, K.; van Duin, A. C. T.; Goddard, W. A. J. Phys. Chem. A 2008, 112, 1040−1053. (31) Hockney, R. W.; Eastwood, J. W. Computer Simulation Using Particles; Adam Hilger, Ltd.: New York, 1989. (32) Plimpton, S. J. Comput. Phys. 1995, 117, 1−19. (33) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (34) Lee, C.; Yang, W.; Parr, R. G. Phys. Rev. B: Condens. Matter Mater. Phys. 1988, 37, 785. (35) Weigend, F.; Ahlrichs, R. Phys. Chem. Chem. Phys. 2005, 7, 3297−305. (36) Valiev, M.; Bylaska, E. J.; Govind, N.; Kowalski, K.; Straatsma, T. P.; van Dam, H. J. J.; Wang, D.; Nieplocha, J.; Apra, E.; Windus, T. L.; de Jong, W. A. Comput. Phys. Commun. 2010, 181, 1477.
Figure 6. Evolution of the number of carbon belonging to C20+ molecules (C@C20+, colored lines) and number of carbon atoms belonging to the heaviest molecule (C@Cmax, colored markers) for 1MNaph compounds at T = 3000 (blue), 3500 (green), 3800 (red), and 4000 K (grey).
First, the activation reaction, which consists of the removal of a hydrogen atom from the initial naphthalene derivative, was considered. Results from simulations show that the position and number of methyl substituents not only influence the kinetics but also the thermodynamics of the reaction. Indeed, the number of methyl substituents influences the amount of naphathalene derivatives, which decompose, and, therefore, the thermodynamics of the reaction, whereas the position of the methyl substituent influences the kinetics of that reaction by making the decomposition faster if a substituent is present in position 2. These results were corroborated by DFT calculation of the abstraction energy. With regard to the formation of the char molecule, a threestep mechanism was highlighted with (i) the activation of the naphthalene derivative, (ii) the dimerization or trimerization of the activated compound inducing the formation of many lowweight char precursor molecules, and (iii) the condensation step inducing the formation of very few high-molecular-weight char molecules. The same mechanism was observed for all of the different derivatives studied. Finally, for this reaction, the temperature was shown to have a stronger effect on the kinetics of the condensation step compared to the dimerization/ trimerization step. These results show that ReaxFF molecular dynamics simulation is an appropriate tool to study complex reaction mechanisms and, in the present case, to understand the formation of char compounds during pyrolysis or cracking reactions.
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. Notes
The authors declare no competing financial interest.
■
REFERENCES
(1) Poutsma, M. L. Energy Fuels 1990, 4, 113−131. (2) Jess, A. Fuel 1996, 75, 1441−1448. (3) Wagner, H. G. In Particulate Carbon: Formation during Combustion; Siegla, D. C., Smith, G. W., Eds.; Plenum Press: New York, 1981; p 1−29. (4) Vandenbroucke, M.; Behar, F.; Rudkiewicz, J. L. Org. Geochem. 1999, 30, 1105−1125. (5) Graber, W.-D.; Huttinger, K. J. Fuel 1982, 61, 505−509. E
DOI: 10.1021/acs.energyfuels.6b01062 Energy Fuels XXXX, XXX, XXX−XXX
