Modeling and simulation of the debonding process of composite solid propellants

In order to study the damage evolution law of composite solid propellants, the molecular dynamics particle filled algorithm was used to establish the mesoscopic structure model of HTPB(Hydroxyl-terminated polybutadiene) propellants. The cohesive element method was employed for the adhesion interface between AP(Ammonium perchlorate) particle and HTPB matrix and the bilinear cohesive zone model was used to describe the mechanical response of the interface elements. The inversion analysis method based on Hooke-Jeeves optimization algorithm was employed to identify the parameters of cohesive zone model(CZM) of the particle/binder interface. Then, the optimized parameters were applied to the commercial finite element software ABAQUS to simulate the damage evolution process for AP particle and HTPB matrix, including the initiation, development, gathering and macroscopic crack. Finally, the stress-strain simulation curve was compared with the experiment curves. The result shows that the bilinear cohesive zone model can accurately describe the debonding and fracture process between the AP particles and HTPB matrix under the uniaxial tension loading.


Introduction
The composite solid propellant as a kind of energetic materials are widely used in the rockets and various strategic and tactical missile, its mechanical performance has a great influence on the survivability and operational capability of the missile. Many experiments have confirmed that the debonding between the particle and matrix interface is a key factor that affects the propellant macro mechanical properties [1]. Therefore, the study of the mesoscopic damage initiation, evolution, gathering and its effects on the mechanical properties of composite solid propellants is quite important and meaningful.
HTPB(Hydroxy-terminated polybutadiene) propellant is a three-phase composites, consisting of HTPB matrix, filled solid particles and the particle/binder interface. The nonlinear macroscopic mechanical response of the propellants dependents on its mesoscopic structure and the multi-scale physical process under loading. At present the meso-experiment of the propellant is difficult to carry out due to the limitations of method, so the numerical simulation method is widely applied to study the influence of mesoscopic structure on macroscopic mechanical behaviour of the HTPB propellants. Matous [2] generated the representative volume element of the propellants' mesoscopic model and used the cohesive elements to simulate the initiation and evolution of the interface debonding between the particles and the matrix. Li [3] adopted the cohesive element model to simulate propellant/linear interface and found that the failure mode experiences damage initiation, damage evolution and failure. Zhi [4][5] investigated the nonlinear mechanical behavior of solid propellants through mesomechanical methods and simulated and discussed the effects of the packing volume fraction and particle size and 2 1234567890 2nd International Conference on Design, Materials, and Manufacturing IOP Publishing IOP Conf. Series: Materials Science and Engineering 220 (2017) 012020 doi:10.1088/1757-899X/220/1/012020 distribution on the debonding process of solid propellants. Han [6] employed the inversion analysis method to obtain the optimized interface parameters and then studied the change rule of the mesoscopic interface performance of HTPB composite solid propellant with different loading rate.
Thus, this paper carried out the uniaxial tension experiment to study the mechanical properties of composite solid propellants. The mesoscopic model was generated by the molecular dynamics method and the cohesive zone model was also employed to simulate the damage evolution process of the interface. Then the Hooke-Jeeves optimization algorithm was used to identify the parameters of cohesive zone model and the optimized parameters can accurately describe the debonding process between the AP particles and HTPB matrix under the uniaxial tension loading.

Packing model and boundary
The components of solid propellants which are studied in this paper are listed in Table 1. In order to simulate the mesoscopic damage evolution process of the propellant interface, with the particle size distribution shown in Fig1, the molecular dynamics method was employed to generate the representative volume element(RVE) of the propellants.
Because of the very high Young's modulus, the AP particles' deformation is quite small, so the mesh type of AP particle must be CPE4. Meanwhile, as the incompressible hypothesis, the mesh type of the matrix should be CPE4H to model the mechanical response of HTPB with large deformation. According to the periodic boundary conditions and the mutual constraints of the RVE deformation, the boundary of the RVE need to keep straight when loading. Due to this fact, the left border and the lower boundary of the RVE are fixed and the right boundary is coupled the displacement in the x axis. The upper boundary, as the load recipient, undertakes the constant displacement load. In this study, the given speed of displacement load is 0.1667 mm/s and the corresponding strain rate is 0.003333 /s.

Cohesive zone model(CZM)
Cohesive technology is widely deployed in simulating the damage between particles and matrix which was proposed by Dugdale [7] and Barenblatt [8].
The interface elements open when damage occurs and lose their stiffness at failure so that the continuum elements are disconnected, and in this paper, a damage model called cohesive zone model that consists of two ingredients a damage initiation criterion and a damage evolution law is accepted. So far, several types of CZM models have been proposed, they can be distinguished by their own cohesive low [9]. The typical bilinear cohesive zone model is given in Fig 3,  traction, t ,to displacement jump,  ,at the interface where a debonding may occur. The interface begins to damage when the maximum contact stress is equal to the maximum stress, and complete failure occurs at the maximum displacement f  . The slope of the traction-displacement curve before damage initial, K , is refereed as the interface stiffness. K , 0 t and f  are the basic parameters to confirm a bilinear CZM model.
where the symbol is the Macaulay bracket, the subscript "n" and "s" represents the normal and the shear direction.
The traction-separation law can be expressed by: where, D is the damage factor, can be expressed by:

Material properties
The components of solid propellants are listed in Table 1 Some researchers considered the matrix of solid propellants as a kind of elastic material. This simplification is inaccurate. In fact, the material usually shows a viscoelastic behavior just like rubber. In this paper, the universal stress relaxation texts on HTPB solid propellants matrix samples are conducted to reveal the viscoelastic response. The applied load is 100% strain for 1200 s. it is assumed that the matrix is a linear viscoelastic material which is defined by a Prony series [10]: where the E  ,

Mechanical parameters
The interface meso-mechanical parameter values have a vital influence on the accuracy of the numerical calculation results. At present, the approximate cohesive parameters are usually estimate according to the experience and literature and it is often inaccurate. In order to solve the disadvantages, many researchers raised the inversion optimization method based on Hooke-Jeeves algorithm to acquire the more exact parameters. It's necessary to conduct the uniaxial tensile test of composite solid propellants for the inversion analysis of the interface parameters. Taking no account of the temperature and strain rate dependency, the applied load is 10 mm/min and the corresponding strain rate is 0.003333 /s at room temperature. The engineering stress-train curve of the propellant is given in Fig.5 There are three basic parameters in a typical bilinear cohesive zone model and they have different influence on the stress-stain curve's tendency respectively. In view of this feature, the step-by-step inversion method was determined and the objective function that characterize the coincidence degree of the test and the simulation curve was build.
When the objective function value is less than the selected tolerance lim R , the current interface parameters are believed to be accurate. The initial and the ultimate value of the cohesive interface parameters are listed in Table 2.

Simulation result
The mesoscopic model generated in the previous section was calculated by the finite element software ABAQUS to simulate the uniaxial tension of HTPB propellants. The Mises stress nephograms of the propellant at 5%, 15%, and 50% strain respectively are given in Fig. 6.
2nd International Conference on Design, Materials, and Manufacturing IOP Publishing IOP Conf. Series: Materials Science and Engineering 220 (2017) 012020 doi:10.1088/1757-899X/220/1/012020 Figure 6. Mises stress distribution in the propellant when at 5%、15% and 50% strain respectively. Figure 6 shows that when the strain is 5%, the stress and strain distribution in solid propellant is uneven, it is caused by the inconsistent deformation among the particles, matrix and interface. As an elastomer, the modulus of particles is quite larger than the matrix, so the deformation of the matrix is much bigger than the particles and this certainly makes the particle/matrix interface the weakness in the propellants. With the load increases to 15% strain, the local strain near the large particles region is very huge and this leads to a very obvious stress concentrain, the stress in the larger particle/matrix interface reaches the cohesive strength fist and as a result, these interfaces begin to debond. When the strain reaches 50%, more and more interfaces have debonded and the micro holes develop, gather and the macro fracture form, eventually, the composite solid propellants break.

Results analysis
In simulation results, the pulling force on the upper boundary of mesoscopic model can be determined by summating all the constraint reaction of each node, and then the emulational stress-strain curve of the solid propellants can be obtained through corresponding calculation. The emulational stress-strain curve and the mean experiment curve that acquired in the previous section were given in Fig. 7.   Figure 7. Engineering stress-strain curve of HTPB propellant. Combining the stress nephograms in Fig.6 and the engineering stress-strain curve in Fig.7, it's found that in stage Ⅰ( 5%   ), the particles and matrix bond in a good condition and there is no damage inside the mesoscopic structure of the representative volume element, so the propellants show the elastic properties and the stress grows linearly with the strain increasing. Then coming to the stage Ⅱ, in this stage the large particles and matrix interfaces begin to debond and it leads to a decrease of the load bearing capability, this phenomenon can also be observed in Fig 7 that the stress-strain curve begins to appear the nonlinearity and enter a "platform area". Finally in the stage Ⅲ， more and more interfaces undertake the serious debonding and the micro holes appear, then the holes develop and gather resulting in macro crack and the stress of the propellants drops rapidly.
Comparing the simulation curve and experiment curve in Fig.7, it's clear that the simulation curve obtained by the optimized inversion method matches the actual experiment curve well, but there also exist a little gap between them. In can be thought that the main causes of the difference are: (1) the mesoscopic model established in this paper is two-dimensional and there is a certain deviation with the three-dimensional state of the actual solid propellants; (2) there are many kinds of initial defects such as micro-hole and micro-crack in the propellants and it's ignored in the mesoscopic model; (3) all the solid particles are considered as regular circular, however there may be irregular particles in the propellants and the stress concentration and interface debonding are more likely to occur near these particles.
From the above analysis, we can see that the interface debonding process is an important impact factor in the propellants' nonlinear property. The mesoscopic model and the inversion optimization interface parameters are reliable and can be used to simulate the debonding process and the mechanical response of the solid propellants.

Conclusion
(1) Based on the periodicity assumption, the molecular dynamics method was employed to establish the representative volume element of the composite solid propellants mesoscopic filling model. Meanwhile, the cohesive element method was employed to simulate the mesoscopic damage initiation, evolution, gathering of the particle/matrix interface and its effects on the propellants' mechanical properties. (2) In order to obtain the exact parameters of cohesive zone model which can characterize the mechanical properties of the particle/matrix adhesive interface, the uniaxial tensile test of composite solid propellants was performed. Combining the parameter optimization inversion method based on "Hooke-Jeeves" algorithm with the experimental data, the relatively accurate parameters are acquired. (3) The interface debonding process between particles and matrix in the solid propellants was simulated by means of the finite element platform ABAQUS. It's found that while the interface stress reaches the cohesive strength, the debonding occurs and with the increase of load, the debonding damage develops, gathers and macro fracture eventually forms. It was proved that the particle/matrix interface debonding is an important cause to the nonlinearity of the propellant stress-strain curve. Besides, the mesoscopic model and simulation method taken in this paper can be used to accurately simulate the macroscopic mechanical behavior of the composite solid propellants.