Molecular Models for Solid Propellant

Introduction and Project Overview

Chemical rocket propellant is a class of materials which provides thrust through the combustion of oxidizer and fuel components. The two most industrially-significant types of rocket propellant are liquid and solid propellants; the oxidizer and fuel chemicals are both in the liquid phase for liquid propellants, while for solid propellants
both components are in the solid phase. Solid phase propellants are utilized due to ease of storage, near instant readiness, and relatively low cost. Liquid phase propellants can have better performance, are throttleable, and have the ability to restart, but can suffer from complications due to handling and storage as well as requiring preparation, such as cryogenic cooling, which introduces significant latency between storage and utilization. For military applications, the advantages of solid propellent often out weigh the better performance of liquid propellant.

Solid energetic compounds used in solid propellants are usually brittle, and thus unable to handle the shock waves present during typical applications. To improve the material properties, the solid energetic compounds, in the form of small crystalline particles, are typically mixed into a viscous pre-polymer which is subsequently cured into a rubbery elastomer possessing the desired material properties. Industrially, a standard solid composite propellant consists of 60–72% ammonium perchlorate (AP) as the oxidizer, up to 22% aluminum powder as metal fuel, and 8–16% polymeric
binder such as hydroxyl-terminated polybutadiene (HTPB)-based polyurethanes, all by percent mass. Due to the unsaturated characteristic of HTPB, over time the cured polymer will undergo thermo-oxidative degradation, continue to cross-link, become brittle, and accordingly lose the desired mechanical properties, leading to potential failure upon use. Even though solid propellants have long shelf lives, their material properties can degrade over the course of a couple decades and the rate of this decay is dependent upon the storage conditions. Once sufficiently aged, there is a non-insignificant risk of device failure upon use due to the formation of defects related to the increase of cross-links in the polymeric binder over time. In particular, microvoids and microcracks associated with embrittlement and debonding of the polymeric binder from the particulate crystalline oxidizer promote uneven burn. Thus, after some empirically determined shelf life devices containing the AP-HTPB propellant are decommissioned.

The unit cell of ammonium perchlorate. Element color representations: red for oxygen, black for chlorine, dark blue for nitrogen, pale blue for hydrogen.
Isophorone diisocyante (IPDI, top) curing agent for hydroxyl-terminated polybutadiene (HTPB, bottom).

Despite the abundance of knowledge on these materials, the relationship between particular degradation mechanisms and the resulting material properties are not well understood. Experimentally, evaluating an aging process that spans decades in a feasible manner relies upon accelerated aging techniques. These techniques involve exacerbating the conditions favorable to increased cross-linking such as higher-than ambient temperature and higher concentration of oxygen, water, or any other contaminant suspected of inducing cross-linking. The material properties of the propellant can then be sampled over the course of months or a couple years as opposed to decades. Unfortunately, accelerated aging studies inherently do not mimic ambient conditions and thus could favor reaction pathways that would otherwise be relatively uncommon, thereby reducing confidence in the extrapolated aging process. Computer simulations are a potential way to overcome these experimental limitations. Thus far, there has been little to no molecular modeling research on the aging of polymeric binders typically used in solid propellants, and what has been reported is mostly bond dissociation energies used to speculate on decomposition mechanisms. To circumvent temporal constraints, a series of multiscale simulations are envisioned to study the aging reaction mechanisms and rates, the role of interfaces in the aging process, and the mechanical properties of the propellant with different extents of cross-linking.

A quantum–classical framework is intended for these studies to inspect the thermodynamic and kinetic quantities associated with the degradation mechanisms as well as the resulting macroscopically observable material properties obtained as a function of the quantity and type of different degradation products present in the propellant. Quantum mechanical simulations by Mark Gordon’s group at Iowa State University and the Ames Laboratory are being designed to calculate pre-exponential factors and rate constants for a representative set of degradation mechanisms pertaining to both inert and oxidative atmospheres. Edward Maginn’s group at the University of Notre Dame is using classical molecular dynamics to evaluate material properties of degraded propellant via responses to uniaxial tensile and compressive deformation. Eventually, the parameters describing the thermodynamics and kinetics of individual degradation mechanisms will be incorporated into a kinetic Monte Carlo scheme that will be used to provide a model that estimates the composition and material properties of the AP–HTPB propellant as a function of time. It is likely that the atomistic models will need to be mapped onto coarser-grained representations in order to access the long timescales that are characteristic of the multimodal mechanical responses that occur during the deformation and relaxation of elastomeric materials. The first step in realizing this multiscale framework is to identify or develop suitable atomistic models for the species constituting conventional solid rocket propellant. My work to date has been wholly focused on the development of classical models for the AP and HTPB-based polymeric binder.

A Fully-Atomistic Model of HTPB

The OPLS-AA force field was chosen for the fully-atomistic model of HTPB. To validate the OPLS-AA force field for adequately modeling the polymer, comparisons of cis-1,4-polybutadiene (PB) using OPLS-AA to a well-established united-atom model and experimental results were made.

Data comparing volumetric and single-chain structural quantities between our OPLS-AA based force field and the well-established united-atom force field of Tsolou et al. for linear cis-1,4-polybutadiene. From left to right: (left) specific volume of our model shows comparable performance to Tsolou et al. when comparing to the horizontal dashed line representing an extrapolation to infinite chain length via experimental data; (middle) the mean-square end-to-end distance (<R2>) as a function of chain size agrees decently with Tsolou et al., but shows somewhat increased scaling of <R2>; (right) the extrapolated characteristic ratio for an infinite chain length (C) for our model is 5.5±0.3 compared to 4.8±0.5 for Tsolou et al. and 4.9±0.2 from experiment, which shows that our model is close to the expected behavior but the chain dimensions are somewhat larger than what is experimentally observed. Future improvements for an atomistic model could attempt to modify the rigidity of the torsional degrees of freedom to achieve smaller single-chain structural characteristics.

Seeing that the OPLS-AA force field adequately reproduced the volumetric and single-chain structural properties of cis-1,4-polybutadiene, a simple model was developed for HTPB. HTPB is industrially produced by free-radical polymerization which results in a wide distribution of chain sizes and number of hydroxyl functional groups per chain. The average chain weight is 2,800 g/mol and the average number of hydroxyl groups per chain is around 2.5. To approximate this behavior, a mostly cis-1,4 microstructure was used to produce difunctional and trifunctional chain models that when used in a 1:1 ratio succeed in reproducing the average molecular weight and hydroxyl functionality.

Model difunctional HTPB structure with a molecular weight of 2251.77 g/mol (a). Model trifunctional HTPB structure with a molecular weight of 3349.60 g/mol (b).

The initial configuration of the packed HTPB polymer melt in the simulation box was produced using the proprietary Amorphous Builder Module in the MAPS software package by Scienomics. Subsequently, the polymer was equilibrated using microsecond-long molecular dynamics runs at constant temperature and pressure (T=333 K, P=1 atm). To evaluate the equilibration progress, the running average of the single-chain structural characteristics was monitored over time.

Mean-square radius of gyration running average calculated separately for difunctional chains and trifunctional HTPB chains.
Mean-square end-to-end distance running average calculated separately for difunctional and trifunctional chains.

Additionally, the end-to-end vector autocorrelation function was monitored to establish the degree to which the end-to-end vectors decorrelated over the microsecond long trajectory.

End-to-end vector autocorrelation function for difunctional and trifunctional HTPB chains.

Yet another metric for improving confidence of equilibration is to assess ensemble average of the center-of-mass displacement of the polymer chains to understand how much phase space sampling is occuring. By moving more than their respective radii of gyration on average, both the difunctional and trifunctional chains are shown to have diffused a significant amount during the simulation.

Average displacement over time for difunctional and trifunctional HTPB chains. Comparisons to ⟨Rg21/2 indicate that the chains have moved at least their own length scale over the course of the simulations, on average.

When visualizing the equilibrated simulation cell, significant hydroxyl association was seen to occur. This behavior was unexpected and is not significantly discussed in the existing experimental literature. The infrared spectroscopy data shown in the literature for HTPB suggests that hydroxyl association does occur, but this feature has not been adequately discussed in prior publications.

Association of hydroxyl groups in the HTPB polymer melt. Polybutadiene backbone not shown for clarity in the righthand image.

To quantify the association of the hydroxyl groups, the hydroxyl radial distribution function was computed and the oxygen-oxygen curve indicates that 3.5 Å is a reasonable cutoff distance between oxygen atoms for defining between hydroxyl groups. Using that definition for association, the statistical distribution of hydroxyl associative cluster sizes was calculated.

The radial distribution function between OH pairs is used to inform a 3.5 Å cutoff between O atoms for defining OH association (left); quantification of the OH cluster size distribution showing that 60%+ of OH groups are associated with another OH group on average (right).

The physical insight gained from these results is that there is a not an even spatial distribution of hydroxyl groups in the polymer melt. This has the potential to influence the distribution of cross-links when the HTPB is cured with the IPDI cross-linking agent. The spatial distribution of cross-links in a cross-linked material impacts the resulting material properties of the cross-linked material. Thus, it is important for a molecular model to adequately represent the spatial distribution of cross-links in the HTPB-IPDI material. Unfortunately, no experimental data exists that characterizes this aspect of the material. Assuming that the hydroxyl groups can serve as an approximate surrogate for the location of cross-links, a cross-linking methodology was developed to maximally retain the OH spatial distribution after cross-linking.

First, hydroxyl pairs had to be identified for cross-linking with an IPDI molecule. To select these pairs, the shortest distance oxygen-oxygen pair was identified, these oxygen atoms removed from the available pool, and then the identification procedure repeated until all 250 hydroxyl groups were paired.

OH pairs are selected such that the shortest distance possible O–O pairs are preferentially selected for cross-links in an attempt to utilize the OH aggregation behavior in the melt as a surrogate for where the IPDI cross-links should be introduced. This is in contrast to a random cross-linking or annealing procedure that would significantly neglect the equilibrium spatial distribution of the reactive OH groups. It is unknown how the spatial distribution of the OH groups would change if equilibrating the unreacted IPDI in the HTPB melt.

This methodology effectively preserves most of the hydroxyl groups’ associative behavior but has the drawback of producing increasingly distant hydroxyl pairs as the extent of conversion increases (when wanting to cross-link all of the hydroxyl groups). When introducing the IPDI cross-links, the IPDI molecules are directly inserted into the simulation box at the midpoint between reacting hydroxyl pairs and the bonding topology is immediately changed to reflect the cross-linked material. All cross-links are added simultaneously. This procedure results in a significant amount of atomic overlap and very energetic bonds. To relax these unstable configurations in a numerically stable manner, the Fast Inertial Relaxation Engine (FIRE) was utilized to minimize the energy of the system. This minimization algorithm succeeded in reducing atomic overlap and producing nominal-length bonds in the relaxed configurations with just 20,000 iterations.

The minimum, maximum, and mean bond lengths of the urethane cross-linkages as a function of the number of minimization iterations. Results shown are for a single system cross-linked using a 1.00 NCO/OH ratio.

Furthermore, the vast majority of the HTPB atoms were relatively unperturbed during the minimization procedure. This was quantified by computing the displacement of the carbon and oxygen atoms before and after the minimization algorithm. Only a handful of HTPB atoms moved more than 4 Å and the atoms that did were likely connected to the few hydroxyl pairs that were initially separated by large distances.

The relative frequency at which unreacted HTPB carbon atoms are displaced by a certain amount during the minimization procedure (a). The relative frequency at which reacted hydroxyl oxygen atoms are displaced by a certain amount during the minimization procedure (b). Both plots are generated by averaging the five configurationally independent systems cross-linked using a 1.00 NCO/OH ratio.

The systems were cross-linked at various isocyanate to hydroxyl ratios (NCO/OH = 0.88, 0.92, 0.96, 1.00) and then subsequently quenched to various temperatures (T in K= 233, 253, 273, 293, 313) and then NPT molecular dynamics was performed for 100 ns to relax stresses remaining after the FIRE minimization algorithm. After the 100 ns of relaxation, the various systems were subjected to uniaxial compressive and tensile deformations to assess the mechanical properties of the HTPB-IPDI cross-linked material. A true strain rate of 107 s-1 was utilized until a 3% strain was achieved at which point the strained dimension was held constant while the lateral dimensions could fluctuate to relax the generated stress. Using this procedure, the temperature dependence of the stress–strain behavior, Young’s modulus, and the stress relaxation behavior was computed.

True stress–time behavior at various temperatures from tensile deformation simulations utilizing a 107 s-1 strain rate. Vertical dashed line denotes where simulation was held constant at 3% strain (a); same as previous plot but utilizing compressive deformations (b).
True stress–strain behavior with linear fits for determining Young’s modulus at various temperatures. Very large moduli are primarily attributed to the extremely fast computation strain rate utilized.
True stress relaxation behavior after tensile deformation fit to a single exponential decay term (a); same as previous plot but utilizing compressive deformations (b). The lack of multimodal decay in the stress over time indicates that only the fastest mode of stress relaxation is being adequately sampling during the time scale of these simulations. Much longer time scales are needed to better characterize the stress relaxation behavior.

Simulations at the various NCO/OH ratios were indistinguishable within the noise of the results; physical experiments contrast starkly with this observation by showing that there is an extreme sensitivity of the mechanical properties to the NCO/OH ratio. The simulated insensitivity coupled with the simulated true stress relaxation behavior indicates that the current time scale for assessing the mechanical properties via simulation is inadequate for probing the multimodal stress relaxation that occurs for the material. Simply put, the atomistic model for the polymeric binder is too computationally expensive to run simulations long enough to adequately gather statistics on the modes of stress relaxation that occur very slowly. Future work is needed to train heavily coarse-grained model of HTPB that reproduces the atomistic details shown here, specifically the associative behavior of the hydroxyl groups.

Evaluating an Atomistic Model for Ammonium Perchlorate

An existing Class II force field was taken from the literature (Zhu et al.) to describe the bonded and non-bonded interactions of ammonium perchlorate. The existing force field was not well-validated in the original publication and what validation was performed was for the material at 10 K. To gain confidence that the force field performs well for describing the physical properties of AP, several material properties were calculated at 10, 78, and 298 K and compared to experimental data. First, the lattice parameters and density of AP were computed and compared to experiment and the original force field publication as a function of the cutoff radius used for interatomic interactions.

Comparison between the experimental and simulated values for lattice constants a, b, and c, and density as a function of VDW cutoff at 10 K and 298 K. Data attributed to Zhu et al. are from 10 K MD simulations.

The results generally show that my calculated values do not exactly match the original force field publication’s values, with the exact reason for this discrepancy being unknown. Compared to experiment, my simulated lattice constant a was ~3% too small, lattice constant b was ~8% too large, and lattice constant c was within ~0.4% of the expected value. When inspecting the hydrogen bonding behavior in the system, the simulated N–H(3)•••O(3) mirror symmetric hydrogen bonds were found to have a much larger angle than what was experimentally observed. It is possible that this larger hydrogen bonding angle and the large error in lattice constant b are both manifestations of an underlying inadequacy of the force field since this hydrogen bond type lies in the plane of the hydrogen bonds.

The two N–H(3)•••O(3) mirror symmetric hydrogen bonds in the 10 K MD simulations. Red is oxygen of perchlorate, blue is nitrogen, and white is hydrogen.

The pressure–volume behavior of the simulated ammonium perchlorate generally reproduces the experimentally observed trend but is overly repulsive resulting in a larger volume for a given pressure.

The pressure–volume data for both simulation and experiment at 298 K.

This behavior is similarly seen when viewing the data in terms of the bulk modulus of the AP crystal as a function of pressure.

The 298K bulk modulus values calculated from simulation compared to the values obtained from the third-order Birch-Murnaghan equation of state fitted to experimental pressure–volume data.

To evaluate the vibrational behavior of the ammonium perchlorate ions, the 298 K infrared spectrum obtained via experiment was compared to the molecular vibrations obtained from simulation. The stretching and scissoring modes in the simulation were identified for both the perchlorate and the ammonium ions by utilizing bond and angle autocorrelation functions. It was generally observed that the simulated model well reproduces the general locations of the peaks, but the N–H stretching mode was somewhat faster in the simulations than what was observed by experiments.

MD simulated and experimentally obtained IR spectra for AP at 298K (a). Stretching frequencies identified for both types of ions in the 298 K MD system (b). Scissoring frequencies identified for both types of ions in the 298 K MD system (c).

To further investigate the vibrational performance of the AP model, the temperature and time step dependence of the N–H stretching mode were assessed. At low temperature, the simulated AP shows three different N–H stretching modes which is expected for the three different hydrogen bonding environments that are experimentally reported. In comparison to Raman spectroscopy data, the simulated N–H stretching mode is significantly blue-shifted and the difference between hydrogen bonding environments is much smaller. The overly similar, simulated N–H stretching modes indicate that the N–H bonds in the force field model are not adequately susceptible to their local chemical environment to respond in a substantially different way, as seen in experiments. When lowering the time step from 1.0 fs to 0.5 fs, the N–H stretching mode is red-shifted, showing that a 1 fs time step is likely too large of a time step if concerned with the exactly vibrational details of the system. Still, the vibrational behavior in the simulated model is qualitatively consistent with experimental observations.

N–H stretching frequencies at 10 K, 78 K, and 298 K calculated from MD with a time step of 1.0 fs compared to Raman data (a). N–H stretching frequencies at 10 K and 298 K calculated from MD with a time step of 0.5 fs are red-shifted compared to a 1.0 fs time step (b).

Experimentally, the perchlorate ion does not rotate between 10–298 K and only vibrates in place while the ammonium ion experiences the onset of nearly free rotation somewhere between 78–298 K. Using a second-order Legendre polynomial to describe the rotational behavior of the ions, the simulated ammonium perchlorate qualitatively agrees with the temperature trend observed in physical experiments.

Second-order Legendre polynomials describing the rotational behavior of the perchlorate anions (a) and the ammonium cations (b) from MD simulations at 10 K, 78 K, and 298 K. Black lines denote an empirical double-exponential fit for each temperature of both data sets.

Holistically, the modified AP model presented here is well-suited for reproducing trends in the material properties of AP between a temperature range of 10–298 K. The error in the quantitative reproduction of experimentally observed values by the model is generally ~5% or greater; however, this is not necessarily a significant issue as classical molecular dynamics simulations are frequently more appropriate for revealing trends and explaining physical phenomena rather than accurately providing quantitative estimates.

Current and Future Work

The development and validation of the models presented in this work are critical first steps in establishing classical molecular models for simulating solid propellant. Currently, validation of an AP-HTPB force field is being performed to match experimentally observed interfacial tension values between AP and HTPB. Once complete, interfacial phenomena important to describing the debonding process of HTPB from AP will enabled. In the future, it is important to develop a model for HTPB that is coarse-grained enough that long time scale simulations are possible for assessing mechanical properties on time scales accessible by experiment so that the simulation results can be corroborated. Eventually, the classical models developed as a result of this work will have degradation products introduced to the polymeric binder to assess how the material properties change as a function of how much thermo-oxidative degradation has occurred.


Support for this work was provided by the Air Force Office of Scientific Research under Contract AFOSR FA9550-18-1-0321. This work was supported in part by high performance computer time and resources from the DOD High Performance Computing Modernization Program as well as the Center for Research Computing (CRC) at the University of Notre Dame. The Materials and Process Simulation (MAPS) code was used for generating initial HTPB configurations. The support of Scienomics is acknowledged (Scienomics SARL, Paris, 2021).