Modelling Microstructural Deformation and the Failure Process of Plastic Bonded Explosives Using the Cohesive Zone Model

A microstructure finite element method combining the cohesive zone model (CZM) is used to simulate the mechanical behavior, deformation, and failure of polymer-bonded explosive (PBX) 9501 under quasi-static loading. PBX 9501 consists of Cyclotetramethylene tetranitramine (HMX) filler particles with a random distribution packaged in a polymeric binder. The particle is treated as elastic and the binder as viscoelastic. Cohesive elements with a bilinear softening law are inserted into the particle/binder interface, the HMX particle, and the binder to study the interface’s debonding and failure evolution. Macroscopic stress–strain curves homogenized across the microstructure under tension and compression with different strain rates are basically consistent with the experimental data. The interface debonding approximately vertical to the loading direction is the primary failure mechanism under tension, while shear failure along the interfaces and particle fracture plays a significant role under compression. The effects of interface strengths and strain rates on the performance of PBX 9501 are also evaluated. The tensile and compressive strengths are dependent on the interface strength and strain rate, but the failure paths are insensitive. This model is shown to accurately predict macroscopic responses and improve our understanding of the relationship between the mechanical behavior and microstructure of PBX 9501.


Introduction
Polymer-bonded explosives (PBXs) are a type of highly filled particulate composite, generally composed of high-energy single-compound explosive particles and a polymeric binder [1,2]. PBXs have been widely used in military and civil products, such as missiles, rockets, and anti-tank bombs. The little insert polymer binder is used to reduce sensitivity and to make the material safer to handle. However, the binder can strongly affect the mechanical and fracture properties of PBXs. Debonding of the particle/binder interface, and even failure of the microstructure can arise as PBXs may be subjected to tension and compression loading. The failure caused by different stimuli influences the mechanical properties and even detonation behavior of PBXs. Understanding the microstructural mechanical properties and failure mechanisms of PBXs are important for their design and safe application.
Many experimental studies of the micromechanical properties and failure mechanisms of PBXs have been conducted in recent years. The Brazilian test with real-time microscopic observations was utilized to investigate the material's microstructural deformation evolution and failure mechanisms under quasi-static tension [3][4][5][6]. Different failure modes, including interface debonding, binder rupture,

HMX Particles
It is generally believed that HMX is brittle at ambient pressure, so it undergoes very little plastic deformation [26]. Hence, we made the simplifying assumption that HMX can be regarded as a linearly isotropic elastic material with Young's modulus of 13.375 Mpa and a Poisson's ratio of 0.32 in this study.

Binder Material
Estane 5703, considered as a polymer binder, is used in PBX 9501. A series of experiments were carried out on Estane 5703 to measure its mechanical properties at strain rates of 10 −4 s −1 to 10 3 s −1 [27]. The material shows viscoelastic properties with a dependence on the strain rate. A 22 element Prony series is used to describe the variation of the shear modulus G with the relaxation time. The Prony series formulation is where G 0 is the instantaneous shear modulus at reference temperature T 0 , G ∞ is the steady-state shear modulus, g i = G i /G 0 is the relative modulus of the ith term, N p is the number of terms in the Prony series, and τ p i is the relaxation times. The terms for the Prony series defining G(t) are listed in Table 1. The stress-strain relationship of the generalized viscoelastic body can be expressed by an integral form where σ is the Cauchy stress, ε D and ε H are the deviatoric and hydrostatic portions of the Eulerian strain tensor, and t and τ are the physical and reduced times, respectively. The bulk modulus K is a constant, as in [24,27].

CZM with a Bilinear Constitutive Law
A cohesive zone model, assuming an imaginary zone in front of the real crack tip, is used to predict the crack in the PBX samples. As shown in Figure 1, the tip of this zone is the boundary point between the intact material and the damaged material. In this zone, the traction linearly increases from zero at the rear crack tip to its maximum at the virtual crack tip. A typical bilinear traction-separation law (as illustrated in Figure 2) is used to characterize the nonlinear relationship between the traction A cohesive zone model, assuming an imaginary zone in front of the real crack tip, is used to predict the crack in the PBX samples. As shown in Figure 1, the tip of this zone is the boundary point between the intact material and the damaged material. In this zone, the traction linearly increases from zero at the rear crack tip to its maximum at the virtual crack tip. A typical bilinear tractionseparation law (as illustrated in Figure 2) is used to characterize the nonlinear relationship between the traction and separation in the cohesive model. Figure 2a shows the curve of the traction's normal component n T with the crack's opening displacement, while Figure 2b shows the curve of the traction's tangential component t T with the crack's sliding displacement. The traction first increases linearly, with separation displacement on the undamaged cohesive crack's surface, and then decreases monotonically after damage and finally reduces to zero after the entire separation. The n T and t T can be expressed as A cohesive zone model, assuming an imaginary zone in front of the real crack tip, is used to predict the crack in the PBX samples. As shown in Figure 1, the tip of this zone is the boundary point between the intact material and the damaged material. In this zone, the traction linearly increases from zero at the rear crack tip to its maximum at the virtual crack tip. A typical bilinear tractionseparation law (as illustrated in Figure 2) is used to characterize the nonlinear relationship between the traction and separation in the cohesive model. Figure 2a shows the curve of the traction's normal component n T with the crack's opening displacement, while Figure 2b shows the curve of the traction's tangential component t T with the crack's sliding displacement. The traction first increases linearly, with separation displacement on the undamaged cohesive crack's surface, and then decreases monotonically after damage and finally reduces to zero after the entire separation. The n T and t T can be expressed as The T n and T t can be expressed as When δ n > 0, When δ n = 0, where n and t are unit vectors along the normal and tangential directions of the cohesive crack surface, and δ is the cohesive crack surface separation displacement vector. δ n = n · δ and δ t = t · δ are the normal and tangential separation displacement components. T nc and T tc , respectively, represent the maximum normal and tangential cohesive strength. δ 0 n and δ 0 t , respectively, represent the normal and tangential separation displacement corresponding to T nc and T tc . δ f n and δ f t are the normal and tangential critical relative separation displacement, respectively. k n = T nc δ 0 n and k t = T tc δ 0 t are the normal and tangential initial stiffness coefficient, respectively. The cohesive zone model involves three parameters: the cohesive strength T, the separation displacement δ, and the initial stiffness coefficient k. The separation displacement δ is also determined by the cohesive fracture energy G, which is equal to the area under the traction-separation curves (see The damage, activated in terms of a maximum stress criterion, can be expressed as where is the Macaulay bracket, defined as The cohesive element is represented by the inserted particles, binders, and particle/binder interface along the element boundaries. The normal and tangential parameters need to be determined for the three types of cohesive elements. The calculation of critical separation displacement caused by different normal and tangential initial stiffness is very complicated. Li et al. proposed an assumption so-called "equal ratio attenuation" to simplify the constitutive model and facilitate the determination of the parameters [28,29]. As shown in Figure 3, as the tangential and the normal displacement are the same, both stiffnesses always maintain the same ratio, therefore the same value of separation displacement is assumed for both the normal direction (mode I) and the tangential direction (mode II).
The cohesive element is represented by the inserted particles, binders, and particle/binder interface along the element boundaries. The normal and tangential parameters need to be determined for the three types of cohesive elements. The calculation of critical separation displacement caused by different normal and tangential initial stiffness is very complicated. Li et al. proposed an assumption so-called "equal ratio attenuation" to simplify the constitutive model and facilitate the determination of the parameters [28,29]. As shown in Figure 3, as the tangential and the normal displacement are the same, both stiffnesses always maintain the same ratio, therefore the same value of separation displacement is assumed for both the normal direction (mode I) and the tangential direction (mode II).  According to Wu [24], the tensile strength of 4.754 MPa is taken as the maximum normal cohesive strength T nc of the cohesive elements in the HMX particles, while the corresponding separation displacement δ 0 n is 2.38 µm, and the normal initial stiffness of k n = 2 GPa/µm. The compression strength σ c of 260 MPa is given for HMX in the literature. The maximum tangential cohesive strength, considered as the shear strength, is 131 MPa, which can be calculated as T tc = σ c +0.3T nc 2 , and then the tangential initial stiffness is k t = 55 GPa/µm. The same critical separation displacement proposed in the literature is used for both the normal direction (mode I) and tangential direction (mode II) (i.e., δ f n = δ f t = 5 µm). Figure 4 shows the calculated and measured stress-strain curves for HMX crystal at uniaxial tensile strain rate of 10 −3 s −1 . It is possible to accurately characterize the HMX crystal behavior with the parameters from literature.  The literature values for the cohesive elements in the binder are used [22,23]. A cohesive strength of 5 MPa and an initial stiffness of 1 GPa/μm are given for the normal direction (mode I), while a cohesive strength of 38.4 MPa and initial stiffness of 7.68 GPa/μm are given for the tangential direction (mode II). The normal and tangential critical separation displacements have the same values of  The literature values for the cohesive elements in the binder are used [22,23]. A cohesive strength of 5 MPa and an initial stiffness of 1 GPa/µm are given for the normal direction (mode I), while a cohesive strength of 38.4 MPa and initial stiffness of 7.68 GPa/µm are given for the tangential direction (mode II). The normal and tangential critical separation displacements have the same values of δ f n = δ f t = 5 mm. Figure 5 presents the calculated and measured compressive stress-strain curves of the binder material at strain rate of 10 −3 s −1 . From the figure, the simulation result is slightly higher than the experimental value. More experimental data of the Estane 5703 under quasi-static loading at the different strain rates are needed to validate the rate sensitivity.
The values for the cohesive elements in the particle-binder interface are taken from the work of Tan et al. [31], who used an experimental method to study a cohesive law for the particle/binder interface in PBX 9501. The normal cohesive strength is T nc = 1.66 MPa and the normal initial stiffness is k n = 1.55 GPa/µm. The normal and tangential separation displacements reach the critical value of δ f n = δ f t = 0.11 mm. Dienes and Kershneer measured the energy consumption per unit area to be about 195 J/m 2 during the binder tearing process [32]. The experimental results show that the mode II crack on the interface is always accompanied by a tearing of the binder. Herein, we take the tangential cohesive fracture energy G tc = 195 J/m 2 , so the calculated value of T tc = 3.55 MPa, and the tangential initial stiffness of k t is taken as 3.3 GPa/µm. Figure 6 shows the normal traction-separation curves of the particle/binder interface at different loading rates [24]. The cohesive strength and fracture energy increased with an increase in the loading rate. Therefore, we should consider incorporating the rate-dependence of the interface into the model. Table 2 summarizes the parameters used for the three types of cohesive elements. Further calibration for the CZ parameters will be performed using experimental data for PBX 9501 in Section 3.3.
The literature values for the cohesive elements in the binder are used [22,23]. A cohesive strength of 5 MPa and an initial stiffness of 1 GPa/μm are given for the normal direction (mode I), while a cohesive strength of 38.4 MPa and initial stiffness of 7.68 GPa/μm are given for the tangential direction (mode II). The normal and tangential critical separation displacements have the same values of 5mm f f n t δ δ = = . Figure 5 presents the calculated and measured compressive stress-strain curves of the binder material at strain rate of 10 −3 s −1 . From the figure, the simulation result is slightly higher than the experimental value. More experimental data of the Estane 5703 under quasi-static loading at the different strain rates are needed to validate the rate sensitivity. The values for the cohesive elements in the particle-binder interface are taken from the work of Tan et al. [31], who used an experimental method to study a cohesive law for the particle/binder interface in PBX 9501. The normal cohesive strength is 1.66MPa nc T = and the normal initial stiffness is The normal and tangential separation displacements reach the critical value of Dienes and Kershneer measured the energy consumption per unit area to be about 195 J/m 2 during the binder tearing process [32]. The experimental results show that the mode II crack on the interface is always accompanied by a tearing of the binder. Herein, we take the tangential cohesive fracture energy , and the tangential initial stiffness of t k is taken as 3.3 GPa/μm. Figure 6 shows the normal traction-separation curves of the particle/binder interface at different loading rates [24]. The cohesive strength and fracture energy increased with an increase in the loading rate. Therefore, we should consider incorporating the rate-dependence of the interface into the model. Table 2 summarizes the parameters used for the three types of cohesive elements. Further calibration for the CZ parameters will be performed using experimental data for PBX 9501 in Section 3.3.
As mentioned in the literature, initially inserted bilinear CZ models reduce cohesive-surfaceinduced stiffness when a finite initial stiffness is used in the cohesive law [33][34][35]. A sufficiently large initial stiffness (k = 1-55 GPa/μm) and a finite element size of 8 μm are used to solve this problem in our work.   As mentioned in the literature, initially inserted bilinear CZ models reduce cohesive-surfaceinduced stiffness when a finite initial stiffness is used in the cohesive law [33][34][35]. A sufficiently large initial stiffness (k = 1-55 GPa/µm) and a finite element size of 8 µm are used to solve this problem in our work.

Microstructure Model of PBX
The PBX 9501 microstructure selected here consists of an HMX particle (approximately 95% by mass) and a typical binder material-Estane 5703. Voronoi tessellation is applied to generate random distributions for the HMX particle polygons of various shapes and sizes, as shown in Figure 7. The PBX microstructure is a 1 mm × 0.8 mm 2D model consisting of 98 HMX particles with a size distribution from 25 µm to 98 µm and an Estane 5703 binder with a width of 8 µm. The densities of HMX and Estane 5703 are 1.58 g/cm 3 and 0.9 g/cm 3 , respectively. The model contains 91% HMX by volume fraction, which is equivalent to a 94.7% mass fraction. On account of the simulation results of the model with 98 and 200 particles, we verified that 98 particles can provide reasonable simulation results, so we did not consider using more particles.

Cohesive element type
Normal Tangential

Microstructure Model of PBX
The PBX 9501 microstructure selected here consists of an HMX particle (approximately 95% by mass) and a typical binder material-Estane 5703. Voronoi tessellation is applied to generate random distributions for the HMX particle polygons of various shapes and sizes, as shown in Figure 7. The PBX microstructure is a 1 mm × 0.8 mm 2D model consisting of 98 HMX particles with a size distribution from 25 μm to 98 μm and an Estane 5703 binder with a width of 8 μm. The densities of HMX and Estane 5703 are 1.58 g/cm 3 and 0.9 g/cm 3 , respectively. The model contains 91% HMX by volume fraction, which is equivalent to a 94.7% mass fraction. On account of the simulation results of the model with 98 and 200 particles, we verified that 98 particles can provide reasonable simulation results, so we did not consider using more particles. The obtained microstructure model is introduced into ABAQUS/Explicit finite element code (ABAQUS, 2016) to partition the CPS3 elements (the triangular plane stress elements) with a mesh size of 8 μm [36]. As mentioned in the literature [17,22,35], cohesive elements are embedded along all finite element boundaries as part of the physical model to predict fracture patterns and the fracture outcome in our analysis. We referred to the insertion procedure proposed by Yang et al. By using our program, COH2D4 elements (a 4-node cohesive element with zero in-plane thickness) are embedded into the model along the CPS3 element boundaries [37]. A node shared by six triangular elements is used as an example to illustrate the embedding process, as shown in Figure 8a. This node is first copied and reset by six separate nodes at the same location (see Figure 8b), and then the reset nodes are connected along the triangular element sides to generate six zero-thickness cohesive elements. The obtained microstructure model is introduced into ABAQUS/Explicit finite element code (ABAQUS, 2016) to partition the CPS3 elements (the triangular plane stress elements) with a mesh size of 8 µm [36]. As mentioned in the literature [17,22,35], cohesive elements are embedded along all finite element boundaries as part of the physical model to predict fracture patterns and the fracture outcome in our analysis. We referred to the insertion procedure proposed by Yang et al. By using our program, COH2D4 elements (a 4-node cohesive element with zero in-plane thickness) are embedded into the model along the CPS3 element boundaries [37]. A node shared by six triangular elements is used as an example to illustrate the embedding process, as shown in Figure 8a. This node is first copied and reset by six separate nodes at the same location (see Figure 8b), and then the reset nodes are connected along the triangular element sides to generate six zero-thickness cohesive elements. This model required a total of 118,377 elements. The three types of cohesive elements are highlighted in red, as shown in Figure 9.

Loading and Boundary Conditions
The macroscopic response is strongly affected by the boundary conditions. Three different cases, including (a) traction free vertical edges,

Loading and Boundary Conditions
The macroscopic response is strongly affected by the boundary conditions. Three different cases, including (a) traction free vertical edges,

Loading and Boundary Conditions
The macroscopic response is strongly affected by the boundary conditions. Three different cases, including (a) traction free vertical edges, (b) straight vertical edges, u(0, y) = u(0, L y ), u(L x , y) = u(L x , L y ), and (c) periodic vertical edges u(L x , y) = u(0, y) + u(L x , 0), where u is displacement and L is length of model, were investigated to simulate the tensile and compressive loading of the model. The loading and boundary conditions are shown in Figure 10. The vertical displacement loading control is applied at the top surface to simulate a uniaxial tension and compression condition. The bottom surface is constrained in the vertical direction. The results presented in Figures 11 and 12 show that the choice between cases (a) and (b) of boundary conditions does not affect the outcome of a given simulation. Case (c) shows some difference due to periodic boundary conditions are generally applied to repeatable units of materials, and the geometry in our study although representative of the material, is random and not repeatable. The same calculated results have been also proposed by Arora et al. [20]. Therefore, all further simulations are performed using case (b) boundary condition. The results presented in Figures 11 and 12 show that the choice between cases (a) and (b) of boundary conditions does not affect the outcome of a given simulation. Case (c) shows some difference due to periodic boundary conditions are generally applied to repeatable units of materials, and the geometry in our study although representative of the material, is random and not repeatable. The same calculated results have been also proposed by Arora et al. [20]. Therefore, all further simulations are performed using case (b) boundary condition.

Model Calibration
The simulations were performed using ABAQUS/Explicit. A massing scale ensured that the kinetic energy was negligible compared to the internal energy ensuring the simulation was quasistatic. As shown in Figure 13, the ratio of kinetic energy to internal energy is less than 5% during the loading process. The effect of element size is analyzed by three different element sizes with 5 μm, 8 μm and 10 μm. The calculated peak stress and failure strain at a tensile strain rate of 10 −3 s −1 and compressive strain rate of 10 −2 s −1 are listed in Table 3. The relative error in the peak stress and failure strain is less than 5.13% of the value for 5 μm element size. Considering the balance between efficiency and accuracy, the element size of 8 μm is chosen for further simulations.

Model Calibration
The simulations were performed using ABAQUS/Explicit. A massing scale ensured that the kinetic energy was negligible compared to the internal energy ensuring the simulation was quasistatic. As shown in Figure 13, the ratio of kinetic energy to internal energy is less than 5% during the loading process. The effect of element size is analyzed by three different element sizes with 5 μm, 8 μm and 10 μm. The calculated peak stress and failure strain at a tensile strain rate of 10 −3 s −1 and compressive strain rate of 10 −2 s −1 are listed in Table 3. The relative error in the peak stress and failure strain is less than 5.13% of the value for 5 μm element size. Considering the balance between efficiency and accuracy, the element size of 8 μm is chosen for further simulations. The average stress σ y is the average value of the stress of all the elements except the cohesive elements in the loading direction, σ y = ( i σ yi S i )/S. The strain is calculated by the displacement u y on the top surface, ε y = u y /L y , L y is the length of model at Y direction.

Model Calibration
The simulations were performed using ABAQUS/Explicit. A massing scale ensured that the kinetic energy was negligible compared to the internal energy ensuring the simulation was quasi-static. As shown in Figure 13, the ratio of kinetic energy to internal energy is less than 5% during the loading process. The effect of element size is analyzed by three different element sizes with 5 µm, 8 µm and 10 µm. The calculated peak stress and failure strain at a tensile strain rate of 10 −3 s −1 and compressive strain rate of 10 −2 s −1 are listed in Table 3. The relative error in the peak stress and failure strain is less than 5.13% of the value for 5 µm element size. Considering the balance between efficiency and accuracy, the element size of 8 µm is chosen for further simulations.  Figure 13. The ratio of kinetic energy to internal energy.  [24,38]. Since the data obtained from the literature only include the stress-strain curves before the damage of PBX 9501, we compared the curves of the strain from 0 to the failure strain. Figure 14 shows a comparison of the measured and simulated stress-strain curves of PBX 9501 at three strain rates of 10 −4 s −1 , 10 −3 s −1 , and 10 −2 s −1 . The measured and calculated responses under uniaxial tension (in Figure  14a) and compression (Figure 14b [24,38]. Since the data obtained from the literature only include the stress-strain curves before the damage of PBX 9501, we compared the curves of the strain from 0 to the failure strain. Figure 14 shows a comparison of the measured and simulated stress-strain curves of PBX 9501 at three strain rates of 10 −4 s −1 , 10 −3 s −1 , and 10 −2 s −1 . The measured and calculated responses under uniaxial tension (in Figure 14a) and compression (Figure 14b Figure 15 shows the macroscopic stress-strain curve calculated at a strain rate of 10 −3 s −1 . The stress increases monotonously and nonlinearly up to the peak stress, and then sharply drops due to the failure from the interface debonding. We noted that the failure strain in uniaxial quasi-static tension is only about 0.0041, revealing the brittle, fractural nature of PBX 9501.  Figure 15 shows the macroscopic stress-strain curve calculated at a strain rate of 10 −3 s −1 . The stress increases monotonously and nonlinearly up to the peak stress, and then sharply drops due to the failure from the interface debonding. We noted that the failure strain in uniaxial quasi-static tension is only about 0.0041, revealing the brittle, fractural nature of PBX 9501. The true strain contours of y-direction y ε and the failure evolution at the macroscopic strains of 0.16%, 0.41%, 0.45%, and 0.5% (marked by dash line in Figure 15) are shown in Figure 16. The soft binder accommodates a great amount of the deformation throughout the process. The strain distribution is heterogeneous and results from the microstructure's features (i.e., the random particle size shapes and interface debonding behavior). Localized strain concentrations labeled in red can be found in some regions around the particle sides (Figure 16a), yet no microcracks are formed at a macroscopic strain of 0.16%. At a macroscopic strain of 0.41% (shown in Figure 16b), a microcrack, as indicated by the arrows, is observed along the interface with a weaker strength as a result of interface debonding, thereby lowering the material's load-carry capacity and contributing to material The true strain contours of y-direction ε y and the failure evolution at the macroscopic strains of 0.16%, 0.41%, 0.45%, and 0.5% (marked by dash line in Figure 15) are shown in Figure 16. The soft binder accommodates a great amount of the deformation throughout the process. The strain distribution is heterogeneous and results from the microstructure's features (i.e., the random particle size shapes and interface debonding behavior). Localized strain concentrations labeled in red can be found in some regions around the particle sides (Figure 16a), yet no microcracks are formed at a macroscopic strain of 0.16%. At a macroscopic strain of 0.41% (shown in Figure 16b), a microcrack, as indicated by the arrows, is observed along the interface with a weaker strength as a result of interface debonding, thereby lowering the material's load-carry capacity and contributing to material softening. The form of the microcrack reduces the effective mechanical properties of the microstructure [39]. With the increasing of loading, greater interface debonding makes the microcrack grow continuously (indicated by Figure 16c). It can be seen that the strain inside the binder on both sides of the crack reduces and redistributes itself. At a macroscopic strain of 0.5%, a single main macrocrack is formed roughly vertical to the loading direction (Figure 16d). However, it is interesting to note that the local maximum strain inside the binder on the crack path reaches 0.588. Residual binders, as marked by the arrows, allow the structure to continue to withstand certain stresses. This also explains why the macrostress of the stress-strain curve in Figure 9 is maintained at about 0.1 MPa. From Figure 16, we can see that the failure mechanism of PBX 9501 in uniaxial quasi-static tension is predominantly interface debonding, which is similar to the phenomenon monitored in the real microscopic examination by Rae et al. and Chen et al. [5,8].

Effect of the Interface Strength
The interface cohesive strength affecting the mechanical properties of PBXs has been reported in some of the literature. Three interface cohesive strengths of 1, 1.66, and 3 MPa are used for the model illustrated in Figure 10b to investigate the mechanical behavior of PBX 9501 under uniaxial tension. The calculated stress-strain curves for the different interface cohesive strength models are shown in Figure 17. As the interface cohesive strength increases from 1 to 3 MPa, the tensile strength increases from 0.78 to 2.41 MPa, and the failure strain increases from 0.0035 to 0.0086. It is remarkable that the predicted crack path is almost identical, as shown in Figure 18a-c. This result suggests that PBX 9501′s tensile mechanical properties always are related to its interface cohesive strength, while the failure path is independent.

Effect of the Interface Strength
The interface cohesive strength affecting the mechanical properties of PBXs has been reported in some of the literature. Three interface cohesive strengths of 1, 1.66, and 3 MPa are used for the model illustrated in Figure 10b to investigate the mechanical behavior of PBX 9501 under uniaxial tension. The calculated stress-strain curves for the different interface cohesive strength models are shown in Figure 17. As the interface cohesive strength increases from 1 to 3 MPa, the tensile strength increases from 0.78 to 2.41 MPa, and the failure strain increases from 0.0035 to 0.0086. It is remarkable that the predicted crack path is almost identical, as shown in Figure 18a-c. This result suggests that PBX 9501's tensile mechanical properties always are related to its interface cohesive strength, while the failure path is independent.
The interface cohesive strength affecting the mechanical properties of PBXs has been reported in some of the literature. Three interface cohesive strengths of 1, 1.66, and 3 MPa are used for the model illustrated in Figure 10b to investigate the mechanical behavior of PBX 9501 under uniaxial tension. The calculated stress-strain curves for the different interface cohesive strength models are shown in Figure 17. As the interface cohesive strength increases from 1 to 3 MPa, the tensile strength increases from 0.78 to 2.41 MPa, and the failure strain increases from 0.0035 to 0.0086. It is remarkable that the predicted crack path is almost identical, as shown in Figure 18a-c. This result suggests that PBX 9501′s tensile mechanical properties always are related to its interface cohesive strength, while the failure path is independent.

Mechanical Response and Failure Mechanism
The mechanical behavior of PBX 9501 under uniaxial compression is obviously different from that under uniaxial tension. Figure 19 plots the stress-strain response of PBX 9501 at a compressive strain rate of 10 −3 s −1 . The similar response features with those under uniaxial tension (such as initial rapid growth, peak stress, and softening behaviors) can also be observed. While some notable differences still exist, as shown in Figure 20, the softening rate of PBX 9501 under uniaxial compression is lower than that under uniaxial tension. The peak strength and failure strain are approximately 9.36 MPa and 0.0118, respectively, which are higher than the 1.16 MPa and 0.0041 under uniaxial tension. This reveals that the compressive properties of PBX 9501 are greater than its tensile properties.

Mechanical Response and Failure Mechanism
The mechanical behavior of PBX 9501 under uniaxial compression is obviously different from that under uniaxial tension. Figure 19 plots the stress-strain response of PBX 9501 at a compressive strain rate of 10 −3 s −1 . The similar response features with those under uniaxial tension (such as initial rapid growth, peak stress, and softening behaviors) can also be observed. While some notable differences still exist, as shown in Figure 20, the softening rate of PBX 9501 under uniaxial compression is lower than that under uniaxial tension. The peak strength and failure strain are approximately 9.36 MPa and 0.0118, respectively, which are higher than the 1.16 MPa and 0.0041 under uniaxial tension. This reveals that the compressive properties of PBX 9501 are greater than its tensile properties.
compression is lower than that under uniaxial tension. The peak strength and failure strain are approximately 9.36 MPa and 0.0118, respectively, which are higher than the 1.16 MPa and 0.0041 under uniaxial tension. This reveals that the compressive properties of PBX 9501 are greater than its tensile properties.     Figure 21 shows the failure evolution of PBX 9501 at a macroscopic strain of 0.96%, 1.2%, 1.77%, and 2.88%, for a 10 −3 s −1 strain rate. At the early stage of loading, the strain increases from 0% to 0.96%, and the soft binder sustains primarily deformation. Debonding occurs along the relatively weak particle/binder interface due to its shearing slide, as shown in Figure 21a. Thus, the tangential stress of the interface first reaches its critical value under compression loading. With an increase in loading, the compression strength reaches its maximum, and more interface microcracks and transgranular fractures of the particle occur (see Figure 21b,c) due to the interface slide and complicated interactions between the particles. Indeed, this interaction affects the failure mechanism of PBX 9501. Next, at a macroscopic strain of 2.88%, the macroscopic failure paths labeled by the red line spread over the whole model, along with the growth and connection of numerous microcracks, so that relatively shear failure modes are exhibited. In Figure 21, we can see that interface debonding and particle fractures (marked by arrows) are the primary modes of failure under uniaxial compression. These results were also observed in an experiment by Gray et al. and Zhou et al. [40,41].
weak particle/binder interface due to its shearing slide, as shown in Figure 21a. Thus, the tangential stress of the interface first reaches its critical value under compression loading. With an increase in loading, the compression strength reaches its maximum, and more interface microcracks and transgranular fractures of the particle occur (see Figure 21b,c) due to the interface slide and complicated interactions between the particles. Indeed, this interaction affects the failure mechanism of PBX 9501. Next, at a macroscopic strain of 2.88%, the macroscopic failure paths labeled by the red line spread over the whole model, along with the growth and connection of numerous microcracks, so that relatively shear failure modes are exhibited. In Figure 21, we can see that interface debonding and particle fractures (marked by arrows) are the primary modes of failure under uniaxial compression. These results were also observed in an experiment by Gray et al. and Zhou et al. [40,41].

Effect of Particle/Binder Interface Strength
The effect of the interface strength on the mechanical response of PBX 9501 at a compressive strain rate of 10 −2 s −1 was also investigated. Four different tangential strengths of the particle/binder interface, including 2 MPa, 4.2 MPa, 5 MPa, and 7 MPa, were selected for the model in Figure 10b. Figure 22 plots the compressive stress-strain curves for a strain rate of 10 −2 s −1 at different interface strengths. As the interface strength becomes stronger, the compressive strength becomes higher. Thus, the failure strain is relatively independent of the interface strength.

Effect of Particle/Binder Interface Strength
The effect of the interface strength on the mechanical response of PBX 9501 at a compressive strain rate of 10 −2 s −1 was also investigated. Four different tangential strengths of the particle/binder interface, including 2 MPa, 4.2 MPa, 5 MPa, and 7 MPa, were selected for the model in Figure 10b. Figure 22 plots the compressive stress-strain curves for a strain rate of 10 −2 s −1 at different interface strengths. As the interface strength becomes stronger, the compressive strength becomes higher. Thus, the failure strain is relatively independent of the interface strength. Figure 23 presents the failure paths for different interface strengths under compression loading. As seen in Figure 23a-d, a large number of microcracks along the interface reveals that interface debonding is the primary failure mode for the weaker interface strength. With an increase of the interface strength, the particle fracture becomes more obvious, especially when the interface strength is 7 MPa. This suggests that a stronger interface strength prevents the relative motion of particle/binder interfaces, thereby enhancing the interactions between particles and increasing the possibility of a particle's fracture. The simulation results agree with Wu's predicted results [24].  Figure 23 presents the failure paths for different interface strengths under compression loading. As seen in Figure 23a-d, a large number of microcracks along the interface reveals that interface debonding is the primary failure mode for the weaker interface strength. With an increase of the interface strength, the particle fracture becomes more obvious, especially when the interface strength is 7 MPa. This suggests that a stronger interface strength prevents the relative motion of particle/binder interfaces, thereby enhancing the interactions between particles and increasing the possibility of a particle's fracture. The simulation results agree with Wu's predicted results [24].

Effect of the Strain Rate
The compressive properties of PBX show obvious strain rate dependence, which has also been reported by Gray et al. and Chen et al. [8,40]. We carried out a series of simulations under three different strain rates by controlling the displacement loading rate. Figure 24 shows the stress-strain curves for strain rates of 10 −4 s −1 , 10 −3 s −1 , and 10 −2 s −1 . As can be seen from the figure, the compressive strength increases as the strain rate increases, while the failure strain has almost no significant change compared with the compressive strength. Based on this phenomenon, Wiegand et al. proposed a constant critical strain failure criterion for PBXs [42].

Effect of the Strain Rate
The compressive properties of PBX show obvious strain rate dependence, which has also been reported by Gray et al. and Chen et al. [8,40]. We carried out a series of simulations under three different strain rates by controlling the displacement loading rate. Figure 24 shows the stress-strain curves for strain rates of 10 −4 s −1 , 10 −3 s −1 , and 10 −2 s −1 . As can be seen from the figure, the compressive strength increases as the strain rate increases, while the failure strain has almost no significant change compared with the compressive strength. Based on this phenomenon, Wiegand et al. proposed a constant critical strain failure criterion for PBXs [42].  Figure 25 shows the failure path of PBX 9501 under compression loading at three strain rates. Clearly, the failure path is similar, as indicated by the red line, which also illustrates that the failure path is insensitive to strain rate. The higher the strain rate, the more serious the particle failure (marked by arrows) is. This may explain why the stiffening in the binder at a high strain rate strengthens particle interactions in compression. This seems to imply that a stronger binder can increase the possibility of the particle's fracture.  Figure 25 shows the failure path of PBX 9501 under compression loading at three strain rates. Clearly, the failure path is similar, as indicated by the red line, which also illustrates that the failure path is insensitive to strain rate. The higher the strain rate, the more serious the particle failure (marked by arrows) is. This may explain why the stiffening in the binder at a high strain rate strengthens particle interactions in compression. This seems to imply that a stronger binder can increase the possibility of the particle's fracture.

Effect of the Strain Rate
The compressive properties of PBX show obvious strain rate dependence, which has also been reported by Gray et al. and Chen et al. [8,40]. We carried out a series of simulations under three different strain rates by controlling the displacement loading rate. Figure 24 shows the stress-strain curves for strain rates of 10 −4 s −1 , 10 −3 s −1 , and 10 −2 s −1 . As can be seen from the figure, the compressive strength increases as the strain rate increases, while the failure strain has almost no significant change compared with the compressive strength. Based on this phenomenon, Wiegand et al. proposed a constant critical strain failure criterion for PBXs [42].  Figure 25 shows the failure path of PBX 9501 under compression loading at three strain rates. Clearly, the failure path is similar, as indicated by the red line, which also illustrates that the failure path is insensitive to strain rate. The higher the strain rate, the more serious the particle failure (marked by arrows) is. This may explain why the stiffening in the binder at a high strain rate strengthens particle interactions in compression. This seems to imply that a stronger binder can increase the possibility of the particle's fracture.

Conclusions
We have presented a method that embeds cohesive elements throughout a microstructure model obtained by Voronoi tessellation, to study the response of PBX 9501 subjected to quasi-static loading. The calculated stress-strain curves agree well with the experimental results. PBX 9501 can be considered a brittle material with a low failure strain. The compressive strength of PBX 9501 is much higher than its tensile strength. The random particle size shape and interface debonding behavior of PBX 9501 produce local heterogeneous strain and failure distribution. Different failure mechanisms are observed. Interface debonding is the primary failure mode of the model under tensile loading. Under compressive loading, however, interface debonding and particle fracture play important roles.
We also used this model to study the effects of interface strength and strain rate on the mechanical properties and failure evolution of PBX 9501. The simulation results show that the higher strength of the particle/binder interface corresponds to greater tensile and compressive strength. The interface strength only affects the tensile failure strain but not the compressive failure strain. The compressive strength of PBX 9501 is rate dependent, while its failure strain is insensitive to the strain rate. The failure evolution path is relatively independent of the interface strength and strain rate. At a high interface strength or strain rate, a stronger interface force or binder can increase the probability of particle failure. These obtained results are important for a good understanding of the material's failure mechanisms and also useful for designing and improving the properties of advanced PBX composite materials by establishing the relationship between their meso-structures and performance.

Conclusions
We have presented a method that embeds cohesive elements throughout a microstructure model obtained by Voronoi tessellation, to study the response of PBX 9501 subjected to quasi-static loading. The calculated stress-strain curves agree well with the experimental results. PBX 9501 can be considered a brittle material with a low failure strain. The compressive strength of PBX 9501 is much higher than its tensile strength. The random particle size shape and interface debonding behavior of PBX 9501 produce local heterogeneous strain and failure distribution. Different failure mechanisms are observed. Interface debonding is the primary failure mode of the model under tensile loading. Under compressive loading, however, interface debonding and particle fracture play important roles.
We also used this model to study the effects of interface strength and strain rate on the mechanical properties and failure evolution of PBX 9501. The simulation results show that the higher strength of the particle/binder interface corresponds to greater tensile and compressive strength. The interface strength only affects the tensile failure strain but not the compressive failure strain. The compressive strength of PBX 9501 is rate dependent, while its failure strain is insensitive to the strain rate. The failure evolution path is relatively independent of the interface strength and strain rate. At a high interface strength or strain rate, a stronger interface force or binder can increase the probability of particle failure. These obtained results are important for a good understanding of the material's failure mechanisms and also useful for designing and improving the properties of advanced PBX composite materials by establishing the relationship between their meso-structures and performance.