# A Step-by-Step Equivalent Microprediction Method for the Mechanical Properties of Composite Solid Propellants considering Dewetting Damage

Authors: Huiru Cui; Xuan Lv; Yurong Xu; Zhiwen Zhong; Zixiang Zhou; Weili Ma
Journal: International Journal of Aerospace Engineering
 (2022)

---

## Abstract

Reliable prediction of the macromechanical properties of composite solid propellants in the microscale can accelerate the development of propellants with high mechanical properties. According to the characteristics of the composition ratio of a four-component hydroxyl-terminated polybutadiene (HTPB) propellant with the component ammonium perchlorate (AP), hydroxyl-terminated polybutadiene, aluminum powder (AL), and cyclotrimethylenetrinitramine (or RDX for short), an improved random delivery algorithm was developed to separately model filler particles with the different sizes. A step-by-step equivalent representative volume element (RVE) model was generated to reflect the microstructures of the propellant. The isotropy and uniformity of the RVE model were also tested using a two-point probability function. The Park-Paulino-Roesler (PPR) cohesive model was introduced to simulate the particle debonding (or dewetting) in solid propellant. The stress-strain curves of the propellant were obtained by the macroscopic test with the extension rate 200 mm/min at different temperatures. Based on these experimental data, the 8 characteristic parameters suitable for the microinterface of the propellant were obtained by using an inversion optimization method. A microscale finite element prediction model of the propellant considering dewetting damage was constructed to study the evolution process of the microdamage of the propellant. The predicted stress-strain curves of the propellant under different loading conditions were in good agreement with the test results.


---

## Body

## 1. Introduction

The development of propellant formulation is a prerequisite to significantly improve the design level of the solid rocket motor. Generally, the development of composite solid propellants mainly relies on empirical or semiempirical methods. It is necessary to conduct a huge number of experimental studies to determine the mechanical properties of the propellant. A composite solid propellant is a kind of particle-reinforced composite material with a high filling ratio, whose macroscopic mechanical properties strongly depend on its microstructure characteristics [1, 2]. In addition, composite solid propellants are manufactured by blending filler particles with polymer. Thus, there exist a huge number of debonding interfaces in the propellant. Studies have shown that debonding at the microinterface of composite solid propellant is an important reason for the nonlinearity of its macroscopic mechanical properties [3, 4]. Therefore, it is urgent to use the microprediction method for studying the mechanical properties of the composite solid propellant with the dewetting damage. The optimization design of microcharacteristic parameters of the propellant can be realized to effectively guide the adjustment and control of the mechanical properties of the propellant, accelerate the development of high-mechanical-performance composite solid propellants, and promote the preparation of high-loading fraction, long-life, and wide-temperature adaptable solid motors. This study can further realize the charge design of the solid rocket motor in the engineering application [5–7].Microprediction methods of composite material mechanical properties mainly have analytical and numerical methods. Analytical methods, such as the Eshelby equivalent inclusion method [8], self-consistent method [9], and Mori-Tanaka method [10], have low calculation accuracy and application limitations, in which it is difficult to describe the composite solid propellants with complex microstructures. Numerical methods have been widely used in the study of micromechanics of composite materials because they can directly construct the microstructure of materials. Inglis [11] combined the homogenization method with the finite element method to study the relationship between the microstructure characteristics of composite materials and the macroscopic equivalent mechanical properties. Ma et al. [12] carried out a multistep homogenization modeling of filler particles of different sizes and types of the matrix. It is proved that the multistep finite element method has a high accuracy and can replace the full-scale finite element method. The equivalent elastic modulus is predicted by the multistep prediction method.Due to the limitations of traditional fracture mechanics models in the study of interface debonding problems, damage mechanic models are commonly used to characterize the cracking and failure processes of the bonding interface. The most popular model is the cohesive model based on damage theory [13, 14]. By introducing the bilinear cohesion model, Chang et al. [15] simulated the damage and failure process of the microadhesive interface and studied the influence of different interface characteristics on the mechanical behavior. Based on the bilinear cohesion model, Zhi et al. [16] analyzes the process of microscopic damage, evolution, polymerization, and macroscopic crack formation of composite solid propellants when the solid content is different, and the effect of this process on the nonlinear mechanical properties of the propellant. As the most typical cohesive model, the bilinear cohesive model is only suitable for the ideal situation of debonding along the normal or tangential direction of the interface, which cannot distinguish the two directions when dealing with mixed interface debonding [17]. Compared with the bilinear cohesion model, the Park-Paulino-Roesler (PPR) cohesive model considers the influence of fracture energy in different directions on the fracture mode, whose constitutive relationship is more complete. The PPR cohesive model can accurately characterize the damage evolution of the interface under complex conditions.According to the characteristics of the composition ratio of a four-component hydroxyl-terminated polybutadiene (HTPB) propellant and the real debonding interface in the microstructure of the propellant, this paper proposes a step-by-step equivalent microprediction method for the mechanical properties of composite solid propellants considering dewetting damage: Based on the improved random delivery algorithm, the propellant representative volume element (RVE) model is generated using a stepwise equivalent strategy. The two-point probability function is used to test the isotropy and uniformity of the stepwise equivalent RVE model, based on the stress-strain curve of the propellant when the extension rate is 200 mm/min at different temperatures. By introducing a PPR cohesive model for the debonding interface, an experiment-based inversion optimization method is used to obtain characteristic parameters for the propellant microinterface. The microscopic finite element prediction models considering dewetting damage are constructed. The damage evolution process of the propellant microinterface is analyzed. The stress-strain curves of the propellant at different temperatures and tensile rates of 2 mm/min, 100 mm/min, and 500 mm/min are predicted. A reliable prediction of the macroscopic mechanical properties of the propellant on a microscale is realized.


## 2. Microgeometric Model of Propellant

### 2.1. Propellant Geometry Modeling

A four-component hydroxyl-terminated polybutadiene (HTPB) propellant is composed of ammonium perchlorate (AP), hydroxyl-terminated polybutadiene, aluminum powder (AL), and cyclotrimethylenetrinitramine (or RDX for short). This article only considers the distribution of filler particles in the propellant microstructure and studies the auxiliary additives into the polymer (HTPB). To establish the microgeometric model of the four-component HTPB propellant, the mass fraction of each component in the propellant formula information needs to be converted into the corresponding volume fraction according to Equation (1). The distribution ratio of each group is shown in Table 1. Among them, RDX and Al are filled with particles with a single grade, while AP is filled with particles with three grades.
(1)vi=mi/ρi∑i=1nmi/ρi,where vi is the volume fraction of each component, mi is the mass fraction of each component, ρi is the density of each component, and n is the total number of types of propellant components.Table 1
The distribution ratio of each group of a four-component HTPB propellant.
ComponentRDXAP1#AP2#AP3#AlPolymer (HTPB)Mass fraction (%)1542771811Density (g/cm3)1.821.951.951.952.70.908Particle size (μm)1063001001024—Volume fraction (%)15.3440.106.686.6812.4118.79Considering the randomness and arbitrariness of the components of the composite solid propellant in the microscale, this paper adopts an improved random delivery algorithm to generate an RVE model which can effectively reflect the microstructure of the four-component HTPB propellant. This algorithm adds a multidirectional settlement algorithm to the random placement algorithm, which greatly improves the applicability and efficiency of the random placement algorithm [18].Because the four-component HTPB propellant has a relatively high particle content and a large particle size span, where the particle content of AP1# with the largest particle size in the RVE model being only 2%, the distribution area of the filled particles corresponds to 593.76μm×593.76μm. At this time, according to the proportion of each component of the propellant, the particle content of RDX and AP2# in the model can be determined to be 6% and 3%, respectively. For AP3# and Al with smaller particle sizes, the particle content in the model will be more, 300 μm and 97%, respectively, as shown in Figure 1.Figure 1
RVE model of a four-component HTPB propellant (unit:μm).Figure1 shows the RVE model of a four-component HTPB propellant established directly based on the composition ratio. In this model, AP1# has fewer particles, while AP3# has more particles. If the established model has representative, it should include enough AP1#. Thus, the size of the entire RVE model will be too large. When the model size is large, the structure of the RVE model becomes quite complicated due to the inclusion of too many AP3#, which makes the subsequent meshing and numerical calculation extremely difficult.To solve this problem, a step-by-step equivalent strategy is proposed to generate the RVE model of a four-component HTPB propellant. This method is mainly to carry out different microgeometric modeling of propellant particles with different sizes, as shown in Figure2. Figure 2(a) shows the RVE model of the propellant containing only AP3# and Al particles of two small sizes. The distribution area which is filled with particles is 342.89μm×342.89μm, where AP3# has a particle content of 100 μm and Al has a particle content of 32%. Figure 2(b) shows that the RVE model of the propellant contains only three large-size particles of RDX, AP1#, and AP2#. The distribution area of the filler particles is 1327.68μm×1327.68μm, where the particle content of RDX is 31%, the particle content of AP1# is 10%, and the particle content of AP2# is 15%.Figure 2
Step-by-step equivalent RVE model of a four-component HTPB propellant (unit:μm). (a) Contains small-sized particles. (b) Contains large-sized particles.
(a)(b)


### 2.2. Statistical Characteristics Analysis of Geometric Model

The microstructure model of composite materials can describe the relationship between the phases in the model throughn-point probability functions [19]. Define the s-phase indicator function as
(2)χsx,α=1ifxinphases0otherwise,where x is the placement point and α is the sample of the placement point.The volume average of thes-phase indicator function is
(3)χsx,α¯=∫χsx,αραdα,where ρα is the probability density.The probability function ofn-point is
(4)Ssp⋯qx1,x2,⋯xn=χsx1χpx2⋯χqxn¯.Therefore, then-point probability function Ssp⋯qx1,x2,⋯xn represents the probability of simultaneously finding the s,p,⋯q phase in the model at the placement x1,x2,⋯xn point.For the two-dimensional microstructure model of composite materials, a two-point probability function is generally used to analyze the statistical characteristics of the model. The two-point probability function is expressed as(5)Sspx,x′=Sspx−x′=Sspr,where r is the distance between the placement points x and x′.The step-by-step equivalent RVE model for a certain four-component HTPB propellant is shown in Figure2. Place 10,000 sample points on a circle with different radial dimensions centered on a random point in the model. The statistical distribution characteristics of the stepwise equivalent RVE model are analyzed further to check the isotropy and uniformity of the model.


#### 2.2.1. Isotropic

The two-point probability function circumferential distribution standard deviation curve of the stepwise equivalent RVE model is shown in Figure3. In order to facilitate identification, the legend subscripts R, 1, 2, 3, A, and M, respectively, represent the RDX, AP1#, AP2#, AP3#, Al, and polymer component phases of the four-component HTPB propellant. It can be seen that the six standard deviation curves of the RVE model with small particles are all within 2%, while the 10 standard deviation curves of the RVE model with large particles are all within 2.5%. This shows that the two-point probability function of the step-by-step equivalent RVE model is basically independent of the direction, indicating that the established step-by-step equivalent RVE model meets the requirements of isotropy.Figure 3
The two-point probability function circumferential distribution standard deviation curve of the stepwise equivalent RVE model. (a) Containing small-sized particles. (b) Containing large-sized particles.
(a)(b)


#### 2.2.2. Uniformity

Figure4 shows the two-point probability function curve of the stepwise equivalent RVE model. When the two components are the same, the two-point probability function of the model gradually decreases from a certain initial value as the distance between the two points increases. When the two components are different, the two-point probability function of the model gradually increases from zero with the increase of distance between the two points. This shows that for two points in the stepwise equivalent RVE model, as the distance increases, the probability of that in the same component phase decreases, while the probability of that in a different component phase increases.Figure 4
Two-point probability function curve of the stepwise equivalent RVE model. (a) Containing small-sized particles. (b) Containing large-sized particles.
(a)(b)Figure5 shows the two-point probability function gradient curves of the stepwise equivalent RVE model. When the distance between two points in the model increases to a certain value, the two-point probability function gradient curves of the model all approach zero. At this time, the components of the four-component HTPB propellant are evenly distributed in the model, and the correlation between the components has nothing to do with the distance. These phenomena indicate that the established stepwise equivalent RVE model can satisfy the uniformity requirements.Figure 5
Two-point probability function gradient curve of the stepwise equivalent RVE model. (a) Containing small-sized particles. (b) Containing large-sized particles.
(a)(b)


## 2.1. Propellant Geometry Modeling

A four-component hydroxyl-terminated polybutadiene (HTPB) propellant is composed of ammonium perchlorate (AP), hydroxyl-terminated polybutadiene, aluminum powder (AL), and cyclotrimethylenetrinitramine (or RDX for short). This article only considers the distribution of filler particles in the propellant microstructure and studies the auxiliary additives into the polymer (HTPB). To establish the microgeometric model of the four-component HTPB propellant, the mass fraction of each component in the propellant formula information needs to be converted into the corresponding volume fraction according to Equation (1). The distribution ratio of each group is shown in Table 1. Among them, RDX and Al are filled with particles with a single grade, while AP is filled with particles with three grades.
(1)vi=mi/ρi∑i=1nmi/ρi,where vi is the volume fraction of each component, mi is the mass fraction of each component, ρi is the density of each component, and n is the total number of types of propellant components.Table 1
The distribution ratio of each group of a four-component HTPB propellant.
ComponentRDXAP1#AP2#AP3#AlPolymer (HTPB)Mass fraction (%)1542771811Density (g/cm3)1.821.951.951.952.70.908Particle size (μm)1063001001024—Volume fraction (%)15.3440.106.686.6812.4118.79Considering the randomness and arbitrariness of the components of the composite solid propellant in the microscale, this paper adopts an improved random delivery algorithm to generate an RVE model which can effectively reflect the microstructure of the four-component HTPB propellant. This algorithm adds a multidirectional settlement algorithm to the random placement algorithm, which greatly improves the applicability and efficiency of the random placement algorithm [18].Because the four-component HTPB propellant has a relatively high particle content and a large particle size span, where the particle content of AP1# with the largest particle size in the RVE model being only 2%, the distribution area of the filled particles corresponds to 593.76μm×593.76μm. At this time, according to the proportion of each component of the propellant, the particle content of RDX and AP2# in the model can be determined to be 6% and 3%, respectively. For AP3# and Al with smaller particle sizes, the particle content in the model will be more, 300 μm and 97%, respectively, as shown in Figure 1.Figure 1
RVE model of a four-component HTPB propellant (unit:μm).Figure1 shows the RVE model of a four-component HTPB propellant established directly based on the composition ratio. In this model, AP1# has fewer particles, while AP3# has more particles. If the established model has representative, it should include enough AP1#. Thus, the size of the entire RVE model will be too large. When the model size is large, the structure of the RVE model becomes quite complicated due to the inclusion of too many AP3#, which makes the subsequent meshing and numerical calculation extremely difficult.To solve this problem, a step-by-step equivalent strategy is proposed to generate the RVE model of a four-component HTPB propellant. This method is mainly to carry out different microgeometric modeling of propellant particles with different sizes, as shown in Figure2. Figure 2(a) shows the RVE model of the propellant containing only AP3# and Al particles of two small sizes. The distribution area which is filled with particles is 342.89μm×342.89μm, where AP3# has a particle content of 100 μm and Al has a particle content of 32%. Figure 2(b) shows that the RVE model of the propellant contains only three large-size particles of RDX, AP1#, and AP2#. The distribution area of the filler particles is 1327.68μm×1327.68μm, where the particle content of RDX is 31%, the particle content of AP1# is 10%, and the particle content of AP2# is 15%.Figure 2
Step-by-step equivalent RVE model of a four-component HTPB propellant (unit:μm). (a) Contains small-sized particles. (b) Contains large-sized particles.
(a)(b)


## 2.2. Statistical Characteristics Analysis of Geometric Model

The microstructure model of composite materials can describe the relationship between the phases in the model throughn-point probability functions [19]. Define the s-phase indicator function as
(2)χsx,α=1ifxinphases0otherwise,where x is the placement point and α is the sample of the placement point.The volume average of thes-phase indicator function is
(3)χsx,α¯=∫χsx,αραdα,where ρα is the probability density.The probability function ofn-point is
(4)Ssp⋯qx1,x2,⋯xn=χsx1χpx2⋯χqxn¯.Therefore, then-point probability function Ssp⋯qx1,x2,⋯xn represents the probability of simultaneously finding the s,p,⋯q phase in the model at the placement x1,x2,⋯xn point.For the two-dimensional microstructure model of composite materials, a two-point probability function is generally used to analyze the statistical characteristics of the model. The two-point probability function is expressed as(5)Sspx,x′=Sspx−x′=Sspr,where r is the distance between the placement points x and x′.The step-by-step equivalent RVE model for a certain four-component HTPB propellant is shown in Figure2. Place 10,000 sample points on a circle with different radial dimensions centered on a random point in the model. The statistical distribution characteristics of the stepwise equivalent RVE model are analyzed further to check the isotropy and uniformity of the model.


### 2.2.1. Isotropic

The two-point probability function circumferential distribution standard deviation curve of the stepwise equivalent RVE model is shown in Figure3. In order to facilitate identification, the legend subscripts R, 1, 2, 3, A, and M, respectively, represent the RDX, AP1#, AP2#, AP3#, Al, and polymer component phases of the four-component HTPB propellant. It can be seen that the six standard deviation curves of the RVE model with small particles are all within 2%, while the 10 standard deviation curves of the RVE model with large particles are all within 2.5%. This shows that the two-point probability function of the step-by-step equivalent RVE model is basically independent of the direction, indicating that the established step-by-step equivalent RVE model meets the requirements of isotropy.Figure 3
The two-point probability function circumferential distribution standard deviation curve of the stepwise equivalent RVE model. (a) Containing small-sized particles. (b) Containing large-sized particles.
(a)(b)


### 2.2.2. Uniformity

Figure4 shows the two-point probability function curve of the stepwise equivalent RVE model. When the two components are the same, the two-point probability function of the model gradually decreases from a certain initial value as the distance between the two points increases. When the two components are different, the two-point probability function of the model gradually increases from zero with the increase of distance between the two points. This shows that for two points in the stepwise equivalent RVE model, as the distance increases, the probability of that in the same component phase decreases, while the probability of that in a different component phase increases.Figure 4
Two-point probability function curve of the stepwise equivalent RVE model. (a) Containing small-sized particles. (b) Containing large-sized particles.
(a)(b)Figure5 shows the two-point probability function gradient curves of the stepwise equivalent RVE model. When the distance between two points in the model increases to a certain value, the two-point probability function gradient curves of the model all approach zero. At this time, the components of the four-component HTPB propellant are evenly distributed in the model, and the correlation between the components has nothing to do with the distance. These phenomena indicate that the established stepwise equivalent RVE model can satisfy the uniformity requirements.Figure 5
Two-point probability function gradient curve of the stepwise equivalent RVE model. (a) Containing small-sized particles. (b) Containing large-sized particles.
(a)(b)


## 2.2.1. Isotropic

The two-point probability function circumferential distribution standard deviation curve of the stepwise equivalent RVE model is shown in Figure3. In order to facilitate identification, the legend subscripts R, 1, 2, 3, A, and M, respectively, represent the RDX, AP1#, AP2#, AP3#, Al, and polymer component phases of the four-component HTPB propellant. It can be seen that the six standard deviation curves of the RVE model with small particles are all within 2%, while the 10 standard deviation curves of the RVE model with large particles are all within 2.5%. This shows that the two-point probability function of the step-by-step equivalent RVE model is basically independent of the direction, indicating that the established step-by-step equivalent RVE model meets the requirements of isotropy.Figure 3
The two-point probability function circumferential distribution standard deviation curve of the stepwise equivalent RVE model. (a) Containing small-sized particles. (b) Containing large-sized particles.
(a)(b)


## 2.2.2. Uniformity

Figure4 shows the two-point probability function curve of the stepwise equivalent RVE model. When the two components are the same, the two-point probability function of the model gradually decreases from a certain initial value as the distance between the two points increases. When the two components are different, the two-point probability function of the model gradually increases from zero with the increase of distance between the two points. This shows that for two points in the stepwise equivalent RVE model, as the distance increases, the probability of that in the same component phase decreases, while the probability of that in a different component phase increases.Figure 4
Two-point probability function curve of the stepwise equivalent RVE model. (a) Containing small-sized particles. (b) Containing large-sized particles.
(a)(b)Figure5 shows the two-point probability function gradient curves of the stepwise equivalent RVE model. When the distance between two points in the model increases to a certain value, the two-point probability function gradient curves of the model all approach zero. At this time, the components of the four-component HTPB propellant are evenly distributed in the model, and the correlation between the components has nothing to do with the distance. These phenomena indicate that the established stepwise equivalent RVE model can satisfy the uniformity requirements.Figure 5
Two-point probability function gradient curve of the stepwise equivalent RVE model. (a) Containing small-sized particles. (b) Containing large-sized particles.
(a)(b)


## 3. Propellant Microscopic Finite Element Prediction Models considering Dewetting Damage

### 3.1. Material Parameters

To construct microscopic finite element prediction models of a four-component HTPB propellant, the material parameters of each component of the propellant need to be known, including the elastic modulus, Poisson’s ratio of the filler particles, and the relaxation modulus of the polymer. Among them, the material parameters of the filler particles are selected according to the literature [4]. The relaxation modulus of the polymer is based on the literature [20] to design the stress relaxation test plan of the matrix. The test machine will stretch the specimen to the initial strain of 5% at an instantaneous speed of 500 mm/min, maintain this state for at least 1000 s, and automatically record the relaxation force curve of the polymer over time during the stress relaxation test. To prevent matrix sample slippage from fixture during sudden loading, the test used a special custom fixture as shown in Figure 6(a). Stress relaxation tests at five temperature points of 60°C, 40°C, 25°C, -20°C, and -30°C were carried out, respectively. The relaxation force of the polymer over time at different temperatures was conducted. Figure 6(b) shows the specimen dimensions in detail. Figure 7 is the master curve of the polymer relaxation modulus when the reference temperature is 25°C fitted by the Prony series.Figure 6
(a) Stress relaxation test of polymer. (b) Specimen dimensions (unit: mm).
(a)(b)Figure 7
The master curve of the polymer relaxation modulus at 25°C fitted by the Prony series.


### 3.2. Interface Parameters

The cohesive model traction-separation displacement relationship is simple and clear, which is conveniently realized and developed by using finite element method, which is the mainstream choice for numerical analysis of bonding interfaces. The PPR cohesive force model was jointly proposed by Park et al. [21]. Each fracture mode contains 8 fracture parameters: fracture energy ϕn and ϕt, cohesive strengths σmax and τmax, shape index parameters α and β, and initial stiffness indicators λn and λt [22–24]. In this paper, based on the UEL interface of ABAQUS, the PPR cohesion unit is developed [25]. The PPR cohesive model is introduced at the debonding interface of the propellant microscopic finite element prediction models.To obtain the characteristic parameters suitable for the microinterface of the four-component HTPB propellant, an inversion optimization method is adopted based on experiments. According to the propellant stress-strain macroscopic results obtained in the experiment, the interface parameters are continuously modified to make the error between the numerical simulation results, and the macroscopic test results meet the requirements, so as to realize the inversion optimization of the interface parameters. The experiment-based inversion optimization method is mainly composed of three parts, including macroscopic experiment, numerical simulation, and optimization algorithm. The macrotest is a standard unidirectional constant-speed tensile test of the propellant to obtain the stress-strain curve of the propellant macrotest. Numerical simulation is to introduce the optimized interface parameters into the microscopic finite element prediction models of the propellant to obtain the stress-strain curve of the propellant numerical simulation. The optimization algorithm is to make the error between the numerical simulation curve and the macrotest curve meet the research requirements, which is used to obtain optimized interface parameters. Because the results obtained by numerical simulation and macroscopic experiments are some discrete data points, the Hooke-Jeeves algorithm [26] can be used to directly obtain interface parameters. The objective function f of the inversion optimization as the mean square error can be defined as
(6)minf=1N∑i=1NPΔisim−PΔiexp2,where N is the total number of discrete data points and PΔisim and PΔiexp are the stresses corresponding to numerical simulation and macroscopic experiment, respectively.The specific implementation process of experiment-based inversion optimization is shown in Figure8. Since this process involves a large amount of data transmission, manual operation is obviously not advisable, and the automated realization process can ensure the efficiency and accuracy of the inversion optimization. For the PPR cohesion model suitable for the propellant micro interface, iSight can be used to establish an automatic inversion optimization task process, which is integrated by three task components, namely Simcode, MATLAB, and Optimization. As a powerful computer-aided optimization platform, iSight has various CAD/CAE integrated interfaces. In iSight, users can create complex task processes by dragging and dropping components to achieve parameterized integration of applications.Figure 8
Flow chart for inversion optimization of propellant microinterface parameters.


### 3.3. Boundary Conditions

It is necessary to apply fixed constraints at the left and lower boundaries of the two-dimensional microscopic finite element prediction models and apply a constant displacement load or constant velocity displacement load at the bearing end of the upper boundary as shown in Figure9.Figure 9
Step-by-step equivalent microscopic finite element prediction model. (a) Containing small-sized particles. (b) Containing large-sized particles.
(a)(b)By constructing a finite element prediction model containing only small-sized particles, the smaller-sized filler particles and the polymer of the propellant are homogenized and equivalent to a viscoelastic material. The viscoelastic parameters of the equivalent matrix can refer to the parameters of the polymer. The initial strain of the finite element prediction model containing small particles is selected as 5%. The constant displacement load applied to the upper boundary of the model is 17.1445μm. Four kinds of extension rates are considered for the propellant constant-rate tensile test program. A constant velocity displacement load of 2 mm/min, 100 mm/min, 200 mm/min, and 500 mm/min will be applied to the upper boundary of the finite element prediction model containing large particles. The strain rate of the model corresponds to 0.0005 s-1, 0.0238 s-1, 0.0476 s-1, and 0.1190 s-1.


## 3.1. Material Parameters

To construct microscopic finite element prediction models of a four-component HTPB propellant, the material parameters of each component of the propellant need to be known, including the elastic modulus, Poisson’s ratio of the filler particles, and the relaxation modulus of the polymer. Among them, the material parameters of the filler particles are selected according to the literature [4]. The relaxation modulus of the polymer is based on the literature [20] to design the stress relaxation test plan of the matrix. The test machine will stretch the specimen to the initial strain of 5% at an instantaneous speed of 500 mm/min, maintain this state for at least 1000 s, and automatically record the relaxation force curve of the polymer over time during the stress relaxation test. To prevent matrix sample slippage from fixture during sudden loading, the test used a special custom fixture as shown in Figure 6(a). Stress relaxation tests at five temperature points of 60°C, 40°C, 25°C, -20°C, and -30°C were carried out, respectively. The relaxation force of the polymer over time at different temperatures was conducted. Figure 6(b) shows the specimen dimensions in detail. Figure 7 is the master curve of the polymer relaxation modulus when the reference temperature is 25°C fitted by the Prony series.Figure 6
(a) Stress relaxation test of polymer. (b) Specimen dimensions (unit: mm).
(a)(b)Figure 7
The master curve of the polymer relaxation modulus at 25°C fitted by the Prony series.


## 3.2. Interface Parameters

The cohesive model traction-separation displacement relationship is simple and clear, which is conveniently realized and developed by using finite element method, which is the mainstream choice for numerical analysis of bonding interfaces. The PPR cohesive force model was jointly proposed by Park et al. [21]. Each fracture mode contains 8 fracture parameters: fracture energy ϕn and ϕt, cohesive strengths σmax and τmax, shape index parameters α and β, and initial stiffness indicators λn and λt [22–24]. In this paper, based on the UEL interface of ABAQUS, the PPR cohesion unit is developed [25]. The PPR cohesive model is introduced at the debonding interface of the propellant microscopic finite element prediction models.To obtain the characteristic parameters suitable for the microinterface of the four-component HTPB propellant, an inversion optimization method is adopted based on experiments. According to the propellant stress-strain macroscopic results obtained in the experiment, the interface parameters are continuously modified to make the error between the numerical simulation results, and the macroscopic test results meet the requirements, so as to realize the inversion optimization of the interface parameters. The experiment-based inversion optimization method is mainly composed of three parts, including macroscopic experiment, numerical simulation, and optimization algorithm. The macrotest is a standard unidirectional constant-speed tensile test of the propellant to obtain the stress-strain curve of the propellant macrotest. Numerical simulation is to introduce the optimized interface parameters into the microscopic finite element prediction models of the propellant to obtain the stress-strain curve of the propellant numerical simulation. The optimization algorithm is to make the error between the numerical simulation curve and the macrotest curve meet the research requirements, which is used to obtain optimized interface parameters. Because the results obtained by numerical simulation and macroscopic experiments are some discrete data points, the Hooke-Jeeves algorithm [26] can be used to directly obtain interface parameters. The objective function f of the inversion optimization as the mean square error can be defined as
(6)minf=1N∑i=1NPΔisim−PΔiexp2,where N is the total number of discrete data points and PΔisim and PΔiexp are the stresses corresponding to numerical simulation and macroscopic experiment, respectively.The specific implementation process of experiment-based inversion optimization is shown in Figure8. Since this process involves a large amount of data transmission, manual operation is obviously not advisable, and the automated realization process can ensure the efficiency and accuracy of the inversion optimization. For the PPR cohesion model suitable for the propellant micro interface, iSight can be used to establish an automatic inversion optimization task process, which is integrated by three task components, namely Simcode, MATLAB, and Optimization. As a powerful computer-aided optimization platform, iSight has various CAD/CAE integrated interfaces. In iSight, users can create complex task processes by dragging and dropping components to achieve parameterized integration of applications.Figure 8
Flow chart for inversion optimization of propellant microinterface parameters.


## 3.3. Boundary Conditions

It is necessary to apply fixed constraints at the left and lower boundaries of the two-dimensional microscopic finite element prediction models and apply a constant displacement load or constant velocity displacement load at the bearing end of the upper boundary as shown in Figure9.Figure 9
Step-by-step equivalent microscopic finite element prediction model. (a) Containing small-sized particles. (b) Containing large-sized particles.
(a)(b)By constructing a finite element prediction model containing only small-sized particles, the smaller-sized filler particles and the polymer of the propellant are homogenized and equivalent to a viscoelastic material. The viscoelastic parameters of the equivalent matrix can refer to the parameters of the polymer. The initial strain of the finite element prediction model containing small particles is selected as 5%. The constant displacement load applied to the upper boundary of the model is 17.1445μm. Four kinds of extension rates are considered for the propellant constant-rate tensile test program. A constant velocity displacement load of 2 mm/min, 100 mm/min, 200 mm/min, and 500 mm/min will be applied to the upper boundary of the finite element prediction model containing large particles. The strain rate of the model corresponds to 0.0005 s-1, 0.0238 s-1, 0.0476 s-1, and 0.1190 s-1.


## 4. Predictive Model Results and Analysis

Based on the step-by-step equivalent RVE model of the four-component HTPB propellant, the microscopic finite element prediction model of the propellant is constructed twice. Combined with the homogenization theory [27], the macromechanical properties of the propellant can be reliably predicted in a micro scale. An initial strain of 5% is applied. The relaxation modulus of the equivalent matrix with time at different temperatures is obtained through numerical simulation as shown in Figure 10. Figure 11 is the master curve of the equivalent matrix relaxation modulus when the reference temperature is 25°C fitted by the Prony series.Figure 10
The simulated relaxation modulus-time curve of the equivalent matrix at different temperatures.Figure 11
The master curve of the equivalent matrix relaxation modulus at 25°C fitted by the Prony series.


### 4.1. Microdamage Development of the Propellant

According to the separation state of the filler particles and the equivalent matrix, the microdamage evolution of the four-component HTPB propellant can be divided into four stages, including unseparated stage, small separation stage, large separation stage, and separated stage.Figure12(a) is the stress distribution contour of the propellant microstructure in the unseparated stage. In the unseparated stage, the microscopic finite element indicates that the deformation of the model is small. Currently, the filler particles do not separate from the equivalent matrix. The adhesion performance at the microinterface is good. In addition, since the modulus of the equivalent matrix and the filler particles are too different, the stress value inside the filler particles is significantly higher than the stress value inside the equivalent matrix. This shows that the filler particles have a strengthening effect on the overall structure. The phenomenon of stress bridging between the filler particles indicates that the interaction between the particles is the key factor affecting the stress distribution of the microstructure.Figure 12
Stress cloud diagram of propellant microstructure. (a) Unseparated stage. (b) Small separation stage. (c) Large separation stage. (d) Separated stage.
(a)(b)(c)(d)Figure12(b) is the stress distribution contour of the propellant microstructure in the small separation stage. In the small separation stage, separation occurs between the filler particles and the equivalent matrix. When the separation displacement of some interfaces reaches the critical fracture displacement, the microinterface starts to be damaged. However, due to the small number of damaged interfaces, the overall structure still has a strong bearing capacity. In addition, the separation displacement between the large-sized particles and the equivalent matrix is significantly larger than the separation displacement between the small-sized particles and the equivalent matrix. The interface between the large-sized particles and the equivalent matrix is the main damage interface at this stage.Figure12(c) is the stress distribution contour of the propellant microstructure in the large separation stage. In the large separation stage, as the deformation of the microscopic finite element prediction models increases, the separation between the filler particles and the equivalent matrix gradually increases, which can further deepen the damage of the microinterface. This leads to a significant reduction in the carrying capacity of the overall structure. At this time, the separation displacement of most interfaces has reached the critical fracture displacement. The separation displacement of some of the interfaces reaches the final failure displacement. The filler particles corresponding to these interfaces are in a dewetting state.Figure12(d) is the stress distribution contour of the propellant microstructure in the separated stage. In the separated stage, because the deformation of the microscopic finite element prediction model is too large, the separation displacement of most interfaces reaches the final failure displacement. The interfaces between the corresponding filler particles and the matrix completely failed. The filler particles are bare and unwrapped.


### 4.2. Predicted Stress-Strain Curve Results and Verification

By constructing the microscopic finite element prediction model of the propellant considering dewetting damage, the stress-strain curve of the four-component HTPB propellant under different loading conditions is predicted. Figure13 shows the comparison between the stress-strain prediction results and the test results of the propellant under different constant-speed tensile loads. To obtain the characteristic parameters suitable for the propellant microinterface, this paper uses the stress-strain curve of the propellant when the extension rate is 200 mm/min at different temperatures to determine the interface parameters. The obtained interface parameters are shown in Table 2.Figure 13
The comparison between the predicted results of the propellant stress-strain curve and the test results under different constant-speed tensile loads. (a) 2 mm/min. (b) 100 mm/min. (c) 200 mm/min. (d) 500 mm/min.
(a)(b)(c)(d)Table 2
PPR cohesion model parameters obtained by inversion optimization.
Characteristic parameters of PPR cohesion modelϕn (N/mm)ϕt (N/mm)σmax (MPa)τmax (MPa)αβλnλt60°C0.5231.3623.531.165.755.410.1240.10825°C0.4281.1073.221.085.885.640.1220.131-30°C0.2840.9542.940.875.765.380.1250.145For the extension rates of 2 mm/min, 100 mm/min, and 200 mm/min, when the propellant is not damaged, the predicted results of the microscopic finite element prediction models considering dewetting damage are in good agreement with the test results. However, when the propellant has been damaged, the predicted results of this model are all too large. This is because only the dewetting of the filler particles is considered in the modeling, but the damage of the matrix is not considered. In addition, for a stretching rate of 500 mm/min, a higher strain rate will make the heat generated by the propellant under load unable to be released in time after a short time accumulation. This will cause the temperature inside the propellant to rise, which will lead to a large deviation between the predicted results and the test results. In general, the prediction results are basically the same as the experimental results, which can validate the prediction model and the accuracy of the interface parameter inversion.


## 4.1. Microdamage Development of the Propellant

According to the separation state of the filler particles and the equivalent matrix, the microdamage evolution of the four-component HTPB propellant can be divided into four stages, including unseparated stage, small separation stage, large separation stage, and separated stage.Figure12(a) is the stress distribution contour of the propellant microstructure in the unseparated stage. In the unseparated stage, the microscopic finite element indicates that the deformation of the model is small. Currently, the filler particles do not separate from the equivalent matrix. The adhesion performance at the microinterface is good. In addition, since the modulus of the equivalent matrix and the filler particles are too different, the stress value inside the filler particles is significantly higher than the stress value inside the equivalent matrix. This shows that the filler particles have a strengthening effect on the overall structure. The phenomenon of stress bridging between the filler particles indicates that the interaction between the particles is the key factor affecting the stress distribution of the microstructure.Figure 12
Stress cloud diagram of propellant microstructure. (a) Unseparated stage. (b) Small separation stage. (c) Large separation stage. (d) Separated stage.
(a)(b)(c)(d)Figure12(b) is the stress distribution contour of the propellant microstructure in the small separation stage. In the small separation stage, separation occurs between the filler particles and the equivalent matrix. When the separation displacement of some interfaces reaches the critical fracture displacement, the microinterface starts to be damaged. However, due to the small number of damaged interfaces, the overall structure still has a strong bearing capacity. In addition, the separation displacement between the large-sized particles and the equivalent matrix is significantly larger than the separation displacement between the small-sized particles and the equivalent matrix. The interface between the large-sized particles and the equivalent matrix is the main damage interface at this stage.Figure12(c) is the stress distribution contour of the propellant microstructure in the large separation stage. In the large separation stage, as the deformation of the microscopic finite element prediction models increases, the separation between the filler particles and the equivalent matrix gradually increases, which can further deepen the damage of the microinterface. This leads to a significant reduction in the carrying capacity of the overall structure. At this time, the separation displacement of most interfaces has reached the critical fracture displacement. The separation displacement of some of the interfaces reaches the final failure displacement. The filler particles corresponding to these interfaces are in a dewetting state.Figure12(d) is the stress distribution contour of the propellant microstructure in the separated stage. In the separated stage, because the deformation of the microscopic finite element prediction model is too large, the separation displacement of most interfaces reaches the final failure displacement. The interfaces between the corresponding filler particles and the matrix completely failed. The filler particles are bare and unwrapped.


## 4.2. Predicted Stress-Strain Curve Results and Verification

By constructing the microscopic finite element prediction model of the propellant considering dewetting damage, the stress-strain curve of the four-component HTPB propellant under different loading conditions is predicted. Figure13 shows the comparison between the stress-strain prediction results and the test results of the propellant under different constant-speed tensile loads. To obtain the characteristic parameters suitable for the propellant microinterface, this paper uses the stress-strain curve of the propellant when the extension rate is 200 mm/min at different temperatures to determine the interface parameters. The obtained interface parameters are shown in Table 2.Figure 13
The comparison between the predicted results of the propellant stress-strain curve and the test results under different constant-speed tensile loads. (a) 2 mm/min. (b) 100 mm/min. (c) 200 mm/min. (d) 500 mm/min.
(a)(b)(c)(d)Table 2
PPR cohesion model parameters obtained by inversion optimization.
Characteristic parameters of PPR cohesion modelϕn (N/mm)ϕt (N/mm)σmax (MPa)τmax (MPa)αβλnλt60°C0.5231.3623.531.165.755.410.1240.10825°C0.4281.1073.221.085.885.640.1220.131-30°C0.2840.9542.940.875.765.380.1250.145For the extension rates of 2 mm/min, 100 mm/min, and 200 mm/min, when the propellant is not damaged, the predicted results of the microscopic finite element prediction models considering dewetting damage are in good agreement with the test results. However, when the propellant has been damaged, the predicted results of this model are all too large. This is because only the dewetting of the filler particles is considered in the modeling, but the damage of the matrix is not considered. In addition, for a stretching rate of 500 mm/min, a higher strain rate will make the heat generated by the propellant under load unable to be released in time after a short time accumulation. This will cause the temperature inside the propellant to rise, which will lead to a large deviation between the predicted results and the test results. In general, the prediction results are basically the same as the experimental results, which can validate the prediction model and the accuracy of the interface parameter inversion.


## 5. Conclusions

According to the composition ratio of a four-component HTPB propellant, a step-by-step equivalent strategy is adopted to generate the RVE model of propellant containing small-sized particles and large-sized particles. The two-point probability function is used to test the isotropy and uniformity of the stepwise equivalent RVE model. A simple and effective implementation method is provided to solve the problem of the composite solid propellant with high filler particle content and large particle size span.Based on the PPR cohesive force model, a microscale finite element prediction model of the mechanical properties of the propellant considering dewetting damage is constructed to predict the stress-strain curve of the four-component HTPB propellant under different constant-speed tensile loads. The predicted results are basically the same as the experimental results, which shows that the constructed microscale finite element prediction model can better simulate the mechanical behavior of the propellant under constant-speed tensile load. A reliable prediction of the macroscopic mechanical properties of the propellant is achieved in a microscale.At present, the micronumerical studies on solid propellants are carried out on a two-dimensional level, and the real microstructure of propellants is composed of three-dimensionally filled particles with complex shapes. However, there are many difficulties in the establishment and calculation of three-dimensional models, especially the modeling of three-dimensional special-shaped filling particles, the calculation of bonding between particles/matrix, and the high hardware requirements of three-dimensional models. If the above-mentioned problems can be effectively solved, a microscale finite element prediction model that is more in line with the actual structure of the propellant can be constructed to accurately simulate the mechanical behavior of the propellant under different loads.



---

*Source: 2427463-2022-02-14.xml*