Establishing Physical and Chemical Mechanisms of Polymerization and Pyrolysis of Phenolic Resins for Carbon-Carbon Composites

The complex structural and chemical changes that occur during polymerization and pyrolysis critically affect material properties but are difficult to characterize in situ . This work presents a novel, experimentally validated methodology for modeling the complete polymerization and pyrolysis processes for phenolic resin using reactive molecular dynamics. The polymerization simulations produced polymerized structures with mass densities of 1.24 ± 0.01 g/cm 3 and Young’s moduli of 3.50 ± 0.64 GPa, which are in good agreement with experimental values. The structural properties of the subsequently pyrolyzed structures were also found to be in good agreement with experimental X-ray data for the phenolic-derived carbon matrices, with interplanar spacings of 3.81 ± 0.06 Å and crystallite heights of 10.94 ± 0.37 Å. The mass densities of the pyrolyzed models, 2.01 ± 0.03 g/cm 3 , correspond to skeletal density values, where the volume of pores is excluded in density calculations for the phenolic resin-based pyrolyzed samples. Young’s moduli are underpredicted at 122.36 ± 16.48 GPa relative to experimental values of 146 – 256 GPa for nanoscale amorphous carbon samples.


Introduction
Structural materials capable of withstanding extreme temperatures are necessary for vehicles used for hypersonic flight and space exploration endeavors.Thermal protection system (TPS) materials have been developed to address this challenge with two main objectives: (1) provide thermal insulation to critical structural components from high temperatures associated with planetary entry and (2) provide structural support while maintaining low weight to reduce payload costs.One class of materials often used for TPS purposes is carbon-carbon (C/C) composites, whose matrix characteristics can vary depending on the polymer precursor and processing conditions.C/C composites can be fabricated via successive polymer-infiltration-pyrolysis (PIP) cycles to counter shrinkage and porosity formation during the polymerization and pyrolysis steps.Phenolic resins are commonly used as C/C composite matrix precursors due to their relatively high char yield, or retained percentage of mass after pyrolysis, of about 55%, which reduces the number PIP cycles required to achieve a dense product acceptable for service [1].Phenolic resins polymerize via a condensation reaction to produce linear phenol chains known as novolacs.Novolacs can be ground into a powder and mixed with a hardener to form a highly branched network of phenols called Bakelite, which is a carbon matrix precursor [2].However, the insoluble and infusible nature of Bakelite and the complex chemical and structural changes that occur during pyrolysis make experimental characterization and process optimization difficult.
Molecular dynamics (MD) is a powerful computational modeling technique that can provide insight into polymerization and pyrolysis conditions and the resulting material structure and properties [1][2][3][4][5].The choice of MD force field, which governs interatomic interactions, is critical to accurately capture the chemistry of high-temperature processes at the atomic scale [6].Reactive force fields have been developed and tested with specific systems to simulate bond formation and dissociation [7][8][9].Although simulations performed with reactive force fields are computationally intensive, one can simulate the bond breaking or forming events which may play an important role in modeling the pyrolysis step in C/C composite manufacturing.Perhaps the most widely-used, bond order-based reactive force field is the ReaxFF force field developed by van Duin et al. [10].ReaxFF has been applied to simulate various materials and processes, including coal combustion [11], protein/DNA interactions in aqueous solutions [12], polymer crosslinking and mechanical properties [8], catalysis [13], and synthesis of 2D materials [14,15].
A number of investigators have used MD simulations to model the polymerization and the early stages of pyrolysis for phenolic resin.One of the earliest MD models was constructed by Izumi et al. [16], which was used to establish methodologies for simulating polymerization.Monk et al. [17] investigated algorithms to generating atomistic models of crosslinked phenolic resins and found that an iterative crosslinking approach, rather than a single step crosslinking, was preferred for predicting thermal properties that match the experimental literature.Shudo et al. [18] developed a probabilistic pseudo-reaction algorithm to establish crosslinked phenolic resin models and characterized the resulting structure based on the location of the methylene linkages.Shudo's study highlighted the importance of the initial system configuration by demonstrating that monomeric phenol produced higher degrees of crosslinking, while longer oligomeric chains were unable to fully crosslink due to structural constraints.In a subsequent study [19], Shudo calculated the Young's modulus and mass density of the pseudo-reaction-based phenolic resin models.Izumi et al. [20] developed a united-atom MD model using the same pseudo-reaction algorithm developed by Shudo and introduced a geometric constraint to avoid highly strained structures.Jiang et al. [21] modeled the initial phase of phenolic resin pyrolysis using the ReaxFF reactive force field and found that H2O was the first gaseous compound to volatilize from the matrix.Similarly, Desai et al. [22] also modeled the initial phase of phenolic resin pyrolysis using a variant of the original ReaxFF and found that H2O and H2 were the primary products of the pyrolysis process and that non-graphitic amorphous carbon was the resulting pyrolyzed structure at higher temperatures (>2000 K).Qi et al. [23] performed a systematic comparison between density functional theory (DFT) methods (DFT and self-consistent charge density-functional tight-binding, or SCC-DFTB) and ReaxFF for modeling the initial stages of phenolic polymer pyrolysis and found that the volatilized products consisted of H2O, CO, CO2, H2, H, CmHn, and CmHnO.With respect to ReaxFF, Qi noted that the carbon structures were relatively inert during pyrolysis at 2500 K.In a subsequent study, Bauschlicher et al. [24] determined that ReaxFF had the highest average absolute error for bond and reaction energies and barrier heights (19.6 kcal/mol) compared to various DFT methods (3.8 kcal/mol for B3LYP, 5.1 kcal/mol for PW91, 17.4 for DFTB(mio), and 13.2 kcal/mol for DFTB(pbc)).Zhong et al. [25] investigated the behavior of phenolic hydroxyl groups during the initial stages of pyrolysis and found that highly active phenolic hydroxyl groups volatilized into phenoxyl and hydroxyl radicals.In their study, it was determined that phenoxyl radicals reduced stability of benzene rings and damaged the backbone of the phenolic resin, while hydroxyl radicals reduced char yield by releasing more CO and CO2 molecules.Harpale et al. [26] used a combination of a fixed-bond force field and a reactive force field to quantify pyrolysis kinetics and found that, for highly-crosslinked phenolic resins, pyrolysis occurred at temperatures greater than 500 K via removal of -OH and -H groups from aromatic carbon rings.Purse et al. [27] recently demonstrated a method for simulating the complete pyrolysis of phenol-formaldehyde resins using reactive MD.The proposed method demonstrated that periodic removal of volatile pyrolysis gases allowed for full carbonization after 2 ns of simulation time.Yan et al. [28] recently employed an approach of switching from a fixed-bond to a reactive force field to investigate the thermostability of crosslinked and non-crosslinked phenolic resins during pyrolysis and determined that highly-crosslinked phenolic resins retained their branched structures and exhibited increased structural stability up to pyrolysis, as compared to non-highly-crosslinked phenolics.Gissinger et al. [29] developed an algorithm for simulating pyrolysis which was able to predict char yields of various polymer systems.These previous works have identified potential methods for modeling pyrolysis within the nanoscale time regime accessible in MD simulations.However, to date there has been no comprehensive modeling approach established for phenolic resins that can capture both the complete polymerization and pyrolysis processes while predicting accurate physical and mechanical properties.
The objective of this work was to establish experimentally validated MD protocols for modeling phenolic resin polymerization and pyrolysis.First, polymerized models were created and validated through comparison of elastic properties with the literature.A parametric study was then conducted to identify optimal MD simulation parameters to produce pyrolyzed models with char yields and graphitic crystallite sizes in agreement with experimental data.The optimal parameters were then used to simulate a larger molecular system to obtain statistically significant results.Developing experimentally validated MD protocols for modeling the complete processing procedure of a carbon matrix precursor can allow for material and process optimization by analyzing material and elastic properties during the polymerization and pyrolysis processes.Furthermore, the established protocols may be applied to other material systems of interest for TPS applications.

Methods
All simulations were performed using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) 2020 software package with two variants of the ReaxFF force field [7] using simulation timesteps of 0.1 fs.The CHON-17_weak (CHON17W) parameter set developed by Zhang et al. [9] was used for the polymerization step due to its ability to simulate polymerization at high temperatures and accurately predict mass densities and mechanical properties.The CHON19 parameter set developed by Kowalik et al. [30] was used for simulating pyrolysis as it was developed specifically for hydrocarbon pyrolysis processes.The complete modeling procedures and corresponding outputs can be divided between the polymerization and pyrolysis processes and are described in terms of four steps (Figure 1): (1) model initialization, (2) polymerization simulation, (3) pyrolysis parameterization, and (4) pyrolysis simulation.The study began with five MD replicates of 4,644 atoms each (hereafter labeled as test systems TS1, TS2, TS3, TS4, and TS5).These smaller, more computationally efficient models were used to investigate the statistical variations that occur in the predicted results due to different initial configurations.A large system of 30,960 atoms was then built to determine how accuracy and precision are affected by model size.Once the optimal simulation parameters were established for one test system (TS1), the other test systems and the large system were subjected to pyrolysis using the same parameters.Molecular visualizations were accomplished using the TopoTools plugin to the Visual Molecular Dynamics (VMD) software [31,32] and OVITO [33].
Figure 1.MD modeling steps and corresponding outputs in bullets.Each color in the densification stage of the initialization panel represents a single cluster of atoms.The green cluster in the polymerization panels indicates the largest cluster of atoms in the system.Outputs are color coded in blue for polymerization and pyrolysis metrics (extent of reaction, mass densities, and char yield), green for mechanical properties (Young's modulus and Poisson's ratio), red for pyrolysis settings (pyrolysis temperature, pyrolysis pressure, and volatile removal interval), and black for crystallographic dimensions (interplanar spacing, stacking height, and crystallite width).

Model initialization
The phenolic monomers were constructed in ChemDraw 3D and consisted of triphenols with single open valencies at chosen reaction sites, as depicted by the blue color-coded functional groups in Figure 2.These initial configurations were chosen such that the polymerized resin would consist of a higher degree of ortho-para linkages relative to ortho-ortho and para-para linkages, resulting in a higher degree of branching, as expected for a Bakelite structure [2].As shown in Figure 2, reactions 1 -4 allow the triphenols to link, while Reactions 5 and 6 allow for further polymerization of the system.The ChemDraw 3D files were converted to ReaxFF data files.The four molecules were then mixed in the constant volume and temperature (NVT) ensemble at 300 K for 100 ps, followed by a replication step of a factor of 27 and 180 to increase the number of atoms to 4,644 and 30,960 for the test systems and the large system, respectively.The systems were then densified to a target mass density of 1.15 g/cm 3 in a two-step process by first decreasing the box size to 40×40×40 Å and then shrinking each axis uniformly at a rate of 15 Å/ps to bring the molecules closer together and increase the probability of the occurrence of the intended reactions.The target density was chosen as a value less than the minimum density for a crosslinked phenolic resin (1.20 g/cm 3 ) [3].

Polymerization
The densified systems were polymerized by heating to 1000 K in the NVT ensemble at a rate of 10 K/ps and holding for 500 ps, followed by cooling to 300 K at 10 K/ps, with 1 ns relaxation at 300 K in the constant pressure and temperature (NPT) ensemble.The degree of polymerization was tracked via the extent of reaction for branching polycondensation reactions [34]: where p is the extent of reaction or conversion, No is the number of initial molecules, N is the number of molecules at any time t, and fav is the average monomer functionality (for example, 3 for phenolic resin).The critical conversion, pc, is the estimated conversion required to achieve gelation and is predicted by the simple relation [35,36]: Since phenolic resin has an average monomer functionality of 3, the estimated pc value is 66.67%.
Another metric used to track polymerization was the size of the largest cluster in terms of percent mass.In MD, a system can be considered to have reached full polymerization when a single cluster comprises most of the periodic simulation box by percent mass.For this work, a fully polymerized system was defined as one in which the largest cluster mass was at least 90% of the system mass and reached an extent of reaction of at least 66%.
The bond order evolution during polymerization was recorded and analyzed using the bond order information outputted via the "fix reaxff/bonds" command in LAMMPS, which writes bond information computed by ReaxFF as specified by Aktulga et al. [37].Bond order information was output every 0.01 ps and averaged over a total time of 0.1 ps to avoid transient bonds from being considered, and then used to identify the existing bonds in the system using a bond order cut-off of 0.3.The extracted bonds were used to compute the total number of clusters, their sizes, and their mass in daltons (Da), which were then used to calculate the extent of reaction and cluster sizes in terms of mass percent.
Once full polymerization was confirmed and equilibrated mass densities verified, uniaxial tension simulations were performed at a strain rate of 2×10 8 s -1 by deforming the simulation boxes along each principal axis up to 10% strain.The uniaxial Young's modulus I for each direction was calculated as the slope of a linear regression along the linear portion of the stress-strain plot as demonstrated by Odegard et al. [38], and the reported values of E are the triaxial averages.
Poisson's ratio () values were calculated through a linear regression of transverse strain plotted against longitudinal strain as specified by Radue et al. [39], and the reported values of  are the triaxial averages.

Parameterization
A parameterization study was performed before the pyrolysis production simulations to find optimal simulation settings to produce pyrolyzed structures in agreement with the reported experimental data [1,[40][41][42][43].The temperatures of the polymerized models were raised to the pyrolysis temperature (Tpyr) at a rate of 10 K/ps in the NVT ensemble.Once the models had reached pyrolysis temperatures, the constant pressure and enthalpy (NPH) and the constant volume and energy (NVE) ensembles were applied to control pressure and temperature, respectively.Hydrostatic pressures (Ppyr) were applied once the systems had reached pyrolysis temperatures to maintain system stability.Due to the number of simulations involved in the parametric study, only a single test system was subjected to the parameterization procedure.Table I shows the parameter values for each study and the range of values used for the variable parameter.Two pressure studies were performed to capture its effect at low and high temperature pyrolysis.The effect of three simulation parameters (pyrolysis temperature, pyrolysis pressure, and the removal interval of volatiles, t_rem) on the pyrolyzed molecular structure was assessed based on graphitic crystallite height (Lc), interplanar spacing of the (002) planes (d002), and char yield (Cy) during a 1 ns pyrolysis simulation.The graphitic crystallite width (La) was not used as a dependent variable due to a significant reduction in simulation box dimensions after pyrolysis.Table II shows the resulting box dimensions after pyrolysis and the number of remaining atoms for all systems.For test system 1, the test system used for the parametric study, all box dimensions are all less than the crystallite width values of 44 Å reported by Vieira et al. [40].Based on the results of the parametric study, the optimal pyrolysis settings established for TS1 were determined as Tpyr = 3200 K, Ppyr = 2 GPa, and t_rem = 0.1 ps.As will be shown, these settings produced char yield and crystallite sizes in reasonable agreement with literature [1,40,44].Once these settings were identified, the remaining four replicates (TS2, TS3, TS4, and TS5) and the larger system (LS) were subjected to the same pyrolysis settings to ensure reproducibility.

Pyrolysis
Several metrics of interest were dynamically tracked during the pyrolysis simulations, including the char yield (Cy), mass of hydrogen (%H), mass of oxygen (%O), hybridization content (sp, sp 2 , and sp 3 carbon), and carbon ring content (5-membered, 6-membered, and 7-membered carbon rings).All metrics were measured through mass percent relative to the total mass of the system.The species volatilized every 0.1 ps were also tracked.Mass density was calculated by averaging the last 100 ps of a 1 ns relaxation at room temperature.
Volatiles were identified using the "fix reaxff/species" command in LAMMPS and removed with the "delete" option developed by Gissinger et al. [29], while also outputting a file with the species deleted at every time step.Species with masses between 2 -50 Da were deleted to simulate the volatilization of low to medium mass compounds.The pyrolyzed samples were then cooled back to 300 K at 10 K/ps and relaxed in the NPT ensemble for 1 ns.Structural analyses were conducted through X-ray diffraction (XRD) via the "compute xrd" command of LAMMPS to obtain the diffuse (002) and (100) peaks characteristic of amorphous carbon, as well as the corresponding diffraction angles (θ002 and θ100, respectively) [45].Diffraction data was output every 0.1 ps and averaged over 1 ps.Depending on the pyrolysis settings used, the resulting XRD spectra was classified as either amorphous or semi-crystalline.Examples of amorphous and semi-crystalline XRD spectra for low-temperature and high-temperature pyrolysis simulations, as well as lowpressure and high-pressure pyrolysis simulations are shown in Figure 3.It is evident that high temperature, high pressure, or a combination of both can cause partial graphitization and produce semi-crystalline XRD spectra, as indicated by the manifestation of a broad peak in the 20 -30 o 2θ range.For spectra in which some crystallinity was observed, peaks were curve-fit to a pseudo-Voigt profile to facilitate identifying the peak position (2θ) and the full-width at half-maximum (β).The pseudo-Voigt profile is commonly used in XRD peak-fitting algorithms and is a convolution of normalized Lorentzian and Gaussian curve fits [46].The interplanar distance (d002) was calculated through Bragg's Law: where λ is the characteristic Cu-Kα radiation wavelength of 1.54 Å.The characteristic graphitic crystallite sizes (stacking height Lc and stacking width La) were calculated through variants of the Scherrer equation as demonstrated by Vieira et al. [40]: and Figure 3. XRD spectra for amorphous and partially graphitized models.The low and high temperature simulations were performed with hydrostatic pressure of 0.5 GPa, while low and high pressure simulations were performed at 2700 K. Low temperatures and low pressures result in amorphous structures (blue and gray curves), while the appearance of a peak around the 20 -30 o 2θ range for higher temperatures and pressures indicate partially graphitic structures (red and black curves).
Uniaxial tension simulations along the three principal axes to obtain E and  were performed as described above for the polymerized samples.

Polymerization
The extent of reaction and the largest cluster concentration are plotted as a function of simulation time in Figures 4 and 5, respectively.Also included in Figure 4 is an inset of the fully polymerized large system.The data in the figures indicate that all MD systems reached a converged polymerization level after 500 ps.All samples had a largest cluster percent mass higher than or very close to 90% and pc values within 1.03% of the theoretical value of 66.67% from Equation 2.
Based on these characteristics, all samples were considered to be fully polymerized.
The data in Figures 4 and 5 also suggest that the smaller test systems and the larger system generally polymerized at the same rate, albeit with some scatter at individual times for the smaller test systems.Specifically, the average values of the extent of reaction and percent mass of the largest cluster of the test systems are also plotted in Figures 4 and 5, respectively, and show agreement with the results of the large system.At every 50 ps interval, the percent mass of the largest cluster for the large system falls within one standard deviation relative to the test systems.The same holds true for the extent of reaction, and the scatter is minimal over the whole 500 ps time range (< 1.16%) for all models, which explains why the standard deviation bars are not shown in Figure 4.The average extent of reaction for the 5 test systems after 500 ps was 65.76 ± 0.31%, and 65.45% for the large system.Similarly, the average percent mass of the largest cluster for the 5 test systems after 500 ps was 92.61 ± 2.00 %, while the larger system's largest cluster comprised 89.50% of the system by mass.This agreement between the small test system averages and the large system suggests that smaller systems should only be simulated if multiple replicates are built.
Otherwise, a single larger system should suffice.
The predicted mass density, E, and  values for all models, which agree with known literature values, are shown in Table III [4].The relatively small standard deviations further reinforce the observation that modeling either several replicates of a relatively small systems (about 5,000 atoms) or a single large atomic system (about 30,000 atoms) will yield similar results.

Pyrolysis
For the parametric study of the pyrolysis process, the effects of pressure, temperature, and volatile removal interval on char yield, mass density, and crystallographic dimensions for TS1 are summarized in Figures 6, 7, and 8, respectively.As shown in Figure 6a, an increase in pyrolysis pressure led to an increase in char yield and denser structures, while stacking height and interplanar spacing decreased as pressure increased, as shown in Figure 6b, indicating that partial graphitization was strongly affected by pressure.Physical evidence for partial graphitization is shown by the formation of parallel graphene planes in the inset in Figure 6a for a model pyrolyzed at a pressure of 10 GPa.Increasing pressure increases char yield by preventing volatilization, while also increasing the mass density.Interplanar spacing of aromatic sheets and stacking heights are decreased with higher pressures.
Increases in pyrolysis temperature resulted in more volatilization and, consequently, a lower char yield, as evident in Figure 7a.Furthermore, the data in Figure 7a indicates that denser structures were formed at higher temperatures.Increasing temperature decreased interplanar spacing and stacking height, as shown in Figure 7b.For pyrolysis temperatures of 2500 K -2700 K, 2900 K, and 3100 K, the resulting XRD spectra were amorphous in nature and no crystallite sizes could be calculated.Both high pressures and temperatures led to a decrease in the interplanar spacing, but only higher temperatures increased stacking height, which was due to the artificially high pressures forcing a decrease in the stacking height by compressing the interplanar spacing between graphene layers, as shown in Figure 6b.
Figure 7. Temperature effect on char yield (a), mass density (a), and crystallite sizes (b).Simulations were performed with hydrostatic pressure of 0.5 GPa while removing volatiles every 0.1 ps.Higher temperatures decrease char yield due to higher volatilization while simultaneously leading to more dense structures with higher density.Interplanar spacing decreases at higher temperatures, while stacking heights increase.
In Figure 8, it can be seen that an increase in the removal interval of volatiles led to higher char yield, possibly due to some carbon-based volatiles recombining with the pyrolyzing polymer within the deletion interval.The mass distribution of deleted carbon-based and non-carbon-based volatiles for simulations with a 0.01 ps and a 0.5 ps removal interval are plotted in Figure 9, including labels of common volatiles found within each mass distribution.Although the distribution shows larger masses of carbon-based volatiles for 0.5 ps removal intervals, the sum of the total mass of deleted carbon-based volatiles was 10.4 kDa, compared to 13.5 kDa for the 0.01 ps removal interval simulation.The larger masses observed for the 0.5 ps removal interval simulation correspond to the slower deletion interval, which allowed for larger carbon-based volatiles to form before the deletion interval was reached.For example, at 0.5 ps removal intervals, more unique compounds were volatilized with masses around 40 Da up to the end of the 1 ns pyrolysis simulation than with 0.1 ps removal intervals.
Figure 8. Volatile removal interval effect on char yield and mass density.Simulations were performed at 2700 K with a hydrostatic pressure of 0.5 GPa.An increase in the removal interval increases char yield while lowering mass density.The mass density decreased with higher removal intervals due to a larger number of volatiles being deleted from the simulation at each deletion interval, as shown in Figure 8.If volatiles were removed at shorter intervals, smaller molecules were deleted from the simulation without affecting the configurations of neighboring structures significantly.If, on the other hand, volatiles were removed at larger intervals, larger molecules were deleted from the simulation in each deletion interval, leading to larger pockets of free volume and less dense structures due to the steric hindrance of neighboring structures.Evidence for this is shown in Figure 9, which shows that deleting volatiles every 0.5 ps led to higher average masses of the deleted species.Finally, the resulting XRD spectra for the removal interval study were amorphous in nature.
Based on the results shown in Figures 6 -8, the pyrolysis simulations were performed at a temperature of 3200 K, with hydrostatic pressure of 2 GPa, and a removal interval of 0.1 ps.While it may seem that a temperature of 3200 K may lead to an underprediction in char yield and an overprediction in stacking height and interplanar spacing (Figure 7), applying a hydrostatic pressure of 2 GPa increased char yield by hindering volatilization while simultaneously reducing the crystallite sizes (Figure 6).The removal interval of 0.1 ps was chosen as it was shown that it maintained char yield within a reasonable range (Figure 8).The average evolutions of all models (TS1-5 and LS), as a function of pyrolysis run time using the parameters from the parametric study, are shown in Figure 10(a-c).The average char yield for all 6 samples was 51.19 ± 1.06 %, which is in agreement with known literature values [1].On average, only 0.29% of hydrogen and 0.02% of oxygen atoms remained after carbonization.All the metrics tracked show similar evolutions between the small test systems and the large system, ensuring reproducibility is feasible with the established protocols.
In terms of structural evolution during pyrolysis, the total ring content continued to increase until leveling off at about 700 ps, as indicated in Figure 10b.The amount of 6-membered carbon rings steadily increased at the expense of 5-membered carbon rings and 7-membered carbon rings.Furthermore, as shown in Figure 10c, the amount of sp 2 hybridized carbon also steadily increased at the expense of sp and sp 3 hybridized carbon.Structural and chemical evolution occurred rapidly up to 500 ps, with more gradual changes occurring past the 500 ps mark, indicating that simulations at the established settings should be taken at least up to the 1 ns mark to ensure the simulated pyrolysis process is complete.The mass distributions of deleted volatiles as a function of simulated pyrolysis run time for each model are plotted in Figure 11, and the chemical species volatilized at each mass distribution layer are labeled.All models had many volatiles deleted upon initiation of pyrolysis due to a lack of volatile removal during the heating up step from 300 K to 3200 K. Within 0.1 ps of pyrolysis, TS1, TS2, TS3, TS4, TS5, and LS each had a total of 3.1, 2.9, 3.4, 2.9, 2.6, and 1.9 kDa deleted in one step (excluded from Figure 11), respectively.As the pyrolysis simulations proceeded, the mass distributions of deleted volatiles decreased until about 200 ps, at which point the mass distributions of deleted volatiles throughout the rest of the simulation remained relatively constant.The volatiles removed within the 200 ps range included small volatiles (small hydrocarbons and H2 gas) and large oxygen-containing hydrocarbons, while the volatiles removed after 200 ps were predominantly H2 gas and small hydrocarbons.After 200 ps, TS1, TS2, TS3, TS4, TS5, and LS had volatiles removed with masses up to 39.02 Da, 39.02 Da, 65.01 Da, 54.05 Da, 39.02 Da, and 63.02 Da, respectively.This indicates that, despite the larger number of atoms in the large system, significant volatilization of larger by-products occurs predominantly within 200 ps of pyrolysis simulation time, as evidenced by the similar deleted volatiles mass ranges after 200 ps between all systems.The predicted XRD spectra for all pyrolyzed and equilibrated models are shown in Figure 12.The cyan lines show the full-width half-maximum (FWHM) of each peak associated with the (002) and (100) planes, and Table IV lists the corresponding 2θ values for each peak, their FWHM (β) values, and the crystallographic information (d002, Lc, and La) computed according to Equations 3-5.The interplanar spacing (d002), stacking height (Lc), and crystallite width (La) are consistent between samples, including the large system.There is a relatively small standard deviation, and values agree reasonably well with the experimental characterization performed by Vieira et al. [40], which used a poly(furfuryl alcohol) derived glassy carbon pyrolyzed at 1000 o C for 30 min.
The lack of scatter in XRD spectra in Figure 12 continues to reinforce a common theme presented in this study: simulating several replicates small atomic systems (about 5,000 atoms) was as effective as simulating a single large atomic system (about 30,000 atoms) in obtaining consistent crystallographic dimensions, as indicated by the lack of scatter in the data.The average mass density evolution of the test systems and the large system during pyrolysis are shown in Figure 13.The starting mass density before pyrolysis corresponds to the equilibrated polymerized mass densities shown in Table II.Throughout the pyrolysis simulation run time, the test systems show consistent values with little scatter (max standard deviation is 0.1 g/cc at 396 ps), and the superimposed large system curve falls within the standard deviation bars.The predicted mechanical properties of the models are considerably higher than the experimental values obtained from macroscopic samples [1,50,51].However, nanoscale thin films of amorphous carbon, whose pristine structures may better correspond to a nanoscale MD model of pyrolyzed phenolic resin, as opposed to macroscale amorphous carbon samples containing a relatively high volume fraction of nano-and micro-level voids, have been shown to yield E values ranging from 146 -256 GPa [47][48][49].Specifically, the overprediction of mechanical properties observed in the MD models relative to macroscale experimental values arises from the volatilization process during laboratory-scale pyrolysis in which volatiles form and either remain trapped within the material or diffuse out, consequently producing porosity and cracks in macroscale specimens that reduce the measured modulus relative to the non-porous and crack-free values.The algorithm used here for volatile removal instantaneously deleted volatiles in the MD models, producing a non-porous carbon network due to a lack of volatile entrapment and diffusion, which allows neighboring structures to recoalesce into dense structures.Thus, a more appropriate mass density comparison for the MD models is the skeletal density, which excludes the volume from accessible pores (as determined via helium pycnometry) in mass density calculations.For glassy carbon, these values have been found to range between 1.95 g/cm 3 -2.10g/cm 3 [41][42][43], which agrees with the MD predictions.
One of the benefits of MD is the ability to examine nanoscale features that are difficult to observe in experiments.For example, even though the average mechanical response of the test systems seems consistent with that of the large system, visualization of the molecular structures of the former reveals a higher degree of anisotropy due to preferential aromatic alignment.Such preferential alignment is expected to increase the mechanical response along that orientation.In Figure 14, snapshots of each pyrolyzed system and the corresponding uniaxial E values are shown, which provides evidence of the anisotropic nature of the test systems.Sheets of aligned aromatic structures in the test systems are highlighted by dotted red lines, and the associated uniaxial E values are shown in red font below the images.Using TS1 as an example, it can be seen that most aromatic alignments are along the Z axis, and Ez is consequently higher than Ex and Ey.The same trends are observed for the other test systems in varying orientations.In contrast, the large system's response along each orientation shows less scatter due to the less ordered nature of the aromatic alignments.The observed anisotropy in smaller test systems emphasizes the need for several replicates if small atomic systems are to be modeled for pyrolysis simulations.

Conclusions
Computational MD models were built to simulate the polymerization and pyrolysis processes for a phenolic resin using reactive force fields.For polymerization, the use of the CHON17W ReaxFF parameter set allowed for the successful production of fully polymerized structures with mass densities and mechanical properties matching known experimental values for all 6 samples.The average density of the polymerized structures was 1.24 ± 0.01 g/cm 3 while the average Young's modulus was 3.50 ± 0.6 GPa.The mechanical response of the polymerized large system matched the average response of the test systems, indicating that modeling a single large system is as efficient as modeling several small system replicates in obtaining models which accurately predict properties.
Pyrolysis was successfully simulated with use of the CHON19 parameter set by volatilizing pure and oxygen-containing hydrocarbons and H2 gas, which were removed using an algorithm to delete volatiles with masses below 50 Da.A parameterization study was first performed to determine the optimal combination of pyrolysis temperature, pyrolysis pressure, and the volatile removal interval to produce structures with char yield and crystallographic dimensions similar to the experimentally reported values.Structural metrics and chemical changes were tracked during pyrolysis and simulated XRD was used to analyze the resulting pyrolyzed structures.It was found that pyrolyzing for 1 ns at 3200 K, 2 GPa, and deleting volatiles every 0.1 ps produced desired structures with an average char yields of 51.2 ± 1.1%, average interplanar spacings of 3.8 ± 0.1 Å, average stacking heights of 10.9 ± 0.3 Å, and average crystallite widths of 35.8 ± 2.6 Å, which match reasonably with known literature values.Mass densities matched known experimental measurements of skeletal densities of glassy carbon, while mechanical properties were underpredicted compared to values reported for thin film amorphous carbon samples.As with the simulated polymerization, the average mechanical response of the pyrolyzed test systems matched the mechanical response of the large system, reinforcing that either approach to modeling phenolic-derived carbon matrices is viable.This study emphasizes the need to have several replicates if the chosen approach is to model small atomic systems, as individual test systems exhibited significant anisotropy stemming from the preferential alignment of aromatic structures along a particular orientation.
The established protocols outlined in this study have demonstrated efficacy in predicting resulting pyrolyzed structures and material properties with the use of atomistic reactive simulations for relatively small models, < 5,000 atoms, and short simulation times, < 1 ns.By extending these protocols to other precursor systems of interest, the material selection process can be expedited to find suitable materials that will produce high char yields and desirable properties.Furthermore, atomistic evolution during the pyrolysis process can elucidate the processing parameters that should be tuned for producing specific pyrolyzed structures.Material properties determined at the MD nanoscale can also serve as inputs into microscale models via an Integrated Computational Materials Engineering (ICME) process modeling approach to optimize processing parameters [52][53][54][55][56][57][58].

Figure 2 .
Figure 2. Initial MD triphenol structures and the corresponding crosslinking reactions.Sites with an open valence are indicated in blue.Reactions 1 -4 link the monomers, while reactions 5 and 6 allow for further gelation.

Figure 4 .
Figure 4. Extent of reaction (p) for all samples.All models achieved the critical extent reaction of 66%.Standard deviation bars are included but are too small to be seen due to minimal scatter.

Figure 4 .
Figure 4. Percent mass of the largest cluster for all samples.All models consistently achieved a largest cluster percent mass of around 90%.

Figure 6 .
Figure 6.Pressure effect on char yield (a), mass density (a), and crystallite sizes (b).Simulations were performed at 2700 K for (a) and at 3200 K for (b) while removing volatiles every 0.1 ps.Increasing pressure increases char yield by preventing volatilization, while also increasing the mass density.Interplanar spacing of aromatic sheets and stacking heights are decreased with higher pressures.

Figure 9 .
Figure 9. Mass distributions of common deleted volatiles for (a) 0.01 ps and (b) 0.5 ps removal intervals.Unlabeled data points represent the sum of masses of volatiles already labeled.Quicker removal intervals (a) remove smaller volatiles than slower removal intervals (b), leading to denser structure as shown in Figure 8. Char yield increases at slower removal intervals, possibly due to some carbon-based volatiles recombining with the main carbon network, as evidenced by the larger total mass of carbon-based volatiles deleted for pyrolysis simulations with 0.01 ps removal intervals (13.5 kDa) compared to 0.5 ps removal intervals (10.4 kDa).

Figure 10 .
Figure 10.Structural and chemical evolution as a function of pyrolysis run time, including (a) char yield and elemental mass percentages, (b) evolution of ring formations, and (c) carbon hybridization.Data points represent the average of the five test systems and the large system, and standard deviation error bars are included but are too small to be visible in some cases.

Figure 11 .
Figure 11.Mass distributions of common deleted volatiles during pyrolysis.Unlabeled data points represent mass sums of volatiles already labeled.After about 200 ps of simulation time, mass distributions of deleted volatiles remained relatively constant for all systems, indicating most volatilization occurred before 200 ps.

Figure 12 .
Figure 12.Curve-fit XRD spectra of pyrolyzed and equilibrated models.Full-width halfmaximum is shown as cyan lines.

Figure 13 .
Figure 13.Mass density evolution during pyrolysis.Error bars (in black) are included every 50 ps and represent standard deviation of the test systems average.

Figure 14 .
Figure 14.Visualization of aromatic alignments and uniaxial mechanical responses of the test systems and the large system.
Gallegos: Conceptualization, Methodology, Software, Validation, Formal Analysis, Investigation, Writing -Original Draft, Visualization.Josh Kemppainen: Conceptualization, Methodology, Software, Data Curation, Writing -Review & Editing.Jacob R. Gissinger: Methodology, Validation, Software, Formal analysis, Writing -Review & Editing.Malgorzata Dr. Małgorzata Kowalik received her PhD in Theoretical Physics from Adam Mickiewicz University in Poznan, Poland.Her background is in Statistical Mechanics and modelling of non-equilibrium systems.Her recent research focuses on the structure and dynamics of carbon-based materials, which she studies using non-reactive Molecular Dynamics and the reactive force fields such as ReaxFF.Adri van Duin obtained his PhD degree in Chemistry from Delft University of Technology.He worked from 1997-2002 at the University of Newcastle upon Tyne on molecular dynamics simulations in organic geochemistry.Between 1999-2001 he invented the ReaxFF method -in collaboration with prof.Goddard's group at Caltech.In 2002 joined the Goddard group.In 2008 he joined the faculty of the Department of Mechanical Engineering at Penn State.In 2013 he co-founded RxFF_Consulting, a company that consults on ReaxFF development and application.He has published over 525 journal papers and has distributed ReaxFF to over 2500 research groups.Dr. Kristopher Wise is a materials researcher working in the area of carbon nanotube reinforced structural composites.He received degrees in chemistry from Drury University (BA) and the University of Oklahoma (MS and PhD).During his time at NASA he has applied his background in computational chemistry and materials modeling to the study of small molecule structure and reactivity, polymer structure and dynamics, and the mechanical properties of carbon nanomaterials.His current work is focused on understanding the elastic and fracture properties of carbon nanotube reinforced composites for application in future space vehicles.After earning an Engineering PhD from Michigan Tech in 2007, Dr. G worked as an application developer at AT&T in Middletown, NJ.He returned to Michigan Tech in 2009 as a post-doctoral fellow in Physics.Since January 2011, has been serving as the Director of Research Computing and a Research Associate Professor with teaching, research, and administrative responsibilities.He oversees just about every aspect of campus HPC program, and represents his University in HPC Advisory Council, NSF ACCESS Campus Champion Program and Supercomputing Conference to develop collaborations and partnerships and discover internship, job and funding opportunities for students and researchers.Dr. Odegard is the John O. Hallquist Endowed Chair in Computational Mechanics in the Department of Mechanical Engineering -Engineering Mechanics at Michigan Tech.He is the Director of the NASA Institute for Ultra-Strong Composites by Computational Design, focused on developing next generation composites materials for manned deep-space missions.Before joining the faculty at Michigan Tech, Greg was a researcher at NASA Langley Research Center.He received his PhD at the University of Denver in 2000.His research is focused on computational modeling of advanced material systems.He is the recipient of the NASA Outstanding Public Leadership Medal, is a Fellow of ASME, and an Associate Fellow of AIAA.

Table I .
Parametric study settings for pyrolysis, which consists of varying temperature, pressure, and the volatile removal interval (variable parameter highlighted in bold).

Table II .
Number of atoms before and after pyrolysis and corresponding simulation box sizes after pyrolysis.

Table III
. Mass density and elastic constants of polymerized models agree with experimental values.Average and standard deviation calculations include all models (test systems and large system).

Table IV .
Calculated XRD data and crystallographic parameters from pyrolyzed models, which are in agreement with experiment.

Table V .
Mass densities and mechanical properties of pyrolyzed and equilibrated models.