Molecular Reactive Force-Field Simulations on the Carbon

Mar 21, 2017 - Yu Ma , Xudong He , Liya Meng , Xianggui Xue , Chaoyang Zhang. Physical Chemistry Chemical Physics 2017 19 (45), 30933-30944 ...
0 downloads 0 Views 6MB Size
Article pubs.acs.org/JPCC

Molecular Reactive Force-Field Simulations on the Carbon Nanocavities from Methane Pyrolysis Xianggui Xue,†,§ Liya Meng,†,‡,§ Yu Ma,† and Chaoyang Zhang*,† †

Institute of Chemical Materials, China Academy of Engineering Physics (CAEP), P.O. Box 919-327, Mianyang, Sichuan 621900, China ‡ College of Material Science and Engineering, Southwest University of Science & Technology, Mianyang, Sichuan 621900, China S Supporting Information *

ABSTRACT: Hydrocarbon pyrolysis is the main way to achieve carbonaceous materials, while most related conversion mechanisms still remain unclear. This work images pyrolysis of methane at various temperatures and densities by molecular reactive force field (ReaxFF) simulations. First, it is interesting to find that the methane decay is dominated by intermolecular collision displacement instead of direct molecular decomposition. Second, a conversion of 1200 methane molecules into a regular carbon nanocavity (CNC) is realized at 3500 K temperature and 0.1 g/cm3 density after a simulation lasting for 10 ns, with 923 carbon atoms and a diameter of 3.4 nm. Such CNC is a perfect precursor of carbon nanotubes, which is confirmed by a sequent simulation on a larger system of 2400 methane molecules and in agreement with several experimental observations. It is found that the CNC growth obeys a polyyne model, without any single aromatic ring formed in the growth. Furthermore, the complex CNC growth appears in some successive stages: primary methane decay, chain elongation and branching, cyclization and condensation, and final sheeting and curling. The regular rearrangement of CNC is thought to be attributed to the limited active centers formed at the initial cyclization and condensation stage; that is, it is a key to control the primary active centers to form regular carbonaceous materials. Polyyne is found in the pyrolysis of both methane and acetylene at high temperatures, suggesting that carbyne, a novel valuable carbonaceous material, may be obtained by hydrocarbon pyrolysis.

1. INTRODUCTION Methane (CH4) is the most stable and simplest hydrocarbon and exists as the maximum component of nature gas. Because of the vast resource of nature gas, most CH4 is directly combusted to supply energy. CH4 is promoted to produce carbonaceous materials,1−4 hydrogen,5−10 and other chemical feedstocks,11,12 such as higher hydrocarbons, synthetic gas, formaldehyde, methanol, and so forth. The chemical conversions from CH4, the same as those from the higher hydrocarbons, are usually implemented by pyrolysis or incomplete combustion, with or without catalyst.13−15 Because of the practical significance of these conversions, many kinetic models have been developed to describe them and set a base for deepening insight and effectively controlling them. In these models, in general, tens or even hundreds of elementary reactions with known kinetic constants were adopted to determine which reaction(s) dominates the total conversions by comparing and fitting with experimental observations.16−29 It implies that these models are too much of a phenomenological sense and are thus insufficient to provide dynamic details of the chemical conversionsit should just be the job of molecular dynamics (MD) simulations. In principle, the MD simulations can provide more details relative to the kinetic models despite much smaller time and space © XXXX American Chemical Society

scales compared with reality; that is, the complex interatomic rearrangements can be carefully checked to clarify the details of a complex chemical process under a given condition.30,31 Moreover, a kinetic model can be deduced from MD simulations if sufficient. The CH4 conversion into a carbonaceous material is a complex chemical process and is strongly dependent on conditions. This is to say, different conditions (e.g., temperature, pressure, catalyst, mixture ratio, etc.) may lead to various polymorphs and morphologies of the C materials. Obviously, it is crucial to reveal the pyrolysis mechanism under various conditions to guide the technological design for manufacturing these carbonaceous materials. In fact, C materials like graphene, C nanotubes and fullerene, and composites based on these C materials have attracted increasing attention due to their excellent potential properties and performances;32−36 however, these materials have not yet been extensively applied, owing to low yields, with a key issue of their controllable and efficient preparation.37−39 As for the conversion mechanism of CH4 by pyrolysis to carbonaceous materials, it has attracted extensive Received: January 11, 2017 Revised: March 18, 2017 Published: March 21, 2017 A

DOI: 10.1021/acs.jpcc.7b00294 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C attention for a long time. Nevertheless, the pyrolysis mechanism, for example, of soot formation, still remains debatable, as two main models of the aromatic hydrocarbons (PAHs) model40 and the polyyne model24 still coexist with each other. With respect to the simulations on pure CH4 pyrolysis, only one study by molecular reactive force field (ReaxFF) MD simulations was reported according to our knowledge. The simulations were conducted on a system of 150 CH4 molecules for 1 ns with various temperatures of 2500 to 3500 K (in 250 K intervals) and densities of 0.05, 0.1, and 0.2 g/cm3. As a result, H2, an extensive range of hydrocarbons, and C dimers were observed in the simulations, and the basic reactions agreed with existing models of thermal CH4 decomposition. Moreover, the formation of the largest molecules was found within the time scale of 1 ns at 3500 K. In principle, the amounts of the C atoms in these largest molecules were adequate to form rings, while ring formation was observed in one case only.41 It seems that the simulations can reproduce some main experimental processes and facilitate to understand the related mechanism. However, it seems that the size and time scales of this simulation are both rather limited. It may cause insufficiency in revealing the underlying mechanism of CH4 pyrolysis. This work also focuses on the CH4 pyrolysis behaviors at high temperatures and pressures by ReaxFF MD simulations, too, with more CH4 molecules and longer time, and an interesting structure, carbon nanocavity (CNC), is obtained. Such CNC is the precursor for forming nanotubes, which have been observed experimentally4 and are confirmed by a subsequent simulation on a system with more CH4. Thus the details of CNC formation by CH4 pyrolysis will be of benefit to carbonaceous material manufacture. In the simulations, we find that the CNC formation obeys a polyyne model, similar to our recent simulation of acetylene pyrolysis to carbon black;42 that is, numerous carbon chains, including polyyne and polyene, appear during the formation. Isolated benzene molecules that are formed in low-temperature CH4 pyrolysis16 are not found in our simulations. Moreover, the finite size effect of simulated systems of CH4 is not found to be significant regarding CH4 decay rates and highest frequency reactions, while the larger simulated system, or the more CH4, facilitates the more sufficient carbon source to form the more perfect carbon nanotubes. This finding significantly supports the polyyne model,24 possibly resulting in a deeper insight into the hydrocarbon pyrolysis at higher temperatures. Moreover, we deem that other hydrocarbon pyrolysis can also be imaged by the methods applied in the present work.

The CH4 pyrolysis can take place at about 1000 K at ambient pressure, while within the time-scale limit of available MD simulations we cannot observe any chemical reaction when the simulation conditions are assigned to be similar to the practical ones. For example, as demonstrated in S1 of the Supporting Information (SI), very few of CH4 molecules are reacted when heated even to 2500 K for 10 ns. Therefore, in contrast with practice, we had to assign three enhanced temperatures 2500, 3000, and 3500 K for simulations to overcome the difficult of simulation time scale of several tens of nanoseconds. As a matter of fact, this strategy of shortening time by increasing temperature has also been adopted for and other hydrocarbons. There is much enhancement of temperatures and pressures applied for simulations.30,42 This enhancement suggests an implication. That is, one would expect recombination to be more active at these increased temperatures and pressures. Meanwhile, the enhancement does not give rise to any artificial feature, as raising temperature does not make any significant change in the results. This will be discussed later. In contrast with reality, it is very common to increase temperature significantly to overcome the time-scale limit for MD simulations based on density function theory (DFT), DFTtight binding, or reactive force-field calculations. Canonical ensemble (NVT) and the molecular reactive force-field ReaxFF combined with some MD methods, like Nosé−Hoover thermostat method, in LAMMPS package43−45 were adopted to simulate the CH4 pyrolysis. Totally, six simulations were individually conducted. Two cells shown in Figure 1 were employed, with each covering three temper-

2. SIMULATION METHODOLOGIES Two factors including structure and condition were concerned for simulating CH4 pyrolysis. In our simulations, first, 1200 CH4 molecules were included in two cubic cells with lengths of 68.38 and 147.3 Å, corresponding to densities (d) of 0.1 and 0.01 g/cm3, respectively. Relative to the density at standard state of 0.0007 g/cm3, these CH4 molecules in the cell of 0.1 g/ cm3 were significantly compressed to a d increased ∼143 times. This compression that causes a higher pressure and a higher CH4 concentration will increase reaction rates and shorten equilibrium time at the same temperature, according to chemical kinetic theory. That is to say, more information will be achieved for these concentrated CH4 molecules within the time scale of our simulations.

atures. The time step and effective relaxation time were set to 0.1 and 20 fs, respectively. Specifically, 200 time steps comprised the coupling frequency of the thermostat to nuclear motion. Atoms in both cells were first relaxed using NVT MD simulations at 300 K for 10 ps (see S2 of the SI, indicating the relaxation equilibrium), and no CH4 decay was found before heating to the given temperatures. Then, six independent simulations each lasted for 10 ns, and the atomic positions and velocities were collected every 1 ps; that is, 10 000 frames were recorded for trajectory analyses. As for the ReaxFF method, its principle can be referred to elsewhere.46−50 ReaxFF is an improved molecular reactive force field and is parametrized to reproduce the density functional theory (DFT) calculation results for given systems and related properties. Compared with many classical force-field methods

Figure 1. Cubic cells containing 1200 CH4 molecules for modeling. (a) d = 0.1 g/cm3 and (b) d = 0.01 g/cm3. C and H atoms are represented in black and green, respectively. These representations are employed in the following figures.

B

DOI: 10.1021/acs.jpcc.7b00294 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C and DFT methods, the prominence of ReaxFF is that it can describe chemical reactions with cheaper calculations. In ReaxFF, bond-order (BO) is regarded as a parameter to determine the bond formation or break. In MD simulations, the BOs of the overall system are updated each step; thereby such bond formation or break can be recognized. Similar to many classical nonreactive force fields, the potential of ReaxFF includes covalent and noncovalent interaction items. The covalent item is a function of bond length, valence angle, torsion angle, lone pair, conjugation, overcoordination, undercoordination, and penalty, and van der Waals and electrostatic interactions are composed of the noncovalent item. Owing to the Pauli principle and the long-range van der Waals attraction (London dispersion), ReaxFF applies a two-body van der Waals to properly describe short-range repulsion. For a special system, the van der Waals interaction can be corrected more exactly by other parametrizations.51 Besides, a shielded Coulomb potential is employed to describe the electrostatic interaction between two charges and makes such electrostatic interaction wellbehaved, even between two bonded atoms. In ReaxFF, the instantaneous atomic charge is determined by the electrostatic field and is described by a simple and second-order of the dependence of the internal energy of the atom on charge.46,47 Three BO cut-offs of 0.55, 0.4, and 0.55 were assigned for simulations in the present work to the atomic pairs of C−C, C−H, and H−H, respectively. That is, the potentials parametrized for CHO systems47 were employed here. ReaxFF has been successfully applied in exploring reaction mechanisms of many hydrocarbons at temperatures ranging from 773.15 to 3500 K,30,41,42,52−57 showing its reliability in revealing underlying mechanism of CH4 pyrolysis by simulations at 2500, 3000, and 3500 K, which is the objective of the this work. A series of FORTRAN scripts were implemented to detailedly describe primary reactions responsible for the earlier decay of CH4 and the subsequent larger molecule formation. All of the identifications (IDs) of the atoms of each reaction can be sought and compared by the script Getreac and can be outputted as a BO file (bonds.reax) through LAMMPS package.43,44 The BO file includes related reaction steps. Besides, all reversible reactions can be deleted by the script Getpath, while keeping the net reactions and their frequencies. Thus we can achieve the primary reactions with the highest net frequencies within given time and preset intervals at various temperatures (see S2 of SI). Besides the above six simulations of 1200 CH4 molecules, an additional three systems of 240, 600, and 2400 CH4 molecules were performed with a same density of 0.1 g/cm3 at 3500 K to clarify the size effect of simulated systems on the pyrolysis mechanism.

Figure 2. Snapshots of CH4 pyrolysis with time proceeding under the condition of the density of 0.1 g/cm3 and the temperature of 3500 K. Only carbon atoms are exhibited in black for clarity.

a stable CNC is formed this time; thus analyses and discussion within this time range are adequate. In comparison, for the other five cases, reactions do not reach completion within 10 ns. In addition, we compare snapshots in Figures S1−S5 of the SI with those in Figure 2 and find a similar evolution mechanism exhibited in all six simulations, only with a difference in reaction rate. That is, to more comprehensively and deeply understand the underlying mechanism of CNC formation, we mainly pay attention to the results corresponding to the highest d of 0.1 g/cm3 and the highest temperature of 3500 K, unless otherwise specified. 3.1. Chemical Species and Their Evolutions. CH4 → C (graphite) + 2H 2 + 75.6 kJ/mol

(1)

The CH4 conversion into perfect graphite and H2 (reaction 1) is thermodynamically disfavored because it is a moderately endothermic process with 75.6 kJ/mol heat absorption in the standard state.58 Reaction 1 cannot be achieved at a common temperature, unless heating to a high temperature. Presumably, heating CH4 will activate them and produce many radicals and molecules. As demonstrated in Figure 3, both the decay of the reactant CH4 and the increase in H2 rapidly proceed within the initial 1 ns. All of these were observed in a previous Reax FF MD simulation.41 As a matter of fact, the pyrolysis of hydrocarbons is a way to H2 production.5−10 This strategy can be evidenced by the rapid H2 increase shown in the Figure. Thereafter, the CH4 decay is complete at ∼1.4 ns, and the number of H2 molecules increases to 2150 and fluctuates around this value, with a difference of about 250 from the theoretical maximum. With respect to atomic H, it reaches a dynamic equilibrium amount of 10 in a short time, only ∼150 ps. Even though these atomic H possess a low occurrence; they play an important role in the initial CH4 decay and subsequent reactions, as they take part in many reactions in combination with other molecules or radicals.

3. RESULTS AND DISCUSSION As pointed out above, six ReaxFF MD simulations were first conducted separately on both CH4 cells (their d are of 0.01 and 0.1 g/cm3, respectively) with 1200 CH4 molecules at 2500, 3000, and 3500 K, and subsequently on three cells of 240, 600, and 2400 CH4 with a density of 0.1 g/cm3 at 3500 K. Because the subsequent simulations were carried out to clarify the size effect, the related results will only be discussed later in Section 3.5. A careful checking shows that, as expected, the CH4 pyrolysis in the cell with the highest d of 0.1 g/cm3 at the highest temperature of 3500 K runs the fastest. From Figure 2, we can find that the CH4 pyrolysis is nearly completed at ∼8 ns; that is, C

DOI: 10.1021/acs.jpcc.7b00294 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

the highest frequency reaction is CH4 + CH4 → CH3 + CH3 + H2 due to very high concentration of CH4. Within initial 0.5 ns, the frequency of eq 2 is 286. As for another main route for the CH4 decay, reaction 6, its frequency is 118, notably below that of reaction 1. Thus the CH4 decay is dominated by the collision with atomic H, followed by the direct decomposition (CH4 → CH3 + H). It also implies an important role of atomic H, that is, promoting the CH4 decay, and the decay is mainly caused by intermolecular collision, instead of direct decomposition. As for the H2 formation, a recent simulation41 exhibited that it proceeds by eq 2 and the combination of two atomic H. This was also proposed in some other experiments of the CH4 pyrolysis.59−61 Nevertheless, the double atomic H combination reaction does not appear in Figure 4 as a high-frequency reaction as well as in Tables S1−S8 of the SI. Even though the combination appears as a high-frequency reaction within 0.2 to 0.5 ns, by careful reaction checking, it still possesses a very low net frequency, attributed to subsequent H2 decomposition, with a close frequency to that of combination. Thus the atomic H combination is not a main source for H2 formation. Combining all reactions related to the H2 formation and dissociation in Table 1 and Tables S1−S8 of the SI, we conclude that H2 is produced by the collisions between H and CH4 (eq 2) and between double CH4 (eq 3) and the direct CH4 decomposition (eq 4). Fitting the CH4 decay data in Figure 3, we obtain a reactant constant k of decomposition of 1.18 × 1012 cm3·mol−1·s−1. The trajectory dynamics study of Marques et al.62 yields k between 7.6 × 1012 cm3·mol−1·s−1 at 3000 K and 2.5 × 1013 cm3·mol−1· s−1 at 3500 K for the reaction Ar + CH4 → Ar + CH3 + H without considering the zero-point energy effect, while k decreases to 2.1 × 1010 cm3·mol−1·s−1 at 3500 K when accounted for. The MD values lie between k obtained by these two different approaches. Lümmen pointed out that the difference is caused by the molecules in the MD simulations colliding with a heat bath in the form of a homogeneous thermostat mechanism instead of Ar atoms.41 Relative to available experimental k, for example, derived from a Arrhenius fitting on shock-wave measurements of dilute CH4−Kr and CH4−Ar systems at 1600−4000 K, 1.7 × 1011 cm3·mol−1·s−1 at 3500 K and 3.5 × 1010 cm3·mol−1·s−1 at 3000 K, and pressures of 30−90% atmospheric pressure,61 our MD result is six times larger. This should be reasonable, as the simulated system is more concentrated. Regarding the C species caused by the CH4 decay at 3500 K, the species possess numerous kinds with varying amounts, and most of these species consist of single to multiple radicals. Conveniently, these species are differentiated in terms of the numbers (n) of the C atoms involved, regardless of the H numbers, and denoted by Cn. At a primary stage, overall, the oligomerization of CH4 takes place with a little decomposition. As illustrated in Figure 5, interestingly, in an n range of 2 to 10,

Figure 3. Evolution of CH4, H2, and H.

By tracing the primary reactions within initial 0.5 ns, we find that, as shown in Figure 4, the CH4 decay mainly proceeds by

Figure 4. Highest frequency reactions within the initial 0.5 ns (the reverse reactions are expunged, and the frequencies are included in brackets). The blue number represents the reaction number.

eqs 2−4, 6, and 9; the atomic H are produced by eq 6, while consumed by eqs 2 and 7; and the H2 production proceeds through the collisions between atomic H with CH4 (eq 2, the highest frequency reaction), between double CH4 molecules (eq 3), and between atomic H and CH3 (eq 7), and the direct decomposition of CH4 (eq 4) and C2H3 (eq 8). Through the frequencies in Figure 4, one can readily understand the change trends of the three chemical species in Figure 3: CH4 is exhausted and H2 is increased, while atomic H reaches equilibrium with a small quantity. From Figure 4 and Tables S1−S8 of the SI we can find that reaction 2, CH4 + H → CH3 + H2, possesses the highest frequency, except from the case of the initial 0.1 ns, in which

Table 1. Frequencies of the Highest Frequency Reactions within the Initial 0.5 ns of Various Systems (The Reverse Reactions Are Expunged) reactions CH4 CH4 CH4 CH3 CH4

+ H → CH3 + H2 + CH4 → CH3 + CH3 + H2 → CH2 + H2 + CH3 → C2H6 → CH3 + H

240 CH4

600 CH4

1200 CH4

2400 CH4

62 37 33 26 16

184 63 70 52 50

286 188 166 129 118

628 379 282 248 274

D

DOI: 10.1021/acs.jpcc.7b00294 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

observations at low temperatures. A careful examination shows C6H6 presents a linear structure only instead of six-membered rings. It should be attributed to the high simulation temperature, disfavoring stabilizing the benzene ring. Our simulation results are consistent with some experiments of thermal CH4 decomposition, in which C2 hydrocarbons (C2H6, C2H4, and C2H2) are formed59,63−66 by the formation of either C2H6 or C2H4 and their stepwise dehydrogenation. The continued dehydrogenation of C2 hydrocarbons can form carbon dimers from C2H2. Nevertheless, they do not aggregate further but are often rehydrogenated to form C2H or C2H2. Arutyunov et al.67 reported a large C2H4 to C2H2 ratio at 1100 and 1300 K. This tendency can also be observed in Figure 6b. Furthermore, Arutyunov et al.66 and Guèret et al.18 found the maxima of the concentration of the main C2 species. This was also confirmed by our simulation, as illustrated in Figure 6b. In addition, the absence of ring formation can be explained with the comparably high temperatures adopted in our simulations, at which the structure of benzene ring is unstable.18 As time proceeds to >1 ns, the polymerization is enhanced (see the structures in Figure 2f−q) with increasing the sizes of certain C species, reducing the total amounts of low C species, and causing the final CNC formation, which will be discussed later. 3.2. Formation Mechanism of CNC. Several typical snapshots in Figure 2 represent different stages of CNC formation. To be clearer, we abstracted these typical structures to Figure 7 to show the stages of the CNC formation as follows. (1) CH4 decay. This has been discussed above. It should be stressed again that the decay mainly undergoes intermolecular collision displacement instead of direct molecular decomposition. (2) Chain elongation and branching. Figure 2a exhibits the primary state of CH4, in which all CH4 are evenly and randomly distributed, like solute molecules dissolved in equilibrated solution. Figure 2b shows several elongated chains. However, the chains found are relatively short this time. As time reached 0.1 ns, the C atoms on the longest chains increase to 4 and to 8 at 0.2 ns, as shown in Figure 7b,c, respectively. As pointed out above, the absence of benzene ring formation can be understood with the comparably high temperatures adopted in our simulations, at which the benzene ring is unstable;18 whereas, the CH4 pyrolysis at 3500 K causes the formation of straight chains, each composed of several tens of C atoms. And few H atoms are found to be linked to the chains. It is just the so-called polyyne and polyene. The polyyne and polyene are formed by intermolecular collisions of low C species. As time resumes, chain length increases and branched chains appear, as exhibited in Figures 2c and 6c. We think that chain branching is an indispensable result when the chains are elongated to a certain extent. In principle, the end−middle fusion of neighboring straight chains is responsible for such branching. It will be enhanced when the chains are elongated, as the middle sites become more and more kinetically dominant during addition than the ends, despite the higher activity of the ends of chains than the middle sites. Figure 2d exhibits the further chain elongation and branching at 0.5 ns. It is similar to the case of chain elongation and branching during acetylene pyrolysis to carbon black.42 (3) Cyclization and cycle condensation. The branched chains are subsequently cyclized, with the fusion of neighboring C chains, instead of the curvature of any straight chain alone. At this point, the cyclization cannot be explained by the PAH

Figure 5. Evolution of C2 to C10 species within the initial 1 ns.

the amounts of Cn reach orderly decreasing maxima with increasing n. In combination with the decrease in CH4, these results suggest that n increases by consuming CH4 or other C1 species, as most n increases with a step of 1 within a time of initial 1 ns and an n range of 2 to 10. In our recent simulation of C2H2 pyrolysis,42 most n increases with a step of 2 at the primary stage. It shows that in the cases of CH4 and C2H2 pyrolysis C1 and C2 are stable elementary units for n increasing, respectively. Because it is crucial to clarify the primary structures and their amounts to reveal the underlying mechanism for the final CNC formation, we carefully examined all Cn (n = 1 to 6) species within the first 1 ns. C1 exists in four types, namely, CH, CH2, CH3, and CH4, and no isolated C atom is found. This result agrees with the recent MD simulations by ReaxFF41 and many experiments.59−61 As illustrated in Figure 6a, the maximal amounts of these C1 species decrease with the decreasing H amounts in C1. From Figure 6a we can find that the CH2 and CH3 increase to maximums at 0.1 ns, and thus the main related primary reactions were checked in this time range in Table S1 of the SI. As a result, CH3 is formed by multiple routes such as CH4 + CH4 → CH3 + CH3 + H2, CH4 + H → CH3 + H2, CH4 → CH3 + H, and CH2 + CH4 → CH3 + CH3, while the CH2 formation proceeds mainly by CH4 → CH2 + H2. These highfrequency primary reactions also appear in Figure 4, showing their dominance within a long period. C2 is immediately produced after CH4 molecules were heated to 3500 K for ∼15 ps. C2H6 appears first to increase to reach a local maximum mainly through the collision between double CH3 (CH3 + CH3 → C2H6), followed by other C2 including C2Hm (m = 0−5), as shown in Figure 6b. A careful examination of primary reactions shows that, similar to C1 formation, most of these C2 are dominantly formed by intermolecular collision displacement instead of direct partition, as C2H + H → C2 + H2, C2H2 + H → C2H + H2, C2H3 + H → C2H2 + H2, C2H4 + CH3 → C2H3 + CH4, C2H5 + CH4→ C2H4 + CH3 + H2, and C2H6 + CH4 → C2H5 + CH3 + H2 for C2, C2H, C2H2, C2H3, C2H4, and C2H5, respectively. As time proceeds above 250 ps, C2H and C2H2 dominate the C2 subsequently. C2H appears as frequently as C2H2, favoring subsequent chain elongation by its polymerization. The cases of C3 to C6 are similar to those of C1 and C2; that is, their formations are mainly attributed to intermolecular collisions (Figure 6c−f). It should be noted that no isolated benzene molecule or ring is found all along, significantly varying from E

DOI: 10.1021/acs.jpcc.7b00294 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 6. Evolution of the amounts of various C species at T = 3500 K and d = 0.1 g/cm3.

model.40 The cyclization is the additional result of the active sites in the middle and terminal ends of C chains. As indicated in the bottom-right section of Figures 2e and 7e, an evident polyring structure was found. The single-ring structure has not been checked. Afterward, the cycle condensation takes place and seems abrupt. As a matter of fact, we find some five- to eight-membered rings, which were observed in a previous quantum-chemical molecular simulation of C deposition.68 Cyclization, in particular, to six-membered rings, is energetically favored. Thus further cycle condensation in Figures 2f−h and 7f is understandable as nucleation exhibits in the CNC formation. (4) Sheeting and curling of condensed rings and final CNC formation. As demonstrated in Figure 2h, at 3 ns, H2, straight and branched C chains, and condensed C cycles with various quantities are contained in the cell. A rudiment of CNC (Figure

Figure 7. Typical structures showing the CNC formation stages within initial 4 ns: (a) initial CH4 molecule, (b) primary chain elongation, (c,d) more elongation and branching, (e) cyclization and condensation, (f,g) more condensation, and (e) sheeting and curling.

F

DOI: 10.1021/acs.jpcc.7b00294 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

checked their distribution against time. It can readily be concluded from Figure 10 that Cn are developed from the

7g) is formed by sheeting and curling the condensed rings. As time elapses, most C species are assembled into a regular CNC, as shown by Figure 2i−q. Figures 8 and 9 clearly show the

Figure 10. Evolution of the C number of the maximal C species.

primary small sizes and large quantities to final large sizes and small quantities. During the growth of Cn into the final CNC, fusion among big Cn plays an important rule. For instance, fusion leads to a rapid size increase in the maximal Cn, as shown by a break within 3.5 to 4 ns in Figure 10; that is, the C number of the maximal Cn does not increase gently, just showing the fusion. Understandably, the rapid growth should not occur through step-by-step additions of C1 or other small Cn, which are the primary pyrolysis product of CH4. The similar case is also exhibited in C2H2 pyrolysis.42 Meanwhile, H partition from C atoms proceeds successively. An initial increase, a succedent fluctuation, and a final decrease in H atoms in the largest Cn are exhibited in Figure S8 of the SI. First, by both the addition and the fusion, forming larger Cn will primarily enhance the H abstraction. Thereafter, fluctuation is generated in a manner of combination−dissociation− combination of the large Cn. This manner can be deduced from Figure 10, representing a fluctuation of C numbers of the largest Cn, as well. Around 3.5 ns, the abrupt H increase, the same as the above-mentioned abrupt C increase, is attributed to the fusion. The final decrease is a result of the successive H partition from the stable maximal Cn. Namely, as shown in Figure S9 of the SI, H atoms are successively partitioned from Cn to form final CNC within 10 ns, resulting in a continuous reduction in H/C ratios, consistent with the increase in the abstracted H atoms (Figure 3), including atomic and molecular ones. Next, we pay attention the C−C bonds, because their types reflect the related chemical structures. As stated above, this system involves complex chemical reactions, and thus a statistical method of radial distribution function (RDF, g(r)) was adopted for analyses. RDF can provide assigned interatomic distances in any complex system within a cutoff.72 We obtain C−C RDFs in Figure 11 by analyzing the structures in Figure 2. The dominant C−C distances of the first 0.5 ps (most evident peaks in Figure 11a−d) are 1.23 and 1.33/1.37 Å, representing CC and CC bonds, respectively. This finding indicates the dominance of the triple and double states of C atoms at this stage. At 1 ns, the g(r) of 1.37 Å is faded and that of 1.42 Å appears (Figure 11e), showing the beginning of

Figure 8. Typical structures showing the final stages in CNC formation: top view (top) and side view (bottom).

Figure 9. CNC structure at 10 ns along various directions. (a−e) Side views with clockwise rotation of 45° in turn and (f) coronal structure of the CNC.

regular structures of CNC: a shallow cavity with a compact wall, composed of five- to eight-membered rings. On the edge of the CNC wall, there are some hanging chains, in expectation of growth to a perfect nanotube. Moreover, The simulated CNC possesses 923 C atoms, a ratio to the total C atoms of about 3/4 (923/1200), and a diameter of about 3.4 nm (Figure 9). Although, the CNC is not perfect, attributed to the deficient C source, presumably, it will become more and more perfect if C source is sufficient enough. Thus we think that our simulations verify the carbon nanotube formation by CH4 pyrolysis. Furthermore, we find that polyyne and polyene chains are chemically bonded in or with these structures. These findings suggest that carbyne69−71 may be manufactured by both CH4 and C2H2 pyrolysis (carbyne was also be observed in the simulations of C2H2 pyrolysis)42 once we have controlled the temperature and concentration perfectly to prevent branching and cyclization because of a higher stability of carbyne at high temperatures, in contrast with other C polymorphs.69 Obviously, the final CNC with a large number of C atoms is formed, at the cost of the great consumption of small Cn. Figure S7 of the SI illustrates the evolution of the total amount of all Cn upon time: a sharp amount decrease in the first 500 ps, followed by a gentle decrease to a final equilibrium amount of ∼100 at 8 ns thereafter. To detail the Cn evolution, we G

DOI: 10.1021/acs.jpcc.7b00294 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 11. Radial distribution functions (RDFs) of C−C at different times.

PAH formation. Subsequently, the g(r) of 1.45 Å is enhanced (Figure 11f−p), showing the CNC growth with a higher and higher degree of perfect cycle condensation. Besides, the dominance of g(r) of 1.33/1.37 Å lasts from 0.1 to 3.5 ns, and the ever-existing of g(r) of 1.23 Å suggest a polyyne model24 for CNC formation. 3.3. Temperature and Density-Dependent C Condensation. As pointed out above, a similar evolution mechanism is undergone for CH4 pyrolysis under various conditions within the simulation time limit of 10 ns, with a difference in rate. In this subsection, we concerned ourselves with the effects of temperature and d on the CH4 decay and the size increase in the largest Cn. In total, Figure 12 indicates that both the higher temperature and the higher d cause the faster CH4 decay, and both temperature and d play important roles in the decay. As stressed before, the CH4 consumption is mainly attributed to the intermolecular collisions, in fact, the intermolecular reactions. It implies that d increases accelerate the decay. High temperature weakens chemical bonds and facilitates reaction increases, too. That is to say, both high temperature and high d aid the CH4 decay. Moreover, comparing the decay gradients of the cases of 0.1 g/cm3 at 3000 K and 0.01 g/cm3 at 3500 K in Figure 12, it seems that increasing temperature is more efficient to increase the CH4 decay than increasing d, as the former exhibits the faster decay. This result accords with the previous MD simulation.41 As for the sizes of the largest Cn, they change widely. In the case of d of 0.01 g/cm3, Figure 13a shows the evolution of the C numbers of the largest Cn at the three temperatures of 2500, 3000, and 3500 K within 10 ns. The difference in size evolution can be explainable, as a higher temperature always promotes

Figure 12. CH4 decay under different conditions.

reactions while it does not always enhance C condensation. Prior to 2 ns, a higher temperature causes a faster size increasing. This should be attributed to the fact that a higher temperature leads to more intermolecular collisions to decompose CH4 and more collision chances among Cn, that is, the faster C condensation. However, for another thing, a higher temperature can also promote the break of condensed Cn. Thus a size is equilibrated at 3500 K after 2 ns. Regarding this, a smooth size increase to a top and then a fluctuation at 3000 K are observed, remaining a larger equilibrium size, relative to 3500 K. At 2500 K, almost no CH4 decay occurs because of the too-low temperature to activate CH4. In the case H

DOI: 10.1021/acs.jpcc.7b00294 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 13. Evolution of C numbers of the maximal C species at densities of (a) 0.01 and (b) 0.1 g/cm3.

of d of 0.1 g/cm3, it exhibits similar tendencies due to the same reason (Figure 13b). Combining Figure 13 and Figures S1−S5 of the SI, we can conclude that CNC formation does not serve as a single result within the 10 ns of the time limit. It suggests that other structures like polyyne may be formed under appropriate conditions. 3.4. Pyrolysis Comparison of Methane and Acetylene. Combined with our recent study of C2H2 pyrolysis to carbon black42 and this study of CH4 pyrolysis to CNC, we have imaged two series of complex processes of conversion of simple C1 to complicated high C particles. Despite a same given condition of 3500 K temperature and 0.01 g/cm3 density for the two separate simulations, significantly different final structures are obtained, suggesting that different C sources may lead to different carbonaceous materials. This difference requires us carry out insight into a comparison of detailed formation mechanism in Figure 14a1−a6,b1−b5. In the fast stage, the primary decay mechanisms of C2H2 and CH4 are much different. For C2H2, its decay is mainly by H abstraction, while for CH4, its decomposition is dominated by intermolecular replacement instead of H abstraction. Sub-

sequently, at the elongation stage of short C chains, additions of C2 species dominate the C 2 H 2 pyrolysis, while the intermolecular replacement controls the C chain elongation as the CH4 decay. This is the reason for C-chain elongation through even-numbered and one-by-one C additions for C2H2 and CH4, respectively. In our simulations, we find that the C lengths of straight and branched chains of C2H2 pyrolysis are much larger than those of CH4 pyrolysis, as chains shown by Figure 14a3 versus 14b3. We think this difference in chain length is stemmed from the difference in the elongation mechanism. When the pyrolysis proceeds to a certain extent, the addition is, in general, easier than the intermolecular replacement because the addition possesses many more active sites than replacement. It is just this chain length difference that makes the difference in subsequent cyclization and condensation. In both pyrolysis processes, the cyclization takes place abruptly without single aromatic rings. However, because of longer chains and therefore more active sites the cyclization in C2H2 pyrolysis is much more multiply localized relative to CH4 pyrolysis (Figure 14a4 and 14a5 vs 14b4 and 14b5). As time elapses, only one C sheet center is formed and sets a base for following sheet enlargement for CH4 pyrolysis, while such centers are found on many places for C2H2 pyrolysis, as Figure 14a6 and 14a7 versus 14b6. As time proceeds further, the final CNC and carbon black particles are formed, with many differences in polymorph and morphology (Figure 14a8 vs 14b7), even though the crystallization zones of both particles are composed of graphene sheets; that is, CNC looks regular, while carbon black appears amorphously. Overall, the different C sources under a same assigned condition may lead to various carbonaceous materials with significant differences in polymorph and morphology. From the above discussion, it seems that it would be an efficient approach to obtain regular structures by controlling the active sites. As a matter of fact, the applications of catalysts are a way to control the active sites to get required carbonaceous materials.13−15 3.5. Size Effect of Simulated Systems on the CNC Formation. To clarify the size effect of simulated systems on the pyrolysis mechanism and product (C particles) formation, three systems with 240, 600, and 2400 CH4 molecules and a same density of 0.1 g/cm3 were applied to conduct NVT MD simulations at 3500 K, and 10 ns lasted for both simulations of 240 and 600 CH4 molecules and 5.6 ns for 2400 CH4 molecules. Thereby, we compare the simulated results of four CH4 systems, composed of 240, 600, 1200, and 2400

Figure 14. Comparison of (a) CH4 pyrolysis to CNC and (b) C2H2 pyrolysis to carbon black. I

DOI: 10.1021/acs.jpcc.7b00294 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

of 10 ns. As a result, apart from the difference in evolution rate, it is found that the pyrolysis mechanism is less dependent on the assigned CH4 d and temperature within the time limit; that is, both the CH4 decay and H2 production are mainly the result of intermolecular replacement reactions involving CH4 instead of H abstraction from CH4. This case is much different from the C2H2 pyrolysis, which is initiated from H abstraction and C2 addition. In all six separate simulations with a same time limit of 10 ns, only one corresponding to a condition of 3500 K and 0.1 g/cm3 reached equilibrium with a stable CNC formed. Careful checking showed the CNC formation undergoes four stages including initial CH4 decay, subsequent chain elongation and branching, cyclization and cycle condensation, and final sheeting and curling. After these stages, a stable CNC with a diameter of 3.4 nm, 923 C atoms, and a very low H content was obtained. Such CNC was a precursor that is structurally close to a C nanotube. Thus a perfect C nanotube was expected with a sufficient C source, as a rather perfect one was confirmed by a subsequent simulation of a larger system of 2400 CH4. Moreover, the comparison of C2H2 pyrolysis to carbon black and CH4 pyrolysis to CNC was carried out. It was thought that the different C sources leading to final C materials should be attributed to the differences in primary decay, chain elongation, and cyclization mechanisms, and it was believed that it could be an efficient way to limit active site amount to obtain regular C materials. The chain elongation, chain branching, chain fusion to cycles, and condensed cycles that serve as the intermediate processes have been observed in our simulations. They are widely different from the reported observations at relatively low temperatures. Besides, the polyyne chains are found to be stable at high temperature. Thereby, the polyyne model should be more appropriate for hydrocarbon pyrolysis at higher temperature, in combination with the previous study of C2H242 and this work of CH4.

molecules, respectively. Regarding the CH4 decay, as shown in Figure S10 of the SI, any evident difference in decay rate is not found. That is, for all four systems, the CH4 decay undergoes almost the same: Half of CH4 molecules are decayed at 0.07 ns, and the decay is finished at ∼1 ns. To detail the initial reaction steps, the highest frequency reactions were checked. As a result, we find that reactions CH4 + H → CH3 + H2, CH4 + CH4 → CH3 + CH3 + H2 and CH4 → CH2 + H2 possess the top three of highest frequencies for all cases; reactions CH3 + CH3 → C2H6 and CH4 → CH3 + H possess, respectively, the third and fourth highest frequencies for the systems of 240, 600, and 1200 CH4 molecules, and this order reverses for the system of 2400 CH4 molecules. With consideration of a small difference in frequency of reactions CH3 + CH3 → C2H6 and CH4 → CH3 + H, we think that there is no difference in the initial steps for all four systems with various sizes. That is to say, it shows no size effect on initial CH4 decay mechanism based on the almost same decay rate and initial steps in Table 1. With respect to the shape of the final carbon cluster structure of each simulation, it depends on the carbon source or the size of the simulated system. As demonstrated in Figure 15 (1−4),



ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.7b00294. Comparison of pyrolysis evolution under various conditions, evidence of the complete relaxation of the simulated system after 10 ps at 300 K, high-frequency primary reactions within various periods, formation mechanism of CNC, and size effect of simulated systems on the CNC formation. (PDF)

Figure 15. CNC structures at 10 ns of 240 (1), 600 (2), and 1200 (3) CH4 systems and carbon a nanotube structure at 5.6 ns of a 2400 CH4 system (4) along various directions. (a−e) Side reviews with clockwise rotation of 45° in turn. (f) Top views.



the carbon clusters vary from a shallow CNC to a carbon nanotube when carbon atoms become increasingly sufficient. Checking the growth details of CNC and carbon nanotubes in Figure 2 and Figures S11−S13 of SI, we find that the growth mechanism are the same, that is, by the processes of abovementioned chain elongation and branching, cyclization and condensation, sheeting, and curling. The final structure of a smaller simulated system can be regarded as an intermediate of a larger one, suggesting the similar growth mechanism. Although both the carbon source and simulation time of this study are limited, we believe that final carbon structure tends to a perfect carbon nanotube by increasing the source and time.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel: 86-816-2493506. ORCID

Chaoyang Zhang: 0000-0003-3634-7324 Author Contributions §

X.X. and L.M. contribute equally.

Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS We are grateful for the financial support from the National Natural Science Foundation of China (21673210 and U1530262).

4. CONCLUSIONS We simulated CH4 pyrolysis at temperatures of 2500, 3000, and 3500 K, and densities of 0.1 and 0.01 g/cm3, with a time limit J

DOI: 10.1021/acs.jpcc.7b00294 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C



(23) Dors, M.; Nowakowska, H.; Jasiński, M.; Mizeraczyk, J. Chemical Kinetics of Methane Pyrolysis in Microwave Plasma at Atmospheric Pressure. Plasma Chem. Plasma Process. 2014, 34, 313− 326. (24) Krestinin, A. V. Detailed Modeling of Soot Formation in Hydrocarbon Pyrolysis. Combust. Flame 2000, 121, 513−524. (25) Billaud, F.; Gueret, C.; Weill, J. Thermal Decomposition of Pure Methane at 1263 K. Experiments and Mechanistic Modeling. Thermochim. Acta 1992, 211, 303−322. (26) Murphy, D. B.; Carroll, R. W.; Klonowski, J. E. Analysis of Products of High-temperature Pyrolysisi of various Hydrocarbons. Carbon 1997, 35, 1819−1823. (27) Arutyunov, V. S.; Vedeneev, V. I. Pyrolysis of Methane in the Temperature Range 1000−1700 K. Russ. Chem. Rev. 1991, 60, 1384− 1397 and references therein.. (28) Dean, A. M. Detailed Kinetic Modeling of Autocatalysis in Methane Pyrolysis. J. Phys. Chem. 1990, 94, 1432−1439. (29) Grenda, J. M.; Androulakis, I. P.; Dean, A. M., Jr; Green, W. H. Application of Computational Kinetic Mechanism Generation to Model the Autocatalytic Pyrolysis of Methane. Ind. Eng. Chem. Res. 2003, 42, 1000−1010. (30) Ding, J.; Zhang, L.; Zhang, Y.; Han, K. A Reactive Molecular Dynamics Study of n-heptane Pyrolysis at High Temperature. J. Phys. Chem. A 2013, 117, 3266−3278. (31) Zhang, C.; Wen, Y.; Xue, X. Self-enhanced Catalytic Activities of Functionalized Graphene Sheets in the Combustion of Nitromethane: Molecular Dynamic Simulations by Molecular Reactive Force Field. ACS Appl. Mater. Interfaces 2014, 6, 12235−12244. (32) Anantram, M. P.; Leonard, F. Physics of Carbon Nanotube Electronic Devices. Rep. Prog. Phys. 2006, 69, 507−561. (33) Harris, P. J. F. Carbon Nanotube Science: Synthesis, Properties and Application; Cambrige University Press, 2009. (34) Wong, H. S. P.; Akinwande, D. Carbon Nanotube and Graphene Device Physics; Cambrige University Press, 2011. (35) Stankovich, S.; Dikin, D. A.; Dommett, G. H. B.; Kohlhaas, K. M.; Zimney, E. J.; Stach, E. A.; Piner, R. D.; Nguyen, S. T.; Ruoff, R. S. Graphene-based Composite Materials. Nature 2006, 442, 282−286. (36) Zhang, J.; Liu, X.; Blume, R.; Zhang, A.; Schlögl, R.; Su, D. S. Surface-modified Carbon Nanotubes Catalyze Oxidative Dehydrogenation of n-butane. Science 2008, 322, 73−77. (37) Bagri, A.; Mattevi, C.; Acik, M.; Chabal, Y. J.; Chhowalla, M.; Shenoy, V. B. Structural Evolution during the Reduction of Chemically Derived Graphene Oxide. Nat. Chem. 2010, 2, 581−587. (38) Stankovich, S.; Dikin, D. A.; Piner, R. D.; Kohlhaas, K. A.; Kleinhammes, A.; Jia, Y.; Wu, S. T.; Nguyen, R. S.; Ruoff, R. S. Synthesis of Graphene-based Nanosheets via Chemical Reduction of Exfoliated Graphite Oxide. Carbon 2007, 45, 1558−1565. (39) Ravagnan, L.; Siviero, F.; Lenardi, C.; Piseri, P.; Barborini, E.; Milani, P.; Casari, C. S.; Li Bassi, A.; Bottani, C. E. Cluster-beam Deposition and in situ Characterization of Carbyne-rich Carbon Films. Phys. Rev. Lett. 2002, 89, 285506. (40) Harris, S. J.; Weiner, A. M. Twenty-Second Symposium (International) on Combustion; The Combustion Institute: Pittsburgh, PA, 1988; p 333. (41) Lümmen, N. ReaxFF-molecular dynamics simulations of nonoxidative and non-catalyzed thermal decomposition of methane at high temperatures. Phys. Chem. Chem. Phys. 2010, 12, 7883−7893. (42) Zhang, C.; Zhang, C.; Ma, Y.; Xue, X. Imaging the C Black Formation by Acetylene Pyrolysis with Molecular Reactive Force Field Simulations. Phys. Chem. Chem. Phys. 2015, 17, 11469−11480. (43) Plimpton, S. J. Fast Parallel Algorithms for Short-range Molecular Dynamics. J. Comput. Phys. 1995, 117, 1−19. (44) Aktulga, H. M.; Fogarty, J. C.; Pandit, S. A.; Grama, A. Y. Parallel Reactive Molecular Dynamics: Numerical Methods and Algorithmic Techniques. Parallel Comput. 2012, 38, 245−259. (45) Hoover, W. Canonical dynamics: Equilibrium Phase-space Distributions. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 1695−1697.

REFERENCES

(1) Barnett, A. Graphite Formation from Low Temperature Pyrolysis of Methane over Some Transition Metal Surfaces. Nature 2008, 221, 1044−1046. (2) Ermakova, M. A.; Ermakov, D. Y.; Kuvshinov, G. G.; Plyasova, L. M. New Nickel Catalysts for the Formation of Filamentous Carbon in the Reaction of Methane Decomposition. J. Catal. 1999, 187, 77−84. (3) Saenko, N. S.; Ziatdinov, A. M. Multi-walled Carbon Nanotubes Synthesized by Methane Pyrolysis: Structure and Magnetic Properties. Solid State Phenom. 2014, 213, 60−64. (4) Kiselev, N. A.; Krestinin, A. V.; Raevskii, A. V.; Zhigalina, O. M.; Zvereva, G. I.; Kislov, M. B.; Artemov, V. V.; Grigoriev, Yu. V.; Hutchison, J. L. Extreme-length Carbon Nanofilaments with Singlewalled Nanotube Cores Grown by Pyrolysis of Methane or Acetylene. Carbon 2006, 44, 2289−2300. (5) Fincke, J. R.; Anderson, R. P.; Hyde, T. A.; Detering, B. A. Plasma Pyrolysis of Methane to Hydrogen and Carbon Black. Ind. Eng. Chem. Res. 2002, 41, 1425−1435. (6) Muradov, N.; Smith, F.; Huang, C.; T-Raissi, A. Autothermal Catalytic Pyrolysis of Methane as A New Route to Hydrogen Production with Reduced CO2 Emissions. Catal. Today 2006, 116, 281−288. (7) Abanades, S.; Kimura, H.; Otsuka, H. Hydrogen Production from CO2-free Thermal Decomposition of Methane: Design and On-sun Testing of a Tube-type Solar Thermochemical Reactor. Fuel Process. Technol. 2014, 122, 153−162. (8) Ibrahim, A. A.; Fakeeha, A. H.; Al-Fatesh, A. S.; Abasaeed, A. E.; Khan, W. U. Methane Decomposition over Iron Catalyst for Hydrogen Production. Int. J. Hydrogen Energy 2015, 40, 7593−7600. (9) Choudhary, T. V.; Sivadinarayana, C.; Chusuei, C. C.; Klinghoffer, A.; Goodman, D. W. Hydrogen Production via Catalytic Decomposition of Methane. J. Catal. 2001, 199, 9−18. (10) Muradov, N. Hydrogen via methane Decomposition: An Application for Decarbonization of Fossil Fuels. Int. J. Hydrogen Energy 2001, 26, 1165−1175. (11) Savinov, S. Y.; Lee, H.; Song, H. K.; Na, B. K. A Kinetic Study on the Conversion of Methane to Higher Hydrocarbons in a RadioFrequency Discharge. Korean J. Chem. Eng. 2004, 21, 601−610. (12) Bhavsar, S.; Veser, G. Chemical Looping Beyond Combustion: Production of Synthesis Gas via Chemical Looping Partial Oxidation of Methane. RSC Adv. 2014, 4, 47254−47267. (13) Robertson, S. D. Carbon Formation from Methane Pyrolysis over Some Transition Metal Surfaces-I. Nature and Properties of the Carbons Formed. Carbon 1970, 8, 365−374. (14) Muradov, N. Catalysis of methane decomposition over elemental carton. Catal. Commun. 2001, 2, 89−94. (15) Horn, R.; Schlögl, R. Methane Activation by Heterogeneous Catalysis. Catal. Lett. 2015, 145, 23−39. (16) Billaud, F.; Gueret, C.; Weill, J. Thermal Decomposition of Pure Methane at 1263 K. Experiments and Mechanistic Modeling. Thermochim. Acta 1992, 211, 303−322. (17) Olsvik, O.; Billaud, F. Thermal Coupling of Methane. A Comparison between Kinetic Model Data and Experimental Data. Thermochim. Acta 1994, 232, 155−169. (18) Guèret, C.; Daroux, M.; Billaud, F. Methane Pyrolysis: Thermodynamics. Chem. Eng. Sci. 1997, 52, 815−827. (19) Hidaka, Y.; Sato, K.; Henmi, Y.; Tanaka, H.; Inami, K. ShockTube and Modeling Study of Methane Pyrolysis and Oxidation. Combust. Flame 1999, 118, 340−358. (20) Slovetskii, D. I.; Mankelevich, Y. A.; Slovetskii, S. D.; Rakhimova, T. V. Mathematical Modeling of the Plasma-Chemical Pyrolysis of Methane. High Energy Chem. 2002, 36, 44−52. (21) Fau, G.; Gascoin, N.; Gillard, P.; Steelant, J. Methane Pyrolysis: Literature Survey and Comparisons of Availabledata for Use in Numerical Simulations. J. Anal. Appl. Pyrolysis 2013, 104, 1−9. (22) Kahle, L. C. S.; Roussiére, T.; Maier, L.; Delgado, K. H.; Wasserschaff, G.; Schunk, S. A.; Deutschmann, O. Methane Dry Reforming at High Temperature and Elevated Pressure: Impact of Gas-Phase Reactions. Ind. Eng. Chem. Res. 2013, 52, 11920−11930. K

DOI: 10.1021/acs.jpcc.7b00294 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C (46) Van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard, W. A., III ReaxFF: A Reactive Force Field for Hydrocarbons. J. Phys. Chem. A 2001, 105, 9396−9409. (47) Chenoweth, K.; van Duin, A.; Goddard, W. ReaxFF Reactive Force Field for Molecular Dynamics Simulations of Hydrocarbon Oxidation. J. Phys. Chem. A 2008, 112, 1040−1053. (48) Naserifar, S.; Liu, L.; Goddard, W. A.; Tsotsis, T. T.; Sahimi, M. Toward a Process-Based Molecular Model of SiC Membranes. 1. Development of a Reactive Force Field. J. Phys. Chem. C 2013, 117, 3308−3319. (49) Naserifar, S.; Goddard, W. A.; Liu, L.; Tsotsis, T. T.; Sahimi, M. Toward a Process-Based Molecular Model of SiC Membranes. 2. Reactive Dynamics Simulation of the Pyrolysis of Polymer Precursor To Form Amorphous SiC. J. Phys. Chem. C 2013, 117, 3320−3329. (50) Naserifar, S.; Goddard, W. A.; Tsotsis, T. T.; Sahimi, M. First Principles-based Multiparadigm, Multiscale Strategy for Simulating Complex Materials Processes with Applications to Amorphous SiC Films. J. Chem. Phys. 2015, 142, 174703. (51) Liu, L. C.; Liu, Y.; Zybin, S. V.; Sun, H.; Goddard, W. A., III ReaxFF-lg: Correction of the ReaxFF Reactive Force Field for London Dispersion, with Applications to the Equations of State for Energetic Materials. J. Phys. Chem. A 2011, 115, 11016−11022. (52) Castro-Marcano, F.; van Duin, A. C. T. Comparison of Thermal and Catalytic Cracking of 1-Heptene from ReaxFF Reactive Molecular Dynamics Simulations. Combust. Flame 2013, 160, 766−775. (53) Page, A. J.; Moghtaderi, B. Molecular Dynamics Simulation of the Low-Temperature Partial Oxidation of CH4. J. Phys. Chem. A 2009, 113, 1539−1547. (54) Cheng, X.; Wang, Q.; Li, J.; Wang, J.; Li, X. ReaxFF Molecular Dynamics Simulations of Oxidation of Toluene at High Temperatures. J. Phys. Chem. A 2012, 116, 9811−9818. (55) Wang, Q.; Hua, X.; Cheng, X.; Li, J.; Li, X. Effects of Fuel Additives on the Thermal Cracking of n-Decane from Reactive Molecular Dynamics. J. Phys. Chem. A 2012, 116, 3794−3801. (56) Wang, Q.; Wang, J.; Li, J.; Tan, N.; Li, X. Reactive Molecular Dynamics Simulation and Chemical Kinetic Modeling of Pyrolysis and Combustion of n-dodecane. Combust. Flame 2011, 158, 217−216. (57) Castro-Marcano, F.; Kamat, A. M.; Russo, M. F.; van Duin, A. C. T.; Mathews, J. P. Combustion of An Illinois No. 6 Coal Char Simulated Using An Atomistic Char Representation and the ReaxFF Reactive Force Field. Combust. Flame 2012, 159, 1272−1285. (58) Luo, Y. Handbook of Bond Dissociation Energies; Science Press: Beijing, 2005. (59) Olsvik, O.; Billaud, F. Modelling of the Decomposition of Methane at 1273 K in a Plug Flow Reactor at Low Conversion. J. Anal. Appl. Pyrolysis 1993, 25, 395−405. (60) Dean, A. J.; Hanson, R. K. CH and C-atom Time Histories in Dilute Hydrocarbon Pyrolysis: Measurements and Kinetics Calculations. Int. J. Chem. Kinet. 1992, 24, 517−523. (61) Kiefer, J. H.; Kumaran, S. S. Rate of Methane Dissociation over 2800−4300 K: the Low-pressure-limit Rate Constant. J. Phys. Chem. 1993, 97, 414−420. (62) Marques, J. M. C.; Martínez-Núñez, E.; Fernández-Ramos, A.; Vázquez, S. A. Trajectory Dynamics Study of the Ar + CH4 Dissociation Reaction at High Temperatures: the Importance of Zero-Point-Energy Effects. J. Phys. Chem. A 2005, 109, 5415−5423. (63) Choudhary, T. V.; Aksoylu, E.; Goodman, D. W. Nonoxidative Activation of Methane. Catal. Rev.: Sci. Eng. 2003, 45, 151−203. (64) Lampe, F. W.; Field, F. H. Reaction of Gaseous Ions. VII. Methane-Hydrogen and the Proton Affinities of Methane and Ethane. J. Am. Chem. Soc. 1959, 81, 3242−3244. (65) Tse, J. S.; Klug, D. D.; Laasonen, K. Structural Dynamics of Protonated Methane and Acetylene. Phys. Rev. Lett. 1995, 74, 876. (66) Kumar, P.; Marx, D. Understanding Hydrogen Scrambling and Infrared Spectrum of Bare CH5+ Based on Ab Initio Simulations. Phys. Chem. Chem. Phys. 2006, 8, 573−586. (67) Arutyunov, V. S.; Vedeneev, V. I.; Moshkina, R. I.; Ushakov, V. A. Pyrolysis of Methane under Static Conditions at 1100−1400K. Kinet. Catal. 1991, 32, 234−240.

(68) Wang, Y.; Page, A. J.; Nishimoto, Y.; Qian, H. J.; Morokuma, K.; Irle, S. Template Effect in the Competition between Haeckelite and Graphene Growth on Ni(111): Quantum Chemical Molecular Dynamics Simulations. J. Am. Chem. Soc. 2011, 133, 18837−18842. (69) Whittaker, A. G. Carbon: A New View of Its High-Temperature Behavior. Science 1978, 200, 763−764. (70) Liu, M.; Artyukhov, V.; Lee, H.; Xu, F.; Yakobson, B. Carbyne from First Principles: Chain of C Atoms, a Nanorod or a Nanorope. ACS Nano 2013, 7, 10075−10082. (71) Chalifoux, W. A.; Tykwinski, R. R. Synthesis of Polyynes to Model the sp-carbon Allotrope Carbyne. Nat. Chem. 2010, 2, 967− 971. (72) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon Press, Oxford Science Publications, 1987.

L

DOI: 10.1021/acs.jpcc.7b00294 J. Phys. Chem. C XXXX, XXX, XXX−XXX