Simulating the complete pyrolysis and charring process of phenol–formaldehyde resins using reactive molecular dynamics

We examine the mechanism of pyrolysis and charring of large (> 10,000 atom) phenol–formaldehyde resin structures produced using pseudo-reaction curing techniques with formaldehyde/phenol ratios of 1.0, 1.5 and 2.0. We utilise Reactive Molecular Dynamics (RMD) with a hydrocarbon oxidation parameter set to simulate the high-temperature thermal decomposition of these resins at 1500, 2500 and 3500 K. Our results demonstrate that the periodic removal of volatile pyrolysis gasses from the simulation box allows us to achieve near complete carbonisation after only 2 ns of simulation time. The RMD simulations show that ring openings play a significantly larger role in thermal decomposition than has previously been reported. We also identify the major phases of phenolic pyrolysis and elucidate some of the possible mechanisms of fragment formation and graphitisation from the RMD trajectories and compute the thermal and mechanical properties of the final pyrolysed structures.


GRAPHICAL ABSTRACT Introduction
The thermal protection system (TPS) is a vital component of any space-faring re-entry vehicle, protecting onboard crew and equipment from the extreme temperatures experienced during passage through the upper atmosphere at hypersonic speeds [1]. Typically, a TPS consists of an arrangement of materials chosen for their exceptional insulating properties. The means of thermal protection may arise either from these materials' ability to re-radiate heat (non-ablators) or to absorb heat, fragment and off-gas during thermal decomposition (ablators) [2]. Ablative materials limit heat transfer via two distinct mechanisms. Firstly, heat is absorbed by the ablator during the endothermic chemical reactions taking place as it pyrolyses [3] and, secondly, that any gasses and particulates formed during pyrolysis act as a buffering layer, limiting thermal contact between the spacecraft and the superheated gas around it while carrying excess heat away [4]. Ablating materials may be further subdivided into ''charring'' and ''noncharring'', with charring materials having the added advantage that the char itself can act as an insulating layer even after ablation has stopped. Figure 1 shows a schematic representation of the ablation process for a charring ablator. Charring ablator TPSs have been of interest to aerospace engineers since the 1960s [5] and remain a popular choice for spacecraft in the modern day [6].
Testing the performance of such materials during re-entry presents numerous technical challenges, principally that the time and costs involved with launching any spacecraft are vast and that recovery of the ablator is often not possible, making in situ testing mostly infeasible until considerable development and refinement of any prototype TPS has occurred. Alternative methods have been employed such as directing acetylene torches [7], rocket exhaust streams and other high temperature gas jets in an attempt to reproduce re-entry conditions [8]. Plasma arc jets have been used extensively for testing TPSs in ionising environments at very high temperatures and pressures [9]. These arc jets consist of an intensely heated gas or gas mixture being expanded through a convergent-divergent nozzle into a vacuum chamber with the test material located in the path of the plasma jet. The gasses reach supersonic velocities and can subject the sample to considerable shear forces. Heat flux can be up to 80 MWÁm -2 , and pressures can exceed 1.4 MPa [10]. Still, arc jet measurements are far more practical than in-flight testing and the conditions that the ablator is subjected to are a good approximation of the hyperthermal environment experienced during re-entry. Despite this, ablation testing using arc jets still comes with considerable cost and requires the preparation of a physical sample of any candidate material. Furthermore, insights J Mater Sci (2022) 57:7600-7620 into the nature of the chemical reactions that are taking place as the ablator is charring still cannot be confirmed since pyrolysis gases are difficult to collect in these setups. Other methods for investigating the thermal decomposition of the ablator are available such as performing thermogravimetric analysis (TGA) [11] or running pyrolysis gas chromatography mass spectrometry (Py-GCMS) [12], but these cannot reproduce the temperatures of an arc jet gas stream and still also require a physical sample to be synthesised. For all the reasons described above, computational chemistry methods offer an alternative in silico approach to test the performance of ablator materials without the tremendous costs and development times needed to produce, characterise, and investigate experimentally the suitability of candidate ablators. Furthermore, computational methods may be successfully employed in the early stages of the development process to screen materials for those possessing suitable properties that may present opportunities for further evaluation as well as providing complementary information to experimentally collected data.
While computational chemistry methods are less costly, they are not without their own restrictions. Chiefly, any model is limited in the size of the structure being modelled, in the length of time that may be simulated and ultimately in the processing power available. Ab initio molecular dynamics (MD) methods using Density Functional Theory (DFT) can provide an accurate representation of materials by calculating the potentials acting on all atoms from first principles [13], but consume vast amounts of computational resources, limiting their usefulness for predicting the behaviour of bulk materials.
Classical MD approximates interatomic potentials using forcefields which have been parameterised using ab initio calculations and empirical data to approximate the behaviour of the system as accurately as possible, with bonds, angles, dihedrals and impropers explicitly defined [14]. ReaxFF is a reactive implementation of MD where bonding is implicit and bond orders between atom pairs are a function of the interatomic distance [15]. With a parameter set for hydrocarbon combustion, reactive MD has been applied to coal combustion and pyrolysis [16][17][18][19], pyrolysis of polyethylene [20], polycarbonate [21,22] and cellulose fibres [23]. Further to this, modelling of the thermal decomposition of cured epoxy resins [24][25][26] and their mechanical response [27] has elucidated unique insights into decomposition steps and failure modes that would have been challenging to determine experimentally.
Studies of the pyrolysis of phenol-formaldehyde resins (PFR), a popular choice of ablator [28] (structure shown in Fig. 2) using ReaxFF, are very limited. Figure 1 Schematic representation of the ablation process for a polymer-based charring ablator. As heat travels through the material from the outside inward, a 'charring front' forms as the polymer reaches sufficient temperature to thermally decompose. Pyrolysis gasses are evolved as thermal decomposition takes place and these gases escape out of the material and away from the charring front through the porous carbon char layer formed by near complete carbonisation of the polymer ablator. Figure 2 The chemical structure of a fully cured phenolformaldehyde resin showing phenol units connected by methylene bridges at the orthoor parapositions. Each phenol unit may form up to 3 methylene bridges with neighbouring phenols, resulting in network structure.
Jiang et al. [29] simulated the initial stages of PFR pyrolysis and determined that dihydrogen, water, carbon monoxide and acetylene were the major volatile pyrolysis gases formed during temperature ramping of the structure. Several mechanisms for formation of water were also identified but due to computational power limitations, the system size in their simulations was limited to 80 repeat units with a formaldehyde/phenol ratio (F/P) of 1.0 and a simulation time of 20 ps.
Qi et al. [30] performed a comparative study of the early stages of pyrolysis of PFR using DFT, densityfunctional tight-binding (DFTB) and ReaxFF but again the system size of just 132 atoms and simulation time of 1 ps were limited to allow for comparison with the quantum methods. The mechanisms of water formation were found to be the same as those reported by Jiang et al. [29]. Mechanisms of carbon monoxide formation and details of the structural changes taking place as charring began to occur were also revealed [31]. Finally, Desai et al. [32] studied larger composite systems and demonstrated that graphitic precursors began to form as thermal decomposition of the system progressed.
It is apparent that the scaling of computational cost as a function of the system size (* o(N) 3 for DFT methods and * o(N) for ReaxFF) has limited the possibility of reactive QM or ReaxFF simulations being run for long enough to produce large, carbonised structures. Since the progression of ablation is driven by formation and expulsion of volatile gasses, any simulation where no atoms are allowed to leave the simulation box may potentially never produce a fully carbonised solid structure.
In this paper, we successfully demonstrate a novel computational method that increases the performance of ReaxFF simulations of ablation and allow us to obtain multiple trajectories revealing the charring process in detail via targeted removal of light pyrolysis gasses. This method enables us to obtain fully charred structures ([ 98% carbon) after only 2 ns of simulation time for systems sizes up to * 10,000 atoms as opposed to control simulations run without pyrolysis gas removal which did not form any significant carbon clusters after 6 ns. Furthermore, we investigated the reactions occurring within the char using trajectory data and a simple kinetic model and correlated the rate of charring to the production of pyrolysis gasses and other fragments during the simulation. We also report the evolution of the structure of the resin and char, characterise the steps of pyrolysis and highlight key structure-property relationships.
The two properties of particular interest to us are the thermal conductivities and elastic constants of the pyrolysed PFRs. In previous studies computing thermal conductivities of graphene and carbon nanotubes, Tersoff and AIREBO potentials have been used as the potentials but are limited to pure C or CH systems [33][34][35]. Diao et al. have shown the ReaxFF potential to be effective in predicting the thermal conductivities of lower dimensional carbon allotropes such as carbon nanotubes (CNTs) and graphene layers when compared to the prediction produced by using the Tersoff and AIREBO potentials particularly in its ability to characterise the ballistic-diffusive behaviours in agreement with experimental data [36]. ReaxFF potentials also have the advantage of being able to model systems containing oxygen and other elements depending upon the parameter set used in the simulation.
Thermal conductivity may be obtained from non-Equilibrium Molecular Dynamics (NEMD) simulations of the system by artificially applying a temperature gradient across the system as proposed by Mü ller-Plathe [37]. The one-dimensional form of Fourier's Law of Heat Conduction (Eq. 1) relates the heat flux, Q/DT, resulting from a temperature gradient, DT, in a single direction to the thermal conductivity, j, and cross-sectional area, A, of a material [38].
By adding and subtracting known amounts of kinetic energy from two regions within the system and measuring the resultant temperature gradient, we can then calculate the thermal conductivity of the system along the axis of thermal transport. NEMD methods are advantageous over equilibrium MD methods from an efficiency standpoint since the temperature gradient converges much more rapidly than heat flux [37]. We utilised a NEMD method and chose to measure the thermal conductivity along the plane of the graphitic sheets due to the structure of the chars we obtained.
The elastic constants may be obtained using MD by following a method first proposed by Sprik et al. where a set of uniaxial strains are applied to the J Mater Sci (2022) 57:7600-7620 structure and the resulting stresses are measured [39].
For a parallelepiped simulation box, the six unique components of the stress vector, r, may be obtained from the so-called virial expression (Eq. 2) where I and J represent the x, y or z dimension, N 0 is the number of interacting atoms when using periodic boundary conditions, m k is the mass of the kth atom, m kI and m kJ are the components of the kth atom's velocity in the Ith and Jth axis, r kI is the positional vector of the kth atom along the Ith axis and f kI is the force vector of the kth atom along the Jth axis.
By sequentially setting one of the components of strain vector, e, to a preselected value while leaving all others at zero and evaluating the new stress tensor, all the stiffness coefficients in the stiffness matrix may be obtained [39].
For an isotropic material, the stress and strain tensors are related by a stiffness matrix (Eq. 3) which only requires the specification of two independent Lamé coefficients: k and l [40].
It is important to carefully select the magnitude of the applied strain so that the material remains in its linear elastic region. Permanent plastic deformation of the material would lead to unrealistic evaluation of the elastic constants.
The elastic constants that may be obtained in this manner are the bulk modulus, B, (Eq. 4) which is the inverse of compressibility and quantifies a materials resistance to volumetric stress; the shear modulus, S, (Eq. 5) which quantifies the materials resistance to shear stress; the Young's Modules, E M , (Eq. 6) which quantifies a materials resistance to extension in a single direction and Poisson's Ratio, t, (Eq. 7) which quantifies the contraction of a material occurring perpendicular to extension along a single axis.

Methodology
The PFR structures were generated in the commercial software package Materials Studio [41] using a modification of the automated curing method proposed in Hall et al. [42]. To produce each structure, a unit cell was first packed with a mixture of phenol and formaldehyde molecules at a predetermined density and then, distance testing between functional groups was used to periodically initiate crosslinking pseudo-reactions. Additional methods such as reaction-specific probabilities, angle limited interactions and interference from neighbouring molecules are incorporated to model the true nature of the curing reactions more accurately. The structure was relaxed with a short period of MD between reaction steps. After many iterations, a cured polymer structure was obtained [43] (Fig. 3 is an example of the structures generated). We investigated PFR structures with F/P ratios of 1.0, 1.5 and 2.0 to gain an insight into if and how crosslinking may impact thermal decomposition. The parameters of the 3 models used for the simulations are listed in Table 1.
We primarily investigated pyrolysis at temperatures of 1500, 2500 and 3500 K. While some of these temperatures are higher than the 2000 K typically experienced during atmospheric re-entry [44], higher temperatures increase the rate of decomposition sufficiently to make carbonised structures accessible on reasonable timeframes (\ 4000 core hours). We also recorded trajectories and bonding information to examine how heating rate and temperature impacted the progression of the pyrolysis process and the final structure of the char.
RMD was carried out in the open-source Largescale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [45] MD suite using the USER-REAXC package, a parallelised C ? ? implementation of the original FORTRAN ReaxFF code [46]. The parameter set from Chenoweth et al. [47] for combustion of small molecule hydrocarbons was used for the RMD since no specific parametrisations exist for PFR and it was deemed to be suitably general for modelling the pyrolysis process. Additionally, the LAMMPS Python package was used to implement in situ detection of species based on bonding data extracted from LAMMPS during simulation runs.
A timestep of 0.1 fs was chosen due to the high energy nature of the system as tests using larger timesteps triggered issues with lost atoms in LAMMPS and led to the emergence of unphysical behaviours. The decision to use NVT dynamics was taken to better represent the rigidity of a bulk material on the charring front. The caveat to this choice of ensemble is that volume expansion during charring is not possible, so the density of the resultant carbon chars may conceivably be higher than those obtained experimentally.
Relaxed structures were first raised to the target temperature using NVT dynamics over 100 ps with a Berendsen thermostat with a coupling parameter of 100 fs to initiate fragmentation while avoiding the large temperature oscillations associated with equilibrating using the Nose-Hoover thermostat. Our testing determined that even though higher heating rates occur when raising to higher temperatures, the fragmentation of the structure converges independently of heating rate and is related to the temperature for this specific case if the equilibration time is sufficient.
Subsequently the simulation was progressed in 2000 individual cycles. In each cycle, NVT dynamics was run for 1 ps using the Nose-Hoover thermostat with a temperature damping parameter of 100 fs (See Online Resource 1 for a justification of the choice of cycle length). This resulted in a total simulation time of 2 ns. Bonding information was periodically recorded ten times at intervals of 25 timesteps at the end of each run and averaged across all ten intervals. Bonds were defined as pairs of atoms whose average computed ReaxFF bond order was greater than 0.3. Unique molecule IDs were assigned to groups of atoms with contiguous bonds. Grouping and Figure 3 An example of a small, fully cured structure containing * 1000 atoms with phenol units joined by methylene bridges. Both branched (a) and linear (b) polymer chains are visible as well as water molecules eliminated during the condensation pseudo-reactions. The structures are generated with periodic boundary conditions on all axes to represent bulk material. Carbon atoms are grey, hydrogen atoms are light grey, and oxygen atoms are red. determination of species were handled by an external Python script that parsed the bond tables output by LAMMPS and returned molecular formulae for all the fragments present at the current timestep. Each molecule's formula was checked against a list of small molecule fragments to remove, and then, all atoms in that molecule were subsequently removed from the simulation if the formula was present in the list. Finally, the structures were slowly cooled to 298 K at the end of the simulation. The list of small molecular fragments selected for removal were chosen to represent the mixture of common fragments found in pyrolysis off-gasses and to sufficiently accelerate the charring process. Dihydrogen, water, carbon monoxide, carbon dioxide and acetylene were noted products in previous experimental and computational studies [29][30][31][32][48][49][50].
Trajectories and ReaxFF connection tables containing all the atomic positions and all bond orders between atom pairs were output from each simulation allowing for detailed analysis of the pyrolysis process.
The elemental composition of the structure was recorded at the end of each cycle to provide an insight into the degree of carbonisation and the losses of oxygen and hydrogen during the various phases of pyrolysis. The elemental fractions at the beginning and end of each simulation were also noted to compare the relative progression of pyrolysis and carbonisation with different models and pyrolysis temperatures.
Ring counts for small heterocycles and homocycles (n = 3-7) were obtained from the smallest set of smallest rings (SSSR) using an implementation of the algorithm proposed by Paton [51] at 1 ps intervals.
The relative amounts of pyrolysis gases removed from the simulation were recorded at the end of each cycle and summed for each simulation to highlight relative formation of each gas at different temperatures and for structures with different F/P rations providing an insight into the thermal decomposition of the PFR structure.
To compute the thermal conductivity, the model was first oriented so that the graphitic sheets lie on the xz plane and was subsequently replicated 4 times along the z-axis, increasing the number of atoms in each layer to * 200. Replicating the model in the z direction mitigates some of the size effects. Prior testing determined that * 300 atoms per region was sufficient to get a smooth temperature gradient after relaxation and application of thermostats.
Each system was first relaxed and equilibrated to 300 K with a short, 10 ps NVT MD run using a Berendsen thermostat followed by a second 10 ps NVT MD run with the Nose-Hoover thermostat. Following this, the simulation box is divided into 20 layers and non-translational kinetic energy is then added to layer 1 and subtracted from layer 11 using the enhanced heat exchange (eHEX) algorithm [52]. The eHEX algorithm retains the advantages of the original heat-exchange algorithm [53] while reducing energy drift. Kinetic energy is exchanged between the two layers at every timestep at a rate of 0.5 kcal mol -1 fs -1 and the temperature of each layer, T L , is evaluated as the kinetic temperature of all N L atoms in that layer. By calculating the kinetic energy contribution of each atom (Eq. 8), the instantaneous temperature of each of the regions, T L , may be evaluated (Eq. 9) by summing the kinetic energy contributions of each of the N L atoms calculating the proportional instantaneous temperature.
A second relaxation period of 50 ps with the thermostats applied was used to obtain a smooth and stable temperature gradient across the system before a 20 ps production run. By computing and averaging the temperature difference between the regions that have thermostats applied, DT may be obtained. Since Dx, A, Q and Dt are known, the thermal conductivity, j, can then be obtained using Fourier's law as discussed previously.
To compute the elastic moduli of the starting materials and the chars, we employ the all-atom Condensed-phase Optimised Molecular Potentials for Atomistic Simulation Studies (COMPASS)-II forcefield which was chosen due to its extensive parametrisation and validation against ab initio and empirical data obtained from studies of both condensed-phase organic molecules and polymers. COMPASS has a broad range of atom types as well as parametrised interactions for non-bonded forces, bonds, angles, torsions, and hydrogen bonds. Unlike ReaxFF, COMPASS is a non-reactive forcefield but for low temperature, low strain applications, consideration of bond-breaking should not be necessary and the increased number of forcefield types should be advantages in accurately evaluating the potential energy of each perturbation. The commercial software package Materials Studio [41] was used for computation of mechanical properties. We selected a small strain magnitude of 3 9 10 -3 m m -1 which was deemed to be small enough to present plastic deformation of the structures.
Six unique patterns of deformation were independently applied to the undeformed structure in both the positive and negative direction as listed below. Linear scaling of the simulation box along the X, Y, and Z axes, respectively, shearing the box in the XZ plane along the X axis, shearing the box in the XY plane along the X axis and shearing of the box in the XY plane along the Y axis. By recording the change in stress occurring due to the applied strain, the elastic constants were obtained as described previously.
The elastic constants returned use the Voigt-Reuss-Hill approximation. This refers to taking the average of the Voigt bounds which assume that the strain remains constant in a material when subjected to some arbitrary amount of stress and the Reuss bounds which assume that the stress remains constant in a material when subjected to some arbitrary amount of strain. These bounds have been shown to give the upper and lower bounds of the elastic moduli, respectively, and therefore, for an isotropic material Hill proposed averaging the two values to obtain the true elastic constants [54].

Results and discussion
Phases of pyrolysis Analysis of the MD trajectories has shown that thermal decomposition of the structures appears to take place in a regular series of steps (Fig. 4) and result in visually evident structural changes (Fig. 5), with a small number of simple reactions such a ring openings and methylene bridge scissions occurring initially. As the number of fragments increases, the number of distinct reactions occurring rises significantly. It is at this point that the system enters the main off-gassing portion of pyrolysis.
The distinct phases of pyrolysis observed across the course of the simulations are as follows: 1. Firstly, water formed during curing that was already present in the structure is removed postheating during the initial fragment removal cycle as the temperature of the system is now high. Therefore, the rate of pyrolysis of the wet variants of the structure is not any higher than the dry variants in these simulations. In reality, the evaporation of water from the structure may slightly slow the initial pyrolysis as the water absorbs heat before vaporising. 2. Thermal decomposition occurs as the temperature continues to rise with methylene bridge breaking, ring opening and reactions between phenolic hydroxyls and neighbouring structures dominating during the early stages, in line with experimental observations [45]. 3. At this point, the reaction profile becomes more complex as ring openings lead to a large increase in fragmentations. Formation of linear carbon chains leads to structures that tend to lose small molecular fragments (CH, CH 2 CH 3 , C 2 H, C 2 , C 4 , Figure 4 A graph showing the 6 phases of pyrolysis demonstrated by changes in elemental composition and ring formation for the PFR structure with an F/P ratio of 1.5 at 3500 K, each stage is labelled and described in this section of the manuscript. Phase 1 is very short with oxygen and hydrogen decreasing significantly, phase 2 shows high three-membered ring formation, and phases 3 and 4 show increasing six-membered ring formation, decreasing three-membered ring formation and continuing dehydrogenation, loss of carbon has stopped by phase 4. Phase 5 shows increasing numbers of six-membered rings and decreasing amounts of hydrogen, and finally, phase 6 shows mostly stable numbers of six-membered rings and almost no oxygen or hydrogen. C 2 O, etc.) as well as the formation of unsaturated carbons chains through loss of hydrogen. 4. Secondary reactions begin to occur with more frequency between the various fragments. Small molecule products form and are removed in large numbers at this time in the simulation. Due to the reduction in mass, voids in the structure begin to form (Fig. 5) and these lower pressure regions further propagate the formation and reaction of fragments to form more products for removal. This would suggest that while the pyrolysis process is initially thermodynamically controlled, in the later stages it is being kinetically driven.
Transient oxygen-containing heterocyclic compounds with furan and pyran groups form at this stage of thermal decomposition in line with experimental observations from lower temperature pyrolysis [42]. 5. At higher temperatures where most of the oxygen have been consumed, smaller hydrocarbons (C n-H m , n \ * 20, 0 B m \ * 5) and C n carbon clusters combine to form large complex structures with greater degrees of graphitisation. The linear carbon chains present during the main offgassing phase tend to become incorporated into graphitised sheets at this stage (Fig. 5).

Figure 5
Examples of structures formed during RMD for the PFR structure with an F/P ratio of 1.5 at 3500 K displaying the various changes throughout the process. Snapshots are taken from the trajectory at six different times to illustrate the six phases of pyrolysis. At 0 ps, the structure is unpyrolysed. At 50 ps, the structure has lost all its water and rings openings and methylene bridge scissions have occurred, three-membered rings are present in significant quantities at this point and are highlighted in the magnified portion of the image. At 150 ps, much greater fragmentation of the structure has occurred with many small molecules present and at 250 ps continuing dehydrogenation of the structure can be seen along with the development of six-membered rings highlighted in the magnified portion of the image. At 350 ps, the agglomeration of carbon clusters has begun to occur along with the beginning of graphitic sheet formation. At 2000 ps, the structure is extensively graphitised and what remains at the end of the simulation is a carbonised char.
6. In the last phase of charring, dehydrogenation of the remaining structure occurs leaving only carbon sheets composed primarily of six-membered rings ( Fig. 4 and 5) that may form stacked graphene sheets or fullerene like structures via radical polycyclic reactions [46].

Mass loss and elemental composition
The fraction of mass lost and the fraction of carbon in the final char were correlated with pyrolysis temperature for all structures ( Fig. 6 and Table 2). Loss of mass occurred rapidly in the initial 100 ps at and above 2500 K before slowing drastically as oxygen was depleted within the structure. The chars obtained from pyrolysis at 3500 K were [ 99% carbon, dropping to [ 88% at 2500 K and [ 63% at 1500 K. Oxygen was almost completely depleted at 2500 K (\ 1%) and fully removed at 3500 K. At 1500 K, the oxygen content plateaued at [ 4%, suggesting that insufficient energy was present in the system at this temperature for rapid formation of oxygen-containing small molecular products. At all temperatures, the hydrogen content decreased Figure 6 Stacked plots showing the atom composition of the structure throughout the simulation with each colour representing a particular element. The mass curves at 2500 and 3500 K show the characteristic inverted S-shape observed experimentally during pyrolysis under inert atmospheres (See Fig. 6 in Cui et al. [40]) due to the rate of ablation initially rising rapidly before falling as oxygen and hydrogen are incorporated into pyrolysis off-gasses. In all cases, oxygen is depleted more rapidly than hydrogen. At 1500 K, complete loss of oxygen is much less rapid and the S curve is not apparent. throughout the simulations, though only pyrolysis at 3500 K showed near total dehydrogenation (\ 1%). For pyrolysis at 1500 K, the fraction of oxygen lost during the simulation was notably higher as the F/P ratio of the structure decreased. The fraction of hydrogen lost increased for the structures with a lower F/P ratio. Convergence toward a fully carbonised structure does still appear to be occurring at lower temperatures but at a very slow rate, particularly at 1500 K with F/P [ 1.
Char yields were lower at higher temperatures (0.801 at 1500 K, 0.575 at 2500 K and 0.428 at 3500 K) and lower F/P (0.550 at F/P = 1, 0.624 at F/P = 1.5 and 0.630 at F/P = 2) indicating pyrolysis progresses more rapidly for less crosslinked structures at high temperatures. Highly cross-linked structures have been experimentally shown to be more stable at higher temperatures under inert atmospheres [55].

Ring formation
Generally, heterocycles occurred infrequently relative to homocyclic rings, particularly at and above 2500 K. Initially, only six-membered rings are present originating from the phenol units in the resin. During equilibration at the target temperatures, fragmentation occurs predominantly via scission of the methylene bridges between rings, preserving many six-membered rings until product removal is initiated. Then, as pyrolysis begins to progress during the main part of the run, ring openings occur frequently and new rings in a range of sizes emerge. The number of six-membered homocycles serves as a useful indicator of the level of graphitisation occurring during pyrolysis and are the most abundant type of ring present at the end of each simulation. Four, five and seven-membered homocycles are also present but in much lower quantities (Fig. 7).
Ring formation is less apparent at lower temperatures with less carbonised structures displaying a greater propensity to form linear carbon chains with Figure 7 Plots showing counts of n = 3-7 rings sampled at the end of each cycle. In all cases, six-membered rings are initially destroyed by ring opening reactions in the early stages of thermal decomposition. At temperatures of 2500 K and above, formation of six-membered homocyclic rings dominates the later stages of pyrolysis due to the formation of polycyclic fragments and later, graphitic sheets. Five-membered rings are also formed, particularly at 1500 K and when the F/P ratio of the structure is greater than 1.0. There are notable spikes in the number of three-membered heterocyclic rings in the first 500 ps of pyrolysis at 3500 K. some crosslinking, usually resulting in much larger rings. The number of small three-membered rings stabilises at * 75 at 1500 K (compared to * 40 at 2500 K and * 25 at 3500 K), and the majority of these are transient cyclopropane rings that decompose as the simulation progresses. Five and sixmembered cyclic ethers are somewhat stable at 1500 K, generally forming after disassociation of a terminal hydroxyl group at the end of a dangling carbon chain and subsequent reaction of the terminal oxygen with an unsaturated carbon. The cyclic ethers occur more frequently and tend to remain more stable in structures with higher F/P ratios, with a few remaining intact for over a nanosecond in structures with an F/P of 2. Qi et al. [30] also reported the presence of five-and six-membered oxygen containing heterocycles, but these were more transient at the higher temperatures used for their MD.

Fragmentation and agglomeration
The trajectories generated by LAMMPS included the bond order calculated by the reax/c package [46] and were output at the end of each cycle. This allowed to track the fragments present at each cycle and analyse the changes in the sizes of fragment present (Fig. 8). We were also able to record cumulative sums of some of the most common fragments which were relevant to pyrolysis gas formation ( Table 3). Some of the observed fragments were small stable molecules such as formaldehyde, dihydrogen, methane, ethene, while others were highly reactive transient species like hydrogen and hydroxyl. The treatment of atoms as particles coupled with the high energy environment allows such small species to exist briefly after fragmentation even before they encounter other fragments and recombine.
At lower temperatures, decomposition initially results in many fragments with similar numbers of atoms clustered into discrete groups. This is likely the Figure 8 Heatmaps showing the distribution of fragment sizes during the pyrolysis processes, molecular data were recorded at 1 ps intervals with the number of atoms in each fragment indicated on the y-axis. Brighter pixels signify a larger number of fragments of that size at that specific timestep in the simulation as shown by the colour scale to the right. Structures with lower F/P ratios are notable for the discrete grouping of fragment sizes because methylene bridge scission is more likely to produce isolated molecules. Notable peaks cantered around cluster sizes of 14, 28, 42, 56 and 70 atoms equate to one, two, three, four and five phenol-formaldehyde units, respectively, as well as variants with different numbers of methylene groups and hydrogens. Higher temperatures tended to produce large numbers of small fragments very rapidly with a broad range of cluster sizes.
result of the scission of methylene bridges since it is understood to be a major mechanism of initial fragmentation. Also, the ReaxFF bond disassociation for a bond order 1 carbon-carbon bond is 173.67 kcalÁmol -1 vs. 228.03 kcalÁmol -1 for a bond order 2 carbon-carbon bond.
Splitting at the methylene bridge gives rise to small fragments such as C 7 H 6 O (a 14 atom unit is equivalent to one phenol-formaldehyde unit), C 14 H 12 O 2 (a 28 atom unit is equivalent to two phenol-formaldehyde units) and C 21 H 18 O 3 (a 42 atom unit is equivalent to three phenol-formaldehyde units) as well as variants with added or lost methylene groups or hydrogens.
This fragmentation appears far more prominent for the PFR structure with an F/P of 1.0 since intuitively; the likelihood of a single bond breaking event leading to the formation of an isolated unit is more probable for polymer chains with fewer crosslinks.
In contrast to Qi et al. [30], significant amounts of ring openings are observed at all temperatures, though higher temperatures promote a greater frequency of ring openings. They appear to be one of the central reactions which ultimately lead to the formation of smaller fragments by further decomposition of the linear molecules formed. At 3500 K, the ring opening reaction is the predominant decomposition mechanism and appears to be hastened by the removal of pyrolysis gases, indicating that higher density structures slow thermal decomposition initially (Fig. 9).

Pyrolysis gas removal
Three initial fragmentations appear to be the main reactions driving the production of all pyrolysis gases in our simulations, and these fragmentations occur readily at all pyrolysis temperatures. Calculations would suggest these fragmentations have the lowest apparent activation energies in our model. The removal of hydrogen from aryl carbons, methylene bridges and hydroxyl groups is one of the earliest decompositions observed in all simulations (Fig. 10).
During the initial thermal decomposition phase, hydrogen can be lost from methylene bridges, hydroxyls, or aromatic rings (Eq. 10-12), but methylene bridges tend to be the most common site while loss from aryl carbons is the least common, reflecting the trend of bond energies. Bond disassociation energies were calculated using ReaxFF for comparison.
ArH Ar þ H DH ¼ 132:15 kcal mol À1 ð11Þ Dihydrogen is the most abundant pyrolysis product at and below 2500 K, and it is still produced abundantly at 3500 K. This is likely due to hydrogen loss occurring much more readily at elevated temperatures along with a kinetic element whereby hydrogens are more likely to migrate and interact at elevated temperatures. Simple combination of  Figure 9 Stacked plots showing the production of pyrolysis gases targeted for removal throughout the pyrolysis process. High temperatures cause much of the off-gassing to occur very rapidly within 100 ps of equilibration, whereas at 1500 K off-gassing occurs very slowly and at a steadier rate. Structures where the F/P is greater than 1 are noticeably more inert throughout pyrolysis, especially at 1500 K. The volume of gas produced is noticeably lower at 1500 K than for higher temperatures and decreases as the F/P ratio of the structure increases at 1500 K. This correlation does not occur at higher temperatures.
disassociated monatomic hydrogen is chiefly responsible for dihydrogen generation.
Water was produced at greater rates during high temperature pyrolysis, likely due to an increased availability of oxygen. As previously described, oxygen tends to be more likely to stay in heterocycles at 1500 K. Peak formation of water occurs earlier than other fragments, in line with both experimental and computational observations [56,57].
Hydroxyl fragments formed due to carbon-oxygen bond breaking in the phenolic hydroxyl groups (Eq. 13) can combine with hydrogen to form water (Eq. 14).
Hydrogenation of hydroxyl group from a b-hydrogen or free hydrogen fragment followed by carbon-oxygen bond breaking can lead to production of water (Eq. 15), but this was observed less frequently.
Carbon monoxide formation is significantly greater at higher temperatures but interestingly does not correlate with F/P, possibly because it is thermodynamically driven and limited by the number of hydroxyl groups present. Formation occurs via one of two major pathways. The first is a step that involves disassociation of the hydrogen from a phenolic hydroxyl group (Eq. 16) followed by a ring opening directly adjacent to the hydroxy group. The terminal carbonyl group is then cleaved from the structure (Eq. 17).
The second is via formation of a formyl fragment (Eq. 18) and subsequent dehydrogenation (Eq. 19), and these fragments are particularly abundant at 3500 K. Their formation usually involves a ring opening adjacent to a hydroxyl and a subsequent migration of hydrogen.
The amount of carbon dioxide produced in all simulations is considerably lower than the amount of carbon monoxide and does not correlate with temperature or F/P. Hydrocarboxyl fragments occur much less, but loss of hydrogen can lead to carbon dioxide formation (Eq. 20-21) analogously to loss of hydrogen from a formyl fragment.
Formation by ring opening analogous to carbon monoxide is less plausible as it would require subsequent combination with another fragment or reaction with a neighbouring hydroxyl group. Free oxygen is rare in the structure, so it is unlikely that a mechanism involving a combination of oxygen and carbon monoxide is present, but it is possible that the combination of carbon monoxide and a hydroxyl fragment can lead to the formation of a hydrocarboxyl fragment. It has been suggested that oxidation of hydroxyl groups to carboxylic acids is a major mechanism of carbon dioxide formation during hydrocarbon oxidation and that this may occur with phenolics as well [58]. Carbon dioxide formation pathways require further investigation given that it is observed in significant quantities in experimental studies.
Acetylene formation is markedly higher at F/P = 1 and increases markedly with pyrolysis temperature. Sivaramakrishnan et al. [59] studied the high temperature thermal decomposition of benzene and suggested that formation of acetylene from benzene occurred in a single step following the loss of an aryl hydrogen (Eq. 22).
It is therefore possible that an analogous mechanism is responsible for acetylene production from the thermal decomposition of the phenol unit.
Observations of the trajectories also suggest that acetylene fragments may be cleaved from the end of a linear chain with a terminal vinyl group and dehydrogenated (Eq. 23).
This observation would tally with the observation that ring openings occur more frequently for structures with less crosslinking as, intuitively, ring openings of aromatic rings with fewer attached groups would be more likely to result in the formation of a linear chain from which fragments could be cleaved. Acetylene production does correlate with ring opening at 1500 K, but at and above 2500 K where near complete breakage of all initial rings occurs, it is correlated with temperature. Large quantities of ethynyl and methylidyne fragments are present along with dicarbon. Dicarbon and ethynyl may combine with hydrogen to form acetylene (Eq. 24-25).
Methylidyne fragments may also combine with one another (Eq. 26) to form acetylene via biradical combination. The increased frequency of ring opening at F/P = 1 may either directly increase the available linear hydrocarbon fragments and/or indirectly accelerate fragmentation.
Thermal conductivity of phenolic chars The NEMD method of computing thermal conductivity proved successful in being able to predict the thermal conductivity of the carbonised chars obtained from pyrolysis at 3500 K along the plane of the graphitic sheets. The thermal conductivity values we obtained are notably lower than those reported for phenolic chars (1.153 W m -1 K -1 ) [60]. Experimental and computational studies of chars obtained from pyrolysis of coal and polymer network materials reveal that the structure of these chars appears to be graphene sheets and small graphite flakes suspended in a complex carbonaceous foam in which other structures such as nanotubes and nanobuds and regions of amorphous carbon present in smaller quantities [61][62][63]. The thermal conductivity of chars is likely reduced by the significant buckling of the graphitic sheets in the chars for the structures with an F/P ratio of 1.0 and 1.5. The significantly lower thermal conductivity of the char obtained from the structure with an F/P ratio of 2.0 is likely due to the much more disordered structure of that char with fewer graphitic sheets present and larger amorphous regions (Table 4). It is likely that larger models with larger, less buckled, and more interconnected regions of graphite would increase the thermal conductivities of the chars (Fig. 11).

Elastic moduli of phenolic resins and phenolic resin chars
In all cases, the shear, bulk and Young's moduli for the un-pyrolysed material are significantly lower than those of the char obtained from pyrolysis, indicating that the starting polymer is significantly softer than any of the resultant chars. The un-pyrolysed structure with an F/P ratio of 1.5 resists deformation that most of all the starting structures, indicating that this ratio results in an optimal number of cross-links to form a rigid structure. The charred structures tend toward the much greater in-plane elastic constant of graphite [64]. For the F/P ratio 1.0 and 1.5 structures, the elastic moduli increase at 1500 K and decrease again at 2500 K before finally increasing at 3500 K.
The variation of elastic moduli as a function of the temperature results from the temperature dependence of the fragmentation rate, which happens less readily at low temperatures, where some cross-linking also occurs. Likewise, at 3500 K polycyclisation reactions between carbons clusters form large, rigid graphitic sheets. For the structure with an F/P ratio of 2.0, the elastic moduli are the highest of the chars and those of the char at 3500 K are lower than those of the char at 1500 K. Visual observation the F/P ratio 2.0 char at 3500 K suggests less ordered layering of the structure than at lower F/P ratios which may account for the lowered resistance to deformation. Further studies with more models and more pyrolysis runs in a greater range of temperatures that could be averaged would be likely to yield clearer trends in mechanical properties (Table 5).

Conclusion
We have successfully demonstrated that the targeted removal of pyrolysis gasses during RMD simulations at 1500, 2500 and 3500 K utilising the ReaxFF forcefield and hydrocarbon oxidation parameter set allows us to simulate the ablation process of PFR ablators for TPS at the atomic scale. We have also been able to obtain chars at all pyrolysis temperatures. Those obtained at 3500 K were found to be almost fully carbonised, heavily graphitised structures and only required 2 ns of simulation time (* 4000 core hours) when 6 ns simulations without product removal yielded no major carbon clusters. Furthermore, we have compared how the degree of crosslinking in the structure affects the processes of thermal decomposition. We determined that ring-openings occurred much more frequently in our larger scale than had previously been reported with small scale RMD simulations of phenolics. We were also able to identify 6 major phases of phenolic pyrolysis from the RMD trajectories. Mechanisms for the formation of water from disassociation of phenolic hydroxyl groups were elucidated from inspection of the trajectories as was the pathway of formation for acetylene fragments via ring-opening and subsequent cleavage of a dangling vinyl group and the mechanism of carbon monoxide formation from dehydrogenation of phenolic hydroxyls followed by ringopening. Carbon dioxide formation was found to be significantly lower than reported in other studies, probably due to removing carbon monoxide from the simulations, since this reduces the rate of formation of carbon dioxide from carbon monoxide. We were able to predict the elastic constants of the un-pyrolysed polymers as well as the polymer chars, determining structure-property relationships and dependencies of pyrolysis temperature. It was discovered that the initial polymerisation and final graphitisation phases increased the stiffness of the material and that the structure with an F/P ratio of 1.5 was the most resistant to deformation. Char structures with less buckled graphitic sheets with more even layering were also determined to be stiffer. Finally, we were able to predict the thermal conductivities of the fully carbonised chars obtained at 3500 K and determine that they were lower than those reported for phenolic ablators chars, possibly owing to the smaller buckled graphitic structures present. In this work, we have sought to determine the utility of this method in simulating the ablation of bulk polymeric materials exposed to high temperatures and obtaining fully carbonised char structures from relatively short MD simulations as well as examining the overall progression of the pyrolysis process. It is important to note that the specific effects of this method on the rates of thermal decomposition processes have not been quantified here and require further study. We have been able to demonstrate effective proof-of-concept for this method of simulated ablation using the open-source LAMMPS package and a general-purpose parameter set for ReaxFF which should allow us to expand to larger structures as well as other polymers and composite materials in future studies.

Declarations
Conflicts of interest The authors declare no conflicts of interest.
Supplementary Information: The online version contains supplementary material available at http s://doi.org/10.1007/s10853-022-07145-4. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licen ses/by/4.0/.