Viscoelastic damage constitutive modeling of the composite propellants considering coupled interfacial debonding and porosity evolution

Abstract

The composite propellant is a type of particle-included composite material and serves as the power source of solid rocket motors. During service, the composite solid propellant sustains complex mechanical loads, which would cause internal damage and ultimately lead to material failure. To investigate the rate-dependent nonlinear mechanical behaviors of nitrate ester plasticized polyether (NEPE) composite propellants under tensile loading, this study focused on the mesoscopic damage mechanisms via observation of mesostructures. The key damage modes were particle-matrix interfacial debonding and expansion of dewetting pores, with porosity evolution showing significant rate dependence. Based on the mesoscopic damage mechanisms and the rate-dependent porosity prediction model, a viscoelastic damage constitutive model was established by introducing an internal state variable (characterizing both damage modes) and a softening function (reflecting modulus degradation). The model effectively predicted the propellant’s nonlinear responses under different strain rates, including the typical three-stage stress trend at high strain rates. Numerical analysis confirmed rate-dependent damage accumulation: higher strain rates delay damage initiation and fracture, with consistent failure modes across strain rates. This work links mesoscopic damage to macroscopic properties, providing a reliable model for propellant design.

Graphical Abstract

Keywords

composite propellant mesoscopic damage mechanism interfacial debonding dewetting pores constitutive model

Introduction

The composite propellant is the solid fuel of the solid rocket motors, and is a kind of polymer -based particles-reinforced composites, which is extensively utilized in high-efficiency propulsion systems.^1,2,3 Due to excellent ductility, stability and aging resistance, the nitrate ester plasticized polyether (NEPE) has been wildly used in aerospace industry.^4,5,6 The propellant grains will sustain mechanical loads for a long time in service. During the cooling stage after fabrication, uniaxial tensile stress is generated by the grain’s shrinkage, while in the ignition stage, the propellant interior is in a biaxial stress state under the action of high-temperature internal combustion gas flow.⁷ As the propellant bears complex mechanical loads, this process is often accompanied by the occurrence of internal mesoscopic damage, thus leading to the degradation of its mechanical properties and the reduction of load-bearing capacity.^8,9 To improve the mechanical properties of propellants, investigating their internal mesoscopic failure mechanisms and establishing the relationship between mesoscopic failure and macroscopic properties has become an effective method.^10,11

Experimental observation is a direct method to obtain the mesoscopic failure process of propellants. By combining optical microscope (OM), scanning electron microscope (SEM), micro computed tomography (micro CT) and synchrotron radiation technology with different loading devices, the mesoscopic damage evolution process inside the propellants can be acquired. Wubuliaisan et al.¹² recorded the evolution of damage morphologies of the NEPE propellant with SEM and OM. Experimental results indicate that the damage distribution inside the propellant exhibits significant rate dependence, manifested as uniform distribution of interfacial damage at low strain rates while damage localization at high strain rates. Meanwhile, particle characteristics also affect the damage evolution. Van Ramshorst et al.¹³ found that the interfacial debonding process tends to occur at the interfaces of large particles, which is also the main failure mode influencing the damage evolution of propellants. Except for OM and SEM which could only acquire the surface morphologies of materials, micro CT and synchrotron radiation technology have also been developed and applied to the observation of mesoscopic damage.^14,15,16,17 Liu et al.¹⁸ and Dong et al.⁷ tracked the morphologies evolution of propellants under tension, and results shown that dewetting pores mainly exhibited hat-like shape and covering particles like ellipsoids. Lai et al.¹⁹ used synchrotron radiation to observe the interfacial debonding process of internal particle-matrix interfaces in a composite propellant under uniaxial tension, and the results indicated that interfacial cracks tend to nucleate at points with maximum particle surface curvature and propagate along the curvature gradient of the particle surface then propagated along the curvature gradient. In addition, Geng et al.²⁰ concluded that the mesoscopic damage modes of propellants under biaxial compression loading are sensitive to temperature, which was observed by micro CT, as transgranular fractures tended to occur at lower temperatures. Owing to the capability of computed tomography to directly capture the internal mesoscopic structure of materials,^21,22,23 it has been widely applied to acquiring the mesoscopic failure mechanisms of propellants under service conditions.^24,25,26,27

The NEPE propellants possess viscoelastic properties,²⁸ which result in distinct strain rate-dependent behaviors when subjected to mechanical loading, in addition, the evolution of internal damage could lead to the nonlinearity of their macroscopic mechanical responses. Therefore, it is a key step for propellant formulation optimization and grain design to establish a constitutive model capable of characterizing the nonlinear mechanical responses of propellants. For the characterization of propellant damage behavior, researchers commonly employ phenomenological property degradation formulas to depict the property degradation of propellants under various loading conditions, primarily incorporating the influences of strain rate,^29,30 aging³¹ and cyclic loading.^32,33 Schapery incorporated the pseudo strain theory based on irreversible thermodynamics and used internal state variables to characterize the degree of internal damage in propellants, thus establishing a damage constitutive model for propellants.^34,35 Based on this work, researchers have subsequently developed a variety of theories for predicting the mechanical behaviors of propellants under the action of complex loads,³⁶ including temperature loading³⁷ and cyclic loading.³⁸ In order to establish a connection between mesoscopic damage evolution and macroscopic mechanical property degradation, numerous researchers have adopted porosity as the damage reference parameter to develop prediction approaches for the performance of materials.^39,40 Xing et al.⁴¹ utilized porosity as the internal damage parameter and formulated a constitutive model for the mechanical behavior of NEPE propellants by random debonding theory. Lei et al.⁸ constructed a constitutive model for propellants based on strength of particle-matrix interfaces, and predicted the macroscopic mechanical responses of a propellant. In this way, porosity has become a critical parameter and connector between mesoscopic damages and macroscopic properties. As propellants evolve towards multi-particle compositions and high particle volume fractions, existing constitutive models are unable to capture the effect of specific internal particle information on the propellant failure process. Therefore, developing a constitutive model that accounts for the propellants’ internal particle characteristics and mesoscopic damage mechanisms has emerged as a reliable approach to achieve accurate prediction of their mechanical properties.

However, existing constitutive models for composite propellants still suffer from several critical limitations. Most models are phenomenological in nature and rely on macroscopic property degradation functions without explicitly accounting for real mesoscopic damage mechanisms. Although some models have introduced porosity as a damage variable, they fail to separately consider the coupled effects of interfacial debonding initiation and subsequent dewetting pore evolution, which are the two dominant mesoscopic failure modes of NEPE propellants. In addition, most constitutive models ignore the real particle size effect and rate dependence of damage nucleation and propagation, leading to insufficient accuracy in predicting the nonlinear mechanical responses. Furthermore, few models can quantitatively link mesoscopic damage evolution (interfacial debonding and pore growth) to macroscopic modulus degradation and rate-dependent fracture behavior.

Since the existing constitutive models are unable to capture the effect of specific internal particle information of propellants during damage analysis, this work developed a particle size-dependent damage prediction equation based on the processes of particle-matrix interfacial debonding and expansion of dewetting pores within propellants. By accounting for the actual particle characteristics and dominant failure mechanisms within composite propellants, the present model demonstrates enhanced predictive accuracy and broader applicability. Subsequently, through the incorporation of a softening function to quantify the impact of internal damage accumulation on macroscopic mechanical properties, a damage constitutive model considering the mesoscopic structural features of propellants was established, which enabled the prediction of nonlinear mechanical responses and numerical analysis of damage behaviors of propellants under various strain rates.

Experiments section

Material and test plan

NEPE composite propellant was used in this work, which consists of polymer matrix and inclusion particles. The primary formulation of the propellant matrix consists of polyethylene glycol (PEG) and nitroglycerin (NG), together with plasticizers, anti-aging agents, and other optimizing components. And inclusions are aluminium (Al) particles, ammonium perchlorate (AP) particles, and high-energy Octogen (HMX) particles.

The propellant is a kind of particulate inclusion polymer matrix composite, so mesoscopic damage evolution within the propellant under loading could result in a nonlinear viscoelastic response. This study conducted macroscopic mechanical property tests and in-situ mesoscopic damage observation of the propellant under uniaxial tension, to characterize its nonlinear mechanical response and the effect of mesoscopic damage on macroscopic properties. Macroscopic tests cover uniaxial tension and stress relaxation, providing stress-strain curves at various strain rates and relaxation modulus. In-situ observation is implemented via sequential CT scanning at the identical propellant location under different strain conditions, to capture the evolution of internal mesostructural damage.

Mechanical property tests

The Zwick/Roell Z005 testing machine (as shown in Figure 1(a)) was utilized to apply tensile loading, and the NEPE propellant was processed into dumbbell specimens (details of the size are presented in Figure 1(b)), in this way, the uniaxial displacement loading can be applied to the propellant specimens (as shown in Figure 1(c)). Both of the uniaxial tension tests and stress relaxation test were conducted at the room temperature (about 25°C), and the nominal stress and strain were adopted to characterize the stress and deformation states of the specimens. In compliance with prevailing industrial standards, uniaxial tensile experiments were conducted on macroscopic dumbbell specimens with a 70-mm effective stretching zone at loading rates of 1, 2, 10, 20, 50 and 100 mm/min. By converting the displacement rates to strain rates, six strain rates incrementally varied from 2.38 × 10⁻⁴ s⁻¹ to 2.38 × 10⁻² s⁻¹ were set. And the propellant specimens were loaded at a strain rate of 2.38 × 10⁻² s⁻¹ until the strain attained 5%, followed by a displacement-holding period of 2000 s, to obtain the stress relaxation curves of the propellant specimens for the stress relaxation tests.

Figure 1.

Mechanical testing machine and specimen. (a) Zwick/Roell Z005 testing machine. (b) A diagram of dumbbell-shaped specimen. (c) The comparison of propellant specimens before and after tension.

Nominal stress-strain curves at different strain rates were plotted in Figure 2. The results indicated that the propellant material exhibited obvious nonlinear viscoelastic characteristics under uniaxial tensile loading. With the increase in strain, the stress first increased linearly, followed by a deviation of the stress–strain curve, which demonstrated a softening behavior. Subsequently, as the tensile load continued to increase, the stress began to rise slowly, resulting in a distinct plateau stage in the latter segment of the stress-strain curve. In addition, the critical strains corresponding to the initial deviation of the stress-strain curves differed with varying strain rates, suggesting that strain rate affects the internal damage evolution of the propellant and results in strain rate-dependent nonlinearity in its macroscopic mechanical response. Notably, the stress-strain curves of the propellant at high strain rates exhibited a characteristic of first increasing, then decreasing, and finally rising slowly, whereas no obvious stress decline stage existed under low strain rates.

Figure 2.

Stress–strain curves at various strain rates of the NEPE propellant.

The polymer matrix of the propellant is a typical viscoelastic material, which renders the overall propellant exhibit viscoelastic characteristics. The relaxation modulus $E (t)$ can be expressed by the Prony series as follows:

E (t) = E_{\infty} + \sum_{j = 1}^{k} [E_{j} \cdot \exp (- \frac{t}{λ_{j}})]

(1)

where

t

means time,

E_{\infty}

represents the equilibrium modulus,

E_{j}

and

λ_{j}

are jth Young’s modulus and relaxation time respectively. Then curve fitting of the relaxation data was conducted to extract parameters of the six-term Prony series representation, with fitting results presented in Figure 3 and parameters numerically summarized in Table 1.

Figure 3.

Prony series fitting of stress relaxation curve.

Table 1.

Fitted parameters of six-term Prony series.

$E_{\infty}$ (MPa)	$E_{1}$ (MPa)	$E_{2}$ (MPa)	$E_{3}$ (MPa)	$E_{4}$ (MPa)	$E_{5}$ (MPa)	$E_{6}$ (MPa)
1.066	0.157	0.190	0.140	0.107	0.130	0.256
--	$τ_{1}$ (s)	$τ_{2}$ (s)	$τ_{3}$ (s)	$τ_{4}$ (s)	$τ_{5}$ (s)	$τ_{6}$ (s)
--	0.332	2.237	12.056	48.008	185.22	1409.0

Mesoscopic damage evolution observation

Variations in the internal mesostructure of the propellant are the intrinsic factor inducing variations in its macroscopic mechanical properties, thus, acquiring the mesoscopic damage evolution process of the propellant constitutes an effective approach for analyzing its nonlinear mechanical behavior. As reported in our previous work,⁷ a high-resolution observation method for propellant mesoscopic damage based on micro-CT technology has been developed, which enables the acquisition of high-resolution grayscale image sets of the propellant’s internal mesostructure under different strain rates and strain states, as shown in Figure 4(a). Subsequently, a mesostructural model of the propellant can be obtained through component segmentation, noise removal and three-dimensional reconstruction, as shown in Figure 4(b). By stretching the specimens from the initial state to a 160% strain state while performing continuous scanning at the same location, the evolution process of mesoscopic damage of the propellant under uniaxial tensile loading can be acquired, as illustrated in Figure 4(c). In previous work, experimental observations of uniaxial tensile tests on the propellant were conducted under four strain rates as 4.76 × 10⁻⁴/s, 1.19 × 10⁻³/s, 2.38 × 10⁻³/s, and 3.57 × 10⁻³/s respectively.⁷ The results indicated that the mesoscopic damage evolution modes of this propellant are generally consistent under different stretching rates. With the increase in strain loading, interfacial debonding first occurred on the surfaces of large particles along the loading direction, generating interfacial cracks. Subsequently, the cracks further propagated and deformed into dewetting pores. When the strain increased to a critical level, most of the large particles were surrounded by dewetting pores, which began to interconnect and coalesce, eventually leading to the occurrence of fracture.

Figure 4.

In-situ mesoscopic damage observation of the NEPE propellant. (a) In-situ tensile tests instruments and procedure. (b) Reconstructions of the propellant components. (c) Evolution process of the mesoscopic damage.

Internal damage analysis

Interface debonding and porosity evolution

According to the in-situ tensile test, the mesoscopic damage evolution of the propellant is dominated by the interface debonding process between particles and the matrix. Subsequently, the interfacial debonding cracks propagate along the interfacial direction under loading, forming damage pores of increasing size. The formation of interfacial debonding cracks causes the particle regions with interfacial debonding to lose their load-bearing capacity, resulting in the degradation of the overall mechanical properties of the propellant. Likewise, the growth of the debonding pores could reduce the solid region inside the propellant, which further lowers the load-bearing capacity of the material. Therefore, this work analyzed the debonding process at the particle-matrix interface and the evolution process of internal damage pores within the propellant, so as to reveal the influence of mesoscopic damage evolution on the macroscopic mechanical properties of the propellant.

As previously reported,^19,42,43 the particle-matrix interface debonding process of propellant is governed by the critical stress threshold, by regarding the particles as spheres with radii r, the critical debonding stress $σ_{d}$ can be expressed as

σ_{d} = {\dot{ε}}^{θ} \sqrt{\frac{Υ E_{e} f_{t}}{r}} = {\dot{ε}}^{θ} \sqrt{\frac{α}{r}}

(2)

where

γ

represents the interfacial fracture energy density,

E_{e}

denotes the equivalent modulus of the propellant,

f_{t}

is the total particle volume fraction. By integrating the parameters determined by the material properties, the critical debonding stress

σ_{d}

is determined by particle radius r, strain rate

\dot{ε}

, and material-dependent constants

α, θ

. It is worth noting that the critical debonding stress denotes that interfacial debonding will occur in all particles with a radius greater than or equal to r when the nominal stress of the bulk propellant material reaches the corresponding level, rather than that fracture takes place once the stress at the particle-matrix interface attains this value. Equation (2) indicates that the stress required to induce interfacial debonding decreases as the particle radius increases. Accordingly, debonding cracks and pores tend to initiate preferentially around large particles during the tensile process.

Porosity is defined as the percentage of the pore volume inside a material relative to its total volume, which could reflect the cumulative degree of mesoscopic damage within the propellant to a certain extent, and thus is commonly utilized for the quantitative analysis of damage. In previous studies, the evolution processes of internal damage in the propellant under four strain rates, obtained from in-situ observation tests, were plotted as shown in Figure 5.⁷ The results indicate that the initial porosity of propellant is about 1%, at the strain stage below 40%, the porosity increased slowly, and the dominant damage mode inside the propellant was particle-matrix interfacial debonding at this stage. Although the formation of debonding cracks had a slight effect on the increase in porosity, it led to a significant degradation of the macroscopic mechanical properties of the propellant, which is manifested by the deviation and softening of the stress-strain curve in Figure 2 when the strain is less than 40%. When the strain exceeded 40%, the porosity exhibited a distinct increasing trend, and the primary damage mode inside the propellant at this stage was the transformation of debonding cracks into dewetting pores, as well as the propagation and coalescence of cracks, resulting in the continuous accumulation of internal damage in the propellant during this stage. Furthermore, the porosity of the propellant at the same strain decreased with an increase in strain rates. This is because the strength of the particle-matrix interface rises with increasing strain rate, and different strain rates also lead to differences in the stress of the propellant at the same strain state, thus resulting in an obvious rate-dependent effect in the growth of porosity.

Figure 5.

Porosity evolution curves at different strain rates of NEPE propellant.

Porosity prediction model

To predict the mesoscopic damage accumulation of the propellant at various strain rates, a porosity prediction model based on particle-matrix interfacial debonding and dewetting pore deformation was developed based on the modified Farris model in prior work.⁷ The model postulates that particles are spheres without deformation under loading (due to the significant mismatch in Young’s modulus between the particles and the matrix), and interfacial debonding initiates once the stress attains the critical value. Under tensile loading, interfacial cracks then develop into ellipsoidal pores that encapsulate the particles (as shown in Figure 6), where the minor axis of the ellipsoid matches the particle radius, and the major axis growth rate with strain has a linear correlation with the particle radius:

\frac{d b}{d ε} = β a

(3)

where

β

is a constant determined by the material properties of the propellant.

Figure 6.

The sketch of dewetting pores.

The total volume of the pores $V_{p o r e}$ can be calculated as:

V_{p o r e} = \sum_{i = 1}^{N} V_{i} = \sum_{i}^{N} [\frac{4}{3} π (a_{i}^{2} b_{i} - a_{i}^{3})] = \sum_{i = 1}^{N} (\frac{4}{3} π a_{i}^{3} \frac{b_{i} - a_{i}}{a_{i}}) = \sum_{i = 1}^{N} (V_{i}^{p} \frac{b_{i} - a_{i}}{a_{i}})

(4)

where

V_{i}

and

V_{i}^{p}

represent the volume of ith pore and particle respectively, N is the total number of the particles, and

a_{i}, b_{i}

are the axis length of ith pore. By classifying the particles into n particles with dewetting pores and N-n particles without pores, growth of the pores’ volume

V_{p o r e}

as the strain increases can be expressed as:

\frac{d V_{p o r e}}{d ε} = \sum_{i = 1}^{n} \frac{d V_{i}}{d ε} + \sum_{i = n + 1}^{N} 0 = β \sum_{i = 1}^{n} V_{i}^{p}

(5)

The porosity $ω$ is calculated by:

ω = \frac{V_{p o r e}}{V_{0} + V_{p o r e}}

(6)

where

V_{0}

is a constant, representing the volume of bulk materials including particles and matrix. Combined with equation (5) the growth rate of porosity with the increase of strain is expressed as

\frac{d ω}{d ε} = {(\frac{V_{0}}{V_{0} + V_{p o r e}})}^{2} (\frac{1}{V_{0}} \frac{d V_{p o r e}}{d ε}) = β {(1 - ω)}^{2} (\frac{1}{V_{0}} \sum_{i = 1}^{n} V_{i}^{p}) = β {(1 - ω)}^{2} f_{d}

(7)

where

f_{d}

is the total volume fraction of all debonded particles. The equation indicates that the porosity growth rate is determined by the current porosity

ω

, the total volume fraction of debonded particles

f_{d}

, and the material constant

β

. For the viscoelastic materials, the stress (without damage)

σ^{R}

can be calculated by the relaxation modulus:

σ^{R} (t) = \int_{0}^{t} E (t - τ) \frac{\partial ε}{\partial τ} d τ

(8)

where

τ

is introduced for integrated. Then the total volume fraction of debonded particles

f_{d}

can be determined by critical debonding stress (equation (2)) and radii distribution of all internal particles (acquired by mesoscopic observation). equations (2), (7), and (8) form the porosity prediction model. The prediction procedure for porosity is as follows:

(1) The evolution of the internal stress $σ^{R}$ in the propellant is calculated using equation (8) based on the prescribed strain rate $\dot{ε}$ .

(2) The particle size range subject to interfacial debonding under different strain states is determined from the internal stress evolution $σ^{R}$ and equation (2).

(3) The total volume of debonded particles is computed using the complete particle information obtained from mesoscopic observations, and the debonded particle volume fraction $f_{d}$ is obtained by dividing the total volume of debonded particles by the initial undeformed volume of the propellant.

(4) Finally, the porosity evolution is determined using equation (7).

Subsequently, the values of the undetermined parameters were obtained by performing an inversion on the experimental results derived from mesoscopic damage evolution observations ( $α = 190$ , $β = 1.2058$ , $θ = 0.28$ ). The comparison between the predicted results and the experimental results is shown in Figure 7, which indicates that the porosity prediction model could accurately predict the mesoscopic damage evolution process inside the propellant under different strain rates, which was illustrated in previous study.⁷

Figure 7.

Comparison porosity evolution cures between experiments and prediction.

Constitutive model

Damage constitutive model based on mesostructures

For viscoelastic materials, Schapery proposed a viscoelastic constitutive model by introducing the pseudo-strain theory,^34,35 which replaces the strain in the constitutive relation of elastic materials with the pseudo-strain, thus correlating the viscoelastic stress response with the pseudo-strain. In this work, the calculation formula for the pseudo-strain $ε^{R}$ is given as follows:

ε^{R} = \frac{1}{E_{i nitial}} \int_{0}^{t} E (t - τ) \frac{d ε}{d τ} d τ

(9)

where

E_{initial} = E (0) = 2.045 MPa

means the initial relaxation modulus of the propellant. Based on equation (8),

ε^{R} = σ^{R} / E_{initial}

, which means that pseudo-strain is a dimensionless quantity, defined as the ratio of the stress without damage to the initial relaxation modulus. Therefore, the correlation between stress and pseudo-strain is characterized by the damage degree of the material at the same strain state.

Section Interface debonding and porosity evolution describes that the evolution of internal damage and the variation of mesostructures in the propellant during loading are the key factors inducing the nonlinear mechanical response of the material. Park and Schapery⁴⁴ derived the irreversible process of damage evolution using thermodynamics and characterized the degree of mesostructural damage of the material with internal state variables (ISVs) S (S₁, S₂, S₃,…,S_m), in which case the stress can be expressed as:

σ = \frac{\partial W^{R}}{\partial ε^{R}} = C (S) E_{0} ε^{R}

(10)

where

W^{R}

means the pseudo-strain energy density,

C (S)

is the softening function of internal state variables S,

E_{0}

denotes the modulus of the materials without any damage. By introducing the effective modulus

E^{e}

at the damage state S, which makes

σ = E^{e} ε^{R}

, then the softening function means the ratio of the effective modulus

E^{e}

to modulus of the materials without any damage

E_{0}

C = E^{e} / E_{0}

. In this case, the softening function

C (S)

represents the remaining fraction of the material modulus in the current damage state relative to that in the undamaged state.

Based on the in-situ observations of mesoscopic damage evolution of the propellant under tensile loading in Section Internal damage analysis, the main damage modes inside the propellant are the interfacial debonding at the particle-matrix interface and the propagation process of dewetting pores. Therefore, the internal damage state variables of the propellant should be determined by these two damage modes. Porosity evolution can characterize the volume growth of dewetting pores, while the contribution of the initial stage of interfacial crack nucleation to the internal damage accumulation of the material cannot be reflected by it (as shown in Figure 8). Therefore, porosity $ω$ cannot be directly used to replace the internal state variables S, instead, the growth trend of S should have the same form as that of porosity $ω$ , yet it should start to accumulate at the onset of interfacial debonding.

Figure 8.

Schematic diagram for evolutions of damage degree and porosity during interfacial debonding and expansion of dewetting pores.

Therefore, the internal state variable S in this work has a single component and adopts the same governing equation with the porosity growth process, with only the critical stress value for the initiation of its growth requiring modification. Combining with porosity prediction model as illustrated in Section Porosity prediction model, the calculation formula for the internal state variable is given as follows:

{\begin{cases} \frac{d S}{d ε} = 1.2058 {(1 - S)}^{2} f_{d} \\ σ_{d} = a_{1} {\dot{ε}}^{0.28} \sqrt{\frac{190}{r}} \\ σ^{R} (t) = \int_{0}^{t} E (t - τ) \frac{\partial ε}{\partial τ} d τ \end{cases}

(11)

where

a_{1}

is the undetermined parameter that modifies the initiation position of damage evolution, and initial value of S is given as 0.01 which equals to the initial porosity.

For the softening function $C (S)$ , (a) value of 0 for the internal state variable indicates the absence of internal damage in the material. At this point, the material modulus does not degrade, and thus $C (0) = 1$ . Then C decreases gradually with the increase in S, reflecting the gradual degradation of the propellant’s overall modulus with damage accumulation. Based on the rate-dependent porosity evolution law identified by in-situ micro-CT tests and the research on the softening behavior of the polymer-based composites,^45,46,47 the form of softening function should be prescribed to capture the slow initial damage accumulation stage and the rapid damage growth stage, which exponential and polynomial functions fail to describe precisely. The functional relationship between C and S in this work is expressed as follows

C (S) = {(1 - a_{2} S^{a_{3}})}^{a_{4}}

(12)

where

a_{2}

a_{3}

a_{4}

are undetermined parameters. By combing equations (9)–(12), the viscoelastic damage constitutive model considering coupled interface debonding and porosity evolution is summarized as

{\begin{cases} σ = C (S) E_{0} ε^{R} \\ ε^{R} = \frac{1}{E_{initial}} \int_{0}^{t} E (t - τ) \frac{d ε}{d τ} d τ \\ C (S) = {(1 - a_{2} S^{a_{3}})}^{a_{4}} \\ \frac{d S}{d ε} = 1.2058 {(1 - S)}^{2} f_{d} \end{cases}

(13)

This formula contains five undetermined parameters, which are $a_{1}$ , $a_{2}$ , $a_{3}$ , $a_{4}$ , $E_{0}$ respectively. Since the manufacturing process cannot produce materials without any initial damage, $E_{0}$ cannot be experimentally measured and is obtained through parameter fitting.

Parameters inversion and model validation

The results of uniaxial tension tests, stress relaxation tests and in-situ mesoscopic damage evolution observation which were summarized in sections Experiments Section and Internal damage analysis were utilized to inverse fitting for undetermined parameters. Firstly, the experiments data of in situ X-ray tomography of porosity under tensile loading and stress relaxation tests were used to acquire porosity prediction model and calculate pseudo-strain $ε^{R}$ , then the evolution of internal state variable S and the softening function $C (S)$ could be calculated by giving the initial values to $a_{1}$ , $a_{2}$ , $a_{3}$ , $a_{4}$ , subsequently, the stress prediction could be conducted by $σ = C (S) E_{0} ε^{R}$ with initial value of $E_{0}$ , finally, the correlation analysis was made between prediction stress and uniaxial tension tests. As a widely accepted statistical indicator in constitutive modeling studies, $R^{2}$ can comprehensively characterize the overall fitting quality and consistency across the entire strain range, thus $R^{2}$ was employed to evaluate the fitting and predicting performance of the constitutive model, with formula as

R^{2} = 1 - \frac{\sum_{i = 1}^{n} {(σ_{i}^{t e s t} - σ_{i}^{p r e})}^{2}}{\sum_{i = 1}^{n} {(σ_{i}^{t e s t} - {\bar{σ}}^{t e s t})}^{2}}

(14)

where n denotes the number of data points on the stress–strain curve,

σ_{i}^{t e s t}

represents the stress value of ith point acquired by experiments,

σ_{i}^{p r e}

represents the stress value of ith point acquired by constitutive model, and

{\bar{σ}}^{t e s t}

represents the mean value of the experimental (true) stress. If coefficient of determination

R^{2}

is higher than 0.9, which means the prediction results matches well with the experiments, the inverse fitting could be completed, otherwise, the values of undetermined parameters will be adjusted and repeated the fitting procedure until meeting the requirement. The inverse fitting procedure is depicted in Figure 9.

Figure 9.

The flow chart of inverse fitting for the parameters determination of constitutive model.

The results of uniaxial tension tests at strain rates 4.76 × 10⁻⁴ s⁻¹, 2.38 × 10⁻³ s⁻¹, 4.76 × 10⁻³ s⁻¹ were employed for parameters determination, and the comparison between inverse fitting curves and experiment results is plotted in Figure 10, with the parameters listed in Table 2. The fitting curves are basically consistent with the experimental curves, with all R-squared values exceeding 0.97, which verifies the accuracy of the obtained parameters. It is worth noting that $E_{initial} < E_{0}$ , indicating that the mechanical properties of the propellant have already degraded after fabrication due to the existence of initial defects in the material.

Figure 10.

Comparison between inverse fitting curves and experiment results.

Table 2.

Fitted parameters of constitutive model.

$E_{0}$ (MPa)	$a_{1}$	$a_{2}$	$a_{3}$	$a_{4}$
3.1786	0.2950	0.0301	0.4417	127.46

The established constitutive model was employed to predict the stress response under different strain rates 2.38 × 10⁻⁴ s⁻¹, 1.19 × 10⁻² s⁻¹, 2.38 × 10⁻² s⁻¹, as shown in Figure 11. The results indicate that the predicted results are in good agreement with the experimental results. Furthermore, the parameters in the model were obtained from the stress-strain curves corresponding to three intermediate strain rates, and the accurate prediction of the stress-strain curves under both higher and lower strain rate could validate the accuracy and generality of the constitutive model. Since the growth form of the internal state variable in the present model is governed by the mesoscopic damage evolution modes and particle information inside the propellant, it can accurately characterize the influence of the propellant’s internal mesoscopic damage state on its macroscopic mechanical properties under different loading conditions, thereby ensuring the model’s accuracy. This also enables accurate prediction of the characteristic that the propellant’s stress-strain curve under high strain rates exhibits an initial rise, followed by a decline and then a gradual increase after the onset of softening deviation.

Figure 11.

Strain-stress curves of tests and predicted results.

Numerical analysis of damage behavior

To analyze the damage evolution process of the NEPE propellant during loading, the variation trend of the softening function C with strain under six tensile rates was calculated according to equation (12), as shown in Figure 12. The results indicate that the modulus attenuation degrees of the propellant at different tensile rates are identical at the initial tensile stage due to the same initial damage represented by porosity. With increasing strain, the damage accumulation gradually rises, leading to a continuous decrease in the softening function, which in turn causes the macroscopic modulus of the propellant to exhibit softening as strain increases. Figure 12 shows that there is a gentle change stage in the softening function curves before the onset of softening evolution, implying a slow increment in the damage degree. Combined with the porosity damage evolution process obtained from in-situ observation tests (as shown in Figure 5), it reveals that the particle–matrix interfaces remain nearly intact and mesoscopic damage accumulation is negligible at the initial loading stage, resulting in an obvious plateau in the softening curve of C. Subsequently, with the occurrence of interfacial debonding and the expansion of dewetting pores, mesoscopic damage accumulates rapidly, leading to a pronounced decrease in the softening function C. Moreover, the higher the tensile rate, the longer the gentle change stage of the softening function, indicating that the critical strain for the initiation of damage evolution is higher and the internal damage accumulation of the propellant occurs later.

Figure 12.

C-ε curves under different strain rates.

In addition, the evolution processes of the softening function under different tensile rates are basically consistent. The only difference is that the greater the tensile rate, the larger the softening function value at the same strain, which demonstrates that the failure mode of the propellant from damage accumulation to fracture is basically consistent under different tensile rates. However, the lower the tensile rate, the faster the internal damage accumulation rate of the propellant with increasing strain and the greater the degradation degree of the propellant’s mechanical properties. This conclusion is consistent with the observation results of the propellant’s mesoscopic damage evolution process in section Internal damage analysis.

For the evolution process of the softening function at each tensile rate, this study adoped the fracture strain obtained from macroscopic tensile tests of the propellant to calculate the softening function value $C_{m}$ at fracture under each tensile rate. The damage accumulation degree parameter D was then defined as $D = C / C_{m}$ , and the propellant fractures when $D = 1$ , indicating that the internal damage accumulation of the propellant reaches the critical threshold for fracture. To analyze the damage accumulation evolution of the propellant at each tensile rate, the evolution of damage degree parameter D from its initial value to the fracture state was calculated. The variation of D is plotted in conjunction with the stress evolution trend, as shown in Figure 13.

Figure 13.

Evolutions of damage degree D and stress under tensile loading at different strain rates. (a) $\dot{ε}$ = 2.38 $\times$ 10^-4s^-1, (b) $\dot{ε}$ = 4.76 $\times$ 10^-4s^-1, (c) $\dot{ε}$ = 2.38 $\times$ 10^-3s^-1, (d) $\dot{ε}$ = 4.76 $\times$ 10^-3s^-1, (e) $\dot{ε}$ = 1.19 $\times$ 10^-2s^-1, (f) $\dot{ε} = 2.38 \times 10^{- 2} s^{- 1}$ .

The results demonstrate that in the initial tensile stage, the damage accumulation parameter D increases extremely slowly, and the propellant exhibits an overall linear viscoelastic behavior. When the strain reaches the dewetting point, damage accumulates sharply, causing continuous softening of the propellant and constant deviation of the stress-strain curve, until the damage accumulation degree parameter reaches 1 and the propellant fractures.

Analysis of the variation trend of $C_{m}$ at six strain rates reveals that softening function value at fracture $C_{m}$ increases with the rise in strain rate and exhibits an exponential relationship with strain rate (as shown in Figure 14), which can be expressed as

C_{m = -} 0.0067 e x p (- 398.0 \dot{ε}) + 0.1061 \exp (0.3111 \dot{ε})

(15)

Figure 14.

Variation trend of softening function value at fracture $C_{m}$ of propellant with strain rate.

Based on this functional relationship, the softening function value of the propellant at fracture under different strain rates can be predicted. Furthermore, the maximum elongation and tensile strength of the propellant can be predicted by equations (11), (12), and (15).

Conclusions

In this work, the mesoscopic damage modes and nonlinear viscoelastic behaviors of NEPE propellant under tensile loading were characterized. Based on the mechanism of interface debonding and expansion of dewetting pores, the governing equations of internal state variable was constructed to represent the damage state in the propellant. Then the softening function was introduced to characterize the influence of mesoscopic damage evolution on degradation of macroscopic mechanical properties. And a viscoelastic damage constitutive model considering coupled interface debonding and porosity evolution was subsequently established. Subsequently, the numerical analysis of damage behavior of propellant was conducted. The following conclusions can be drawn:

The internal damage evolution of the propellant under tensile loading is mainly dominated by particle-matrix interfacial debonding and the expansion of dewetting pores. Interfacial debonding causes the covered particle regions to lose their load-bearing capacity, and the expansion of dewetting pores reduces the solid region inside the propellant. These two factors are the primary causes for the nonlinear characteristics of the propellant’s macroscopic mechanical response.

The viscoelastic damage constitutive model established based on the mesoscopic damage evolution mode of the propellant can accurately predict the nonlinear mechanical behaviors of the propellant under different strain rates, and precisely characterize the feature that the propellant’s stress first increases, then decreases, and finally rises gently with increasing strain at high strain rates. Under tensile loading, the propellant’s internal damage state initially remains unchanged, then damage starts to accumulate once the strain exceeds a certain critical value, and this critical strain increases with increasing strain rate. This indicates that the damage accumulation process inside the propellant is rate-dependent. At the same strain state, a higher strain rate corresponds to a smaller damage accumulation amount and a later occurrence of fracture.

Footnotes

Acknowledgments

The authors sincerely thank the technical support by the Core Facilities and Experiment Center of Xi’an Jiaotong University. And authors specially thank the advices for color blending of graphs by Zishi Wang.

ORCID iD

Zhelin Dong

Author contributions

Zhelin Dong: Methodology, Writing-original draft, Conceptualization, Validation, Investigation, Data curation. Kaining Zhang: Formal analysis, Methodology, Resources. Chunguang Wang: Writing-review and editing, Resources, Supervision, Conceptualization. Hanlin Wang: Investigation, Formal analysis. Feifei Zhu: Investigation, Formal analysis.

Funding

The authors received no financial support for the research, authorship, and/or publication of this article.

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Data Availability Statement

Data will be made available on request.*

References

Meda

Marra

Galfetti

, et al. Nano-composites for rocket solid propellants. Compos Sci Technol 2005; 65: 769–773. https://doi.org/10.1016/j.compscitech.2004.10.016

Deng

Shen

, et al. Long-term stability analysis and failure mode of composite solid propellant: experiment and modeling. Compos Struct 2025; 356: 118878. https://doi.org/10.1016/j.compstruct.2025.118878

Chen

Hart

Romero-Calvo

. Thermal propellant management devices for CubeSat cold gas propulsion. Acta Astronaut 2026; 243: 291–304. https://doi.org/10.1016/j.actaastro.2026.01.057

Zhang

Lei

, et al. Research on the low-frequency fatigue behavior of NEPE solid composite propellant based on fractional derivative constitutive model. Int J Solid Struct 2024; 300: 112931. https://doi.org/10.1016/j.ijsolstr.2024.112931

Duan

Yang

, et al. Experimental and numerical investigation on dynamic deformation and ignition response of NEPE-based propellants. Polym Test 2022; 116: 107814. https://doi.org/10.1016/j.polymertesting.2022.107814

Chen

Xue

, et al. Research on microstructure-driven heat transfer and thermal instability in heterogeneous NEPE propellants. Int Commun Heat Mass Tran 2026; 172: 110471. https://doi.org/10.1016/j.icheatmasstransfer.2025.110471

Dong

Wang

Zhang

, et al. Mesoscopic damage evolution and rate-dependent porosity prediction model of a composite propellant by in situ X-ray tomography. Compos Sci Technol 2025; 270: 111267. https://doi.org/10.1016/j.compscitech.2025.111267

Lei

Wang

Cheng

, et al. A constitutive model of the solid propellants considering the interface strength and dewetting. Compos Sci Technol 2020; 185: 107893. https://doi.org/10.1016/j.compscitech.2019.107893

Hou

Zhou

, et al. Microstructural simulations of debonding, nucleation, and crack propagation in an HMX-MDB propellant. Mater Des 2021; 207: 109854. https://doi.org/10.1016/j.matdes.2021.109854

10.

Nie

Liang

Zhang

, et al. Evolution of structural damage of solid composite propellants under slow heating and effect on combustion characteristics. J Mater Res Technol 2023; 25: 5021–5037. https://doi.org/10.1016/j.jmrt.2023.07.012

11.

Wang

Qiang

Wang

, et al. Macro-microscopic study on the damage threshold strain of particle-filled polymer composites. Polym Test 2025; 145: 108762. https://doi.org/10.1016/j.polymertesting.2025.108762

12.

Wubuliaisan

Hou

, et al. Viscoelastic debonding criterion-based interface for modeling the mechanical behavior of solid propellants subjected to large deformation. Eur J Mech A/solids 2023; 98: 104873. https://doi.org/10.1016/j.euromechsol.2022.104873

13.

van Ramshorst

MCJ

Di Benedetto

Duvalois

, et al. Investigation of the failure mechanism of HTPB/AP/Al propellant by In-situ uniaxial tensile experimentation in SEM. Propellants, Explos Pyrotech 2016; 41: 700–708. https://doi.org/10.1002/prep.201500264

14.

Buffiere

Maire

Adrien

, et al. In situ experiments with X ray tomography: an attractive tool for experimental mechanics. Exp Mech 2010; 50: 289–305. https://doi.org/10.1007/s11340-010-9333-7

15.

. Time-resolved 3D CT with failure physics (5D CT) for quantifying internal defect evolution behaviors of criticality safety AM components. Eng Fract Mech 2025; 319: 111036. https://doi.org/10.1016/j.engfracmech.2025.111036

16.

Xiao

Withers

. The imaging of failure in structural materials by synchrotron radiation X-ray microtomography. Eng Fract Mech 2017; 182: 127–156. https://doi.org/10.1016/j.engfracmech.2017.07.027

17.

Withers

Preuss

. Fatigue and damage in structural materials studied by X-Ray tomography. Annu Rev Mater Res 2012; 42: 81–103. https://doi.org/10.1146/annurev-matsci-070511-155111

18.

Liu

Qian

Wang

, et al. In situ X-ray tomography study on internal damage evolution of solid propellant for carrier rockets. Mater Sci Eng, A 2023; 882: 145451. https://doi.org/10.1016/j.msea.2023.145451

19.

Lai

Sang

Bian

, et al. Interfacial debonding and cracking in a solid propellant composite under uniaxial tension: an in situ synchrotron X-ray tomography study. Compos Sci Technol 2024; 256: 110743. https://doi.org/10.1016/j.compscitech.2024.110743

20.

Geng

Qiang

Wang

, et al. Macroscopic and mesoscopic properties of HTPB propellant under low temperature dynamic biaxial compression loading. Polym Test 2023; 119: 107922. https://doi.org/10.1016/j.polymertesting.2023.107922

21.

Liao

Zhao

, et al. Investigations of the failure behaviors of open-cell copper foam based on in-situ X-ray tomography compression experiments and image reconstructed finite element modeling. Eng Fract Mech 2022; 263: 108323. https://doi.org/10.1016/j.engfracmech.2022.108323

22.

Drodge

Williamson

. Understanding damage in polymer-bonded explosive composites. J Mater Sci 2016; 51: 668–679. https://doi.org/10.1007/s10853-013-7378-6

23.

Tan

Liu

Huang

, et al. The cohesive law for the particle/matrix interfaces in high explosives. J Mech Phys Solid 2005; 53: 1892–1917. https://doi.org/10.1016/j.jmps.2005.01.009

24.

Pei

Qiang

Wang

, et al. Mesoscopic failure behavior of HTPB propellant bonding interface under multi-angle pull-and-shear loading. Polym Test 2024; 132: 108365. https://doi.org/10.1016/j.polymertesting.2024.108365

25.

Kerschen

Sorensen

Guo

, et al. X-Ray phase contrast imaging of the impact of a single HMX particle in a polymeric matrix. Propellants, Explos Pyrotech 2019; 44: 447–454. https://doi.org/10.1002/prep.201800002

26.

Lei

Ren

Chen

, et al. Correlation between solid propellant failure and interface debonding in solid rocket motors. Polym Test 2022; 115: 107755. https://doi.org/10.1016/j.polymertesting.2022.107755

27.

Gallier

Hiernard

. Microstructure of composite propellants using simulated packings and X-Ray tomography. J Propul Power 2008; 24: 154–157. https://doi.org/10.2514/1.30454

28.

Ranjan

Murthy

. A review on material characterization of composite solid propellant. J Energetic Mater 2023; 43: 1–35. https://doi.org/10.1080/07370652.2023.2287466

29.

Chen

Wang

Zhang

, et al. A nonlinear viscoelastic constitutive model for solid propellant with rate-dependent cumulative damage. Materials 2022; 15: 5834. https://doi.org/10.3390/ma15175834

30.

Gao

Tang

Yang

, et al. Parameter determination and verification of ZWT viscoelastic dynamic constitutive model of Al/Ep/W material considering strain rate effect. Int J Impact Eng 2024; 184: 104816. https://doi.org/10.1016/j.ijimpeng.2023.104816

31.

Wang

Cheng

Lei

, et al. Mechanical properties and constitutive model of a composite solid propellant under the synergistic effects of accelerated aging time, pre-strain, and damage growth. Mech Mater 2020; 148: 103453. https://doi.org/10.1016/j.mechmat.2020.103453

32.

Tunç

Özüpek

. Constitutive modeling of solid propellants for three dimensional nonlinear finite element analysis. Aero Sci Technol 2017; 69: 290–297. https://doi.org/10.1016/j.ast.2017.06.025

33.

Tong

Doghri

, et al. A nonlinear viscoelastic constitutive model for cyclically loaded solid composite propellant. Int J Solid Struct 2020; 198: 126–135. https://doi.org/10.1016/j.ijsolstr.2020.04.011

34.

Schapery

. Correspondence principles and a generalizedJ integral for large deformation and fracture analysis of viscoelastic media. Int J Fract 1984; 25: 195–223. https://doi.org/10.1007/BF01140837

35.

Schapery

. Analysis of damage growth in particulate composites using a work potential. Compos Eng 1991; 1: 167–182. https://doi.org/10.1016/0961-9526(91)90017-M

36.

Fazekas

Goda

. New numerical stress solutions to calibrate hyper-visco-pseudo-elastic material models effectively. Mater Des 2020; 194: 108861. https://doi.org/10.1016/j.matdes.2020.108861

37.

Chen

Wang

, et al. Thermo-damage-viscoelastic constitutive model of HTPB composite propellant. Int J Solid Struct 2014; 51: 3209–3217. https://doi.org/10.1016/j.ijsolstr.2014.05.024

38.

Yun

K-S

Park

J-B

Jung

G-D

, et al. Viscoelastic constitutive modeling of solid propellant with damage. Int J Solid Struct 2016; 80: 118–127. https://doi.org/10.1016/j.ijsolstr.2015.10.028

39.

Jung

G-D

Youn

S-K

Kim

B-K

. A three-dimensional nonlinear viscoelastic constitutive model of solid propellant. Int J Solid Struct 2000; 37: 4715–4732. https://doi.org/10.1016/S0020-7683(99)00180-8

40.

Özüpek

Becker

. Constitutive equations for solid propellants. J Eng Mater Technol 1997; 119: 125–132. https://doi.org/10.1115/1.2805983

41.

Xing

Wang

Zhang

, et al. Mechanical behavior and constitutive model of NEPE solid propellant in finite deformation. Mech Mater 2022; 172: 104383. https://doi.org/10.1016/j.mechmat.2022.104383

42.

Liu

Wang

Liu

, et al. Dewetting stress of solid propellant under tensile loading In: 2018 9th international conference on mechanical and aerospace engineering ICMAE, Budapest, Hungary, 2018. https://doi.org/10.1109/ICMAE.2018.8467648

43.

Zhang

Dong

Zhai

, et al. Damage analysis of solid propellants with default defects based on macro-microscopic approach. Mater Des 2024; 245: 113261. https://doi.org/10.1016/j.matdes.2024.113261

44.

Park

Schapery

. A viscoelastic constitutive model for particulate composites with growing damage. Int J Solid Struct 1997; 34: 931–947. https://doi.org/10.1016/S0020-7683(96)00066-2

45.

Toulemonde

P-A

Diani

Gilormini

, et al. On the account of a cohesive interface for modeling the behavior until break of highly filled elastomers. Mech Mater 2016; 93: 124–133. https://doi.org/10.1016/j.mechmat.2015.09.014

46.

Jabalquinto

Mai

Mohammed

, et al. Influence of humidity on the decomposition process of double base rocket propellant. Polym Degrad Stabil 2023; 218: 110546. https://doi.org/10.1016/j.polymdegradstab.2023.110546

47.

Liu

Chen

, et al. A one-dimensional model for the aging and damage effects on the response of composite double-base propellants. Mech Time-Dependent Mater 2023; 27: 1007–1023. https://doi.org/10.1007/s11043-022-09554-z