Comparison of ReaxFF, DFTB, and DFT for ... - ACS Publications

Article pubs.acs.org/JPCA

Comparison of ReaxFF, DFTB, and DFT for Phenolic Pyrolysis. 1. Molecular Dynamics Simulations Tingting Qi,∥ Charles W. Bauschlicher, Jr.,† John W. Lawson,*,‡ Tapan G. Desai,§ and Evan J. Reed*,∥ ∥

Department of Materials Science and Engineering, Stanford University, Stanford, California 94305, United States Mail Stop 230, Entry Systems and Technology Division, NASA Ames Research Center, Moffett Field, California 94035, United States ‡ Mail Stop 234, Thermal Protection Materials Branch, NASA Ames Research Center, Moffett Field, California 94035, United States § Advanced Cooling Technologies, Inc., Lancaster, Pennsylvania 17601, United States †

ABSTRACT: A systematic comparison of atomistic modeling methods including density functional theory (DFT), the selfconsistent charge density-functional tight-binding (SCC-DFTB), and ReaxFF is presented for simulating the initial stages of phenolic polymer pyrolysis. A phenolic polymer system is simulated for several hundred picoseconds within a temperature range of 2500 to 3500 K. The time evolution of major pyrolysis products including small-molecule species and char is examined. Two temperature zones are observed which demark cross-linking versus fragmentation. The dominant chemical products for all methods are similar, but the yields for each product differ. At 3500 K, DFTB overestimates CO production (300−400%) and underestimates free H (∼30%) and small CmHnO molecules (∼70%) compared with DFT. At 3500 K, ReaxFF underestimates free H (∼60%) and fused carbon rings (∼70%) relative to DFT. Heterocyclic oxygencontaining five- and six-membered carbon rings are observed at 2500 K. Formation mechanisms for H2O, CO, and char are discussed. Additional calculations using a semiclassical method for incorporating quantum nuclear energies of molecules were also performed. These results suggest that chemical equilibrium can be affected by quantum nuclear effects at temperatures of 2500 K and below. Pyrolysis reaction mechanisms and energetics are examined in detail in a companion manuscript.

I. INTRODUCTION Pyrolysis is the thermochemical decomposition of organic materials at elevated temperatures in the absence of oxygen. This fundamental chemical process is important for ablative heat shields for space vehicles which experience extreme, high temperature, chemically reactive environments during atmospheric entry.1 Phenolic resins are a primary constituent of ablative thermal protection systems. During atmospheric entry, endothermic pyrolytic reactions of the phenolic resin absorb heat, produce gaseous products that block further heating, and produce a solid carbon residue (char) on the surface. The highly emissive char layer efficiently radiates heat and thermally insulates the interior of the vehicle. Improving the performance of ablative materials will require understanding the gaseous product distribution, increasing the char yield, and producing char with improved properties, such as superior mechanical strength to prevent spallation of the surface char layers. For these purposes, it is important to obtain a fundamental understanding of phenolic pyrolysis reaction mechanisms, the process of char formation, and the resulting structural, mechanical, and thermal properties of char.2−6 Atomic level information about the complex series of reactions involved in phenolic pyrolysis is limited. This is due in part to the lack of detailed experimental data for these complicated systems. © 2013 American Chemical Society

The principal source of experimental information is thermogravimetric analysis (TGA) which can study the kinetics of solidstate reactions4 by tracking the mass loss of a pyrolyzing system as a function of temperature. Trick et al., for example, analyzed four TGA (300−800 °C) mass-loss peaks for phenolic/carbon pyrolysis and estimated the activation energies to be in the range of 17−47 kcal/mol4 with average pre-exponential factors, A, having values from 1.1 × 10−11 ps−1 to 6.7 × 10−6 ps−1. However, there is not a clear relationship between TGA curves and detailed decomposition mechanisms. Other methods to study pyrolysis reactions are gas chromatography and mass spectrometry (GC/MS) of the volatile small-molecule pyrolysis products.1 These techniques yield important information about reaction products, but the actual reaction pathways are not obtained. Molecular dynamics (MD) simulations offer another route toward understanding the detailed chemistry of phenolic pyrolysis. Recently, the reactive force field ReaxFF7 has been used for MD simulations of phenolic pyrolysis.6,8−10 Jiang et al. simulated the initial stage of pyrolysis.8 They found that H2O is the first reaction Received: August 13, 2013 Revised: October 3, 2013 Published: October 4, 2013 11115

dx.doi.org/10.1021/jp4081096 | J. Phys. Chem. A 2013, 117, 11115−11125

The Journal of Physical Chemistry A

Article

with DFT and DFTB in the context of pyrolysis, with attention to reaction details at the atomistic level.

product observed when temperatures reach approximately 2900 K. This is followed by H2 and C2H2 and finally by CO. Multiple dehydration channels were reported: β-hydrogen elimination, interchain dehydration, and dissociated hydrogen extraction. An activation energy for water production of 34 ± 7 kcal/mol was derived from their simulation data which is within the range inferred from the TGA experiments.4 Desai et al. also performed reactive MD simulations investigating the initial stage of pyrolysis.9 They found that water is the primary product in the temperature range of 2750 to 3250 K with an activation energy of 30 ± 5 kcal/mol, consistent with the results obtained by Jiang et al. Despite these interesting results, MD simulations suffer from two significant drawbacks for examining chemical processes such as phenolic pyrolysis. First, MD simulations typically are limited to relatively short time scales, usually on the order of nanoseconds, whereas realistic chemical processes may unfold on time scales many orders of magnitude longer. Thus, MD simulations can only consider the initial stages of a complex process such as pyrolysis. The time scale problem may have a solution by using accelerated molecular dynamics or kinetic Monte Carlo methods. Second, the fidelity of MD simulations depends critically on the quality of the interatomic potential used. Recently, reactive force fields such as ReaxFF have been developed that have the capability of describing processes involving the creation and breaking of chemical bonds. ReaxFF for example is a bond order based potential whose parameters have been fit to ab initio data. While ReaxFF can be approximately 2 orders of magnitude more costly than more simple force fields, for many applications, it produces highly accurate results at a fraction of the computational cost of ab initio methods. Nanosecond time scales are possible with ReaxFF. On the other hand, ReaxFF parameter sets are obtained by fitting to a limited set of reactions. Whether a given parameter set will perform well for reactions outside its training set must be verified explicitly. In this paper, we perform MD simulations for phenolic pyrolysis using three state-of-the-art methods to compare their accuracy for this application. In particular, we consider density functional theory (DFT), density functional tight binding (DFTB), and ReaxFF. Compared to ReaxFF, self-consistent charge (SCC) DFTB is, in principle, a more transferable method. DFTB is a tight-binding method where Hamiltonian integrals, overlap integrals, and pairwise repulsive nuclear interactions are fit to DFT calculations. The self-consistent charge variant is based on a second-order expansion of the Kohn−Sham total energy with respect to charge density and spin density fluctuations to enable some degree of charge transfer between atoms and representation of charged chemical species.11,12 Molecular dynamics simulations with a DFTB potential for a few hundred atoms are achievable on 1 ns time scales. These time scales are comparable to ReaxFF; however, fewer atoms can be simulated with DFTB due to the currently superior parallel implementations of ReaxFF. DFT is expected to provide the highest accuracy of the methods employed for MD in this work since the other methods are fit to DFT calculations. DFT computes the most accurate potential energy surface, so we expect the most accurate kinetics and thermodynamics to be extracted from the DFT calculations. However, DFT is considerably more expensive computationally. Clearly, only methods with computational cost comparable to or less than ReaxFF will be appropriate for simulating macroscopic polymers with many thousands of atoms on nanosecond or longer time scales. Our goal is to compare the performance of ReaxFF

II. METHOD We use the parameter set of Chenoweth et al.13 and the MD package LAMMPS14 for our ReaxFF MD simulations. The Chenoweth parameter set was originally developed for hydrocarbon combustion. We additionally employ the energy penalty EC2 described in ref 15 to prevent a large unphysical binding energy for molecular C2. As in the original Chenoweth paper,13 an extra triple bond energy correction for stabilizing CO is also activated. We use the C implementation of ReaxFF in LAMMPS. The DFTB calculations are performed with the selfconsistent charge density functional tight binding method (SCC-DFTB)12 using the DFTB+ code package16 and both the mio-1-1 and pbc-0-3 parameter sets. The mio set was developed for organic molecules with atoms including O, N, C, and H.12 The pbc set is developed with a focus on solid state systems.17 DFT trajectories are generated using the planewave code VASP18,19 and the PW91 exchange-correlation functional.20 An energy cutoff of 400 eV is used, and the projector augmentedwave (PAW) potentials developed for VASP are applied.21 The k-space sampling is restricted to the Γ point only. Our computational cell for all three energy methods contains 132 atoms where periodic boundary conditions have been used as shown in Figure 1. We choose this size to facilitate a

Figure 1. The noncross-linked phenol formaldehyde resin molecule (a) used in the initial 132-atom periodic computational cell (b), where carbon, oxygen, and hydrogen atoms are displayed in black, red, and silver, respectively.

computationally tractable direct comparison between methods. Our computational system consists of four phenol formaldehyde chains with two repeating units each in the ortho−ortho sequence and a methyl group termination. This molecule has a stoichiometry of 11116



Article

C15H16O2. The simulation box size is set to 1.065 × 1.065 × 1.065 nm3 yielding a typical experimental room-temperature resin density of 1.25 g/cm3. We investigate pyrolysis at temperatures of 2500 K, 3000 K, and 3500 K. These temperatures are chosen higher than those of typical pyrolysis processes so that the simulations exhibit chemical decomposition on the time scales accessible with these MD methods (hundreds of ps to 1 ns). We use enhanced temperatures to increase the kinetics of these reactions. The use of increased temperatures to accelerate kinetics can alter the relative weights of reaction channels and also alter the chemical equilibrium toward which the system evolves. We employ high temperatures to enable a systematic comparison between energy models for conditions that approach those of lower temperature pyrolysis. Activation barriers for reaction channels applicable to lower temperature pyrolysis are identified in our companion manuscript. Both DFTB-MD and ReaxFF-MD simulations are performed using a computer cluster comprised of AMD Opteron processors with a clock speed of 2.3 GHz. On this machine, it typically takes about 2 h (1 core) and 7 min (8 cores) to propagate 1 ps of simulation time for DFTB and ReaxFF, respectively. All DFT trajectories are computed on 32 Intel Xeon 3 GHz processors, requiring 8.5 h per picosecond of simulation time. The system is initially equilibrated at a constant temperature of 2000 K and using a constant volume NVT canonical ensemble for 5 ps. A Nosé-Hoover thermostat was employed along with a 0.25 fs time step. No obvious reactions are observed during the 5 ps equilibration with any of the methods considered here. The system is then heated from 2000 K to the target temperature within 1 ps using a constant rate T increase. After the temperature ramp, we perform NVT simulations for 500 ps for the ReaxFF and DFTB, or 100 ps for the DFT simulations. To identify molecular species during the simulations, we use the following maximum bond length criteria: O−H = 1.2 Å, C−O = 1.85 Å, H−H = 0.9 Å, C−H = 1.5 Å, O−O = 1.6 Å, and C−C = 1.85 Å. We have set these bond length cutoffs at the first minimum of the corresponding pair radial distribution function.22 In addition to the length criterion, atoms must be within the specified distance for 50 or more femtoseconds to be considered bonded. In addition to MD trajectories, we also compute certain bond energies to compare different computational methods. In addition to the ReaxFF, DFTB, and DFT approaches for these bond energies, we employ the coupled cluster singles and double approach CCSD(T) including the effect of connected triples determined using perturbation theory. The spin unrestricted approach is used for open shell systems, and the carbon and oxygen 1s electrons are not correlated. The correlation-consistent polarized valence triple-ζ (cc-pVTZ) set23 is applied in CCSD(T) calculations.

Table 1. C6H5OH Bond Energies Calculated Using DFT, CCSD(T), DFTB (mio-1-1 and pbc-0-3 Parameter Sets), and ReaxFFa

DFTB bond

DFT

CCSD(T) cc-pVTZ

O−H C−OH C−H(1) C−H(2) C−H(3) C−H(4) C−H(5)

91.1 124.0 116.7 116.4 117.8 116.5 118.4

95.3 115.9 121.0 120.3 121.6 120.3 122.5

mio

pbc

ReaxFF

125.4 148.4 150.8 150.7 151.9 151.1 151.7

124.9 154.0 145.3 145.3 146.5 145.6 146.2

110.3 115.4 156.5 142.7 145.8 143.5 158.2

a Energy units are kcal/mol. C−H bonds are labeled as in the figure. The color codes for atoms are the same as in Figure 1.

of C−OH, O−H, and various C−H bonds are displayed in Table 1. The DFT values are in closest agreement with the CCSD(T) values. DFTB and ReaxFF overestimate the CCSD(T) numbers but have similar energetic trends: the O−H bond is the weakest with all others showing similar bond energies. ReaxFF overestimates O−H, C−H, and C−C bond energies but yields a C− OH bond energy that agrees well with the CCSD(T) result. The bond energy difference between C−OH and O−H is relatively small within ReaxFF (5 kcal/mol), compared with the other energy models (20−30 kcal/mol). In our 3500 K trajectories, we found that the ratio of OH vs H population is roughly 10% for ReaxFF and 2.5% for DFT, which is consistent with what one might expect based on the bond energies in Table 1. This suggests that the OH population is overestimated within ReaxFF. Since H and OH are both typical polymer-reaction radicals, overestimation of OH could lead to artificially enhanced decomposition of phenolic resins. Several possible C−C bond breaking reactions are depicted in Figure 2, and energies for the various methods are given in Table 2. The right side of Figure 2 shows energy minimized structures using DFT. All calculations in Table 2 are the energy difference between products and reactant of five phenolic ring breaking reactions. Table 2 shows considerable variation in the C−C bond energy when hydrogen transfer occurs. For example, the blue box in Figure 2 shows the same C−C bond broken in reaction ii and reaction iii, but in reaction iii the H in the neighboring −OH group transfers to the H-free carbon with C−C bond breaking. The H transfer in reaction iii is found to significantly lower the bond energy by over 66 kcal/mol for all energy methods. A similar effect is found in comparing reaction iv and reaction v. For all five C−C bond breaking reactions, ReaxFF systematically overestimates the bond energy. If the energy barrier for these processes correlates positively with the bond energies, this suggests that breaking six-membered carbon rings may be less likely in ReaxFF than the other methods. This possibility is reflected in the population of six-membered carbon ring structures observed in our MD simulations, as shown in Figure 3. At 2500 K, ReaxFF exhibits slower decomposition of

III. RESULTS AND DISCUSSION A. Bond Energies. The pyrolysis process consists of breaking and forming chemical bonds. Therefore, we begin our method comparison by computing bond energies for some representative species. We choose the isolated (gas phase) molecules C6H5OH, C13H12O2, and C13H13O2, which permit us to compare ReaxFF, DFTB, and DFT to the more accurate CCSD(T) approach. Although the surrounding chemical environment in the pyrolyzing phenolic resin might play a role in the bond energies, the comparison of gas phase molecules can be useful for assessing the accuracy of these methods. The energies 11117



Article

five-membered and six-membered ring populations at temperatures of 3000 and 3500 K relative to the other energy methods. This may point to differences in the structure of the resulting char that may impact material properties as well. Small fused-ring patches serve as seeds for the formation of large graphitic structures (e.g., graphene sheets and nanotubes), which could play a substantial role in determining thermalprotection properties. Figure 4 shows the number of carbon atoms in each bonded cluster as a function of time, providing complementary information to the ring populations in Figure 3. Figure 4 reveals whether carbon atoms are clustering (rings, polymers, etc.) or fragmenting into small molecular products. Each of these simulations starts with a peak at 15, associated with the C15H16O2 molecules used for the initial configuration. Fragments of 7 and 8 are also observed early in the simulations and are associated with the transient breaking of a C15H16O2 into two separate fragments. All three of the DFT, DFTB, and ReaxFF simulations exhibit two distinct types of carbon chemistry that have been reported previously.9 Polymerization occurs at lower temperatures (2500 and 3000 K), forming fragments with larger numbers of carbon atoms. At higher temperatures (3500 K), the polymers dissociate forming smaller polymer fragments/molecules. The ReaxFF trajectory at 2500 K appears to be less reactive than DFT and DFTB, consistent with the bond energies in Table 2 suggesting that ring breaking may require more energy with ReaxFF. The formation of small carbon clusters at 3500 K is qualitatively different from the lower temperature pyrolysis process that occurs during high speed atmospheric flight, indicating that temperatures of 3500 K may have more impact than only to accelerate the pyrolysis process. The O atoms contained in the phenolic resin generally form small molecules, such as H2O, CO, and CO2. Interestingly, we also find O atoms contained in both five-membered (four carbon members plus one oxygen member) and six-membered (five carbon members plus one oxygen member) ring structures. Such heterocyclic oxygen-containing ring structures are stable for several hundred femtoseconds or several picoseconds within DFT-MD. They are more often found at the lower temperature of 2500 K. Some of these heterocyclic oxygencontaining ring structures can be seen in Figure 5, where the snapshots of the largest fragment at the end of simulations are displayed. Because oxygen atoms form fewer covalent bonds than carbon atoms, oxygen-containing ring structures may prevent the continued growth of fused rings around the oxygen site. Some oxygen-containing structures in graphene oxide have been reported to exhibit high thermodynamic stability, hindering the production of pristine graphene during annealing.25 One might expect these oxygen atoms to affect the growth rate and the resulting structure of the char, having potential impact on the mechanical and thermal properties of the char. C. Small Molecular Product Evolution. During the pyrolytic decomposition of phenolic resin, a variety of smaller products are produced. Figure 6 presents the time evolution of some of the principal pyrolysis products: H, C2H2, H2O, and CO. Each data point is averaged over 5 ps. In total, there are sixty-four hydrogen atoms and eight oxygen atoms in our computational cell. For all computational methods we considered, the first reaction we observe is always an H atom loss from the C15H16O2 reactant molecules. This suggests that the loss of an H atom has one of the lower, if not the lowest, reaction barriers. H is also found to be the most abundant product species.

Figure 2. Table 2 lists the C−C bond energies for the molecules C13H12O2 and C13H13O2 shown in this figure. We evaluate bond energies calculated by applying DFT, DFTB (mio-1-1 and pbc-0-3 parameter sets) and ReaxFF to compare the methods for chemical reactions that might occur during pyrolysis. Five different C−C bond breaking reactions are labeled (i, ii, iii, iv, and v) as in the figure. Products in the same colored box are derived from breaking the same C−C bond.

Table 2. C−C Bond Energies for the Isolated Molecules C13H12O2 and C13H13O2 Shown in Figure 2 Using DFT, DFTB (mio-1-1 and pbc-0-3 Parameter Sets), and ReaxFFa DFTB reaction label

DFT

mio

pbc

ReaxFF

i ii iii iv v

90.6 141.8 47.9 143.9 57.0

114.1 164.0 57.3 173.1 72.9

111.8 162.1 49.2 170.9 70.6

117.6 166.0 99.3 171.3 75.7

a

Five different C−C bond breaking reactions that might occur in pyrolysis are labeled (i, ii, iii, iv, and v) as in Figure 2. The energy unit is kcal/mol. DFT energies are lowest, suggesting the least stable carbon rings. A difference in C−C bond energy (e.g. between ii and iii) is observed when a nearby hydrogen atom is allowed to transfer in addition to C−C bond breakage.

the six-membered rings than the other energy methods. DFT sixmembered rings show the most reactivity at 2500 K in Figure 3, consistent with the bond energy trends in Table 2. B. Char Formation. Char is produced from phenolic resin after the loss of hydrogen- and oxygen-containing compounds. Over the durations of our simulations, we observe many fusedring fragments consisting of mainly five- and six-membered rings. Figure 3 shows the carbon ring statistics for fivemembered, six-membered, and all carbon rings (summation of carbon rings containing eight and fewer carbon members) at different temperatures. Ring statistics are characterized using a shortest-path criterion.24 All counted rings are primitive and not composed of smaller rings. Each data point in Figure 3 is averaged over 5 ps. All energy methods exhibit qualitatively similar behavior, but Figure 3 shows that ReaxFF appears to underestimate the 11118



Article

Figure 3. The number of five- and six-membered carbon rings and the total number of carbon rings (summation of three-, four-, five-, six-, seven-, and eight-membered) as a function of time at temperatures of 2500 K, 3000 K, and 3500 K. ReaxFF underestimates ring populations at high temperatures relative to the other energy methods. ReaxFF also exhibits slower ring breaking kinetics at 2500 K.

Figure 4. The fraction of C atoms contained in N-carbon clusters, where N is shown on the y axis. The total number of carbon atoms in the simulation is 60. Fractions equal to or greater than 0.3 are plotted with the same color. At lower temperatures, an evolution to larger carbon clusters or polymers is observed, while higher temperatures are marked by fragmenting of carbon structures into small molecules for all energy methods. Note that the time scale on the DFT plots is shorter than those of the other energy methods.

formation energies in Table 3 are in the closest agreement with experimental measurements. To better understand the complex pyrolysis reactions, we trace the major consumption channels of hydrogen and oxygen, as shown in Figure 7. Each data point is averaged over 5 ps. The top three rows of Figure 7 show that hydrogen atoms tend to go into small-molecule hydrocarbons CmHn (m ≤ 5), H, OH, H2, and H2O. We notice however that ReaxFF does not

Of the four chemical species we present, the most significant difference between methods is the overestimation of CO by DFTB. To better understand this disagreement, we computed the formation energies for these species (Table 3) and found a positive correlation between the value of the formation energies and the size of the populations. Compared with experimental data, the DFTB/mio and DFTB/pbc parameter sets bias the CO formation by −22 kcal/mol and −13 kcal/mol in the favored direction. DFT 11119



Article

Figure 5. Snapshots of the largest fragment in the computational cell at the end of the simulations. At the lower temperature of 2500 K, heterocyclic oxygen-containing ring structures are stable for several hundred femtoseconds up to several picoseconds. Slightly more fused ring structures are observed at 3000 K than at lower (2500 K) and higher (3500 K) temperatures.

breaks from the phenolic chain and then forms water by combining with a free H atom. Free H generation occurs first, and therefore it will not be the rate-determining step. Additionally, the combination step H + OH → H2O does not have a reaction barrier at all. As a result, the rate of this particular water formation mechanism is determined by the H and OH populations. We find that the OH dissociation does not involve any transition state, making the corresponding reaction barrier height equivalent to the bond energy of breaking the OH group from the phenolic ring. Unlike H2O, CO formation mechanisms during pyrolysis of phenolics have not been previously reported. Even with our minimal 132-atom computational cell, a detailed analysis of full reaction path is not possible. An analysis of the DFT-MD trajectories shows that the carbon comes from the benzene rings while the oxygen originates from the OH side group. This suggests that CO formation requires a phenolic ring to open first. We are only able to identify the final steps of CO formation as in Figure 8 (e) and (f): the carbonyl group directly detaches (called “one-step”), or an aldehyde group detaches first and then splits into one H atom and one CO molecule (called “two-step”). The companion manuscript contains additional detail regarding possible reaction pathways and their energy barriers. We have attempted to infer some approximate kinetics for the initial stage of pyrolysis from our simulations. These results may depend on the system size that we employ and the associated statistical noise, which can be large. If we assume all reactions are first-order, with a single step from reactant A to product P (A → P), the population of P can be expressed as [P] = [A](1-e−kt), where

produce small hydrocarbons at 2500 K. When the temperature is elevated to 3000 K, the quantity of small hydrocarbons produced by ReaxFF is significantly lower than the amount produced by the other methods. The oxygen atoms are usually found in CmHnO (m ≤ 5), CO, CO2, OH, H2O, and O. At all three temperatures, the overestimation of CO by DFTB leads to the underestimation of oxygen contained in CmHnO in DFTB calculations. ReaxFF shows better agreement with DFT in predicting oxygen consumption channels. D. Reaction Mechanisms and Barriers. The companion paper to this manuscript focuses in detail on reaction mechanisms and barriers for phenolic pyrolysis, with some additional detail from the MD simulations presented here. In this manuscript, we provide a summary description and refer the reader to the companion manuscript for additional detail. Jiang et al.8 have reported three dehydration channels which are shown in Figure 8 (a)-(c): β-hydrogen elimination, the interchain dehydration, and dissociated hydrogen extraction. All mechanisms involve an intermediate step where a hydroxyl group combines with a nearby H. The H could be a free, dissociated atom, or it could come from a functional group (CH2 group, benzene ring, or OH group). All mechanisms share a similar transition state, which resembles a loosely attached H2O where one of the O−H bond lengths is normal, and the other is slightly stretched. The C−OH bond becomes weakened and therefore also appears longer than normal. From our DFT trajectories, we identified an additional fourth dehydration process, also shown in Figure 8 (d). This process is only observed at higher temperatures: a hydroxyl group 11120



Article

Figure 6. The numbers of major pyrolysis products H, C2H2, H2O, and CO as a function of time at temperatures of 2500 K, 3000 K, and 3500 K. In total, there are sixty-four hydrogen atoms and eight oxygen atoms in our computational cell. The dashed lines are fits to the function [P] = [A](1-e−kt) used to determine the reaction barriers in Table 4. DFTB gives larger populations of CO than the other energy methods. ReaxFF gives lower free H populations than the other methods at high temperatures.

molecules are considered to be “char”. In Figure 9, we show the percentage of all carbon atoms contained in so-defined char. For all computational methods, we obtain char formation barrier values which range from 43 to 56 kcal/mol shown in Table 4. While the primary value of Table 4 is to compare the methods, it is possible to make some comparisons to other values in the literature with the caveat that the temperatures and phenolic polymer lengths considered in this work are higher and shorter, respectively. A crude comparison can be made with the analysis by Trick et al. of four TGA (300−800 °C) mass-loss peaks for phenolic/carbon pyrolysis and estimated the activation energies to be in the range of 17 to 47 kcal/mol4 with average pre-exponential factors, A, having values from 1.1 × 10−11 ps−1 to 6.7 × 10−6 ps−1. The higher temperatures and noncross-linked polymer model employed in this work may play roles in differences observed from the values reported by Trick et al. We speculate that many reactions could be multistep and therefore not well described by the single-step rate fits described above. Multistep reactions may exhibit induction periods before the generation of some products. Possible induction periods can be seen in some of our CO and H2O time-dependent population plots at 2500 K. By allowing an extra parameter describing the induction period in fitting [P] = [Ai](1-e−ki(t−ti)), we find some resulting energy barriers are lower than without the induction parameter. The energy barriers of CO and H2O are lowered by 1−5%, and the energy barrier for char production is lowered by 20%. The precomponent factor [Ai] of H2O production is reduced by less than 5% from the single-step [A] and is still unphysical. E. Quantum Nuclear Effects. The molecular dynamics simulations employed in this work are expected to be appropriate when the temperature is significantly above the highest vibration

Table 3. Formation Energies for Several Pyrolysis Products (DFT Results Are without Zero-Point Energy) Calculated Relative to the Compound’s Elements in Their Standard Statesa DFTB chemical

expt

DFT

mio

pbc

ReaxFF

H C2H2 H2O CO

52 54 −58 −26

52 51 −74 −29

62 25 −103 −48

62 39 −93 −39

55 12 −66 −24

a

In units of kcal/mol. The experimental standard formation enthalpy values are collected from NIST Chemistry WebBook. The relatively low, negative DFTB CO formation energy correlates with the relatively high populations observed in the DFTB MD in Figure 6.

[A] is the final population of product P, k is the first-order reaction rate, and t is time. By performing a least-squares fit to this function, as shown in Figure 6, we can extract reaction rates at different temperatures for each chemical species. With an assumption of Arrhenius kinetics, we fit the slope of ln(k) versus −1/T to determine Eb/kB, where T is the temperature and kB is the Boltzmann constant. All computed reaction barriers are listed in Table 4. We find larger error bars for reaction barriers obtained from DFT-MD which is due in part to the short simulation runs used with this method. It is important to note that all error bars listed here are determined from a single trajectory for each method; other error sources such as system size are not reflected and may be comparable or larger. Values omitted (including all fits for H2O) exhibited unphysical prefactors >1015 s−1, suggesting the presence of a large error in these fits. Table 4 shows that the DFTB CO barriers are comparable. Molecules with carbon content higher than the initial phenolic 11121



Article

Figure 7. The fractions of H and O contained in different pyrolysis products as a function of time at temperatures of 2500 K, 3000 K, and 3500 K. The top three rows show H channels and the bottom two rows show O channels. At 2500 and 3000 K, the quantity of small hydrocarbons CmHn (m ≤ 5) produced by ReaxFF is significantly lower than the amount produced by the other methods. ReaxFF is generally in better agreement with DFT on the oxygen production channels due to the enhanced CO production in DFTB. DFT and ReaxFF show more oxygen contained in small hydrocarbons, CmHnO (m ≤ 4).

Figure 8. Several reaction mechanisms of H2O and CO formation, as observed in the DFT trajectories. See the companion paper for additional details regarding reaction pathways for pyrolysis products.

approximation to vibrational energies on the fly during the course of the simulation. Langevin-type (classical) equations of motion are applied to the atoms to force the vibrational energies to obey a Bose-Einstein spectrum associated with a quantum harmonic oscillator. Instead of computing the random forces from a white-noise spectrum employed in the usual Langevin thermostat, the random force obeys a "colored” Bose-Einstein spectral distribution. This semiclassical approach has the potential to describe energetic contributions to pyrolysis product free-energies that are associated with quantum nuclear effects, e.g. zero point energies. Shifts in product energies can potentially change the chemical

frequencies in the system, where quantum nuclear effects can be neglected. C−H bonds have stretch frequencies in the vicinity of 3000 cm−1 with quantum ground state energies corresponding to a temperature T = ℏω/2kB of 2200 K. Therefore, these high frequency stretching modes are primarily in the ground state even at the high temperatures considered in this work, and one might expect that quantum nuclear effects could play some role, particularly for the 2500 K case. To explore the potential impact of quantum nuclear effects in phenolic pyrolysis, we adopt a semiclassical approach using a so-called colored thermostat for quantum nuclear effects by Dammak et al.26 This approach applies a quasi-harmonic 11122



Article

quantum and classical cases exhibit some differences at 2500 K, while the 3000 and 3500 K cases appear to be similar. Figures 10 and 11 show the carbon ring statistics and carbon cluster sizes, indicating increased reactivity and polymerization at 2500 K in the quantum case. Nearly half of the original six-membered phenolic carbon rings are observed to break, forming either larger carbon clusters or small molecule products in the quantum case. Figure 12 shows H and H2O production is enhanced in the quantum case at 2500 K. The water yield is more than doubled at the end of 500 ps of simulation time. The comparison of Figure 12 to Figure 6 is particularly noteworthy and indicates that the quantum bath results in larger changes to the H and H2O populations at 2500 K than the variations between the three energy methods considered in this paper. This suggests that while the energy methods may differ, the quantum nature of the atoms may play as important or perhaps a more important role at temperatures of 2500 K and below which may approach realistic pyrolysis temperatures. The semiclassical thermostat that we have employed here is subject to the usual caveats of semiclassical approaches, and further study with more accurate quantum methods is warranted.

Table 4. Reaction Barriers for CO and Char As Obtained by Fitting an Arrhenius Law to MD Simulation Data in Figure 9 and Figure 6a DFTB chemical

DFT

mio

pbc

ReaxFF

CO char

56 ± 11

55 ± 3 -

58 ± 8 45 ± 3

43 ± 2

a

Energy units are kcal/mol. All error bars listed here are determined from a single trajectory for each method; other error sources such as system size are not reflected and may exceed the listed error bars. Values omitted (including all fits for H2O) exhibited unphysical prefactors >1015 s−1.

composition of the chemical equilibrium. This approach has been employed recently to study shock wave induced chemistry in liquid methane, where significant changes in temperatures, pressures, and chemical reactivity were observed along the shock Hugoniot27 Using ReaxFF, we ran MD simulations with a quantum thermal bath at three target temperatures 2500 K, 3000 K, and 3500 K, and compared the results to the classical MD calculation in Figures 10, 11, and 12. We find that the

Figure 9. The percentage of carbon contained in char with respect to the total carbon in our computational cell as a function of time. We define molecules with carbon content higher than the initial phenolic molecules to be “char”. The dashed lines are fit to function [P] = [A](1-e−kt) used to determined the reaction barriers in Table 4.

Figure 10. Carbon ring statistics extracted from ReaxFF-based MD pyrolysis simulations. The black curve is classical MD (corresponding to the data in Figure 3), and the orange curve is produced using a Langevin quantum thermal bath thermostat for semiclassical quantum nuclear effects. Some reduction in carbon rings is observed at 2500 K in the quantum case. The classical and quantum cases are comparable at higher temperatures. 11123



Article

Figure 11. The evolution of carbon cluster size extracted from ReaxFF-based MD pyrolysis simulations. The top row is classical MD (corresponding to the data in Figure 4), and the bottom row is produced using a Langevin quantum thermal bath thermostat for semiclassical quantum nuclear effects. The quantum case at 2500 K exhibits more polymerization into larger carbon fragments than the classical case. Similar behavior is observed in the quantum and classical cases at higher temperatures.

Figure 12. The evolution of major reaction products extracted from ReaxFF-based MD pyrolysis simulations. The black curve is classical MD (corresponding to the data in Figure 6), and the orange curve is produced using a Langevin quantum thermal bath thermostat for semiclassical quantum nuclear effects. H and H2O populations at 2500 K are higher for quantum the case. Comparison to Figure 6 indicates that the quantum bath results in larger changes to the H and H2O populations at 2500 K than the variations between the three energy methods considered in this paper.

IV. CONCLUSIONS We present a comparison of different energy methods (ReaxFF, two parametrizations of DFTB, and DFT) for atomistic modeling of the initial stage of phenolic resin pyrolysis. Chars are characterized with respect to carbon cluster size and ring structure. For ReaxFF, the carbon structures are relatively inert during pyrolysis at 2500 K, while the other methods exhibit ring-opening and polymerization. ReaxFF underestimates ring populations at high temperatures relative to the other energy methods. Heterocyclic oxygen-containing ring structures are observed in the

char at 2500 K, which could have implications for the structure and properties of the char. In addition to char, small-molecule products are observed in our DFT-MD simulations, mainly including H2O, CO, CO2, H2, H, CmHn, and CmHnO. H dissociation from reactant phenolic chain molecules is the first reaction to occur. DFTB produces more CO than the other methods. This fact is correlated to more favorable CO formation energies for DFTB than the other methods. Therefore, proper reproduction of the formation energy is important for correct prediction of the product populations at chemical equilibrium. 11124



Article

(21) Kresse, G.; Joubert, D. Phys. Rev. B 1999, 59, 1758. (22) Goldman, N.; Reed, E. J.; Fried, L. E. J. Chem. Phys. 2009, 131, 204103. (23) Dunning, T. H. J. Chem. Phys. 1989, 90, 1007. (24) Franzblau, D. S. Phys. Rev. B 1991, 44, 4925. (25) Bagri, A.; Mattevi, C.; Acik, M.; Chabal, Y. J.; Chhowalla, M.; Shenoy, V. B. Nat. Chem. 2010, 2, 581. (26) Dammak, H.; Chalopin, Y.; Laroche, M.; Hayoun, M.; Greffet, J.-J. Phys. Rev. Lett. 2009, 103, 190601. (27) Qi, T.; Reed, E. J. Phys. Chem. A 2012, 116, 10451.

From the DFT trajectories, we identify pathways for H2O formation including a high temperature mechanism not elucidated previously. Two approximate CO formation channels are also identified. We make a rough extraction of single-step reaction barriers of H2O, CO, and char from the time-dependence of the populations in our MD, subject to statistical fluctuations. We do not find substantial differences between the mio and pbc parameter sets employed for the DFTB calculations. Finally, by coupling our pyrolysis system with a semiclassical quantum thermal bath, we find that quantum nuclear effects may have a significant impact on the chemistry at temperatures of 2500 K and below. Classical MD calculations with ReaxFF indicate relatively inert chemistry: the original phenolic ring structures are well preserved, and only a small amount of free H and water molecules are produced. When quantum nuclear effects are applied at the same temperature, however, we observe half of the phenolic rings break and the system forms large-size carbon clusters. The yield of both free H and water molecules are nearly doubled in the quantum case.

■

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]; [email protected] Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by a NASA Ames Center Innovation Fund grant. We also gratefully acknowledge the computational resources from the National Energy Research Scientific Computing Center (NERSC) and the Stanford NNIN Computing Facility (SNCF).

■

REFERENCES

(1) Bennett, A.; Payne, D. R. American Institute of Aeronautics and Astronautics 2013, 0183. (2) Yamashita, Y.; Ouchi, K. Carbon 1981, 19, 89. (3) Trick, K. A.; Saliba, T. E. Carbon 1995, 33, 1509. (4) Trick, K. A.; Saliba, T. E.; Sandhu, S. S. Carbon 1997, 35, 393. (5) Lawson, J. W.; Srivastava, D. Phys. Rev. B 2008, 77, 144209. (6) Till, S.; Heaton, A.; Payne, D.; Stone, C.; Swan, M. American Institute of Aeronautics and Astronautics 2013, 0182. (7) van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard, W. A., III. J. Phys. Chem. A 2001, 105, 9396. (8) Jiang, D.-e.; van Duin, A. C. T.; Goddard, W. A., III.; Dai, S. J. Phys. Chem. A 2009, 113, 6891. (9) Desai, T. G.; Lawson, J. W.; Keblinski, P. Polymer 2011, 52, 577. (10) Movahed, F. S.; Cheng, G. C.; Venkatachari, B. S.; Cozmuta, I. 42nd AIAA Thermophysics Conference, Honolulu, HI, 2011. (11) Seifert, G. J. Phys. Chem. A 2007, 111, 5609. (12) Elstner, M.; Porezag, D.; Jungnickel, G.; Elsner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Phys. Rev. B 1998, 58, 7260. (13) Chenoweth, K.; van Duin, A. C. T.; Goddard, W. A., III J. Phys. Chem. A 2008, 112, 1040. (14) Plimpton, S. J. Comput. Phys. 1995, 117, 1 http://lammps. sandia.gov. (15) Nielson, K. D.; van Duin, A. C. T.; Oxgaard, J.; Deng, W.-Q.; Goddard, W. A., III. J. Phys. Chem. A 2005, 109, 493. (16) Aradi, B.; Hourahine, B.; Frauenheim, T. J. Phys. Chem. A 2007, 111, 5678. (17) Koehler, C.; Frauenheim, T. Surf. Sci. 2006, 600, 453. (18) Kresse, G.; Furthmueller, J. Phys. Rev. B 1996, 54, 11169. (19) Kresse, G.; Furthmueller, J. Comput. Mater. Sci. 1996, 6, 15. (20) Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. Phys. Rev. B 1992, 46, 6671. 11125


Comparison of ReaxFF, DFTB, and DFT for ... - ACS Publications

Recommend Documents