Mesoscale Damage Evolution, Localization, and Failure in Solid Propellants Under Strain Rate and Temperature Effects

TL;DR

This paper studies how solid rocket propellants break under high-speed stress and temperature changes, using a new modeling approach.

Contribution

A cohesive finite element method framework is developed to quantify thermomechanical responses at the microstructural scale.

Findings

Damage evolution in propellants is strongly influenced by strain rate and temperature.

The proposed CFEM model accurately predicts damage progression and mechanical responses.

Neural network-based methods improve parameter identification for cohesive models.

Abstract

High-energy solid propellants are multiphase engineering materials, whose mechanical behavior is predominantly governed by the characteristics of embedded crystalline particles. While microstructural influences have been extensively examined, quantitative correlations between microstructure and macroscopic mechanical properties remain underexplored. This work develops a cohesive finite element method (CFEM) framework to quantify the thermomechanical response of high-energy solid propellants at the microstructural scale. The analysis focuses on impact loading at strain rates ranging from 103 to 104 s−1, accounting for large deformation, thermomechanical coupling, and microcrack-induced failure. Damage evolution under impact conditions was evaluated using a combined neural network-based inverse identification method and a three-dimensional cohesive finite element model to determine temperature-dependent bilinear-polynomial cohesive parameters. Results demonstrate a strong dependence of the propellant’s mechanical behavior on both strain rate and temperature. Validation against experimental data confirms that the proposed temperature-sensitive CFEM accurately predicts both damage progression and macroscopic mechanical responses.

Figures

Full text

1. Introduction

Solid propellants, such as hydroxyl-terminated polybutadiene (HTPB) and nitrate ester plasticized polyether (NEPE), exemplify highly filled particulate composites. These materials primarily consist of rigid energetic constituents—including ammonium perchlorate (AP), hexogen (RDX), and hexanitrohexaazaisowurtzitane (CL-20)—dispersed in a lightly cross-linked viscoelastic binder [1,2,3]. As the core material in solid rocket motors, the mechanical performance of these propellants is critical to maintaining structural stability under operational loads [4,5,6]. Accurate characterization of their mechanical response under diverse loading conditions remains a key research objective. Despite substantial efforts to develop constitutive models, significant challenges persist due to nonlinearities introduced by evolving damage mechanisms [7,8,9]. Variations in temperature affect the binder’s properties, thereby altering the kinetics and morphology of damage accumulation. This temperature dependence influences both the initiation and evolution of fine-scale damage, impacting the material’s time-dependent mechanical behavior [10,11,12]. Therefore, investigating the temperature-dependent mechanical response of solid propellants is essential for advancing predictive modeling.

Computational methods—particularly finite element and meshfree approaches—have been widely adopted to simulate failure processes in solid propellants. Cohesive zone models (CZMs) have been applied to capture interfacial failure phenomena [13,14,15]. During loading, substantial mismatches in stiffness between particles and the binder matrix result in stress concentrations at interfaces, leading to interfacial debonding and microcrack formation [16,17,18,19]. Increased particle packing density tends to elevate porosity levels introduced during fabrication. Under mechanical loading, these pores often coalesce with microcracks, promoting damage propagation and culminating in macroscopic crack formation, which significantly degrades mechanical performance. While earlier studies primarily emphasized qualitative or phenomenological observations, a more comprehensive understanding necessitates quantitative analysis of damage progression. To this end, fracture mechanics principles have been integrated into numerical frameworks to examine the mechanical and damage characteristics of propellants and particle-reinforced composites [20,21,22]. Recent advances in mesoscale damage modeling have facilitated a deeper exploration of damage mechanisms and their macroscopic implications. These efforts incorporate both surface-level characterization techniques—such as optical microscopy (OM) [23] and scanning electron microscopy (SEM) [24]—and internal damage tracking via high-resolution X-ray computed tomography (CT). Hu et al. [25] analyzed the effects of micropore expansion on bulk mechanical properties, establishing links between microstructural features and macroscopic behavior. Subsequent studies employed CZMs for quantitative simulation of damage evolution, further bridging the gap between microscale phenomena and macroscale mechanical degradation in propellants [15,26,27]. However, existing studies still have two key limitations. First, most cohesive zone models (CZMs) remain primarily two-dimensional [28,29,30]. Such models are only applicable to preliminary analyses involving geometrically simple structures, planarized loads, and low precision requirements. Second, the majority of existing methods for obtaining cohesive zone element parameters rely on the Kalman filtering algorithm [31], which depends on accurate system models and noise assumptions. This renders it constrained in scenarios characterized by complexity, nonlinearity, or unknown models.

To address the aforementioned issues, this study investigates the compressive damage behavior of solid propellants under high strain rates and varying temperature conditions using a framework based on a three-dimensional thermo-mechanical cohesive zone model. An artificial neural network (ANN) model is employed to determine the parameters of the continuum finite element model (CFEM). The performance of the neural network system and the effectiveness of the proposed model are evaluated by comparison with experimental results. Damage mechanisms across a range of strain rates are examined and quantitatively characterized. Compared with existing works, the method presented in this paper enables more accurate prediction of the damage evolution process of propellants under different temperature and strain rate conditions.

2. Model and Methodology

The constitutive model of each component was defined. A linear elastic model was applied to crystalline phases, while the binder was modeled as a linear viscoelastic material. To simulate damage in the crystals, binder and interfaces, zero-thickness cohesive elements governed by a bilinear traction–separation law were incorporated throughout the structure. Parameters for the bilinear model were determined using a hybrid inverse optimization strategy.

2.1. Viscoelastic Model for Polymer Binder

The generalized Maxwell model (GMM) was used to represent the linear viscoelastic behavior of the polymer binder. The model, expressed in integral form, was derived from the Boltzmann superposition principle, as shown in Equation (1) [32]:

$$[eqn]$$

where $[eqn]$ and $[eqn]$ are the instantaneous stress and strain tensors, respectively, $[eqn]$ is the equilibrium relaxation modulus tensor, and $[eqn]$ is the relaxation modulus tensor. $[eqn]$ is expressed using Equation (2):

$$[eqn]$$

where $[eqn]$ and $[eqn]$ are the bulk and shear moduli, respectively.

The relaxation functions were modeled using a Prony series expansion. Both the shear modulus G(t) and the bulk modulus K(t) were expressed as follows:

$$[eqn]$$

where G^(n)^ and K^(n)^ represent the shear and bulk stiffnesses, respectively, $[eqn]$ represents the shear and bulk relaxation times, G^(∞)^ and K^(∞)^ are the long-term elastic shear and bulk moduli, respectively, and N is the number of Maxwell elements.

2.2. Thermomechanical Cohesive Zone Element

The method proposed by Xu and Needleman [33] was employed to insert cohesive elements between physical elements to simulate crack propagation. As illustrated in Figure 1, the cohesive zone comprises virtual cracks and a plastic zone, with real cracks forming beyond the cohesive zone. No stress was assumed on surfaces adjacent to real cracks, and the stress field σ(x) at any point within the plastic zone depended on the distance from the crack tip. Stress evolution followed a traction–separation law, expressed as a function of the cracking displacement S(δ). When the stress exceeded a critical threshold, the stiffness of the virtual crack diminished, eventually forming a real crack.

To analyze damage and failure mechanisms under thermal expansion at fine scales, the thermomechanical response of cohesive elements was investigated. This analysis incorporated thermal–mechanical coupling effects to capture the initiation and evolution of microscopic defects. The failure mechanism progressed through three stages: thermally induced stress, damage initiation, and crack propagation. Thermal mismatch stresses arose at interfaces due to differences in the coefficients of thermal expansion among material phases. The thermal strain for each component was calculated using the following:

$$[eqn]$$

where $[eqn]$ is the thermal expansion coefficient, and $[eqn]$ are the temperature changes.

The interfacial failure strain was defined as follows:

$$[eqn]$$

where $[eqn]$ and $[eqn]$ represent the thermal expansion strains of the matrix and grain phases, respectively. Among them, the grain phase refers to the energetic particles dispersed in the propellant, and the matrix phase refers to the polymer binder in the propellant.

Normal and shear stresses in the cohesive element increased linearly with interfacial separation. The slopes of these relationships defined the stiffness of the cohesive element. Upon reaching the critical normal and shear stress values, S_n_max and S_t_max, damage evolution was initiated according to

$$[eqn]$$

where $[eqn]$ and $[eqn]$ are the temperature-dependent normal and shear stresses, respectively.

Cohesive element separation was controlled using the damage variable, defined as follows:

$$[eqn]$$

where D is damage variable, $[eqn]$ and $[eqn]$ are the normal and shear displacement of the cohesive element, and $[eqn]$ and $[eqn]$ are the normal and shear temperature correction parameter of the cohesive element.

2.3. Fine-Grained Modeling

The investigated HTPB propellant is a particle-reinforced composite comprising 41 wt.% AP, 9.2 wt.% cyclotrimethylene trinitramine (RDX), 48 wt.% HTPB matrix, and additional minor additives. The typical morphology of AP and RDX crystals observed under scanning electron microscopy (SEM) is shown in Figure 2. Most AP particles exhibit a near-spherical shape with diameters ranging from 200 to 300 μm; a minority are smaller than 100 μm and display relatively smooth contours. By contrast, RDX particles are generally finer, with diameters between 40 and 60 μm. Previous studies on the microscale mechanical behavior of solid propellants have predominantly utilized two-dimensional cohesive zone models. However, to more accurately capture the effects of particle morphology and spatial distribution on mechanical responses—particularly in relation to crack initiation, propagation, and damage mechanisms—a three-dimensional (1 × 1 × 1 mm) fine-scale representative volume element was constructed. After conducting a grid convergence analysis on the numerical simulation model, the optimal grid size was determined to be 0.012 mm. At this mesh size, the computer calculation time is 18.2 h, and the peak stress of the model is basically converged to a fixed value. When the mesh size increases, the peak stress of the model is lower; when the mesh size decreases, the calculation time is longer. Therefore, the optimal mesh size for the numerical simulation model in this study is determined to be 0.012 mm. Based on the particularity of the particle shape in this three-dimensional model, we selected the 6-node three-dimensional cohesive element. Compared with the 4-node and 8-node cohesive elements, the 6-node cohesive element can significantly reduce the calculation time, enhance the bending degree of freedom, and accurately predict the interface failure. Zero-thickness cohesive elements were embedded at the interfaces and interiors of AP, RDX, and the HTPB matrix to simulate interfacial interactions and fracture evolution.

3. Model Validation

3.1. Modeling Parameters

The oxidizer crystals in solid propellants, namely ammonium perchlorate (AP) and RDX, are generally classified as quasi-brittle materials, exhibiting negligible plastic deformation upon fracture. Both AP and RDX particles were modeled as linear elastic solids. The Young’s modulus and Poisson’s ratio for AP were set at 32.4 GPa and 0.13, respectively [34] while, for RDX, the corresponding values were 8.6 GPa and 0.3 [35]. HTPB, a polymeric binder, displays significant viscoelastic behavior, including pronounced creep and stress relaxation, as demonstrated in extensive prior studies [32,36]. Based on experimental data from previous relaxation tests on HTPB [37], a master curve of the relaxation modulus was constructed using the time–temperature superposition principle, as shown in Figure 3. The master curve was fitted using an 18th-order Prony series. The fitted parameters are listed in Table 1. In the Prony model, E(t) denotes the time-dependent Young’s modulus. The corresponding bulk modulus K_i_ and shear modulus Gi were derived from the elastic modulus components using the relations $[eqn]$ and $[eqn]$ , where the Poisson’s ratio ( $[eqn]$ ) is 0.45.

3.2. Mechanical Characterization of Solid Propellants

Dynamic compressive properties were evaluated using a split Hopkinson pressure bar (SHPB), a standard apparatus for high strain rate testing of energetic materials [38,39,40]. The SHPB system included an impact bar, incident bar, transmission bar, and absorption bar, as illustrated in Figure 4. Due to the low density and acoustic velocity of the propellant, aluminum was selected for the bars to reduce impedance mismatch. The aluminum bars had a Young’s modulus of 73 GPa and a density of 2700 kg/m^3^. The lengths of the impact, incident, transmission, and absorption bars were 300 mm, 3000 mm, 3000 mm, and 1500 mm, respectively, each with a diameter of 20 mm.

Figure 5 presents the stress–strain curves of the solid propellant at various strain rates and temperatures, demonstrating typical viscoelastic behavior. Because of the high binder content, the material exhibited significant nonlinearity and toughness during deformation. At 20 °C, increasing the strain rate from 1000 s^−1^ to 4000 s^−1^ elevated the peak yield stress from 8.23 MPa to 28.64 MPa, a 3.48-fold increase. At a constant strain rate of 2000 s^−1^, raising the temperature from 20 °C to 80 °C reduced the elastic modulus from 131.82 MPa to 93.65 MPa, indicating a maximum reduction of 14.2%, consistent with temperature-dependent softening behavior.

As shown in Figure 6, SEM imaging revealed microstructural changes in the propellant from 20 °C to 80 °C. Interfacial dehumidification became increasingly prominent with temperature. Initial damage localized at the particle–binder interface. With increasing temperature, cracks initiated at these interfaces propagated into the binder matrix, forming microcracks. Thermal expansion of the binder contributed to volumetric increase and growth of internal voids.

3.3. Inverse Optimization for Cohesive Model Parameters

The bilinear cohesive zone model requires specification of three parameters for each interface: bond strength, critical displacement, and initial stiffness. Accordingly, a total of 12 parameters must be identified for the four considered interfaces. Accurate determination of these cohesive parameters necessitates an adaptive inverse algorithm characterized by rapid convergence and high precision. An artificial neural network (ANN) based on the back-propagation learning algorithm provides an efficient and accurate recursive solution. Figure 7 illustrates the training process of the ANN, which establishes a precise mapping between the stress–strain response and the cohesive parameters.

The ANN architecture comprises an input layer, two hidden layers, and an output layer. The input consists of dynamic compressive stress–strain curves extracted from the simulation database, while the output corresponds to the associated sets of cohesive parameters for each interface. Sigmoid functions serve as the activation functions for the hidden layers, and a Purelin linear function is employed in the output layer. The number of nodes in the hidden layers significantly affects network performance: insufficient nodes lead to large errors, whereas excessive nodes increase computational time. To balance accuracy and efficiency, each hidden layer was configured with 50 neurons. The network error (e) is calculated as follows:

$$[eqn]$$

where O_k_ represents the calculated value, and T_k_ is the desired network output.

The ANN exhibits strong adaptability to complex nonlinear relationships. In this study, mean square error (MSE), R^2^, and root mean square error (RMSE) are selected as the performance evaluation metrics for the model. They are defined as follows:

(1) Mean square error

$$[eqn]$$

where $[eqn]$ and $[eqn]$ represent the actual and predicted values of the ith sample, respectively, and n is the number of specimens.

(2) Coefficient of determination (R^2^)

R^2^ is the ratio of the sum of squares of the regression to the sum of squares of the total deviation in linear regression:

$$[eqn]$$

where $[eqn]$ represents the average of the actual values, and its value equals the square of the correlation coefficient.

(3) Root mean square error

$$[eqn]$$

RMSE is a common metric used to measure the difference between values.

The training process utilized MATLAB2020’s neural network toolbox. Initially, the ANN was trained using sample data from a simulation database. Once training was complete, the inversion model was applied to estimate the cohesive parameters.

The procedure is outlined below:

Step 1: Define material property parameters for the training samples.

Step 2: Perform three-dimensional fine-grained compression simulations in ABAQUS to generate training data for the ANN. Use the resulting data to construct the simulation database.

Step 3: Input experimentally measured compressive stress–strain curves into the trained ANN model (Figure 7) to predict cohesive zone parameters.

Step 4: Estimate bond strength and critical displacement values using the MRS method, (The MRS method is an intelligent inversion strategy that combines response surface models with global optimization algorithms, specifically designed for identifying material constitutive parameters, such as interface strength and fracture energy of cohesive models). Its core idea is to iteratively construct high-precision surrogate models, which reduces computational costs while ensuring the robustness and global convergence of parameter identification) [41,42]. Initial stiffness values are determined through the inverse ANN model. Assign the optimized parameters to the propellant material model and perform dynamic uniaxial compression simulations in ABAQUS. Table 2 presents the inverted cohesive zone parameters for the four interfaces across different temperatures. Validation using SHPB experimental results confirms the reliability of the finite element simulations.

3.4. Experiment-Based Validation

Figure 8 compares the computed and experimentally measured stress–strain curves for the solid propellant. The simulated results show strong qualitative agreement with experimental data, particularly in the elastic deformation region. The maximum deviation is limited to 7.3%, With reference to the error standards in Standards for Verification and Validation in Computational Solid Mechanics and Guide to the Expression of Uncertainty in Measurement, experimental results and numerical simulations with an error within 10% can be regarded as meeting the requirements of international standardization, indicating that the model accurately captures the material response and the identified parameters are valid.

Regarding the training of the artificial neural network (ANN), our ANN model was trained and optimized using a total of 1728 synthetic data samples, with 70% allocated to the training set and 30% to the validation set. Based on the results of our previous material testing experiments, the ranges for relevant cohesive parameters were determined as follows: interface initial stiffness K_n_ ∈ [15, 150] MPa/mm, maximum traction S_n_ ∈ [5, 10] MPa, and initial damage displacement δ_n_ ∈ [0.0005, 0.005] mm. In this study, the ANN model was utilized to accurately predict the stress–strain behavior of the propellant under different temperatures and strain rates. To ensure the model’s reliability, a two-step validation process was implemented. First, the synthetic data (used for training) were split into a 70% training set and a 30% validation set to guide model optimization. Second, independent experimental data were randomly divided into a 70% training subset and a 30% testing subset for further validation; the trained ANN model was then applied to predict results for the testing subset. We conducted an error analysis on the cohesion parameters inverted by the neural network, as shown in Table 3 and Table 4. These tables present the error statistics of the interfacial cohesion parameters at 40 °C with a strain rate of 2000 s^−1^, including the mean squared error (MSE), root mean squared error (RMSE), and coefficient of determination (R^2^) for both the training set (synthetic data) and the testing set (experimental data). From these statistics, it can be calculated that the overfitting coefficient of each parameter is less than 1.2, which indicates that the training result is satisfactory. In addition, we introduced 5% Gaussian white noise into the ANN model and calculated the deviation between the parameters predicted by the ANN and the real parameters under the noisy data. The results show that the relative error is approximately 6.8% under 5% noise, indicating that the model is insensitive to noise and has strong robustness. To further demonstrate the model’s accuracy, we analyzed the prediction errors of the propellant at 60 °C under four different strain rates (Figure 9). The prediction errors decreased as the loading process progressed, highlighting the high precision of the ANN model in predicting the micromechanical properties of HTPB propellants.

4. Analysis and Discussion

4.1. Thermal Damage Behavior of Propellants at Different Temperatures

A temperature–displacement coupling analysis was performed to establish a three-dimensional detailed model of the solid propellant, as depicted in Figure 10, this figure presents a comparison between the numerical simulation results of crack damage in solid propellants under different ambient temperatures and the SEM (Scanning Electron Microscope) images. The bottom surface was fully constrained, while thermal loads were applied to the remaining six surfaces. Based on experimental observations, to enable the propellant to stably heat up to different temperatures, we set a relatively slow heating rate. This ensures that the propellant can reach the various target temperatures for the study within 45 min. Meanwhile, we also set the same heating environment in the numerical simulation to guarantee the accuracy of the simulation. With increasing temperature, interfacial debonding emerged at the interface between the oxidizer particles and the binder. This process initiated the formation and gradual propagation of interfacial cracks. As illustrated in Figure 10, interfacial debonding was the initial damage mechanism. Subsequently, cracks nucleated at the tips of AP particles, evolved into matrix microcracks, and further propagated through the binder matrix.

4.2. Damage Behavior of Propellant at Different Strain Rates

The damage behavior of the solid propellant’s three-dimensional microscopic model under varying strain rates at a constant temperature was investigated. Figure 11 presents deformation cloud diagrams of the model at 20 °C for strain rates of 1000 s^−1^, 2000 s^−1^, and 3300 s^−1^. Black dashed arrows indicate the microcrack evolution trends, while white arrows highlight the formation of macroscopic main cracks caused by the coalescence of microcracks. Figure 12 displays the cross-sectional profile of Figure 10 at a position of 0.42 mm along the x-axis. This view clearly illustrates crack initiation and propagation within the oxidizer particles and binder matrix, revealing the complete evolutionary history of damage under compressive loading.

The damage evolution can be categorized into three distinct stages based on the deformation cloud and cross-sectional profile. In the first stage, elastic deformation occurs within the structure. The binder accommodates most of the compressive deformation due to its relatively low elastic modulus. The crystals exhibit minimal deformation. The second stage is marked by a nonlinear structural response. Strain concentrations develop in the binder due to lateral expansion, while the crystals remain largely undeformed. Once the surface stress of the binder exceeds its interfacial debonding strength, microcracks initiate. Strain energy localizes at the microcrack tips, leading to the development of shear slip zones. In polymeric binders, inter-chain shear slip contributes to plastic deformation and strain softening, resulting in reduced structural stiffness. As the binder’s plastic deformation capacity is exhausted, microcracks propagate, further compromising stiffness.

The third stage represents the post-peak damage evolution, during which microcracks begin to destabilize and propagate. Due to the high deformability of the binder, the cracks remain discontinuous and dispersed within the shear zone, rather than forming a dominant crack path. Significant crack tip blunting and pronounced strain localization are observed, consistent with previously reported findings [43].

Crack propagation functions as a primary mechanism for strain energy release. Accordingly, variations in crack density and propagation rates reflect changes in the energy dissipation process. To quantify damage characteristics under different strain rates in the three-dimensional microscopic model of the solid propellant, the evolution of failed units and maximum crack width was recorded, as illustrated in Figure 13.

Figure 14 illustrates a steady increase in crack volume and average crack propagation rate, indicating that the additional energy introduced by higher strain rates is predominantly dissipated through increased damage and faster crack propagation. Although both the number of failed units and propagation rates rise with strain rate, the overall increments remain moderate. As the strain rate increases, the propagation rate becomes more significant, while the number of failed units increases substantially. Specifically, under a compression strain rate of 3300 s^−1^, the crack propagation rate increases by 2.36 times, the number of failed units rises by 11.24%, and the crack volume expands by approximately 53.7%, relative to the values observed at 1000 s^−1^.

4.3. Damage Behavior of Propellant at Different Temperatures

Figure 14 and Figure 15 present deformation cloud diagrams of the solid propellant at a constant strain rate of 1000 s^−1^ under ambient temperatures of 40 °C, 60 °C, and 80 °C. The deformation sequence under different thermal conditions proceeds through consistent stages: elastic deformation, interfacial displacement between particles and matrix, dewetting at the particle–matrix interface, and irreversible microcracking in the substrate. During the third stage, longitudinal elongation of the substrate increases significantly with temperature. Dewetting becomes more pronounced at the particle–matrix interface, resulting in enlarged void regions and the redistribution of mechanical loads to the matrix. This behavior indicates that, under identical strain rates, elevated temperatures induce enhanced dewetting, more complex microcrack evolution pathways in the matrix, and reduced material toughness. This softening effect promotes increased fragmentation. The observed trend is attributed to the elevated molecular mobility in the viscoelastic matrix at higher temperatures, which leads to shorter macroscopic stress relaxation times.

Figure 16 presents the evolution of total crack length, crack volume, maximum crack width, and average crack tip expansion velocity as functions of impact time for the fine-scale solid propellant model at different temperatures. The term “crack tip” refers to cracks actively propagating at the current time step. The crack expansion rate was determined by differentiating the total crack length with respect to time, representing the cumulative propagation of all cracks in the model. The average crack expansion velocity, shown in Figure 16d, was obtained by dividing the total expansion rate by the number of active crack tips at each time step.

The data in Figure 16 indicate a gradual increase in both the number of crack tips and the average crack propagation velocity as temperature increases from 40 °C to 80 °C. This suggests that, at elevated temperatures, the strain energy introduced by dynamic loading is predominantly dissipated through increased crack tip formation and propagation rates. According to dynamic fracture theory, the Rayleigh wave velocity for the solid propellant is approximately 180 m/s. Despite increasing temperature, the peak crack propagation velocity remains below 200 m/s with limited fluctuation, indicating that propagation is constrained by the Rayleigh velocity. Therefore, the formation of additional crack tips, rather than an increase in individual crack speed, serves as the principal mechanism for strain energy dissipation at higher temperatures. Although the total crack length does not change significantly, the increase in crack tip density contributes to sustained fracture energy accumulation, thereby enhancing the system’s capacity for energy dissipation.

Interfacial debonding is widely recognized as the primary failure mode in solid propellants, with macroscopic strength strongly influenced by interfacial adhesion. Damage initiation predominantly occurs at the interface, while crack growth is largely driven by matrix deformation and fracture processes. The spatial distribution of failed units and the ratio of tensile to shear-induced damage in each component are presented in Figure 17. Across all temperature conditions, the matrix and interface exhibit a substantially higher concentration of damaged units than the particle regions. Specifically, the number of damaged units in the matrix increases with temperature, whereas the corresponding values for the interface and particles decrease. These observations indicate a temperature-driven transition in damage localization from the crystalline regions to the matrix and interface. This shift in failure mode is attributed to increased interfacial bond strength at elevated temperatures, which promotes greater matrix tearing and deformation between crystalline domains. As a result, a larger proportion of strain energy is dissipated through matrix fracture and viscoelastic deformation, thereby reducing stress concentration within the crystals and minimizing crystal rupture. However, as the matrix softens at higher temperatures, crystal regions within shear-dominant zones are subjected to enhanced compressive forces and collisions due to matrix slip. Under such conditions, compressive shear becomes the dominant damage mechanism in the crystals. These findings align with previous observations reported by Wiegand and Pinto [44]. Furthermore, the increased proportion of tensile failure at the interface suggests that tensile deformation predominates in this region. Given that mode I fracture toughness is the lowest among fracture modes, cracks are most likely to initiate at the interface. With increasing temperature, a decline in the number of fractured crystals is observed. Therefore, the post-peak damage behavior of the material exhibits more ductile characteristics, as reflected by a reduced slope in the post-peak segment of the stress–strain curve.

5. Conclusions

This study presents a novel methodology that integrates artificial intelligence-based inverse identification algorithms with three-dimensional cohesive zone modeling to investigate the damage mechanisms in solid propellants under coupled thermal and mechanical loading conditions. The analysis reveals that energetic solid propellants exhibit pronounced sensitivity to both strain rate and temperature. The primary conclusions are as follows:

(a) Increasing temperature significantly increases the number of crack tips, and the fracture energy growth rate exceeds the crack propagation rate limit. This enables efficient dissipation of strain energy at elevated temperatures.

(b) At higher temperatures, dewetting and fragmentation become more severe under equivalent strain levels. These effects are accompanied by more complex microcrack networks in the matrix, indicating a marked reduction in material toughness (softening behavior).

(c) The mode of crack propagation is strongly temperature-dependent. Lower temperatures are associated with higher crack velocities and longer total crack lengths, while higher temperatures promote an increase in crack tip density. This shift enhances peak fracture energy and its growth rate, facilitating effective energy dissipation through this mechanism.