Article pubs.acs.org/EF
Revealing the Initial Chemistry of Soot Nanoparticle Formation by ReaxFF Molecular Dynamics Simulations Song Han,†,‡ Xiaoxia Li,*,†,‡ Fengguang Nie,† Mo Zheng,† Xiaolong Liu,†,‡ and Li Guo†,‡ †
State Key Laboratory of Multiphase Complex Systems, Institute of Process Engineering, Chinese Academy of Sciences, Beijing 100190, People’s Republic of China ‡ University of Chinese Academy of Sciences, Beijing 100049, People’s Republic of China S Supporting Information *
ABSTRACT: This work presents long time ReaxFF MD simulations of fuel-rich combustion for up to 10 ns to explore the initial mechanism of soot nanoparticle formation. A 24-component rocket propellant 1 (RP-1) model based on the major components of RP-1 fuel was employed. Simulations were performed by GPU-accelerated code GMD-Reax, and reactions therein were revealed with the aid of VARxMD. Simulated evolution of physical and chemical properties of the largest molecule exhibits the overall structural transitions of three stages for incipient ring formation, nucleation, and graphitization from fuel molecules to the formation of a single soot nanoparticle. The incipient ring formation takes place in stage 1 by large ring generation from activated aliphatic polyyne-like chains, ring number increase from internal bridging between carbon atoms of large rings, and consequent formation of PAH-like molecules with aliphatic side chains. Nucleation of a nanoparticle in stage 2 is the result of coalescence of PAH-like molecules, accompanied by ring closure reactions that occurred at side chains of the PAHlike molecules, and following formation of internal bridged bonds. Graphitization of the nanoparticle in stage 3 is dominated by the transformation from C5/C7 rings to C6 rings with C3 rings as intermediates. This work demonstrates that ReaxFF MD simulation might be a promising approach for qualitatively characterizing morphological evolution and the underlying chemical complexity of soot formation using multicomponent fuel models.
1. INTRODUCTION Understanding the soot formation mechanism is of great importance for improving the performance of engines and emission control. It is well-accepted that the soot formation process starts with thermal cracking of the fuel and formation of radical precursors, followed by formation of the first aromatic ring and its growth to polycyclic aromatic hydrocarbons (PAHs), then nucleation or inception of particles with enlargement by surface growth of gas phase hydrocarbons, and coalescence by particle−particle collisions, graphitization, or carbonization of particles from amorphous to graphitic, oxidation, and particle aggregation.1−4 Most soot formation mechanisms include detailed chemical reaction pathways with involvement of PAHs.5−8 While recent experimental results9−13 and modeling work14−16 provide some evidence on possible roles of aliphatics in soot formation, which may be important to improve soot models,17 description at an elementary level of the transition from soot precursors to particles is not straightforward due to the limited experimental data and complementary predictive models.18−20 It is still difficult to fully characterize the evolution of soot formation from fuel molecules. With the rapid development and fast growing applications of the reactive force field (ReaxFF) proposed by van Duin and Goddard et al.,21,22 reactive molecular dynamics (ReaxFF MD) simulation shows the capability and great potential in capturing the complex chemical evolution at high temperature and pressure in organic systems. ReaxFF is a potential based on bond order to describe breaking and formation of chemical bonds which remedies the deficiency of the classical MD © 2017 American Chemical Society
methods in describing chemical reactions, with parameters extensively optimized from experimental data and quantum mechanics calculations. It enables fast computation for large models (>1000 atoms) with precision close to the DFT method.23−25 In addition to bonded interaction energy terms, the ReaxFF potential allows for descriptions of conjugated, hydrogen bonding, van der Waals, and Coulomb interactions, which are important for representing chemical formation of aromatic compounds and crystallite of soot particles in the soot formation process.21 Polarization effects are accounted for by dynamic atomic charge updated by electron equilibration method (EEM) at each MD time step.21,26,27 Moreover, no predefined reaction pathways required makes it a prospective method in predicting the complex chemistry in fuel combustion, especially under extreme conditions.23,24,28 Recent applications of ReaxFF MD simulations of coal char24 and carbon black model29 to study the structural evolution have shown the potential of ReaxFF MD to depict some details on soot formation. The coal char combustion simulations by Castro-Marcano et al.24 investigated the char oxidation process and pointed out that the char transition pathways included the conversion of six-membered rings into five-membered and seven-membered rings that would further decompose or react with mostly O and OH radicals. Zhang et al.29 reported simulations of carbon black formation at 3500 K by ReaxFF MD using a model of 1000 acetylene molecules, in which Received: April 26, 2017 Revised: June 21, 2017 Published: June 27, 2017 8434
DOI: 10.1021/acs.energyfuels.7b01194 Energy Fuels 2017, 31, 8434−8444
Article
Energy & Fuels Table 1. Species and Composition of the 24-Component RP-1 Model
morphological evolution of the carbonaceous material and finally the structure of folded graphite layers were observed within 2 ns. The methodology development for ReaxFF MD in the authors’ group further extends its capability for investigating reaction complexity with GMD-Reax30 and VARxMD.31 GMDReax30 is a GPU-accelerated code for ReaxFF MD simulations, which allows significantly improved computational efficiency for simulations of large models within accepted machine hours. Comparing the performance of GMD-Reax running on a single GPU of C2050 and LAMMPS32 Reax codes on eight CPU cores, GMD-Reax allows for speed up as high as 11 times over the LAMMPS FORTRAN codes and 5 times over LAMMPS C codes for single precision, and up to 7.8 times and 4 times for double precision.33 VARxMD31 is a unique code specialized for reaction analysis and visualization of ReaxFF MD simulations, which helps in revealing the complex reaction mechanism that is difficult to access manually. In this study, long time ReaxFF MD simulations of 24component rocket propellant 1 (RP-1) models were performed for 5 ns to study the influences of temperature, equivalence ratio, and concentration of reactants (density) on soot formation. A fuel-rich condition with an equivalence ratio of 5, density of 0.1 g/cm3, and 3000 K was chosen, at which simulations were extended to 10 ns to explore the soot nucleation and graphitization mechanism. Fuel model construction, simulation strategy, and reaction pathway analysis are described in this work. Morphological evolution of a wellshaped fullerenic nanoparticle with a diameter of over 3 nm and
the underlying reaction pathways in its nucleation and graphitization are presented and discussed.
2. METHOD 2.1. Model Construction. Liquid fuels derived from petroleum sources are complex mixtures containing hundreds or even more components. Due to the expensive computational cost, a fuel model is often limited to a surrogate in which only a few components can be directly described.34 In order to explore the detailed chemical mechanism of soot formation in fuel combustion with models as close as possible to real fuel, a RP-1 model with 24 components was constructed using assignment of major components based on GC-MS analysis of RP-1 fuel with the chromatographic peak area counts in excess of 1% in Bruno’s work.35 Weight fractions of the components in Table 1 were estimated based on the uncalibrated peak area, and the number of each species is assigned based on its weight fraction. Table 1 lists structures and compositions of the RP-1 components. All of the 102 fuel molecules (3702 atoms) were used in the simulations. The model contains various paraffins, cycloalkanes, and aromatics to describe the diversity of chemical structures as much as possible in real fuel. Soot formation in fuel oxidation is a complex process that can be affected by temperature, equivalence ratio, and concentration of soot building blocks (related to concentration of fuel species).36 To study the effects of oxygen amount on soot formation, RP-1 models were created with different equivalence ratios of 0.5, 1, 2, 5, and 10 at a density of 0.1 g/cm3 for RP-1 oxidation simulations of 4068−11040 atoms. Details can be found in Table S1 in the Supporting Information. To investigate the effects of density (concentration of fuel species) on soot formation behavior, RP-1 models were constructed with different densities of 0.01, 0.02, 0.05, 0.1, 0.2, and 0.5 g/cm3 at an equivalence ratio of 5 as shown in Table S1. Since the 8435
DOI: 10.1021/acs.energyfuels.7b01194 Energy Fuels 2017, 31, 8434−8444
Article
Energy & Fuels model construction processes are similar, the construction of the RP-1 model with a density of 0.1 g/cm3 and an equivalence ratio of 5 is described here as an example. The model was constructed at 300 K using the amorphous cell module in Materials Studio software.37 The edge length of the cubic simulation box is 78.2 Å. Geometry optimization was carried out with a DREIDING force field38 using the forcite module in Materials Studio. Then dynamics simulation of the canonical ensemble (NVT) was performed for system equilibrium at 300 K to stretch the molecular structures and to get better spatial distributions, followed with the energy minimization to refine the configurations to reduce the effects introduced by temperature jump from the NVT equilibrium to the followed high temperature reactive MD simulations. Both the NVT-MD simulation and energy minimization were performed using nonreactive DREIDING force field. The temperature was controlled by Berendsen thermostat39 with a damping constant of 0.1 ps. 2.2. Simulation Strategies and Details for RP-1 Combustion. To investigate the temperature effects on soot evolution, the ReaxFF MD simulations were performed at 1500−4000 K with an interval of 500 K using the NVT ensemble. To improve the computational efficiency, the GPU-accelerated code30 was employed. A strategy of jumping temperature to the target temperature was employed to avoid the influences on analysis of chemical reactions in the temperatureramping process. The high simulation temperatures favor the observation of soot formation at practical computational costs and are close to or within the range of the combustor temperature of the liquid oxygen/RP-1 rockets.40 Three parallel ReaxFF MD simulations were performed for a simulation condition at an equivalence ratio of 5, a density of 0.1 g/cm3, and 3000 K for 5 ns, which was chosen for feasible observation of soot formation in RP-1 oxidation. One of the simulations was extended to 10 ns to have an overall evolution picture of the formed nanoparticle as much as possible. The temperature was controlled by Berendsen thermostat39 with a damping constant of 0.1 ps. Periodic boundary condition was applied. The time step was set to 0.1 fs to ensure the precision of the simulations. The output interval of the simulation trajectories was set to 1 ps for close observation of morphology and analysis of the chemical reactions. The 10 ns of ReaxFF MD simulation took 43.6 days running on a CentOS 6.5 server with 2* Intel Xeon E5-2650v2 2.6 GHz, 64 GB RAM and an NVIDIA Tesla K20C GPU attached. To facilitate observation of chemical reactions within practical machine time, high temperature and density were employed. The 10 ns of RP-1 combustion simulation in this work was performed with a dense system at 54−384 MPa, which is somewhat close to the condition for RP-1 combustion in rocket engines, though the pressure may be high for the typical flame conditions.28,41 2.3. Analysis of the Reaction Mechanism. Using ReaxFF MD trajectories as input, detailed chemical reactions were obtained with the aid of VARxMD.31 The analysis started with generation of all the detailed chemical reactions in a list. Then the reactions, particularly reactions representing evolution of soot formation, were selected and analyzed by observing the 3D structure of reactants, products, and reaction sites. Combined with species analysis, nucleation and graphitization of a nanoparticle in the soot formation process were identified with the 3D frame views of the whole reaction model.
Figure 1. Evolution of C number and H/C ratio of the largest molecule in parallel fuel-rich combustion simulations of a 24component RP-1 model by ReaxFF MD.
reactivity evolution with time was analyzed by counting the number of bond breaking and formation at an interval of 0.1 ns for 0−2 and 1.0 ns for 2−10 ns, as shown in Figure 2. The
Figure 2. Reactivity involved by the largest molecule at different times.
smaller interval used at the early stage within 0−2 ns was chosen based on the fact of faster growing of the C number, and correspondingly a larger interval was used during the period of 2−10 ns for its slower evolution at the later stage. Number evolution of five-/six-/seven-membered rings in the largest molecule is shown in Figure 3. Density evolution of the largest molecule is shown in Figure 4, which was calculated based on the mass and its van der Waals (VDW) volume.37 The data point intervals both in Figures 3 and 4 are 0.1 ns within 0− 2 ns and 0.5 ns during 2−10 ns. As shown in Figure 1, C number increased to 934 on average at 2 ns and then slowly approached 1087 at 10 ns. The largest molecule observed at 10 ns is the only nanoparticle in the whole system, and details on the remaining small molecules can be found in Table S4. There was a rapid increase of C number in the largest molecule beginning at ∼1 ns, followed by drastic fluctuations until ∼2 ns. After that, its size growing became moderate and tended to be stable. The number of reactions in the largest molecule was found to be much larger during 1−2 ns than in 0−1 and 2−10 ns. The drastic fluctuation of C
3. RESULTS AND DISCUSSION 3.1. Evolution of Physical and Chemical Properties for the Largest Molecule. To characterize the soot formation process in the simulations of the 24-component RP-1 model, evolution of physical and chemical properties for the largest molecule was investigated. The largest molecule identified from the simulation trajectory was considered as a representative for the evolution of the soot formation process, which was defined as the chemical structure with maximal number of carbon atoms formed in the simulation system at any moment during the soot formation process. The evolution trends of C number and H/C ratio of the largest molecule are shown in Figure 1. Its 8436
DOI: 10.1021/acs.energyfuels.7b01194 Energy Fuels 2017, 31, 8434−8444
Article
Energy & Fuels
number during 1−2 ns (Figure 1) corresponds to frequent fluctuations of the bond breaking and formation (Figure 2), which suggests the existence of reversible reactions in the early stage of soot inception process. More details on the fluctuations observed in parallel simulations during this period can be found in Figure S6. This observation supports the proposed reversibility in soot formation process, whose important role was recently addressed by Eaves et al.11,42 Based on the transitions of the C number and the number of reactions (reactivity), evolution of the largest molecule can be divided into three stages of 0−1 ns (stage 1), 1−2 ns (stage 2), and 2− 10 ns (stage 3). Stage 1 is characterized by starting of the molecular size growth as indicated by the carbon number increase in Figure 1, which results in the birth of a larger structure with C5−C7 rings formed as shown in Figure 3. At ∼1 ns, the start of a quick C5−C7 ring increase corresponds to the increase of C number and reactivity as displayed in Figures 1 and 2, which may suggest that coalescence of large molecules was occurring (transition from stage 1 to stage 2 at ∼1 ns). Thus, in stage 2, both the number of C6 rings and the number of C5 and C7 rings increase rapidly. In stage 3, the number of C6 rings increases continuously at a fast rate, while the continuous increase rate of C5 and C7 rings slows down and tends to stabilize with time. The phenomenon can be explained by graphitization starting in the largest molecule with the transformation from C5 and C7 rings into the aromatic C6 ring structures, which would be further explained in section 3.3. Meanwhile, the H/C ratio (Figure 1) decreases continuously in the simulation, with a value of ∼0.14 at 2 ns and ∼0.02 at 10 ns, of which the descending tendency is consistent with the decrease of H/C ratio from 0.3−0.4 to about 0.1 observed in the premixed flame along the soot formation region.43 The relatively small H/C ratio at 10 ns obtained in the simulation was close to that of a fullerenic carbon43 due to the continuous dehydrogenation reactions in the graphitization process (stage 3) without further adding of fuel species in the simulation box. The graphitization process will be discussed in section 3.3.
Figure 3. Evolution of ring number in the largest molecule.
Figure 4. Density evolution of the largest molecule in parallel simulations.
Figure 5. Snapshots of system evolution representing soot nucleation within 2 ns of fuel-rich combustion simulation. Carbon atoms are in black, hydrogen in cyan, and oxygen in red. The structures representing incipient ring formation and coalescence into nucleation are highlighted in red circles. 8437
DOI: 10.1021/acs.energyfuels.7b01194 Energy Fuels 2017, 31, 8434−8444
Article
Energy & Fuels
Figure 6. Example of structural evolution leading to the formation of a PAH-like molecule with aliphatic chains observed in the fuel-rich combustion simulation. Key reaction sites of bond formation are highlighted in red circles.
3.2. Initial Chemistry for Nucleation of Soot Nanoparticle. The molecular structure and morphological evolution obtained in the fuel-rich combustion simulation confirmed the first transition from the stage of incipient ring formation to nucleation of a carbonaceous nanoparticle (from stage 1 to stage 2) within 2 ns as shown in Figure 5. The initial stage for the nucleation of a nanoparticle started with thermal decomposition of fuel species, followed by formation of polyyne-like structures (very low H/C ratio as shown in Table S3) at 0.2 ns. The ring formation was observed at 0.3 ns, followed by formation of PAH-like molecules and nucleated nanoparticle between 0.7 and 2 ns. Detailed structural evolution for the formation of a PAH-like molecule with aliphatic chains is shown in Figure 6. As shown in Figures 5 and 6, a large amount of dehydrogenation reactions of chain species and cycloparaffin components can be observed within 0.2 ns, accompanied by fewer intramolecular H-migration and C−C bond-scission reactions, leading to the generation of aliphatic radicals including polyyne-like radicals. The dominance of bimolecular H-abstraction over unimolecular dehydrogenation, C−C bond scission, and intramolecular H-migration in the initial stage with data are provided in Table S2, which suggests that the very high pressure condition (0.1 g/cm3) favors bimolecular reactions (H-abstraction) more than unimolecular reactions (mainly C−C bond scission). Length growth of a few chains and ring closure reactions of polyyne-like radicals can be found at 0.3 ns in Figure 5, which roughly agree with the description that large polyyne molecules or radicals can react with each other and may contain ring closure reactions.45 Formation of a PAH-like molecule with aliphatic chains can be found at 0.7 ns in Figure 5. The reaction pathways to form the PAH-like molecule are further depicted in Figure 6.
Average relative errors in three stages for the 5 ns parallel simulations are 24.6%, 15.9%, and 5.6% for C number and 9.3%, 5.4%, and 13.2% for H/C ratio of the largest molecule. More details can be found in Figure S8 and Table S6 in the Supporting Information. As shown in Figure 4, the density of the largest molecule increased quickly between 0 and 0.7 ns, which corresponds to the reactions from fuel molecules to the formation of PAH-like molecules. Then the value tended to a plateau at ∼1 ns, meanwhile the C number of the largest molecule experienced a steep jump with fluctuations, indicating nucleation of a nanoparticle from coalescence of the PAH-like molecules. Then the density increased quickly to ∼2.0 g/cm3 at 2 ns. After that, the increase of density became slow and approached a value of ∼2.3 g/cm3 at 10 ns, corresponding to the structural reorganization of the nanoparticle. Evolution of the density is related to the structural transformation of soot nanoparticle and its H/C ratio (Figure 1). Reaction details of structural transformation will be discussed in sections 3.2 and 3.3. The density of ∼2.0 g/cm3 at 2 ns is close to the reported value of 1.8−2.0 g/cm3 for the soot particles,11,44 but the density of ∼2.3 g/cm3 at 10 ns is larger. Continuous dehydrogenation and excessive reorganization of the soot nanoparticle to the graphitic or fullerenic carbon at 10 ns in the high temperature condition are likely responsible for the large density. Relative error of density is within 3% for the 5 ns parallel simulations; more details can be found in Figure S8 and Table S6. Evolution of properties shown in Figures 1−4 suggests three stages for soot formation, stage 1 for incipient ring formation before 1 ns, stage 2 for nucleation of the carbonaceous nanoparticle during 1−2 ns, and stage 3 for the consequent graphitization of the nanoparticle during 2−10 ns. The reaction mechanisms in nucleation and graphitization of soot nanoparticle will be discussed in sections 3.2 and 3.3, respectively. 8438
DOI: 10.1021/acs.energyfuels.7b01194 Energy Fuels 2017, 31, 8434−8444
Article
Energy & Fuels
Figure 7. Schematic representation of a nanoparticle generated from coalescence of two PAH-like molecules. The PAH-like molecules are in simplified representation as PAH-like cores with side chains considering the fact that it is difficult to distinguish the structures and reaction sites of the deformed morphology of the molecules from the snapshots of trajectory.
to the enlargement of the molecule into C920H127 at 2.0 ns. As shown in Figure 5, bending of the large PAH-like core with its increase in ring number was observed; meanwhile the multilayered PAHs tend to evolve into an ellipsoid-shaped nanoparticle. In addition to the carbonaceous nanoparticle, major products obtained in the simulation of RP-1 combustion are H2O and H2 molecules, with the evolution provided in Figure S5. At the early stage of the 10 ns fuel-rich simulation, a large number of small hydrocarbons (e.g., CH3, C2H2, C2H4, and C4H2) and a small number of CO and CO2 can be observed. Much of them consumed later to form larger components. The reaction pathways were provided in Figure S4 in Supporting Information. Reaction pathways of the nanoparticle nucleation revealed in 0−2 ns roughly agree with the four stages of chain elongation, chain branching, cyclization and ring densification, and condensed ring folding for carbon black particle formation obtained from the 2 ns of ReaxFF MD simulations of acetylene pyrolysis by Zhang et al.29 Evolution of various C-containing species and morphology of key structures were characterized in Zhang’s simulations. This work revealed further detailed reaction pathways on the formation of PAH-like molecules and their coalescence into the nucleation of a soot nanoparticle, especially the effects of aliphatic side chains on the coalescence of PAH-like molecules for nucleation of the nanoparticle. Moreover, the long time evolution of the nucleated nanoparticle observed uncovers the consequent graphitization process, which will be discussed in section 3.3. 3.3. Graphitization of the Nucleated Nanoparticle. The 10 ns simulation in the fuel-rich combustion condition
As shown in Figure 6, with the radical generation and dehydrogenation at the initial stage, ring closure can be found resulting in the formation of a six-membered ring with aliphatic chains at 0.513 ns or the formation of a large ring at 0.525 ns, both by cross-link of polyyne-like structures. Internal bridged bond formation of a 14-membered large ring observed at 0.589 ns led to the formation of two interconnected eight-membered rings. Moreover, the frequent internal bond formation of the large ring was responsible for the structural transformation of the interconnected rings during the period of 0.589 and 0.650 ns, followed by the generation of one PAH-like molecule with aliphatic chains at 0.740 ns. Then it became a stage of amount increase for five-/six-/seven-membered rings by coalescence of small PAH-like molecules, as indicated by the quick increase of ring number in Figure 3. Formation of a nanoparticle occurred that was accompanied by surface growth reactions with small hydrocarbon molecules and structural transformation of PAHlike molecules. A schematic representation for formation of carbonaceous nanoparticle from coalescence of large PAH-like molecules is provided in Figure 7. The coalescence of two large molecules of C398H73O and C355H59 was initiated by linking aliphatic side chains to form C753H131O at 1.342 ns, followed by further linking side chains that led to the closure of a large ring in C753H131O at 1.344 ns. Then the condensation of the large ring started through internal bond cross-linking or bridged bond formation at the radical sites, resulting in the increase of ring number of C753H131O at 1.5 ns. The ring condensation of C753H131O continued, accompanied by size growth by side chains linking with other small hydrocarbon molecules, leading 8439
DOI: 10.1021/acs.energyfuels.7b01194 Energy Fuels 2017, 31, 8434−8444
Article
Energy & Fuels
of large rings still exists by cyclization of polyyne-like chains and formation of internal bonds for the generation of PAH structures. Very slow reorganization of the C skeleton driven by the reactions on dehydrogenized sites can be observed. The inner and outer surfaces of the nanoparticle became clearly observable after 2 ns. Surface growth reactions between the nanoparticle and small hydrocarbon molecules occurred on the sites of edge defects or on the outer surface of the nanoparticle, but none can be observed on the inner surface in this work. The process of dehydrogenation would lead to the increase of the aromatic carbon network.47 An example of typical reaction sites of graphitization in the nanoparticle for bond breaking at 2.061 ns and bond formation at 2.062 ns is shown in Figure 9. Bond breaking generated new radical sites; particularly breaking of the bridged bond of a C3 and a C5 ring would convert the C5 ring into the C6 ring structure. C6 ring was generated either through formation of the bridged bond between a C3 and a C6 ring within the C7 ring or through formation of internal bonds within large ring structure that was formed by ring closure of polyyne-like chains. Therefore, three-membered rings are unstable, serving as intermediate structures for the reversible transformation from C5/C7 rings to C6 rings, which play an important role in the slow transformation of the skeleton structure in the nucleated nanoparticle to the fullerenic structures. This mechanism for graphitization revealed agrees with the conclusion that rearrangement of skeleton structures is promoted by reduction of dangling bonds and elimination of edge defects from the experiments.48−50 These reaction pathways can be helpful for explaining the evolutions of C number and C5−C7 rings as shown in Figures 1 and 3 that the number of C atoms and six-membered rings in the largest molecule keep their increase with time, but become slower at the later stage, and the increase of number of five-/sevenmembered rings at an early stage tends to a plateau at the later
allows for investigating the graphitization of the nucleated carbonaceous nanoparticle formed at 2 ns. Figure 8 shows the
Figure 8. Well-shaped fullerenic structure formed in the 10 ns simulation.
morphological evolution of the nanoparticle after its nucleation in stage 3. The nanoparticle formed at 10 ns is a well-shaped, spherical fullerene-like nanoparticle under the condition of fuelrich combustion. Close observation of the morphology of nanoparticle suggests that the incepted nanoparticle evolved from the multidirectional oriented PAH structure at 2 ns into a three-layered onion-like structure with less aliphatic side chains at 4 ns (similarities and differences of the morphology of the nanoparticles observed in parallel simulations are provided in Table S9). Further cyclization of side chains of the nanoparticle and transformation of its ring structures would finally lead to the formation of a well-shaped fullerenic structure with few side chains at 10 ns. Evolution from the carbonaceous nanoparticle in irregular shape at 2 ns to the relatively ordered fullerenic structure at 10 ns is slower compared with the nucleation process, which agrees with the conclusion that fullerenic soot formation under the flame conditions can be much longer than nascent soot formation.46 In this stage, C-containing species of small radicals nearly ran out without further supplement in the simulation box. Closure
Figure 9. Example of typical reaction sites for graphitization of the soot nanoparticle observed in the 10 ns fuel-rich combustion simulation. Bond breaking is highlighted in red, and bond formation in green of RP-1 by ReaxFF MD. 8440
DOI: 10.1021/acs.energyfuels.7b01194 Energy Fuels 2017, 31, 8434−8444
Article
Energy & Fuels stage. Although the increase of C number and six-membered rings of the largest molecule can be explained by the widely accepted hydrogen abstraction acetylene addition (HACA) mechanism,51 the reaction pathways revealed above may suggest other pathways for the phenomenon. C5−C7 rings can be generated from ring closure reactions (Figures 6 and 9). The five-/seven-membered rings will convert into sixmembered rings to reduce the structural defects, which may counteract the amount increase of C5 and C7 rings due to ring closure reactions. Since there is no more feed of fuel molecules into the simulation box, the number of five-/seven-membered rings becomes almost stable in the later stage to maintain the spherical structure of the nanoparticle. In addition, long-range organized multiple-shell hollow structure in spherical morphology was observed in the graphitization of the nanoparticle as shown in Figure 8, which was stable as indicated by its properties approaching steady close to 10 ns, including the evolution of C number and H/C ratio (Figure 1), number of ring structures (Figure 3), and density (Figure 4). The graphitized soot nanoparticle structure can endure high temperature conditions of thermal decomposition. The observed continuous graphitization agrees with the point of view that sufficient time under high temperature conditions is feasible for high energy barrier reactions of soot skeleton rearrangement to fullerenic structures.49,52 The diameter of the fullerenic nanoparticle generated at 10 ns in the combustion simulation is over 3 nm. Given the limited scale of the model (1212 C atoms), the observed nanoparticle is almost the largest nanoparticle that can be formed in the RP1 model. Considering the size of the formed nanoparticle (over 3 nm) and the box (edge length of 7.82 nm), as well as the fact that there is no more feed of fuel molecules into the simulation box, the remaining molecules and space are not favorable to form another nanoparticle. Large models with C atoms more than ∼10,000 should have good opportunity to have more than one nanoparticles generated. Though the description of a primary soot particle generated in the flame conditions with a diameter of 20−60 nm has to use a much larger model of ∼1,000,000 atoms,53 simulation of such a large model is well beyond the current computational capability. With a computational cost of more than 40 days for one condition of 10 ns simulations in this work, the chemical reactions and the structural details observed for a nanoparticle nucleation and graphitization in this work should be useful in understanding the soot formation process. 3.4. Effects of Temperature, Equivalence Ratio, and Density on Soot Formation. Effects of a single factor including temperature, equivalence ratio, and density on the soot particle formation were investigated. Simulations were performed for 5 ns at conditions of 1500−4000 K at an interval of 500 K, equivalence ratios of 0.5, 1, 2, 5 and 10, and densities of 0.01, 0.02, 0.05, 0.1, 0.2, and 0.5 g/cm3. The C number of the largest molecule representing the evolution tendency of soot nanoparticle formation is shown in Figure 10. Morphology of the largest molecules for all these simulations at 5 ns was provided in Table S7 in the Supporting Information. More detailed morphological evolution of the whole systems and the largest molecules at 3500 and 4000 K can be found in Table S8. On the influence of temperature, 3000 K is a condition most favoring formation of soot nanoparticle as observed in Figure 10a (equivalence ratio of 5 and density of 0.1 g/cm3), at which the biggest nanoparticle generated among all simulation temperatures. C number increases slowly when simulation
Figure 10. Effects of equivalence ratio, temperature, and density on C number of the largest molecule formed in 5 ns ReaxFF MD simulations of RP-1 oxidation.
temperature is lower than 3000 K, among which 1500 and 2000 K are not feasible simulation temperatures to have a nanoparticle observed within 5 ns due to its very slow reactions, and 2500 K is relatively low that soot nucleation 8441
DOI: 10.1021/acs.energyfuels.7b01194 Energy Fuels 2017, 31, 8434−8444
Article
Energy & Fuels
membered rings, with the involvement of three-membered ring intermediates. The spherical and hollow-centered morphology of the nanoparticle is a result of continuous dehydrogenation and excessive reorganization of the soot nanoparticle. By combining high performance computing of GPU-enabled program GMD-Reax with unique capability of VARxMD for chemical reaction analysis, long time ReaxFF MD simulations of the multicomponent fuel model may be capable to qualitatively show an overall picture of morphological evolution for the formation of a soot nanoparticle and have its reaction pathways uncovered computationally, though further validations are needed in future work.
occurs later than that of 3000 K. At 3500 K, decomposition of the formed soot nanoparticle can be observed during 3.5−5 ns as shown in Figure 10a, due to the breaking of hydrocarbon side chains from the nanoparticle (Table S8). A temperature of 4000 K is too high that none of a molecule with C number larger than 300 can be observed (Figure 10a), at which only large rings from polyyne-like molecules can be generated, but none of PAH-like molecules or nanoparticles can be observed. The observations provide reasonable descriptions of temperature effects on soot formation tendency that suggests the simulation temperature 3000 K can be accepted for observing chemical events of soot formation process in this work. As shown in Figure 10b, at the conditions of 3000 K and 0.1 g/cm3, the fuel-rich combustion conditions of the equivalence ratios of 2, 5, and 10 are feasible for soot formation, among which the equivalence ratio of 5 is the most favorable condition for observing the formation of a soot nanoparticle. For the stoichiometric (equivalence ratio of 1) and fuel-lean (equivalence ratio of 0.5) conditions, none of the PAH-like molecules or nanoparticles can be observed. The simulated observations support that fuel-rich conditions favor formation of soot nanoparticles. The observations provide reasonable description of the equivalence ratio effects on soot formation tendency, which suggests the rationality of choosing equivalence ratio of 5 at 3000 K for soot formation observation. The effect of density on soot formation tendency at 3000 K and equivalence ratio of 5 is shown in Figure 10c. The density of 0.1 g/cm3 not only allows for generation of a soot nanoparticle within 5 ns, but its mild evolution rate can facilitate the observation of its structural evolution and reaction details. Apparently, soot nanoparticle can be formed within 5 ns for density of 0.02−0.5 g/cm3. The varied density of the simulated system as shown in Figure 10c can give some sense of pressure effects on the evolution tendency of soot nanoparticle formation. The soot nanoparticle forms early, and the size of the soot nanoparticle increases with density increasing from 0.02 to 0.5 g/cm3, which means that high pressure favors the formation of a large soot nanoparticle. When pressure is very high (density from 0.1 g/cm3 increased to 0.2 g/cm3 or even to 0.5 g/cm3), the carbon number of the soot nanoparticle tends to a stable value around 1100, which seems to agree with what has been reported in a recent review, that the maximal soot generation seems to asymptotically reach a plateau when the pressure exceeds the critical value.54
■
ASSOCIATED CONTENT
S Supporting Information *
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.energyfuels.7b01194. Details on the RP-1 oxidation models at different equivalence ratio and density, further discussions on reaction paths and species distributions, simulation results of the three parallel simulations, morphologies of the largest molecules, simulation results using threecomponent and 24-component RP-1 models, discussions on time step and ensemble for the simulations, and ReaxFF force field parameters used in this work (PDF)
■
AUTHOR INFORMATION
Corresponding Author
*Phone: 86-10-82544944. Fax: 86-10-62561822. E-mail: xxia@ ipe.ac.cn. ORCID
Song Han: 0000-0002-4584-6136 Mo Zheng: 0000-0002-4467-5494 Xiaolong Liu: 0000-0003-4970-0406 Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS This work was supported by grants from the National Natural Science Foundation of China (Nos. 91641102, 21373227, and 91434105), from China’s State Key Laboratory of Multiphase Complex Systems (No. COM2015A004).
■
4. CONCLUSIONS Combustion simulations of a 24-component RP-1 model via ReaxFF MD were performed. A well-shaped spherical nanoparticle can be observed at the end of the 10 ns simulation. The underlying reaction mechanism of the morphological transformation of the whole system is revealed. Property evolutions of the largest molecule suggest that the soot formation process can be characterized with three stages. Stage 1 is for incipient ring formation that PAH-like molecules were initiated by ring closure of aliphatic polyyne-like radicals produced from fuel decomposition during the early stage. Stage 2 is for nucleation indicated as coalescence of two PAH-like molecules achieved by cross-link and ring closure of their side chains, and then the internal bonding leads to increase of ring number and formation of a PAH-like core. Aliphatic constituents play important roles as observed in the nucleation process. Stage 3 is for graphitization from PAH-like to fullerenic dominated by transformation from five-/seven-membered rings to six-
REFERENCES
(1) Richter, H.; Howard, J. B. Formation of polycyclic aromatic hydrocarbons and their growth to soota review of chemical reaction pathways. Prog. Energy Combust. Sci. 2000, 26, 565−608. (2) McEnally, C. S.; Pfefferle, L. D.; Atakan, B.; Kohse-Höinghaus, K. Studies of aromatic hydrocarbon formation mechanisms in flames: Progress towards closing the fuel gap. Prog. Energy Combust. Sci. 2006, 32, 247−294. (3) Tree, D. R.; Svensson, K. I. Soot processes in compression ignition engines. Prog. Energy Combust. Sci. 2007, 33, 272−309. (4) Bockhorn, H., Ed. Soot formation in combustion: mechanisms and models; Springer: Berlin, Germany, 1994. (5) Wang, H.; Frenklach, M. A detailed kinetic modeling study of aromatics formation in laminar premixed acetylene and ethylene flames. Combust. Flame 1997, 110, 173−221. (6) Appel, J.; Bockhorn, H.; Frenklach, M. Kinetic modeling of soot formation with detailed chemistry and physics: laminar premixed flames of C2 hydrocarbons. Combust. Flame 2000, 121, 122−136.
8442
DOI: 10.1021/acs.energyfuels.7b01194 Energy Fuels 2017, 31, 8434−8444
Article
Energy & Fuels (7) Richter, H.; Granata, S.; Green, W. H.; Howard, J. B. Detailed modeling of PAH and soot formation in a laminar premixed benzene/ oxygen/argon low-pressure flame. Proc. Combust. Inst. 2005, 30, 1397− 1405. (8) Raj, A.; Prada, I. D. C.; Amer, A. A.; Chung, S. H. A reaction mechanism for gasoline surrogate fuels for large polycyclic aromatic hydrocarbons. Combust. Flame 2012, 159, 500−515. (9) Akhter, M. S.; Chughtai, A. R.; Smith, D. M. Aromaticity of elemental carbon (soot) by 13C CP/MAS and FT-IR spectroscopy. Carbon 1985, 23, 593−594. (10) Cain, J. P.; Camacho, J.; Phares, D. J.; Wang, H.; Laskin, A. Evidence of aliphatics in nascent soot particles in premixed ethylene flames. Proc. Combust. Inst. 2011, 33, 533−540. (11) Wang, H. Formation of nascent soot and other condensedphase materials in flames. Proc. Combust. Inst. 2011, 33, 41−67. (12) Santamaria, A.; Yang, N.; Eddings, E.; Mondragon, F. Chemical and morphological characterization of soot and soot precursors generated in an inverse diffusion flame with aromatic and aliphatic fuels. Combust. Flame 2010, 157, 33−42. (13) Wang, L.; Song, C.; Song, J.; Lv, G.; Pang, H.; Zhang, W. Aliphatic C−H and oxygenated surface functional groups of diesel incylinder soot: Characterizations and impact on soot oxidation behavior. Proc. Combust. Inst. 2013, 34, 3099−3106. (14) Krestinin, A. V. Polyyne model of soot formation process. Symp. (Int.) Combust., [Proc.] 1998, 27, 1557−1563. (15) Chung, S. H.; Violi, A. Peri-condensed aromatics with aliphatic chains as key intermediates for the nucleation of aromatic hydrocarbons. Proc. Combust. Inst. 2011, 33, 693−700. (16) Elvati, P.; Violi, A. Thermodynamics of poly-aromatic hydrocarbon clustering and the effects of substituted aliphatic chains. Proc. Combust. Inst. 2013, 34, 1837−1843. (17) Saffaripour, M.; Veshkini, A.; Kholghy, M.; Thomson, M. J. Experimental investigation and detailed modeling of soot aggregate formation and size distribution in laminar coflow diffusion flames of Jet A-1, a synthetic kerosene, and n-decane. Combust. Flame 2014, 161, 848−863. (18) Commodo, M.; Tessitore, G.; De Falco, G.; Bruno, A.; Minutolo, P.; D’Anna, A. Further details on particle inception and growth in premixed flames. Proc. Combust. Inst. 2015, 35, 1795−1802. (19) Karataş, A. E.; Gülder, Ö . L. Soot formation in high pressure laminar diffusion flames. Prog. Energy Combust. Sci. 2012, 38, 818−845. (20) Desgroux, P.; Mercier, X.; Thomson, K. A. Study of the formation of soot and its precursors in flames using optical diagnostics. Proc. Combust. Inst. 2013, 34, 1713−1738. (21) van Duin, A. C. T; Dasgupta, S.; Lorant, F.; Goddard, W. A. ReaxFF: A Reactive Force Field for Hydrocarbons. J. Phys. Chem. A 2001, 105, 9396−9409. (22) Senftle, T. P.; Hong, S.; Islam, M. M.; Kylasa, S. B.; Zheng, Y.; Shin, Y. K.; Junkermeier, C.; Engel-Herbert, R.; Janik, M. J.; Aktulga, H. M.; Verstraelen, T.; Grama, A.; van Duin, A. C. T. The ReaxFF reactive force-field: development, applications and future directions. Npj Comput. Mater. 2016, 2, 15011. (23) Russo, M. F., Jr; van Duin, A. C. T. Atomistic-scale simulations of chemical reactions: Bridging from quantum chemistry to engineering. Nucl. Instrum. Methods Phys. Res., Sect. B 2011, 269, 1549−1554. (24) Castro-Marcano, F.; Kamat, A. M.; Russo, M. F., Jr; 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. (25) 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. (26) Mortier, W. J.; Ghosh, S. K.; Shankar, S. Electronegativityequalization method for the calculation of atomic charges in molecules. J. Am. Chem. Soc. 1986, 108, 4315−4320. (27) Smirnov, K. S.; van de Graaf, B. Consistent implementation of the electronegativity equalization method in molecular mechanics and molecular dynamics. J. Chem. Soc., Faraday Trans. 1996, 92, 2469− 2474.
(28) Chenoweth, K.; van Duin, A. C. T.; Goddard, W. A. ReaxFF Reactive Force Field for Molecular Dynamics Simulations of Hydrocarbon Oxidation. J. Phys. Chem. A 2008, 112, 1040−1053. (29) 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. (30) Zheng, M.; Li, X.; Guo, L. Algorithms of GPU-enabled reactive force field (ReaxFF) molecular dynamics. J. Mol. Graphics Modell. 2013, 41, 1−11. (31) Liu, J.; Li, X.; Guo, L.; Zheng, M.; Han, J.; Yuan, X.; Nie, F.; Liu, X. Reaction analysis and visualization of ReaxFF molecular dynamics simulations. J. Mol. Graphics Modell. 2014, 53, 13−22. (32) LAMMPS. Available via the Internet at http://lammps.sandia. gov/, accessed May 11, 2016. (33) Li, X.; Mo, Z.; Liu, J.; Guo, L. Revealing chemical reactions of coal pyrolysis with GPU-enabled ReaxFF molecular dynamics and cheminformatics analysis. Mol. Simul. 2015, 41, 13−27. (34) Pitz, W. J.; Mueller, C. J. Recent progress in the development of diesel surrogate fuels. Prog. Energy Combust. Sci. 2011, 37, 330−350. (35) Bruno, T. J.; Smith, B. L. Improvements in the Measurement of Distillation Curves. 2. Application to Aerospace/Aviation Fuels RP-1 and S-8. Ind. Eng. Chem. Res. 2006, 45, 4381−4388. (36) Wang, T. Thermophysics Characterization of Kerosene Combustion. J. Thermophys. Heat Transfer 2001, 15, 140−147. (37) Accelrys.Materials-Studio. Available via the Internet at http:// accelrys.com/products/collaborative-science/biovia-materials-studio/, accessed May 11, 2016. (38) Mayo, S. L.; Olafson, B. D.; Goddard, W. A. DREIDING: a generic force field for molecular simulations. J. Phys. Chem. 1990, 94, 8897−8909. (39) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984, 81, 3684−3690. (40) Yang, W.; Sun, B. Numerical Simulation of Liquid Film in a Liquid Oxygen/Rocket Propellant 1 Liquid Rocket. J. Thermophys. Heat Transfer 2012, 26, 328−336. (41) Weismiller, M. R.; Russo, M. F., Jr; van Duin, A. C. T.; Yetter, R. A. Using molecular dynamics simulations with a ReaxFF reactive force field to develop a kinetic mechanism for ammonia borane oxidation. Proc. Combust. Inst. 2013, 34, 3489−3497. (42) Eaves, N. A.; Dworkin, S. B.; Thomson, M. J. The importance of reversibility in modeling soot nucleation and condensation processes. Proc. Combust. Inst. 2015, 35, 1787−1794. (43) Russo, C.; Tregrossi, A.; Ciajolo, A. Dehydrogenation and growth of soot in premixed flames. Proc. Combust. Inst. 2015, 35, 1803−1809. (44) Haynes, B. S.; Wagner, H.Gg. Soot formation. Prog. Energy Combust. Sci. 1981, 7, 229−273. (45) Vlasov, P. A.; Warnatz, J. Detailed kinetic modeling of soot formation in hydrocarbon pyrolysis behind shock waves. Proc. Combust. Inst. 2002, 29, 2335−2341. (46) Grieco, W. J.; Howard, J. B.; Rainey, L. C.; Vander Sande, J. B. Fullerenic carbon in combustion-generated soot. Carbon 2000, 38, 597−614. (47) Sirignano, M.; Alfè, M.; Tregrossi, A.; Ciajolo, A.; D’Anna, A. Experimental and modeling study on the molecular weight distribution and properties of carbon particles in premixed sooting flames. Proc. Combust. Inst. 2011, 33, 633−640. (48) Kholghy, M. R.; Veshkini, A.; Thomson, M. J. The core−shell internal nanostructure of soot − A criterion to model soot maturity. Carbon 2016, 100, 508−536. (49) Rodgers, R. P.; Reilly, P. T. A.; Whitten, W. B.; Ramsey, J. M. Soot-free synthesis of C60. Carbon 2003, 41, 687−692. (50) Vander Wal, R. L.; Choi, M. Y. Pulsed laser heating of soot: morphological changes. Carbon 1999, 37, 231−239. (51) Frenklach, M. Reaction mechanism of soot formation in flames. Phys. Chem. Chem. Phys. 2002, 4, 2028−2037. 8443
DOI: 10.1021/acs.energyfuels.7b01194 Energy Fuels 2017, 31, 8434−8444
Article
Energy & Fuels (52) Homann, K. H. Fullerenes and Soot Formation New Pathways to Large Particles in Flames. Angew. Chem., Int. Ed. 1998, 37, 2434−2451. (53) Fredrik Ahlström, A.; Ingemar Odenbrand, C. U. Combustion characteristics of soot deposits from diesel engines. Carbon 1989, 27, 475−483. (54) Karataş, A. E.; Gülder, Ö .L. Soot formation in high pressure laminar diffusion flames. Prog. Energy Combust. Sci. 2012, 38, 818−845.
8444
DOI: 10.1021/acs.energyfuels.7b01194 Energy Fuels 2017, 31, 8434−8444