Micro-damage analysis and numerical simulation of composite solid propellant based on in situ tensile test

In order to quantitatively analyze the mesoscopic damage process of hydroxyl-terminated polybutadiene composite solid propellant under external load, periodic boundary conditions were applied to the representative volume element model based on sample composition and morphology, the mixed matrix containing aluminum powder was homogenized, and the hyperelastic matrix damage and bilinear/exponential particle– matrix interface cohesive model with initial damage were compiled through the secondary development of Abaqus. At the same time, a data interaction platform was constructed by means of Python and MATLAB, matrix and cohesion parameters were inverted according to the optimization algorithm and experimental data, and the whole process of propellant damage and fracture was simulated from the mesoscopic perspective. The results show that combining the adaptive particle swarm optimization algorithm and the Hooke–Jeeves algorithm can achieve the global optimal parameter inversion in 102 calculations, compared with the single local search algorithm, which can cut about 11% of the objective function values. Considering the matrix damage and the exponential cohesion model with initial damage, the optimal objective function value is 0.01635, which canmore accurately simulate the propellant damage and fracture process compared with 0.02136 of a bilinear cohesion model.


Introduction
Owing to its appearance, solid propellant has been widely used as the power source of solid rocket motors in the military field because of its stable performance, simple structure, and various types. Different types of equipment formulations will also change due to different emphasis, but the requirements for maintaining stable high energy output have never changed. Therefore, solid propellants are generally composed of oxidizers, combustion agents, and auxiliary materials, which are in a multiphase state geometrically. Under long-term or excessive external load, particle dehumidification, matrix deformation damage, and even fracture will occur in the interior of this multiphase composite material, making its macroscopic mechanical properties change complex. Therefore, it is greatly significant to explore the whole process of propellant damage under load from the microscopic perspective to improve the mechanical properties and maintain the structural integrity of the propellant.
With the development of observation means and computer technology, such as scanning electron microscope and digital image processing technology, the change process of propellant micromorphology under load is established, and the mechanism of propellant micromorphology damage is further clarified. However, due to the limitations of experimental means and observation technology, it is difficult to obtain the interface parameters and damage threshold directly, and it is impossible to achieve quantitative analysis. Therefore, numerical analysis is gradually applied to the mesoscopic damage process of propellant. In order to achieve the accurate acquisition of relevant parameters, in addition to the basic experiments, the construction of geometric models, the application of boundary conditions, and the setting of material constitutive equations become the key to the analysis. Collins et al. [1] and Lee et al. [2] were the first to obtain the three-dimensional micromorphology of solid propellant by using micro-CT and reconstruct the cracks caused by the loading of the propellant. Li et al. [3] discovered the change rule of propellant micro-damage morphology under different tensile strains using high-precision micro-CT and further analyzed the impact of microdamage on macroscopic mechanical properties. Zhang et al. [4] calculated the average stress of the mesoscopic model based on the finite-element theory and homogenization method and predicted the modulus of composite solid propellant. Tunc and Ozupek [5] and Yun et al. [6] applied the constitutive model including initial damage to the study of mechanical properties of solid propellants and further improved the simulation accuracy. Li et al. [7] proposed the superposition model of Young's modulus of propellant, established the four-stage damage constitutive model of propellant, and verified the rationality of the model prediction. Pan and Qu [8] used a nonlinear ultrasonic method to monitor the nonlinear dynamic response of propellant in real time to determine the effect of particle dehumidification on propellant damage. Gao and Li [9] obtained the crack growth trend and stress-strain curve under different tensile rates by prefabricating different crack angles. Hou et al. [10] established a mesoscopic crack model of high-energy propellant by using the uniform displacement boundary condition and inserting the cohesive element globally. Zhou et al. [11] realized the simulation of the propellant matrix cracking process based on the secondary development of Abaqus software. Zhi et al. [12] inversely calculated the micro-particle dehumidification process by establishing the objective function of the bonding interface damage parameter and comparing the stress-strain curve of the simulation experiment. However, the aforementioned methods are limited to the optimization of a certain aspect of the numerical analysis of propellants and focus on exploring the influence of single factors such as model, constitutive, or matrix on the mechanical properties of propellants. The aforementioned main factors are not taken into account in the microscopic analysis at the same time. The use of a single inverse algorithm is also easy to make the numerical analysis fall into local optimization, and the setting of the cohesive model is mainly based on a bilinear model so as to a significant difference between the analysis result and the experiment.
In this study, on the basis of considering the cohesive model of the grain-matrix interface and the hyperelastic matrix damage, the mesoscopic particles with different grading are periodically arranged to generate representative volume element (RVE), and periodic boundary loads are applied. At the same time, zero-thickness bonding elements are inserted into the grain-matrix interface, and the matrix damage and exponential cohesion model subprograms are compiled. A data interaction platform is built through Python and MATLAB, and the global optimization and direct search algorithm are combined and compared with the in situ tensile test data to inverse the cohesion and matrix parameters, so as to realize the simulation and quantitative analysis of the whole process of mesoscopic damage of solid propellants.

Experimental section
Compared with the macrotensile test, the microtensile test regards the propellant sample as a mixture of microparticles and matrices. In addition to focusing on the overall mechanical properties, it focuses more on the changes of matrices, microparticles, and interface in the process of external force loading. In addition, the addition of ammonium perchlorate (AP) particles and other auxiliary agents, while strengthening the matrix, also took place a certain chemical reaction, which made the mechanical changes of the interface between the matrix and the micro-particles more complex. Therefore, the micro-dynamic tensile test was carried out using a scanning electron microscope to further analyze the micro-change of the propellant under load.

Test piece and instrument
Test piece: 20 sets of dumbbell-type test piece molds can be made according to the QJ924-85 standard [13] to ensure that the thickness of the test piece is 5 mm, the width of both sides is 12.0 mm, and the width of the central load-bearing part is 5 mm. At the same time, the raw materials purchased from the 41st Institute of the Sixth Research Academy of China Aerospace Science and Technology Corporation were pre mixed into the required slurry for casting. The formed solid propellant sample uses hydroxyl-terminated polybutadiene (HTPB) as the adhesive, AP as the oxidant, aluminum (Al) powder as the combustion aid, toluene diisocyanate (TDI) as the curing agent, and combined bonding agents (aziridine polyester [PAZ], amine polyester [PAM]) are used. The mass fractions of adhesives, oxidizers, combustion aids, and other components are 11, 18.5, 64, and 6.5%, respectively. The treatment of oxidant includes crushing, batching, drying, analysis and inspection, and other processes. For the acceptance of purity, the heating titration and decomposition absorption methods are mostly used, and for the particle size index, the mechanical vibration standard screening is used. For the determination of the properties of Al powder, the oxidation-reduction titration method and mass screening method are selected.

Instrument:
Sigma300 scanning electron microscope: Gemini, Germany; ETD-2000 ion sputtering instrument: Rigaku, Japan; MICROTEST200 in situ dynamic tensile test platform: Deben Company, UK; The tensile test piece and fixture are shown in Figure 1.

Meso component
The density of Al powder and AP is 2.70 and 1.95 g cm −3 , respectively, and AP particle size is divided into 50-80 and 100-150 μm, which basically conforms to Gaussian distribution. Due to the small particle size of aluminum powder, most of it is embedded in the matrix, so the observed image is mostly AP particle distribution, as shown in Figure 2. Among them, the proportion of Al powder as auxiliary combustion agent is relatively large. However, due to its small particle size, the amount is very large. For the microparticles in the propellant, the size of the particle is proportional to the possibility of particle dehumidification, so before the aluminum powder particle dehumidification, the propellant matrix has broken [14]. Therefore, the matrix and aluminum powder mixture is homogenized according to the Mori-Tanaka model [15] to obtain the mechanical properties of the mixed matrix. According to the tensile test of the matrix film used to prepare the sample and the performance test of the filler particles conducted by the Aerospace Institute, the material parameters used in this article are shown in Table 1.

In situ tensile test
The prepared specimen is placed on the tensile fixture after spraying gold and the fixture is placed in the scanning electron microscope. After vacuuming, in-situ tensile tests will be conducted, and larger AP particles with significant damage will become the focus of observation. In addition, due to the slow imaging of SEM, the tensile rate of 0.05 mm/min is used to observe at 50 times magnification, and the tensile clamp is connected with the mechanical sensor, so as to transmit the tensile load data from time to time. At the mesoscopic level, the different particle distribution in the propellant makes the observed images of each specimen different, but the mechanical property curves and damage process of dozens of specimens in the in situ tensile test are basically consistent. In this article, the representative experiment with better observation effect is selected for mechanical property analysis, and the stress-strain curve and the corresponding micromorphology changes during the in situ tensile process are shown in Figure 3.
According to the figure, the change of mechanical properties of solid propellants can be divided into four stages: linear viscoelastic section, stress softening section, unloading section, and damage fracture section. At the initial stage of loading, the fine particles and the matrix are well bonded, the propellant structure is complete, and the stiffness value fluctuates within a certain range, which is in the linear viscoelastic stage. With the application of the external force load, the strain of the propellant further increases, and some micro-particles and the matrix are dehumidified, resulting in the   appearance of internal pores, the reduction of homogenization modulus, the weakening of the performance of the particle reinforced matrix, and the trend of stress rise has slowed down, gradually entering the stress softening stage. Under a continuous external load, a large number of micro-particles dehumidify, which causes the pores caused by dehumidification to gather and form cracks and continue to expand. Once the strain limit of the matrix is reached, it immediately enters the damage and fracture stage from the unloading section. Compared with the standard dumbbell test piece, the volume of the loaded part of the in situ tensile test piece is smaller. Once the pores caused by particle dehumidification gather, the cracks will have a greater impact on the loaded part; at the same time, in view of the improvement of the mechanical properties of the propellant, the use of the composite bonding agent further strengthens the grain-matrix bonding interface, making the linear viscoelastic section and the stress softening section larger, and the stress unloading and damage fracture stage relatively small.

RVE and periodic boundary conditions (PBC)
For the microscopic damage analysis, how to maximize the reflection of macroscopic phenomena with less calculation is the focus of attention. Therefore, RVE should contain enough information, including the composition, proportion, and distribution of fine particles, and not too large to cause a sudden increase in the calculation [16]. In general, the macropropellant is composed of numerous periodic arrays of mesoscopic parts, so when constructing RVE, the internal particles must be arranged periodically. In addition, because the propellant grain is cast and formed by premixed slurry, the distribution of fine particles is irregular, so with the help of the compiled Python statement, by setting the particle size, grading, spacing, and controlling the order of particle generation and the filling mass fraction, the plane space is filled with particles starting from the randomly generated initial fine particles, and the radial ratio can be set to control the particle shape. The mass fraction and particle size of AP and Al powder particles are basically consistent with the production data of samples provided by the aerospace research institute. In addition, on the premise of ensuring the accuracy of the finite-element model, the mixed phase composed of the propellant matrix and aluminum powder is homogenized and pretreated [17], and its mechanical properties are equivalent to the green area in Figure 4 to reduce the calculation cost. When the RVE size is more than three times the maximum particle size, the representative volume unit is considered to be effective [18]. Therefore, this article comprehensively considers the finite-element calculation cost and calculation accuracy, and the RVE size is set as 900 × 900 based on the particle size distribution rule. The RVE model is shown in Figure 4.
In the process of strengthening the mechanical properties of the matrix, which is closely related to the strength of the bonding interface of the matrix, once the particles are dehumidified, it tends to show an unloading phenomenon on the macrolevel. Moreover, the current research mainly focuses on how to strengthen the bonding interface [19], so that the micro-particles can avoid dehumidification as much as possible under the action of external forces. Therefore, when investigating the mechanical properties of propellants, the bonding interface becomes an inevitable factor. Dehumidification of the interface between the micro-particles and the matrix is actually the crack propagation process of the bonding interface. However, the traditional fracture mechanics research based on linear viscoelasticity not only needs to pre-crack but also often has singularity at the crack tip, which is difficult to reproduce the process of crack initiation and propagation. Compared with the traditional fracture mechanics, the cohesive zone model (CZM), as a phenomenological model, does not need to prefabricate cracks while greatly reducing the singularity of stress, and it can set different tension displacement relationships according to different materials, making it applicable and widely used. For CZM, in addition to the different constitutive design for different materials, it mainly includes two aspects, one is the initial damage criterion, and the other is the damage evolution law. The CZM element can be regarded as two faces separated by a certain thickness. These two faces share nodes with other solid elements. Once the strength limit is reached, it will immediately enter the stiffness-softening stage [20]. Therefore, when building the propellant model, it is necessary to embed the bonding element between the particles and the matrix to simulate the particle dehumidification process. Taking the bilinear cohesion model as an example, the element embedment and damage diagram is shown in Figure 5.
For bilinear cohesion model, its standard form is shown in the following equation: where δ 0 and δ f are initial damage displacement and maximum failure displacement, mm, respectively; T max is the maximum cohesion strength [21], N.  By establishing PBC for the RVE model, that is, creating the equation relationship between the two opposite sides and the reference point, the displacement difference between the two sides of the RVE in the loaded direction is consistent with the applied load displacement, and by pre-calculating the overall Poisson's ratio of the RVE, the shrinkage of the RVE in the direction perpendicular to the loaded direction is realized. The application of PBC ensures the continuity of stress and strain on the opposite sides of the RVE, that is, to a certain extent, the array of multiple RVEs can reflect the loading condition of the macro test piece [22]. Because the set reference point has an equation relationship with both sides of the RVE, the overall tension of the RVE can be achieved by applying displacement load to the reference point and setting the loading rate. At the same time, according to Newton's third law, the reaction force (RF) on the reference point is the load on the RVE, and the displacement is consistent with the RVE, then the forcedisplacement data can be converted into the overall stress-strain data of the RVE through calculation. The schematic diagram is shown in Figure 6.

Secondary development of ABAQUS
With the development of computer technology, finite-element analysis software is gradually applied to research various engineering problems. Because of its low resource cost and high calculation accuracy, its application scope is more and more extensive, and the problems studied are more and more complex. However, many disciplines and complex engineering problems make updating the finiteelement software difficult. Therefore, many software have secondary development ports. Secondary development realized by subroutines can further expand the application field and reduce the workload [23].
Since the cohesive force model was proposed by Barenblatt and Dugdale, for many complex materials, the cohesive force model represented by bilinear, exponential, and polynomial has evolved. With the further improvement of the finite-element software, the bilinear cohesion model has been embedded in the software. However, previous studies [24] have shown that the bilinear cohesion model is applicable to brittle materials, and the exponential cohesion model is widely applied to the interface cracking of composite materials and the crack initiation process of film coating on the ductile matrix. When using the finite-element method to analyze ductile composite materials, the calculation cost is low, and the convergence is better. The fracture energy control equation is shown in the following equation: (2) where Δ n and Δ t are the normal and tangential displacement values on the interface, mm, respectively; ϕ n is the interface fracture energy when the interface is completely cracked under the pure normal cracking state, mJ; and δ n δ t are the characteristic displacement values of normal and tangential interface cracking, that is, the displacement values corresponding to the maximum stress point, mm. The parameters q and r are where ϕ t is the pure tangential interface fracture energy, * Δ n is the normal displacement value when the normal stress is zero and the tangential crack is complete [25].
The normal and tangential stress equations can be obtained from the derivation of (2) as follows:  For the bonding interface of two-phase composite materials, due to the inevitable effect of external loads such as vibration, temperature or humidity during storage or transportation, the irrecoverable cumulative damage may occur, resulting in the weakening of its bearing capacity, so it is necessary to add damage factors and modify the above formula, as follows: Since no corresponding setting of exponential cohesion model exists in the finite-element software [26], this study uses VUMAT subroutine to write tension-displacement control statements. In addition, in the finite-element software, the essence of simulating fracture is to make the damaged element stiffness zero, make it lose the ability to bear stress, and hide it to achieve the effect of fracture. Because the hyperelastic model fits the matrix constitutive model, and its damage will break when it accumulates to a certain extent. Therefore, the VUSDFLD subroutine is used for the secondary development of ABAQUS to realize the setting of fracture strain. In the explicit calculation, the DEPVAR state variable is set for the matrix material, and the field variable is customized to make the material attribute and the state variable linked, and the state variable is used to realize the deletion of the element.

Parameter inversion 4.1 Combinatorial optimization algorithm
Because the cohesion model parameters are difficult to measure directly and can only be optimized with the help of experimental data, the data calculated by ABAQUS under the initial parameters are imported into MATLAB, which is a data interaction platform [27], and compared with the experimental data put in advance, and the cohesion parameters and hyperelastic fracture strain data are adjusted through the optimization algorithm to achieve the optimal simulation of the real experiment. Multiparameter optimization often involves multiple factors and the interaction between parameters. The set objective function may be multimodal, nonlinear, non-continuous, and nondifferentiable. Variables and constraint functions may also be linear, nonlinear, continuous, or discrete variable sets. In the face of these complex situations, traditional numerical optimization and direct search algorithms are difficult to find the required derivative and gradient information or often get the local optimal solution. Therefore, most scholars use the global search method to solve such complex optimization problems.
Although the global search method has strong adaptability and the ability of global search, its application is limited to a certain extent by too much computation. Therefore, the global search and numerical optimization algorithms are combined considering the advantages and disadvantages of global search and direct search algorithms. At first, the global optimization algorithm is used to locate the area where the target extreme value is located in the design space, and then the direct search algorithm is used to accurately optimize the area. This gives full play to the advantages of the global optimization algorithm in the overall design space traversal and can quickly locate the design-sensitive area, while avoiding the efficiency problem of the global algorithm in detail optimization. For the global optimization algorithm, Toushmalani [28] realized the inversion of gravity anomaly parameters generated by faults by virtue of the global optimization ability of the PSO algorithm; Zhang et al. [29] realized the inversion of the mechanical parameters of gravity dams by improving the PSO algorithm and verified it through practical engineering applications; Han et al. [30] reasonably determined the rock mass mechanical parameters of the dam foundation during construction and operation through the joint use of PSO and ABAQUS; Zeng et al. [31] used the parallel PSO algorithm to globally retrieve the optimal value in the parameter space, reducing the relative error to less than 4%. PSO algorithm is one of the heuristic search algorithms [32]. Like a swarm intelligence optimization algorithm, this algorithm starts initialization through a group of random population (particles) and adjusts the fitness of the population through an iterative update to find the global optimal solution.
The updated expression of the improved PSO can be expressed as follows: n n n n n n n n 1 1 2

(10)
where + X n 1 and + V n 1 are the current position and speed of particles, and units are mm and − mm s 1 , respectively, which are the current optimal solution and the direction of the next search; ωV n is the inertia term, which makes the particles have the tendency to expand the search space, where ω is the inertia weight; ( ) − c r P X n n n 1 represents the thinking of the particle itself on the direction of improvement, where c 1 and c 2 are the acceleration factors to adjust the convergence speed of the algorithm, R n is a random number uniformly distributed from 0 to 1, and P n is the individual best position (pbset) determined by the particle from initialization to the nth iteration; ( represents the sharing of optimal information among particles. G n is the position with the best fitness value among all pbset positions in the nth iteration, that is, the global best (gbest). Among them, the location update of pbset depends on the established target function ( ) ⋅ f , namely,

(11)
At the same time, the adaptive method is used to adjust ω to improve its adaptive ability to the nonlinear and complex optimization process and further reduce the possibility of falling into the local optimal solution.
where ω i n is the weight value of the ith particle in the nth iteration, ω max and ω min are the maximum and minimum values of the inertia weight, Δx i is the distance between the ith particle and the adaptive optimal value, k is the iteration coefficient, and N is the maximum number of iterations.
For the local search algorithm, because the Hooke-Jeeves method does not need continuous objective function and line search [33] and can deal with non-continuous parameter space, the Hooke-Jeeves method is used to add a penalty term to the objective function, transform the constrained problem into an unconstrained problem for optimization, and further search the fitness interval roughly located by the PSO algorithm. The Hooke-Jeeves method mainly includes two parts, one is the detection of movement, and the other is the mode movement [34]. After the initial position and target function are determined, the detection movement is carried out along each dimension of the search space in order to determine the new base point and search direction, while the mode movement is carried out along the direction of the adjacent base point in order to make the target function meet the requirements as soon as possible. The main implementation steps are as follows: (1) Initial base point ∈ x E n 1 , n coordinate directions ⋯ e e e , , , n 1 2 , initial step length δ, acceleration factor ≥ α 1, , let's go to step 4; otherwise, go to step 3; ( , let's go to step 4; otherwise, set = + y y j j 1 , go to step 4; (4) If < j n, set = + j j 1, and go to step 2; otherwise, proceed to step 5; , proceed to step 6; otherwise, proceed to step 5;

Data interaction and inversion analysis
Although the numerical calculation and experimental results are constructed by a discrete point fitting curve, considering the calculation cost, the numerical calculation result data are relatively small, and the simulation model has no fixed length unit. Therefore, this study is based on the experimental data, smoothes the numerical calculation results with the help of the smooth function in MATLAB, and uses bilinear interpolation to supplement the missing stress values relative to the experimental data, and defines the objective function as where σ ε s i and σ ε e i , respectively, represent the simulation and experimental stress values under the same strain, and M is the number of sampling points.
At the same time, in order to reduce the amount of calculation, on the one hand, all the calculated parameters are combined and stored to avoid repeated calculation in the process of parameter optimization. On the other hand, when the simulation process is more than halfway through, the target function value is pre-calculated. Once the result is greater than 1/2 of the previous set of data, the calculation will be stopped immediately. In addition, in addition to ABAQUS, the inverse identification of mechanical parameters also involves data interaction platform and data processing analysis, which mainly includes two parts: Python's reading of the result data file odb and the submission of the calculation file inp; MATLAB processes the result data, calculates the objective function, and adjusts the parameters of the optimization algorithm. On the one hand, the optimization algorithm needs to preprocess the simulation data, including singular value elimination and data smoothing; on the other hand, it needs to search and judge the optimal value of the inversion parameters, so as to modify the inp file parameters. The specific implementation process is shown in Figure 7.
Inverse analysis is mainly divided into three steps: the first is to obtain the experimental stress-strain response curve, the second is the finite-element numerical analysis and result output, and the third is the optimal parameter search under the optimization algorithm. The specific analysis process is shown in Figure 8.
For the PSO algorithm, on the basis of comprehensive consideration of the calculation cost and accuracy, it sets the calculation parameters maximum iterations as 15, number of particles as 4, inertia as 0.9, global increment and particle increment as 0.8, maximum velocity as 0.05, and maximum failed runs as 10; for the Hooke-Jeeves algorithm, it sets the calculation parameter max evaluations to 100, relative step sizes to 0.2, and step size reduction factor to 0.5. At the same time, in order to verify the advantages of the combinatorial optimization algorithm, two modes are used for parameter inversion, one is to use only Hooke-Jeeves for local optimal search, and the other is to use the combinatorial optimization algorithm. When the Hooke-Jeeves model is used for parameter inversion alone, it is estimated to calculate 201 times, but only 101 times are actually calculated. When the combination optimization strategy is used, which is estimated that 410 times will be calculated, 61 times will be calculated in the first stage, and 41 times will be calculated in the second stage. For the three types of models, except for the numerical differences, the overall change trend of various parameters is relatively similar.

Inversion result analysis
Although the combined optimization algorithm has more computation than the single local search algorithm in theory, it can effectively reduce the computational cost and obtain relatively superior inversion parameter values through the preprocessing of partial calculation results. Under the two inversion modes, the objective function values are shown in Figure 9.  Table 2.

Mesoscopic simulation analysis
With the help of parameter inversion data, the whole process of propellant tensile damage and fracture is simulated and analyzed. The damage evolution and stress distribution of each simulation model at different stages are shown in Table 3.
Due to the difference between the micro-particles and the matrix material, the internal stress distribution is uneven, but the stress value of the micro-particles is larger than that of the bonding interface, and the interface stress value is larger than that of the matrix part. This phenomenon reflects the strengthening effect of the micro-particles on the mechanical properties of the propellant, on the other hand, it also explains the cause of damage initiation from the interface. Without considering the matrix damage, it can be clearly seen that with the progress of the micro-particle dehumidification, the stress concentration phenomenon occurs at the matrix, which causes it to undergo large deformation. Once the matrix damage is considered, it will lead to pore expansion and tensile fracture. The three types of models can reflect the damage evolution of the propellant during the tensile process to a certain extent, but the difference of stress distribution and damage state reflects the influence of the cohesive model and the setting of matrix damage parameters on the simulation process. Through the set reference points, the stress-strain curves of three models are obtained, as shown in Figure 10. By comparing Table 3 and Figure 10, it can be seen that under the action of periodic boundary load conditions, the stress of mesoscopic propellant shows nonlinear change during the tensile process. At the initial stage, the fine particles and the matrix are well bonded, showing a short linear viscoelastic stage. With the increase of tensile strain, the propellant appears to have softening phenomenon, and the slope of the curve decreases, that is, the equivalent elastic modulus begins to decline, which shows partial micro-particle dehumidification on the simulation cloud diagram. When a large number of particles are dehumidified, the stress value increased due to stretching is less than the stress value reduced due to particle dehumidification, that is, the stress reduction stage is entered. With the further intensification of the tension, the matrix reaches the strain limit, and some of the matrices appears fracture, which makes the stress drop sharply until it returns to zero.
By comparing the simulation and experimental stressstrain curves, it can be seen that the three models can reflect the stress evolution of the propellant during stretching to a certain extent in the linear viscoelastic and stress softening stages, but the linear viscoelastic stage is closer to the experimental data than the stress softening stage, and the maximum stress value is higher than the experimental results. In the stress unloading stage, the deviation of model I is larger than the actual one, and the stress will rise again with the increase of strain. Compared with model III, the  transition from stress softening to stress unloading stage of model II is very rapid, and at the end of the stress softening stage, there is a sudden increase in the slope of the curve. In the stage of damage and fracture, except for model I, model II and model III reflect the process from stress unloading to tensile fracture, but the damage and fracture trend is relatively slow compared with the experimental data. In general, model III with exponential cohesion interface and  4) ( ) δ δ n t is the interface cracking characteristic displacement; 5) λ is the damage factor; 6) LE is the matrix fracture strain.  considering matrix damage is close to the actual data, and its optimal objective function value is 0.01635, which is smaller than the 0.02136 deviation of model II and closer to the real value. Although the simulated stress-strain curve of model III is close to the experimental data, there are still three problems: the curvature of the stress softening stage is too large, the maximum stress value is too high, and the damage and fracture stage is relatively flat. Compared with the experiment and simulation analysis, the error mainly comes from the following aspects: 1) After the fabrication of the microsample, some original defects or damages will inevitably occur, which are not considered in the simulation process, resulting in the high maximum stress value. 2) In the process of tensile fracture of propellant, as a ductile material, it can be considered as a layer by layer fracture. Considering the calculation cost, the simulation model is only based on the plane strain element, which makes the stress softening stage deviate from the actual data curvature. 3) At the microlevel, the simulation calculation is based on the establishment of RVE to characterize the macromodel, which results in the high proportion of the fractured matrix in the RVE. At the macrolevel, because the proportion of cracks is relatively small and the local stress concentration is more obvious, the damage and fracture stage in the simulation calculation is more gentle.

Conclusion
(1) The organic combination of the global optimization of the adaptive particle swarm optimization algorithm and the local search of the Hooke-Jeeves algorithm, the optimal objective function values of model 1, model 2, and model 3 parameter inversion are 0.03033, 0.02136, and 0.01635, respectively. Compared with the single search algorithm, the objective function values are reduced by 10.85%, 10.67%, and 14.01%, respectively, and the actual calculation amount of the two inversion modes is basically the same, which shows the superiority and applicability of the combined optimization algorithm for the inversion of propellant parameters. (2) The optimal objective function values of the three types of RVE models decrease in turn, with model 2 reduced by 29.57% compared with model 1, and model 3 reduced by 23.46% compared with model 2, indicating that the exponential cohesive model can more accurately simulate the mechanical properties of propellants compared with bilinear, while considering the matrix damage makes the simulation data closer to the experimental data. (3) Comparing the simulation nephogram with the stress-strain curve, it can be seen that after the particle dehumidification, the micro-particles will no longer bear the load, the matrix between the particles will have a high-stress zone, and the overall stiffness of the RVE will be reduced, and the mechanical properties will be weakened. With the tensile fracture of the local matrix, the continuous superposition of pores caused by particle dehumidification will lead to the initiation and expansion of cracks, and the overall stress will change from slow growth to rapid decline. Therefore, the weakening of the mechanical properties of the propellant starts from the particle dehumidification and accelerates the matrix damage, and the two have a symbiotic effect.