Modeling initial stage of phenolic pyrolysis: Graphitic precursor formation and interfacial effects
Tapan G. Desai1,, John W. Lawson2 , Pawel Keblinski3
1 Advanced Cooling Technologies, Inc., Lancaster, PA 17601, USA
2 Thermal Protection Materials Branch, NASA Ames Research Center, Moffett Field, California 94035, USA
3 Department of Materials Science and Engineering, Rensselaer Polytechnic Institute, Troy, NY 12180, USA
Received 29 September 2010
Received in revised form 9 November 2010
Accepted 10 November 2010
Available online 18 November 2010
Reactive molecular dynamics simulation
Reactive molecular dynamics simulations are used to study the initial stage of pyrolysis of phenolic polymers with carbon nanotube and carbon fiber. The products formed are characterized and water is found to be the primary product in all cases. The water formation mechanisms are analyzed and the value of the activation energy for water formation is estimated. A detailed study of graphitic precursor formation reveals the presence of two temperature zones. In the lower temperature zone (<2000 K) polymerization occurs resulting in the formation of large, stable graphitic precursors, while in the high temperature zone (>2000 K) polymer scission results in formation of short polymer chains/molecules. Simulations performed in the high temperature zone of the phenolic resin (with carbon nanotubes and carbon fibers) show that the presence of interfaces does not have a substantial effect on the chain scission rate or the activation energy value for water formation.
Ablative materials are thermal insulators used in hypersonic space vehicles. They are typically carbon reinforced composites with a phenolic resin matrix, which absorb heat in part through endothermic, pyrolysis of the matrix. The char produced as a result of these reactive processes yields a thermally insulating and protective layer at the material surface. Mechanical strength of the char layer is required to prevent premature flaking and subsequent exposure of the virgin material underneath. The amount of material should be enough to compensate for any part of the char layer that might be ablated away during the re-entry. However, the weight of the material should be kept to a minimum. Optimization of the thermal protection system requires accurate prediction of the pyrolysis process and the evolution of char morphology. The pyrolytic processes are important for other applications as well, including carbonization of polymers, development of fire retardants, etc.
One recently developed ablative material is “phenolic impregnated carbon ablator” (PICA). This low-density material is a good choice due to its low thermal conductivity and efficient ablation properties. Phenolic resins have a high char yield under pyrolytic conditions. These reactions involve a complex sequence of reactions, and further development of these materials will require fundamental understanding of these molecular processes, and subsequent char evolution under transient conditions. Computational modeling of the resin-to-char process using atomistic-level simulation technique such as Molecular Dynamics (MD) with reactive force field can provide insights into the involved reaction barriers and pathways.
Recently, De-en Jiang et al.,  performed MD simulations with the reactive force field, ReaxFF, to simulate the initial stages of phenolic resin pyrolysis. Results from these simulations are benchmarks for the studies presented in this paper. The time scale accessible for routine MD simulations is limited to the nanosecond range, which allows study of only the initial stage of pyrolysis. Lawson and Srivastava  developed a high temperature, pyrolytic MD simulation method for polyethylene, based on the gradual removal of free Hydrogen atoms from the system. With their method, they were able to run the reactions to completion, however because of the approximations involved, accurate information about reaction kinetics was difficult to obtain.
In an effort to improve the mechanical properties of the char layer, NASA Ames Research Center has recently explored the effect of adding of carbon nanotubes (CNT) to PICA. Tensile test results showed a 35% increase in strength with the addition of 0.4 wt% functionalized CNTs in PICA, without any change in insulation properties . This increase in strength at low concentrations is possible because the percolation threshold for CNTcomposites is very low, on the order of 0.1 vol %. The CNT-polymer interface is expected to play a major role in maintaining the insulation properties by limiting the effective thermal conductivity of CNTs in the composite due to a large interfacial thermal resistance . Thus, without an appreciable increase in the weight or thermal conductivity, CNT addition can provide large mechanical reinforcement to the composite. In addition, the char morphology is known to change due to the addition of CNTs, as the char strongly adheres to the CNT yielding a cigar-like structure .
Furthermore, addition of CNT in flame retardant materials has been recently explored. Research has shown that the addition of only 0.5 wt% concentration of single-walled CNT (SWCNT) resulted in formation of a crack-free network-structured protective residual layer in PMMA composites and also reduced the rate of heat release . Formation of this crack-free protective layer is essential to prevent exposure of virgin ablative material underneath. However, in both the cases (ablation and flame retardant), the role of CNTs in altering the (already) complex pyrolysis reaction is not known.
In this paper, we attempt to understand the role of CNT and carbon fibers on the initial stage of pyrolysis. We performed detailed simulations with the ReaxFF method to characterize the products involved and the kinetics of pyrolysis. For the pristine phenolic resin case, (i.e., without carbon fillers) we also performed simulations to determine the conditions and processes involved in the formation of graphitic precursors. Two temperature zones were found to be important: in the lower temperature zone (<2000 K) polymerization occurred resulting in formation of stable graphitic precursors, and in the high temperature zone (>2000 K), disintegration resulted in the formation of short polymer chains/molecules. The effect of two fillers, carbon fibers and CNT, on the initial stage of pyrolysis was also studied. No substantial effect was found due to filler addition on the product (water) formation rate or polymer disintegration rate in the high temperature zone. The paper is organized as follows: Section 2 describes the simulation method and details; results are presented in Section 3; and Section 4 provides a summary and conclusion of this work.
MD simulations were performed with the LAMMPS molecular dynamics code . The reactive force field, ReaxFF , was employed to perform pyrolysis simulations. ReaxFF uses an empirical relationship between bond orders and bond lengths to allow smooth bond formation and breaking. This method has been quite successfully applied to a variety of chemical problems, including thermal degradation of polymers . Different ReaxFF parameter sets have appeared in the literature and in general we expect our results will depend on the particular choice of potential. To be consistent with previous work , we use the version of ReaxFF given by Chenoweth et al. . Further optimization of the potentials may be required in future work.
Simulations were performed on a non cross-linked phenolic formaldehyde resin (Fig. 1: top panel) in which two, neighboring phenol moieties are connected by a methylene group. The simulation cell, shown in Fig. 1: bottom panel, consisted of 16 phenol formaldehyde chains with 8 repeating units each in the orthoeortho sequence and a methyl group termination. The 16 polymer chains were randomly placed in the simulation cell by using the software Xenoview , which uses a Monte Carlo algorithm to build the amorphous structure. The simulation box size was adjusted to 2.62 nm 2.62 nm 2.62 nm to yield the experimental resin density of 1.25 gm/cc at 300 K. The polymer, simulation cell size and density are similar to the work performed by De-en Jiang et al. . The time step used to perform the simulation was 0.25 femtosecond (fs). Periodic boundary conditions were applied in all the directions of the simulation cell. A Berendsen thermostat with a temperature damping constant of 100 fs was used to maintain constant temperature of the simulation cell.
For composite systems, we chose to model two fillers e CNT and carbon fibers. Two CNT composite systems with two different single walled CNT chirality, (5,5) and (10,10) were studied. The radius of (5,5) and (10,10) CNT was 3.5 and 7.0 Å, respectively. In both systems, the CNT was 27 Å long and uncapped (or open) as shown in the inset in Fig. 2a. The phenolic resin structure was similar to the pristine case. The phenolic resin chains were introduced randomly around the CNT. Periodic boundary conditions were applied in all the directions of the simulation cell. The size of the simulation cell was adjusted such that the resin density was 1.25 gm/cc. The excluded volume due to the van der Waal’s distance between the CNT and resin was also accounted during this adjustment. Our simulations indicate that the van der Waal’s distance between the resin and CNT is 3.4 Å, similar to the work published by Wei and Srivastava for CNT polyethylene composite modeled with a Tersoffe Brenner potential and united atom model potential .
Currently, the ablation materials employed on space vehicle are made of phenolic resin impregnated within carbon fibers. Owing to the nanometer length-scale of the present study, the carbon fibers were simulated as surfaces with zero curvature, i.e., planar graphite sheets. As shown in the snapshot in Fig. 2b, the graphite block had 3 graphine sheets in the AeBeA stacking sequence. The surface area parallel to the sheet had the dimensions, 25.56 24.6 Å2 . The value of the surface area for the graphite block was selected to match the surface area of CNT (10,10). This was intended to provide insights on the effects introduced by the presence of tube curvature. The phenolic resin structure was similar to the pristine and the CNT-filled systems. The phenolic resin chains were introduced randomly on both sides of the graphite sheet. Periodic boundary conditions were applied in all the directions of the simulation cell and the size of the simulation cell was adjusted such that the resin density was 1.25 gm/cc.
The equilibration procedure for all systems (pristine and composite) involved (i) equilibration at 0.1 K with an NVT-MD simulation for 40 ps (ps), followed by, (ii) thermal annealing procedure e heating the system with a ramp rate of 0.02 K/fs to 1000 K, equilibrating at 1000 K for 20 ps and then cooling to 300 K with the same ramp rate. This annealing procedure was performed twice followed by an equilibration at 2000 K for 200 ps. Our in-house code was used to examine and identify the new molecules that were produced during the pyrolysis simulations. Our analysis confirmed that no reactions took place during the 200 ps equilibration run. Pyrolysis simulations were performed at 5 different temperatures e 2750, 2875, 3000, 3125 and 3250 K to obtain reaction rates and activation energy. To improve the statistical accuracy and provide error bars for the reaction rate values, all simulations were performed with 4 different starting structures at each temperature.
3. Results and discussions
This section is divided into two sub-sections e pristine phenolic resin and phenolic resin composites.
3.1. Pristine phenolic resin
Under this subsection, first the product evolution is discussed including the associated reaction rates and activation energies. Then, the details on graphitic precursor formation are presented.
3.1.1. Evolution of products
A detailed structural analysis on the equilibrated phenol formaldehyde resin is performed to analyze the accuracy of the structure yielded by the relatively new ReaxFF potential. NPT simulation performed to estimate the resin density at 300 K indicates that the resin density is only 4.5% lower than the experimental resin density MD simulations at 300 K (see Fig. 12). A good agreement was seen which validates the ReaxFF potential at room temperature. The details on this comparison are provided in Appendix A.
High temperature simulations at 5 different temperatures in the range of 2750e3250 K are performed on 4 different equilibrated structures. Water is the primary product in all these simulations. The second most abundant product is hydrogen gas. The number of water molecules formed (averaged over 4 different starting structures) is plotted as a function of time in Fig. 3a. The rate of water formation (i.e., the slope) increases with temperature, suggesting a thermally activated reaction process. At each temperature, the value of the water formation rate is calculated from a least-square fit of the data sets and plotted in Fig. 3b. The error bars in Fig. 3b indicate the difference between the maximum (minimum) rate from the different starting structures and the average rate. Using the Arrhenius equation, a value of 124 20 Kj/mol is obtained for the activation energy for water formation. This value matches well with the value of 144 28 kJ/mol published by Jiang et al. .
Additionally, simulations are performed for an equilibrated structure at three different temperatures, 2750, 3000 and 3250 K with simulation time extended to 200 ps to reach completion of reactions. Initially the water formation rate depends on the temperature. At all temperatures, the total number of water molecules formed reaches a plateau after a certain time and converges to a value of w100. The simulation cell contains 128 oxygen atoms, which indicates a yield (defined as the ratio of actual product formed to the theoretical) of w78%. The different mechanisms involved in water formation (including the transition states) are analyzed and presented in Appendix B. The remaining oxygen atoms are consumed during formation of CO, CO2, and, small amounts of aldehydes, ketones and alcohols.
3.1.2. Graphitic precursor formation
The ablation process during the re-entry of the space vehicles involves carbonization of the phenolic resin under aninert atmosphere at high temperatures ranging from 1500 to 2000 K that leads toformation of glassy carbon with large graphene fragments. Shown in Fig. 4 is a snapshot of a small graphitic precursormolecule containing 3 fused rings (two 6-membered rings and one 7-membered ring) formed after a 40 ps run at 3250 K. A commonly proposed reaction pathway  is that an open benzene ring and a closed benzene ring combine to form fused ring. However, our simulations illustrate that alkyl-like fragments from four different chains break, rearrange themselves and combine to form the fused ring fragment. Thus, it reveals another complex reaction pathway for graphitic precursor formation. This fused-ring structure is unstable at longer times, possibly due to the extremely high temperatures, and fragments into smaller carbon molecules.
To identify a temperature zone where the graphitic precursors would be stable, detailed chain disintegration analysis is performed. In Fig. 5, the number of chains as a function of chain length (averaged over 4 different starting structures and a simulation time of 2 ps) is plotted on a semi-log scale for two different times e 20 ps (top panel) and 40 ps (bottom panel) and 3 different temperatures e 2750, 3000 and 3250 K. At the start of the simulation run, the equilibrated structures have 16 chains with a carbon chain length of 56. At all temperatures after 20 ps, polymerization occurs between polymer chains resulting in longer chains of carbon chain length up to 120 at 3250 K and 190 at 3000 K. However, at longer times, all the long chains break to form smaller, shorter chains. This fragmentation into smaller chain occurs for other temperatures (2750 and 3000 K) at simulation times longer than 40 ps. The number of C2 molecules (C2 stands for molecule which has two carbon atoms) is seen to increase. Included in the number of C2 molecules are a large amount of C2H2, C2H radicals and C2 gas molecules. This increase in C2 molecules with time is seen at all three temperatures. The fragmentation indicates that the temperatures are too high to polymerize and form large graphitic molecules.
The results from Fig. 5 suggest that pyrolysis kinetics can be divided into two temperature zones. The high temperature (fragmentation) zone involves breaking of polymers into smaller molecules and the low temperature (polymerization) zone promotes bonding between polymers to form large connected graphitic network. In order to estimate a temperature range for these zones, in Fig. 6, the total number of C2 molecules formed after 40 ps is plotted as a function of temperature on a semi-log scale. The slope through the data points is extrapolated to intersect the x-axis (temperatureaxis). The intersection occurs at 2050 K signifying a higher probability of graphitization at temperatures below 2000 K.
To capture large (stable) graphitic precursor formation a simulation was performed at 1750 K. Since, no pyrolysis reactions are observed at 2000 K during the equilibration run for 250 ps, it can be expected that the time scale to observe any chemical reaction at 1750 K would go beyond the time scale accessible by MD simulations. Hence, to capture pyrolysis at 1750 K, a partially pyrolyzed structure from a 2750 K run after 40 ps was chosen as an input structure. In this structure more than 60% of the chains were intact. A reaction generally proceeds in 3 stages: initiation, propagation and termination. By using a partially pyrolyzed structure, the initiation stage is drastically reduced. In addition, after every 1.25 ns, stable molecules such as H2O and CO2 are removed from the simulation cell. These stable gases interfere during the graphitization process preventing formation of large graphene precursors. At the end of a 12 ns simulation run two large graphitic structures are found, one containing 7 fused rings as shown in Fig. 7 and other containing 5 fused rings. The precursor starts with 2 fused rings and grows over the 12 ns simulation run. This indicates that the reaction kinetics is different in the two temperature zones and it is necessary to perform the MD simulations at the temperature of experimental interest for accurate representation of the reaction behavior rather than merely increasing the MD simulation temperature to capture the reactions on the MD time scale.
3.2. Phenolic resin composites
Under this subsection, the effect of addition of CNT and carbon fibers on the initial stage of pyrolysis is discussed.
3.2.1. CNT-Phenolic resin composite
Pyrolysis simulations are performed at three different temperatures 2750, 3000 and 3250 K, on two composites containing CNT with (5,5) and (10,10) chirality. The CNT with two different chirality, is chosen to evaluate the effect of diameter on the pyrolysis simulation. In Fig. 8, the evolution of the primary product, water, is plotted as a function of time at three different temperatures, for CNT composites and compared with the pristine phenolic resin case. For all cases, the data is averaged over 4 different starting structures to improve the statistical accuracy. The addition of fillers did not change the water production rate. The activation energy for water formation for the two types of CNT composites and the pristine case is same, 124 20 kJ/mol. Thus, the chemical reactions in the initial stages of the pyrolysis are not affected by the presence of CNT, even at very high weight loading of 19% and 38% for CNT (5,5) and (10,10), respectively.
Detailed analysis on the structure of CNT (5,5) after a simulation run of 40 ps for 4 different starting structures is performed. As shown in the various snapshots in Fig. 9, different starting structures yield a different number of chemically adsorbed molecules on the CNT surface. Hence, it is important for statistical purposes to perform 4 to 5 different simulation runs at each temperature of interest. Fig. 9a, suggests that a CNT can play host to a wide range of molecules from a single hydrogen atoms to a 5-membered ring to large polymer chains. Even at 3000 K, the CNT is still stable after 40 ps and chemical adsorption of polymers. At 3250 K, in Fig. 9c, CNT is distorted for 2 different starting structures, but it is found to be stable for other two starting structures after 40 ps. Separate simulations were performed on CNT without the phenolic resin at 3250 K to test the stability of CNT at this high temperature with ReaxFF potential. The CNT was stable at 3250 K during a 40 ps simulation run implying that the distortion of CNT is caused due to the physical and chemical interactions with the phenolic resin.
3.2.2. Carbon fiber-phenolic resin composite
As mentioned in Section 2, the carbon fiber is modeled as a block of graphite at the nanometer length-scale. Pyrolysis simulation results are performed at 3250 K and averaged over 4 different starting structures. The number of water molecules formed is plotted against time and compared with the results from the pristine, CNT (5,5) and CNT (10,10) composite case in Fig. 10. A strong overlap of the water formation rate with no appreciable difference is seen in all cases. Similar overlap was seen at lower temperatures (2750 and 3000 K). Thus, it can be concluded that the presence of fillers such as CNT and carbon fibers do not play a catalytic role in the dehydration stage of the pyrolysis of the phenolic resin, even under high loading conditions.
Even at 3250 K, for all the 4 structures studied, the graphite block is found to be stable during the simulation run of 50 ps. At the end of the simulation run, small molecules such as methyl, ethyl group are bounded to the graphite layer. For the lowest pyrolysis simulation temperature studied at 2750 K no molecules were bounded to graphite (in contrast to CNT, see Fig. 9a). For a simulation run at 3000 K, similar to other two temperatures for the fiberresin composite, no large molecules were found bounded to the graphite. Thus, in comparison to carbon fiber, a stronger covalent binding occurs in the case of CNT filler addition.
The dehydration rate does not represent the entire kinetics of the phenolic resin pyrolysis and therefore it is unjustifiable to claim that the filler addition does not have any effect on the pyrolysis process. However, disintegration of polymer chains to form shorter molecules encompasses a large part of the overall pyrolysis process. Hence, in Fig. 11, we track the number of polymer chains with carbon chain length greater than 50 as a function of time for 4 different starting structures and at 3 different temperatures. The original carbon chain length of the phenolic resin in the starting structure is 56. The value of 50 during the disintegration analysis is used for statistical purposes to prevent any chain end effects. At 3000 and 3250 K, the disintegration rate for the pristine case (in blue) overlaps the CNT case (in green) and the graphite case (in red). At 2750 K, there is a slight discrepancy between the disintegration rates for simulation times greater than 20 ps. However, to confirm the effect of fillers, the simulations need to be performed in the low temperature zone (<2000 K). Experimentally, it has been observed the addition of CNT produced cigar-like structures where the char strongly adheres to the CNT.
The primary objective of this paper was to model the effect of filler materials such as CNT and carbon fibers on the initial stage of the phenolic resin pyrolysis. Water was found to be the primary product for all the cases. There was no major effect of the filler addition on the water formation rate or the polymer disintegration rate. The evolution of graphitic precursors was studied rigorously. Two temperature zones were identified: the high temperature zone (>2000 K) involved fragmentation leading to non-graphitic structures and the low temperature zone (<2000 K) involved polymerization leading to predominantly graphitic structures. Formation of stable large (7 fused rings) graphitic precursors was captured at 1750 K during a simulation run of 12 ns. Thus, running the simulations at low temperatures was identified as the key to study the formation and growth of graphitic precursors. In the high temperature zone, no effect of fillers was observed on the graphitic precursor formation. To capture the cigar-like structures observed at 1000 K experimentally for the CNT-resin composites, the simulations need to be performed in the low temperature reaction zone.
The dehydration observed in our simulations is also found to be prominent in experiments during the curing stage of the resin which involves crosslinking of polymer chains to improve the mechanical properties of resin. However, the curing occurs at lower temperatures w450e550 K. Jiang and coworkers showed that the simulation results for activation energy falls within the experimental range of 74e199 kJ/mol for four mass-loss peaks . The experiments were performed by Trick and Saliba on a cured sample of carbon-phenolic resin composite, where the pyrolysis kinetics was calculated from thermo-gravimetric analysis. It is difficult to compare our MD simulation results with such experiments because: (i) experiments are performed on a composite structure with carbon fibers as filler material; (ii) the phenolic resin is cured at 450 K in experiments, where most of the dehydration takes place resulting in crosslinking of polymer chains; (iii) for pyrolysis experiments, the heating rate is 1 C/min in the temperature range of 573e1173 K, whereas simulations are performed with a stepchange in temperature in the range of 2750e3250 K. Thus, in addition to a different starting structure, an additional assumption is needed that the reaction kinetics in the experimental temperature range of 573e1173 K is same as the kinetics in the simulation temperature range of 2750e3250 K. According to our knowledge, there are no experiments to which we can directly compare our simulation results.
The reactive MD simulations are definitely an important tool to provide insights on the involved reaction pathways and energy barriers. However, corresponding experiments are needed to validate these simulations. The experiments should include high heating rates to reach high temperatures (within a nanosecond time scale) such that the kinetics are monitored on the time scale accessible to MD simulations. State-of-the-art laser heating of polymer samples couple with spectroscopy is an option. In addition, the MD time scale needs to be increased from nanosecond to micro-seconds in order to perform simulations at temperature of interest (1500e2000 K) without losing the atomistic information. This increase in MD time scale can be achieved by applying appropriate accelerated MD techniques . Nevertheless, the results serve as a good model to provide insights on reaction mechanisms with caution while direct comparison with experiments.
The work was supported by the NASA Small Business Innovation Research Grant (SBIR), under Contract No. NNX10CC69P. We thank Dr. Howard Pearlman from Advanced Cooling Technologies, Inc. for his input during the project.
Appendix A. Comparison of radial distribution function e ReaxFF vs. ab-initio calculations
Since ReaxFF is a relatively new force field, we tested its reliability for this specific application by comparing the radial distribution function (rdf) obtained from ReaxFF against results from ab-initio based MD simulations. In particular, we performed simulations for a single phenolic chain with four repeating units at T ¼ 300 K and confined in a periodic box of size 20 Å on a side. The ab-initio simulations were performed with the code VASP within the context of density functional theory (DFT). Note that unlike other implements of ab-initio MD, VASP evolves the electronic and ionic degrees of freedom separately: forces obtained from DFT are converged first, and then used to advance the ionic positions. The GGA of Perdew and Wang was used with a cutoff of 300 eV and computations were done using the Projected Augmented Wave (PAW) potentials. Results are shown in Fig. 12. As can be seen the agreement is very good. Various bond lengths identified by peaks on the rdf agree very well although there is some variation in the peak heights. It should not be concluded that this represents a definitive test of ReaxFF. However, DFT is very good at describing bonded interactions, but is well known to have difficulties with weaker, long-ranged interactions that are common to polymer systems and are in principle included in ReaxFF.
Appendix B. Water formation mechanisms
Several dehydration mechanisms were observed during the different simulation runs. Similar to observations by Ref. , the dehydration mechanisms involved were: (i) between two eOH groups on phenols, (ii) a disassociated H and eOH group on a phenol, (iii) H from a benzene ring and eOH group from a phenol, and (iv) H from the methyl group (between two benzene rings) and eOH group on a phenol. For these four mechanisms, both intermolecular and intramolecular reactions took place.
Shown in Fig. 13 is the temporal snapshot of a single polymer chain (for the pristine case) undergoing the fourth mechanism during intramolecular dehydration. The single polymer chain is isolated from the rest of the 15 polymer chains in the simulation cell for clarity. The simulation temperature is 2750 K and the simulation time is in the range of 9.025e9.2 ps. After 9.05 ps, the first transition state involves bonding between the eOH group from a phenol and hydrogen atom in the adjacent methyl group. At 9.075 ps, the eOH group breaks from the phenol, however, at 9.1 ps, thewater molecule reattaches to the benzene ring. After 9.2 ps, the water molecule breaks free of the polymer and diffuses away from the parent phenolic resin chain. During this period, the polymer chain is intact without undergoing any substantial disintegration (except thermal dissociation of couple of hydrogen atoms). Thus, Fig. 13 provides an atomistic picture of the transition states involved during a particular water formation mechanism for a thermally stable polymer chain. This process can be adopted to develop reaction pathways and extract associated chemical kinetic information to build highly detailed atomistically-informed chemical kinetic models.
 Jiang D, van Duin ACT, Goddard WA, Dai S. J Phys Chem A 2009;113 (25):6891e4.
 Lawson JW, Srivastava D. Phys Rev B 2008;77(14):144209-1e144209-6.
 NASA biennial report, http://research.jsc.nasa.gov/BiennialResearchReport/ PDF/Eng-14.pdf; 2007.
 Shenogin S, Xue L, Ozisik R, Keblinski P, Cahill DG. J Appl Phys 2004;95 (12):8136e45.
 Fan W. Private communication; 2010.
 Kashiwagi T, Du F, Douglas JF, Winey KI, Harris RH, Shields JR. Nat Mater 2005;4(12):928e33.
 Plimpton S. J Comput Phys 1995;117(1):1e19, http://lammps.sandia.gov.
 van Duin ACT, Dasgupta S, Lorant F, Goddard WA. J Phys Chem A 2001;105 (41):9396e409.
 Chenoweth K, Cheung S, van Duin ACT, Goddard WA, Kober EM. J Am Chem Soc 2005;127(19):7192e202.
 Chenoweth K, van Duin ACT, Goddard WA. J Phys Chem A 2008;112 (5):1040e53.
 Shenogin S, Ozisik R. Xenoview: visualization for atomistic simulations, http:// www.rpi.edu/wshenos3/xenoview.html; 2007.
 Berendsen HJC, Postma JPM, van Gunsteren WF, Dinola A, Haak JR. J Chem Phys 1984;81(8):3684e91.
 Wei C, Srivastava D. Nano Lett 2004;4(10):1949e52.
 Schürmann BL, Vogel L. J Mater Sci 1996;31(13):3435e40.
 Fyfe CA, McKinnon MS, Rudin A, Tchir WJ. Macromolecules 1983;16 (7):1216e9.
 Trick KA, Saliba TE. Carbon 1995;33(11):1509e15.
 Voter A, Montalenti F, Germann T. Annu Rev Mater Res 2002;32(1):321e46.