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

2.1 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.

Figure 1

In situ stretching fixture and specimen.

2.2 Meso component

The density of Al powder and AP is 2.70 and 1.95 g cm⁻³, 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.

Figure 2

AP particle distribution.

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.

2.3 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 micro-morphology changes during the in situ tensile process are shown in Figure 3.

Figure 3

Micro-tensile damage and fracture of propellant.

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.

3.1 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.

Figure 4

Solid propellant RVE model.

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.

Figure 5

Schematic diagram of CZM tensile damage: (a) schematic diagram of interface CZM, (b) schematic diagram of CZM tensile damage, and (c) bilinear cohesive model of initial damage.

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 force–displacement data can be converted into the overall stress–strain data of the RVE through calculation. The schematic diagram is shown in Figure 6.

Figure 6

Schematic diagram of periodic boundary application.

3.2 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 finite-element 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:

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.

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. Multi-parameter 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:

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 1 r n ( P n − X n ) 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; c 2 R n ( G n − X n ) 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,

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 1 ∈ E n , n coordinate directions e 1 , e 2 , ⋯ , e n , initial step length δ , acceleration factor α ≥ 1 , reduction factor 0 ≤ β ≤ 1 , allowable error ε ≥ 0 , set y 1 = x 1 , k = 1, j = 1;

(2) If f ( y j + δ e j ) < f ( y j ) , set y j + 1 = y j + δ e j , let’s go to step 4; otherwise, go to step 3;

(3) If f ( y j − δ e j ) < f ( y j ) , set y j + 1 = y j − δ e j , let’s go to step 4; otherwise, set y j + 1 = y j , go to step 4;

(4) If j < n , set j = j + 1 , and go to step 2; otherwise, proceed to step 5;

(5) If f ( y n + 1 ) < f ( x k ) , proceed to step 6; otherwise, proceed to step 5;

(6) Set x k + 1 = y n + 1 , y 1 = x k + 1 + α ( x k + 1 − x k ) , k = k + 1 , j = 1 , go to step 2;

(7) If δ ≤ ε , stop the iteration and get the point x k ; otherwise, set δ = β δ , y 1 = x k , x k + 1 = x k , k = k + 1 , j = 1 , and go to step 2.

4.2 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 σ ε i s and σ ε i e , 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.

Figure 7

Data interaction process.

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.

Figure 8

Parameter inversion structure diagram.

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.

5.1 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.

Figure 9

Parameter inversion objective function value. (a) Single local search algorithm, (b) first stage of combinatorial optimization algorithm, and (c) second stage of combinatorial optimization algorithm.

Model 1 and model 2 indicate whether to consider the matrix damage and the RVE constructed by the bilinear cohesion model considering the initial damage at the bonding interface, and model 3 indicates the RVE constructed by the exponential cohesion model considering the matrix damage and the initial damage at the bonding interface. Under the single search algorithm, the minimum objective function values of model 1, model 2, and model 3 are 0.03362, 0.02364, and 0.01864, respectively. Under the combined optimization algorithm, the minimum objective function values of model 1, model 2, and model 3 are 0.03033, 0.02136, and 0.01635, respectively. The optimal values of the inversion parameters of the three models are shown in Table 2.

5.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.

Table 3

Stress nephogram of propellant damage evolution

State evolution	Model 1	Model 2	Model 3
Linear viscoelastic stage
Dehumidification of particles
Bearing substrate
Substrate damage
Tensile failure

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.

Figure 10

Microscopic simulation stress–strain diagram.

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 stress–strain 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 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.