Molecular Simulations of Interface-Driven Crosslinked Network Formation and Mechanical Response in Composite Propellants

TL;DR

This study uses molecular simulations to understand how different components in composite propellants interact and affect mechanical properties, offering insights for better propellant design.

Contribution

The paper introduces a computational framework combining coarse-grained simulations and reactive force fields to study interfacial interactions in composite propellants.

Findings

A two-step reaction mechanism was identified at the AP interface involving self-polymerization and crosslinking.

Optimized parameters led to a ~20% improvement in tensile modulus and strength of the propellants.

Key factors like HTPB chain length and IPDI content significantly influence mechanical performance.

Abstract

Composite solid propellants, which serve as the core energetic materials in aerospace and military propulsion systems, necessitate tailored enhancement of their mechanical properties to ensure operational safety and stability. A critical challenge involves elucidating the interfacial interactions among the multiple propellant components (≥6 components, including the polymer binder HTPB, curing agent IPDI, oxidizer particles AP/Al, bonding agents MAPO/T313, plasticizer DOS, etc.) and their influence on crosslinked network formation. In this study, we propose an integrated computational framework that combines coarse-grained simulations with reactive force fields to investigate these complex interactions at the molecular level. Our approach successfully elucidates the two-step reaction mechanism propagating along the AP interface in multicomponent propellants, comprising interfacial self-polymerization of bonding agents followed by the participation of curing agents in crosslinked network formation. Furthermore, we assess the mechanical performance through tensile simulations, systematically investigating both defect formation near the interface and the influence of key parameters, including the self-polymerization time, HTPB chain length, and IPDI content. Our results indicate that the rational selection of parameters enables the optimization of mechanical properties (e.g., ~20% synchronous improvement in tensile modulus and strength, achieved by selecting a side-chain ratio of 20%, a DOS molar ratio of 2.5%, a MAPO:T313 ratio of 1:2, a self-polymerization MAPO time of 260 ns, etc.). Overall, this study provides molecular-level insights into the structure–property relationships of composite propellants and offers a valuable computational framework for guided formulation optimization in propellant manufacturing.

Figures

Full text

1. Introduction

Composite solid propellants (CSPs) serve as the primary energetic materials in modern aerospace and military propulsion systems due to their exceptional energy density, tunable combustion characteristics, and superior performance [1,2,3,4]. These materials are highly filled elastomeric composites, typically composed of a polymeric binder (fuel component) and a majority of solid fillers, such as oxidizer particles and aluminum powder, uniformly dispersed within the matrix. However, the extreme manufacturing and working conditions impose stringent mechanical demands on CSPs, requiring them to withstand thermal and pressure loads while resisting aging-induced degradation [5,6,7]. Additionally, CSPs must endure complex mechanical stresses arising from shocks, vibrations, and pressure fluctuations. Under such demanding conditions, tensile and shear stresses may induce microscopic defects within CSPs, potentially leading to significant performance degradation or even catastrophic failure [8,9,10]. Therefore, a fundamental understanding of the microscopic structures of CSPs at the molecular level and their correlation with mechanical properties is critical for ensuring the safety and reliability of propulsion systems [11,12].

Usually, CSPs contain many propellant components, such as a polymer binder (hydroxyl-terminated polybutadiene, HTPB), solid oxidizer (Ammonium perchlorate, AP), plasticizer (Dioctyl sebacate, DOS), curing agent (Isophorone diisocyanate, IPDI), and bonding agent (Tris-1-(2-methylaziridinyl) phosphine oxide, MAPO, Boron trifluoride triethanolamine complex, T313), etc. Thus, the formulation selection of polymer binders and additives is crucial for controlling the mechanical performance of CSPs [13,14,15,16,17,18]. The polymeric binders not only provide structural integrity due to its high content and molecular weight but also ensure strong interfacial adhesion with oxidizer particles and other functional additives. The HTPB binder undergoes curing reactions with isocyanate groups, forming a polyurethane elastomer that exhibits optimal stiffness, high elongation at break, excellent thermal stability, and superior water resistance. Bonding agents, such as aziridines and their derivatives (e.g., MAPO), as shown by previous studies [19,20], can adsorb onto AP surfaces through hydrogen bonding and undergo ring-opening polymerization, forming an interfacial layer with high modulus and tear resistance. Furthermore, these bonding agents react with isocyanate-based curing agents, significantly strengthening the interface between AP particles and the polymer matrix. These synergistic reactions contribute to an integrated crosslinked network, enhancing the overall mechanical robustness of CSPs [1,20,21]. On the other hand, interfacial defects—particularly debonding under tensile loading—are among the primary causes of mechanical property degradation in CSPs. Such defects often arise from voids near the interface due to spatial heterogeneity, leading to crack initiation and propagation [14,21,22]. Collectively, these intricate microscopic factors, primarily governed by interfacial interactions and reactions between polymer binders and additives, play a crucial role in determining the mechanical behavior of CSPs [23,24,25,26]. Nevertheless, a comprehensive molecular-level understanding of interfacial properties and their influence on mechanical performance remains a significant challenge [8,27,28].

In this work, an integrated computational framework combining coarse-grained molecular dynamics (CGMD) simulations and reactive force field methods is developed to elucidate the microscopic mechanisms governing the mechanical properties of multicomponent CSPs. This framework reproduces the sequential two-step reaction process during CSP production well: (1) surface-catalyzed self-polymerization of the bonding agent MAPO on AP particles, followed by (2) diverse crosslinking reactions mediated by curing agents. Moreover, through tensile simulations, we find that microscopic defects preferentially nucleate at interfacial regions, while the conformational changes of HTPB chains and interfacial interactions collectively determine the system’s mechanical response. In addition, based on these results, we demonstrate that CSP formulations can be further optimized to achieve an approximately 20% simultaneous improvement in both elastic modulus and maximum tensile strength by properly changing critical factors, such as the self-polymerization time, chain length and side-chain ratio of HTPB, IPDI content, etc. Our study provides molecular-level insights into the interfacial interactions within CSPs and the formation mechanisms of crosslinked networks and offers fundamental guidance for designing propellants with enhanced mechanical performance.

2. Materials and Methods

2.1. Coarse-Grained Models for Simulations

All MD simulations were performed using the Galamost software package version 1.40, (Institute of Chemistry, Chinese Academy of Sciences, Beijing, China), a GPU-accelerated program designed for efficient coarse-grained molecular dynamics simulations, under isothermal–isobaric (NPT) ensemble conditions [29]. The Martyna–Tobias–Klein (MTK) barostat and Nosé–Hoover thermostat were utilized to maintain a constant system pressure of 1 bar and temperature of 300 K [30]. The selected temperature corresponds to typical experimental sample preparation and storage conditions (i.e., room temperature). However, we note that temperature could significantly affect the mechanical properties of analogous material systems [8,19]. Short-range interactions were modeled using the LJCoulombShiftForce module, which ensures a smooth decay to zero at the cutoff distance. Figure 1a illustrates the coarse-grained molecular representations of the simulated system components, including the polymer binder HTPB, curing agent IPDI, bonding agents MAPO/T313, and plasticizer DOS (note that oxidizer particle AP is shown in Figure 1c). The all-atom-to-coarse-grained mapping protocols and corresponding force field parameters are provided in detail in Figure S1 and Table S1. Specifically, the HTPB chains were initially constructed with a polymerization degree of 50, where chain length, functionality, and side-chain reactivity were implemented as adjustable parameters for subsequent investigations into mechanical properties. The IPDI molecule was coarse-grained into a three-bead model with asymmetric arm lengths to reflect the differential reactivity of its benzene ring termini [31]. Both MAPO and T313 were represented by four-bead models with enhanced interaction potentials to account for their strong hydrogen-bonding characteristics with solid particles. Due to the relatively large size of the AP particle, it was simplified as a solid surface in our simulations. Specially, the AP surface was constructed by arranging CG beads along the (001) crystal orientation, with a surface thickness of 2 nm.

In addition, in our work, the component molecules were basically modeled using standard Martini bead types (e.g., P1, P3, C3, C5, and Q0) [32,33,34,35,36]. Moreover, due to the requirements of the reactive force field (which necessitates distinct naming conventions for different molecule types), certain bead designations were modified: P1 beads were re-named as IP or HP1, and C3 beads were uniformly designated as CR (see Table S1). Given the system’s specificity, some molecular structures could only be approximated using Martini beads with similar properties. This approximation might potentially affect component density distributions. To compensate, we proportionally scaled up the σ values for all bead pairs while specifically strengthening the P3-Q0 interactions (unique to the bonding agent and AP, respectively) to preserve experimental density fidelity (see Figure S2). Note that the extreme complexity of our multicomponent system made it challenging to achieve simultaneous convergence for all molecules by using classical coarse-graining approaches. But the emerging machine learning approaches (e.g., machine learning-informed energy renormalization) may be promise for addressing such challenges to pursue more accurate parameterization schemes [37,38].

2.2. Simulations System Setup and Workflow

Based on experimental observations [20], the simulation workflow consists of three key stages (Figure 1b): (1) component mixing, (2) MAPO self-polymerization (Reaction-1), and (3) IPDI-involved network crosslinking reactions (Reaction-2). Experimental evidence suggests that the bonding agents MAPO and T313 can ultimately coat the AP surface and undergo self-crosslinking under appropriate conditions to form a tear-resistant layer, though this process requires prolonged mixing durations. Accordingly, in our simulation protocol, we directly established a three-layer composite model with systematically arranged components as the initial mixing configuration, followed by 100 ns NPT equilibrium simulations of the component mixtures. Figure 1c (visualized using OVITO software, version 2.9.0 [39]) (OVITO GmbH, Darmstadt, Germany) illustrates this hierarchical structure, consisting of (1) a blended binder layer (top) containing IPDI and HTPB for crosslinked network formation, with uniformly dispersed DOS plasticizer; (2) an interfacial reaction layer of MAPO and T313, which serves as the site for subsequent MAPO self-polymerization; and (3) a rigid solid particle layer (bottom).

The subsequent stages involve the formation of a crosslinked network through a two-step reaction process (Figure 1d,e). The first step (Reaction-1, simulation time varies from 20 ns to 340 ns) corresponds to MAPO self-polymerization, where reactions occur between arms of different MAPO molecules (see Figure 1c). Based on experimental observations [20,40], these reactions are spatially confined to the AP interface via a constraint algorithm that deactivates reactive sites beyond a critical distance (~1 nm) from AP particles. The second step (Reaction-2, 600 ns simulations) involves IPDI-mediated network crosslinking through reactions between IPDI, HTPB, MAPO, and T313, ultimately yielding a fully integrated network structure (Figure 1d). Following MAPO polymerization, IPDI initiates reactions with HTPB terminal hydroxyl groups, T313 functional groups, and ring-opened MAPO arms, forming a dense, three-dimensional network. Moreover, in addition to the primary hydroxyl-terminated reactions, the vinyl-1,2 side-chain isomer of HTPB—containing reactive double bonds—also participates in crosslinking, further increasing network connectivity. Notably, previous studies have demonstrated that such side-group incorporation substantially influences both HTPB’s properties and the resulting polyurethane’s mechanical characteristics, including its tensile strength, elongation at break, and hardness [41]. In addition, it should be noted that MAPO participates in both reactions, whereas T313 is only involved in Reaction-2 through its interaction with the curing agent IPDI.

In our simulations, these reactions are modeled using the stochastic polymerization and reaction algorithm developed by Liu et al. [42,43,44]. In this approach, an active center reacts with nearby monomers within a defined reaction radius with probability Pr; once the reaction event is achieved, the active center transfers to the newly incorporated monomer, while the original center becomes deactivated. This “active center transfer” mechanism governs the polymerization process, with each active site permitted to participate in only one reaction event to prevent multiple reactions. In our cases, for optimal simulation efficiency while preserving physical accuracy, all reaction probabilities are standardized at 2%, except for the vinyl-1,2 side reaction, which was set to 0.002% to reflect its high activation energy barrier [45,46].

2.3. Mechanical Characterization Through Tensile Simulations

Uniaxial tensile simulations are conducted along the z-axis (normal to the solid surface) under NVT ensemble conditions. As shown in Figure 2a, the equilibrated crosslinked system was subjected to a constant strain rate of 5 × 10^−7^ s^−1^, allowing for sufficient simulation time for complete elastoplastic deformation and defect formation.

The mechanical response is characterized by monitoring the stress–strain evolution (Figure 2b), from which three key parameters are derived through curve fitting, including the linear strain, maximum tensile strength, and the corresponding strain. Notably, the simulation protocol excludes chemical bond rupture events. System size independence is verified through scale-up tests (Figure S3), demonstrating consistent mechanical properties across different simulation dimensions.

3. Results and Discussion

Our simulation approach provides a molecular-scale resolution of the CSP production process, enabling fundamental understanding of performance-determining mechanisms for targeted mechanical optimization. The formation of crosslinked networks through the two-step reactions is the core of CSP production. In Reaction-1 (MAPO self-polymerization), the temporal evolution of bond formation follows a power-law growth trend, with decreasing monomer concentration leading to progressively slower bond formation kinetics (Figure 3a). This reaction primarily occurs near the AP interface. Figure 3b shows the 2D spatial distribution of MAPO bonding positions, revealing non-uniform patterns where MAPO self-polymerization occurs preferentially in specific localized regions (indicated by purple and red dots in Figure 3b). These reactive zones are interspersed with areas containing T313 and partially unreacted MAPO (white regions in Figure 3b). Note that such irregular distributions of MAPO bonding positions are commonly observed across different systems (Figure S4), probably arising from the stochastic nature of the single-armed polymerization reaction. This heterogeneous spatial arrangement may create localized stress differences that could further initiate mechanical failure during tensile deformation. Notably, while the reaction probability affects the rate of MAPO self-polymerization, it ultimately leads to less than 10% variation in the total number of bonds formed after sufficient reaction time (Figure S5). In Reaction-2, IPDI reacts with HTPB, MAPO, and T313 to form an integrated network structure. Figure 3c reveals that the bond formation kinetics still follow a power-law trend, with most reactions reaching plateaus rapidly (<30 ns). Notably, the IPDI-vinyl (HTPB) reaction, though slower (saturating within 200 ns), ultimately forms more bonds than the IPDI–OH (HTPB) reaction. A critical finding is IPDI’s role as a mobile crosslinker: it diffuses through the upper layer to participate in all secondary reactions, as evidenced by the z-directional distribution of IPDI bonding sites (Figure 3d). The dominant reaction peak occurs near the AP surface (around z = 1.5 nm), corresponding to interfacial crosslinking between IPDI and the pre-formed 2D MAPO/T313 network. This interfacial reaction zone is directly linked to the defect formation mechanisms, as discussed in the subsequent section.

Following crosslinked network formation, we conduct tensile simulations to evaluate the system’s mechanical response (Figure 2). Especially, defects nucleate at sufficiently large strains, ultimately leading to material rupture. As shown in Figure 4a and Figure S6, defect initiation occurs preferentially at the AP particle–polymer interface, corresponding to a characteristic stress drop in the stress–strain profile (Figure 4a). Our analysis reveals a two-stage defect evolution process mediated by interfacial interactions: (1) initial bonding agent (MAPO/T313)-AP contact maintenance through strong adsorption (red dashed line, Figure 4b), followed by (2) interfacial debonding. Figure 4c quantifies this process through two complementary metrics: (i) the mass density near the AP interface shows a linear decrease with strain, while (ii) the bonding agent–AP contact number exhibits a biphasic behavior—a gradual initial decline followed by rapid decrease beyond a critical strain threshold (~5%)—leading to rapid crack propagation. Notably, this transition point coincides precisely with the maximum tensile strength in Figure 4a, indicating the onset of catastrophic interfacial failure (Figure S7). This hierarchical failure mechanism indicates that pre-existing interfacial defects from the crosslinking step may influence the composite’s ultimate mechanical performance.

Our computational framework, which integrates the atomic-resolution capabilities of MD with a parametrically tunable propellant model, provides a way to explore the CSP formulation space for optimizing mechanical properties, including the tensile strength and elastic modulus. Herein, we mainly focus on three classes of molecular-level design parameters: (1) the self-polymerization duration in Reaction-1 governing reaction kinetics, (2) HTPB chain length and side-chain vinyl group reactivity associated with polymer architecture, and (3) compositional factors such as MAPO/T313 ratio and DOS/IPDI numbers. Notably, optimizing these parameters is a non-trivial task due to their complex influence on mechanical properties. As a benchmark approach, we herein employ a one-factor-at-a-time optimization strategy, mirroring experimental protocols by varying individual parameters step by step while maintaining others at baseline values. The optimal value for each parameter is determined by either maximizing key mechanical properties or achieving an optimal balance among them (see Figure S8). Figure 5 illustrates the performance enhancements achieved relative to the initial formulation. Our analysis reveals Reaction-1 duration and IPDI content as primary determinants of the tensile strength, whereas the elastic modulus shows strongest dependence on the HTPB chain length. Generally, these three parameters affect the molecular architecture of the polymer binder and the critical interfacial interactions within the composite system.

As a representative case, we analyze the impact of the MAPO self-polymerization duration in Reaction-1 on the mechanical properties. Figure 6a,b demonstrate the changes of tensile strength and elastic modulus with this duration. Below 140 ns, both properties increase monotonically, correlating with progressive MAPO network formation (as shown in Figure 3a). However, beyond this time point, distinct behaviors emerge: the tensile strength exhibits non-monotonic variations (characterized by rapid decline followed by partial recovery), while the elastic modulus continues increasing until reaching 260 ns. The optimal balance between these two properties occurs at 260 ns, where these competing trends reach an advantageous equilibrium.

Furthermore, the MAPO self-polymerization duration significantly influences HTPB chain conformation, as evidenced by the evolution of the mean end-to-end distance (⟨d_ee_⟩) and radius of gyration (⟨R_g_⟩), shown in Figure 6c. It is found that ⟨d_ee_⟩ initially increases linearly before entering an oscillatory regime, while ⟨R_g_⟩ exhibits monotonic growth until reaching a plateau at ~300 ns. The optimal 260 ns time point coincides with maximum ⟨d_ee_⟩ and elevated ⟨R_g_⟩, indicating enhanced HTPB chain extension and spatial expansion that contribute to improved mechanical performance. Moreover, we observe a structural transition in HTPB from amorphous entanglement to quasi-crystalline alignment with parallel chain orientation during MAPO self-polymerization (Figure 6c insets) [26,47,48]. This transition is further supported by flexion angle distributions (Figure 6d), where the characteristic angle (defined in Figure 6d inset) [49] increases substantially with polymerization time. However, this crystallization phenomenon presents a mechanical paradox: while improving stiffness metrics, it simultaneously reduces material ductility, as evidenced by decreased elongation at break (Figure S8). The competing effects of conformational ordering (enhancing rigidity) and crystallinity (causing embrittlement) highlight the critical need for balanced parameter optimization in propellant design.

In addition, as a key crosslinking agent, IPDI actively participates in network formation. Therefore, the effect of IPDI number on the mechanical properties is also analyzed. Our analysis reveals its non-monotonic concentration-dependent effects on mechanical properties (Figure 7a,b): both the tensile strength and elastic modulus initially increase with IPDI content, peak at N = 1200, then undergo decline with further addition. This behavior correlates directly with interfacial enhancement at the AP particle–polymer interface. As illustrated in Figure 7c, interfacial reactions between IPDI and bonding agents intensify with increasing IPDI content until reaching saturation at N = 1200. Higher IPDI concentrations improve connectivity between the binder and coupling agent layers, thereby enhancing the elastic modulus. Figure 7d compares stress–strain curves in the linear elastic region for systems with N = 800 and 1200. The inset demonstrates that a higher IPDI content (N = 1200) leads to a denser and more uniform distribution of reactive IPDI, which strengthens interfacial resistance to tensile deformation under strain. All these results suggest that optimal mechanical performance requires balanced IPDI content to maximize interfacial connectivity.

In general, our simulations reveal that while individual compositional variables modulate specific aspects of the crosslinked network, their synergistic interactions govern the emergent structural properties, including network topology, connectivity, and stiffness. Critical formulation parameters—particularly MAPO self-polymerization duration, HTPB chain architecture, and IPDI stoichiometry—collectively determine the system’s mechanical performance and load-bearing capacity. Especially, using the optimized parameters such as referred in Figure 6a,b and Figure 7a,b, we develop a hierarchical parameter optimization protocol to systematically improve mechanical performance. While competing effects among intrinsic factors, rational parameter selection nevertheless achieves significant enhancements, as demonstrated in Figure 5. Our optimized formulation (MAPO: T313 ratio = 1:2, MAPO self-polymerization time = 260 ns, DOS molar ratio = 2.5%, HTPB chain length = 30 monomers, vinyl side-chain ratio = 20%, and IPDI molar ratio = 59%) yields simultaneous 20% improvements in both the tensile modulus and strength relative to baseline systems (Figure 5). It is also worth noting that formulation optimization presents a significant challenge for CSPs. Herein, we only demonstrate a preliminary application of our integrated computational framework for quantitatively analyzing the effects of key parameters on mechanical properties that are difficult to characterize experimentally. In the future, the implementation of advanced optimization strategies, such as genetic algorithms, Bayesian optimization, and machine learning approaches, could offer a promising pathway for efficient identification of optimal CSP formulations [50,51].

4. Conclusions

In summary, we have developed an integrated computational framework combining coarse-grained molecular dynamics simulations with reactive force field methods to elucidate the mechanisms of self-polymerization and crosslinking network formation in CSPs. This approach successfully reproduces the two-step reaction process and characterizes the mechanical response under strain conditions. Most importantly, we uncover fundamental structure–property relationships that govern the performance of CSPs. For instance, we demonstrate that the duration of MAPO self-polymerization is a critical control parameter, with an optimal duration of approximately 260 ns enhancing HTPB chain extension while promoting favorable molecular ordering for synergistic mechanical improvement. Additionally, the IPDI content exhibits a non-monotonic influence on mechanical properties, achieving a peak at N = 1200 due to optimal interfacial bridging effects. Through systematic optimization of parameters—including the MAPO/T313 ratio (1:2), DOS content (2.5%), HTPB chain length (30 monomers), vinyl side-chain ratio (20%), and IPDI molar ratio (59%)—we achieve a concurrent 20% enhancement in key mechanical properties compared to conventional formulations. Overall, the developed framework provides unparalleled molecular-level insights into CSP structure–property relationships and establishes a paradigm for the rational design of high-performance propellants through computationally guided formulation optimization.