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.
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. where is the volume fraction of each component, is the mass fraction of each component, is the density of each component, and is the total number of types of propellant components.
Considering 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 . 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 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 Figure 2. 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 , 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 , where the particle content of RDX is 31%, the particle content of AP1# is 10%, and the particle content of AP2# is 15%.
(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 through -point probability functions [19]. Define the -phase indicator function as where is the placement point and is the sample of the placement point.
The volume average of the -phase indicator function is where is the probability density.
The probability function of -point is
Therefore, the -point probability function represents the probability of simultaneously finding the phase in the model at the placement 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 where is the distance between the placement points and .
The step-by-step equivalent RVE model for a certain four-component HTPB propellant is shown in Figure 2. 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 Figure 3. 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.
(a)
(b)
2.2.2. Uniformity
Figure 4 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.
(a)
(b)
Figure 5 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.
(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.
(a)
(b)
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 and , cohesive strengths and , shape index parameters and , and initial stiffness indicators and [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 of the inversion optimization as the mean square error can be defined as where is the total number of discrete data points and and are the stresses corresponding to numerical simulation and macroscopic experiment, respectively.
The specific implementation process of experiment-based inversion optimization is shown in Figure 8. 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.
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 Figure 9.
(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.
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.
Figure 12(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.
(a)
(b)
(c)
(d)
Figure 12(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.
Figure 12(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.
Figure 12(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. Figure 13 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.
(a)
(b)
(c)
(d)
For 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.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work is supported by the Natural Science Foundation of Jiangsu Province (BK20210435).