Mesoscale Simulation to Study Constitutive Properties of TATB/F2314 PBX

Material Point Method (MPM) mesoscale simulation was used to study the constitutive relation of a polymer bonded explosive (PBX) consisting of 1,3,5-triamino-2,4,6-trinitrobenzene (TATB) and a fluorine polymer binder F2314. The stress-strain variations of the PBX were calculated for different temperatures and different porosities, and the results were found to be consistent with experimental observations. The stress-strain relations at different temperatures were used to develop the constitutive equation of the PBX by using numerical data fitting. Stress-strain data for different porosities were used to establish the constitutive equation by fitting the simulation data to an improved Hashion-Shtrikman model. The equation can be used to predict the shear modulus and bulk modulus of the PBX at different densities of the sample. The constitutive equations developed for TATB/F2314 PBX by MPM mesoscale simulation are important equations for the numerical simulations of the PBX at macroscale. The method presented in this study provides an alternative approach for studying the constitutive relations of PBX.


Introduction
The development of constitutive models for polymer bonded explosives (PBXs), especially for load-deformation or stress-strain relationships, has received considerable attention in the energetic material field. The development of the constitutive equations usually starts with a series of experimental measurements with sophisticated equipment and techniques, followed by curve fitting of the experimental data based on a proposed model, and further testing and refinement by predicting mechanical properties [1][2][3][4].
This approach has two shortcomings. First, the design of experimental set up can be challenging, considering the uncertain and extremely sensitive behavior of explosives. Experiments are very expensive and difficult to carry out. Second, the constitutive models obtained in such a way are macroscopic approximate constitutive models [5,6], or empirical constitutive models [7,8]. These models fail to account for complex, multiphase, multiscale structural, dynamical, and chemical properties. For instance, the macroscale deformations applied to PBX composites in experiments are generally not the same as the local deformation fields in a component crystal within the composite. It is important for constitutive modeling to bridge the particulate nature at the mesoscale to the mechanical properties at the macroscale. Computer simulation provides a safer and convenient way to study the mechanical properties of energetic material and has been used to assess and understand the constitutive relations of these materials. Most studies have been performed to understand the constitutive relations of energetic crystal materials or binders at atomic scale using molecular dynamics computation [9][10][11][12]. Only a few studies dealt with the mechanical properties of PBX composition at the microstate level.
The Material Point Method (MPM) [13][14][15] offers some advantages over other numerical techniques for simulations of mechanical properties of composites with complicated geometries found in PBX. It is easy in MPM to discretize complex geometries of PBX composites compared to mesh generation needed for Finite Element Method (FEM) calculations. Investigations reported by Banerjee et al. show that for simulations of large deformation of PBX with a high-strain rate, MPM has an advantage over other numerical methods [16][17][18].
This paper reports the mesoscale study of the constitutive relations of TATB/F 2314 PBX using the MPM approach. We start with the constitutive equations and equations of state of both components, TATB grain, and F 2314 binder and take the grain-grain, binder-binder, and grain-binder interactions into consideration. Based on the simulations results, the constitutive equations of TATB/F 2314 PBX are developed and can be used in macroscale simulations. This study provides a numerical simulation approach to study the constitutive relation of PBX.

Material Point Method
MPM was described by Sulsky et al. [13]. The use of a regular background grid in MPM has a lot of computational advantages: it does not require periodic remeshing steps and remapping of state variables, and is, therefore, more suitable for problems involving contact between materials, cracks, large deformations, and other high-velocity problems [14].
In MPM, materials are discretized into particles on a background grid and each particle is assigned a mass, coordinate, velocity, and other properties. The discrete momentum equation is given as: where m is the mass matrix, a is the acceleration vector, f ext is the external force vector, and f int is the internal force vector resulting from the divergence of the material stresses.
In each MPM step, all particle states are extrapolated to the nodes of the computational grid, then the mass matrix m, the nodal external forces f ext , and velocities v are calculated at individual nodes using the following equations: where represents a summation over all particles, i refers to individual nodes of the grid. m p is the particle mass, v p is the particle velocity, f ext p is the external force on the particle, and S ip is the generalized shape function of the node i evaluated at the particle p.
The nodal internal forces f int i is computed at the nodes as a volume integral of the divergence of the stress on the particles: where G ip is the shape function gradients of the node i evaluated at the particle p, while σ p and V p are stress and volume over the particle domain, respectively. As the momentum equations for grid nodes are solved, then mapped back to the material points, and state variables are calculated at the material points using the updated kinematic variables.

Interface friction Contact Algorithm
PBX consists of TATB grains and F 2314 polymeric binder, while F 2314 is a fluoropolymer binder containing a 1:4 molar ratios of comonomer vinylidene fluoride (VDF) and chlorotrifluoroethylene (CTFE). F 2314 is one of the most commonly used polymer binders for PBX [19]. Two components of PBX have different constitutive models and therefore, interface interaction between two components should be taken into consideration in mesoscale simulation. The interaction of material interfaces can be realized using the Multi-material algorithm [20]. In Multi-material MPM, each material extrapolates to its own velocity field. All nodes with more than one velocity field must determine if the materials are in contact. When the surfaces are determined to be in contact, the task is to modify nodal momenta and forces to reflect the implemented contact mechanics; details of the algorithm are as follows: To implement friction between TATB and F 2314 , the tangential traction experienced by frictional sliding is calculated by: where S slide is sliding traction and g N, A c , ∆v i , .... is any arbitrary friction law which may depend on various parameters such as contact normal pressure N, contact area A c , and relative sliding velocity ∆v i . An improved multi-material method was used to model contact between TATB and binder and the sliding force was calculated by: where S a is a shear adhesion strength of the interface, µ is the coefficient of friction, and d n is the normal momentum changes required to stick. S a takes zero in this calculation because for surfaces adhesion, the calculation is equivalent to MPM treatment of self-contact [20]. The calibration of interfacial parameters is taken from the literature where available.

Simulation Model
In this mesoscale interface model for TATB-based PBX, the analysis focuses on a two-phase meso-structure consisting of TATB grains and F 2314 binder. In MPM code [21], each material has a different constitutive equation and equation of state and two materials are discretized into material points. They carry the history-dependent state variables such as stresses and strains, as well as mass and kinematic variables such as position, velocity, and acceleration. An interface between materials is considered and all those state variables for PBX then are calculated at the material points using the updated kinematic variables of the two materials and their interfaces.
An elastic-plastic model and Mie-Grüneisen equation of state are used for TATB grains. The elastic-plastic material model [22] for deviatoric stress is expressed as: where, σ is stress, σ 0 is initial yield stress, all units are MPa, . ε is the strain rate, and C c and P c are Cowper-Symond strain rate parameters. ε P e f f is the effective plastic strain, E P is plastic hardening modulus, and β is the hardening parameter. The Cowper Symonds strain rate parameters were sent to zero as a constant strain rate was applied and only one yielded stress for PBX. The Mie-Grüneisen equation [22] of state is given by: where, P is pressure; µ = ρ/ρ 0 − 1 = 1 − V 0 /V, ρ 0 is for initial density, V 0 is initial volume; c is the bulk speed of sound at room temperature and pressure; γ 0 is Grüneisen's gamma at the reference state; E is internal energy per unit reference volume; s 1 ,s 2 , and s 3 are linear Hugoniot slope coefficients, and α is a first order volume correction factor. A viscoelastic constitutive model [23] and Tait equation [24] of state are used for the F 2314 binder. The Cauchy stress for F 2314 binder is expressed as: where σ represents the Cauchy stress. G(t) and K(t) are the deviatoric part of the shear relaxation modulus and bulk modulus of the constitutive behavior, respectively. e and ∆ refer to the deviatoric and hydrostatic portions of the Eulerian strain tensor, and t and τ refer to physical and reduced times, respectively. I is the unit tensor. The Tait equation of state is expressed as: where V(T, P) represents the final volume of compression at temperature T and pressure P; V 0 (T) is the temperature dependence of the volume at zero pressure. The parameter B is a function of temperature and depends on the specific system, and B(T) = CK(0, T). The parameter C is independent of structure and temperature, while K(0, T) is the temperature dependence of the bulk modulus at zero pressure. The parameters of TATB grain [25,26] and F 2314 [27] binder used in this study are listed in Table 1. In this study, an idealized meso-structure, with packed circular particles, was chosen for TATB grains. Research reported by Barua et al [28] showed that the particle morphology of explosive particles has some effects on the interface interaction, interface debonding, and energy distribution, while it has little effect on the effective elastic response of PBX, as reported by Dai [29] et al.
All particles are placed randomly within the packing box with a bimodal distribution of the diameters, which is consistent with experimental observations [30]. The mesoscopic calculation model based on the MPM was designed as shown in Figure 1a. The model consists of a rigid mold, rigid punch, TATB grains, F 2314 binder, and void between particles. Furthermore, interfacial properties between particles and binders that describe frictional contact were included in our model. The calculations were performed on a rectangular region of 1.0 × 1.2 mm. Materials are confined between three fixed rigid walls and a fourth wall that moves at a constant velocity of 0.3 m/s, which is about 0.01% of the bulk wave speed for the TATB grains in PBX. The velocity of the top surface yields overall strain rates of 250 s −1 . The diameters of TATB grains are chosen based on experimental data [30]. However, a minimum diameter is not selected, this results in a size range from 0.05 mm to 0.3 mm. F 2314 binder thickness ranging from 0.0013 to 0.0078 mm was selected to cover the outer layer of particles. The mass ratio between TATB and the binder was 95:5. The friction coefficient of 0.3 was used in the calculation. This value was also used for frictional contact between HMX particles and binder in an earlier simulation [20], where the HMX particles and binder were modeled with Mie-Grüneisen pressure response and an elastic shear response. The simulation results were consistent with experimental observation. Figure 1b shows the snapshot of final state of compression. The snapshot is similar to the realistic PBX morphology, and is used for discussions. 1 Parameter values of TATB were taken from reference [25] and reference [26]; 2 Parameter values of F2314 were taken from reference [27].
In this study, an idealized meso-structure, with packed circular particles, was chosen for TATB grains. Research reported by Barua et al [28] showed that the particle morphology of explosive particles has some effects on the interface interaction, interface debonding, and energy distribution, while it has little effect on the effective elastic response of PBX, as reported by Dai [29] et al.
All particles are placed randomly within the packing box with a bimodal distribution of the diameters, which is consistent with experimental observations [30]. The mesoscopic calculation model based on the MPM was designed as shown in Figure 1a. The model consists of a rigid mold, rigid punch, TATB grains, F2314 binder, and void between particles. Furthermore, interfacial properties between particles and binders that describe frictional contact were included in our model. The calculations were performed on a rectangular region of 1.0 × 1.2 mm. Materials are confined between three fixed rigid walls and a fourth wall that moves at a constant velocity of 0.3 m/s, which is about 0.01% of the bulk wave speed for the TATB grains in PBX. The velocity of the top surface yields overall strain rates of 250 s −1 . The diameters of TATB grains are chosen based on experimental data [30]. However, a minimum diameter is not selected, this results in a size range from 0.05 mm to 0.3 mm. F2314 binder thickness ranging from 0.0013 to 0.0078 mm was selected to cover the outer layer of particles. The mass ratio between TATB and the binder was 95:5. The friction coefficient of 0.3 was used in the calculation. This value was also used for frictional contact between HMX particles and binder in an earlier simulation [20], where the HMX particles and binder were modeled with Mie-Grüneisen pressure response and an elastic shear response. The simulation results were consistent with experimental observation. Figure 1b shows the snapshot of final state of compression. The snapshot is similar to the realistic PBX morphology, and is used for discussions. A compression simulation was performed for different temperatures. Temperature-force coupled loading was performed on the simulation model of the same fill rate with an initial porosity of 0.08 and a corresponding density of 1.787 g/cm 3 . The four loading temperatures are 293 K, 333 K, 353 K, and 373 K. The initial temperature of the PBX particles and surrounding environment was set to the corresponding temperature. Simulations using different densities of the PBX sample (1.761 g/cm 3 , 1.787 g/cm 3 , and 1.839 g/cm 3 ) with an initial temperature of 293 K, were performed to study the effect of porosity on mechanical properties as well. A compression simulation was performed for different temperatures. Temperature-force coupled loading was performed on the simulation model of the same fill rate with an initial porosity of 0.08 and a corresponding density of 1.787 g/cm 3 . The four loading temperatures are 293 K, 333 K, 353 K, and 373 K. The initial temperature of the PBX particles and surrounding environment was set to the corresponding temperature. Simulations using different densities of the PBX sample (1.761 g/cm 3 , 1.787 g/cm 3 , and 1.839 g/cm 3 ) with an initial temperature of 293 K, were performed to study the effect of porosity on mechanical properties as well.

Calculation of Mechanical Properties at Different Temperature
Stress and strain plots for compression of the PBX are presented in Figure 2 for the temperatures at 293 K, 333 K, 353 K, and 373 K. The stress-strain variations with different temperatures all show two stages: a fast-linear rise at early time followed by a slow non-linear variation, due to the elastic behavior and plastic deformation behavior of the PBX, respectively. Young's modulus, the slope of the linear changes of stress in the linear stage, clearly decreases with increasing temperatures. When the linear stage terminates at what is known as the yield point, stress changes slowly with strain. The yield stress of PBX estimated from the curve levels of stress-strain variations and the beginning of plastic deformation ranges from 13 MPa to 25 MPa and decreases as temperature increases. The yield stress of PBX is obviously lower than that of TATB, indicating the influence of binder on the properties and response of TATB. The notable temperature dependence of quasi-static compression for PBX comes from the fact that polymeric binder softening at high temperature affects the interface interaction between binders and particles and, therefore, changes the mechanical properties of PBX. The simulation results demonstrate the value of Young's modulus and yield stress decrease as the loading temperature increases, indicating a negative correlation of temperature with mechanical properties. However, the critical strains of PBX transition from linear elasticity to elastic-plastic stage were not affected. All these results are consistent with the experimental results of Guo [31] et al.
at 293 K, 333 K, 353 K, and 373 K. The stress-strain variations with different temperatures all show two stages: a fast-linear rise at early time followed by a slow non-linear variation, due to the elastic behavior and plastic deformation behavior of the PBX, respectively. Young's modulus, the slope of the linear changes of stress in the linear stage, clearly decreases with increasing temperatures. When the linear stage terminates at what is known as the yield point, stress changes slowly with strain. The yield stress of PBX estimated from the curve levels of stress-strain variations and the beginning of plastic deformation ranges from 13 MPa to 25 MPa and decreases as temperature increases. The yield stress of PBX is obviously lower than that of TATB, indicating the influence of binder on the properties and response of TATB. The notable temperature dependence of quasi-static compression for PBX comes from the fact that polymeric binder softening at high temperature affects the interface interaction between binders and particles and, therefore, changes the mechanical properties of PBX. The simulation results demonstrate the value of Young's modulus and yield stress decrease as the loading temperature increases, indicating a negative correlation of temperature with mechanical properties. However, the critical strains of PBX transition from linear elasticity to elastic-plastic stage were not affected. All these results are consistent with the experimental results of Guo [31] et al. The change of bulk modulus, shear modulus, and Young's modulus with temperatures are plotted in Figure 3. The elastic modulus of PBX was negatively correlated with temperature change. Bulk modulus did not change significantly from 293 K to 353 K because of the property of F2314 binder, meaning the binder softening due to temperature can be restored when the loading temperature is not particularly high. The three elastic moduli for TATB-based PBX at different temperatures obtained by fitting the stress-strain curves are compared with values from experiment [32] and molecular dynamics simulation [33] in Table 2. All the moduli decrease with increasing loading temperature because the flexibility of polymeric binder chain segments increases at higher temperature. This leads to an increase in the elasticity of PBX. The bulk modulus and shear modulus for PBX at 293 K are close to those obtained from the molecular dynamic simulation [33], and Young's modulus at 298 K is close to experimental observation [32]. In addition, Young's modulus at 333 K is in consistent with experimental values from 323 to 338 K [32]. The change of bulk modulus, shear modulus, and Young's modulus with temperatures are plotted in Figure 3. The elastic modulus of PBX was negatively correlated with temperature change. Bulk modulus did not change significantly from 293 K to 353 K because of the property of F 2314 binder, meaning the binder softening due to temperature can be restored when the loading temperature is not particularly high. The three elastic moduli for TATB-based PBX at different temperatures obtained by fitting the stress-strain curves are compared with values from experiment [32] and molecular dynamics simulation [33] in Table 2. All the moduli decrease with increasing loading temperature because the flexibility of polymeric binder chain segments increases at higher temperature. This leads to an increase in the elasticity of PBX. The bulk modulus and shear modulus for PBX at 293 K are close to those obtained from the molecular dynamic simulation [33], and Young's modulus at 298 K is close to experimental observation [32]. In addition, Young's modulus at 333 K is in consistent with experimental values from 323 to 338 K [32].    The compression strength of PBX at different temperatures obtained from stress-strain variations is plotted in Figure 4. The compression strength decreases with temperature and a rapid decrease above 353 K results from binder softening. Experimental dynamic mechanical thermal analysis of F2314 [34] have shown that the behavior of F2314 binder with temperature falls into following steps: glassy state, the transition zone between glassy state and rubber state, and rubber state, followed by the change to the viscous flow state. In the viscous flow state, the binder loses more adhesion to explosive particles and the compression strength of PBX decreased significantly. Experimental results [32,[35][36][37] of PBX compression strength are shown in Figure 4 for comparison. Different experimental data come from either different techniques or different temperatures. The simulation results of this study are within the range of these data.  The compression strength of PBX at different temperatures obtained from stress-strain variations is plotted in Figure 4. The compression strength decreases with temperature and a rapid decrease above 353 K results from binder softening. Experimental dynamic mechanical thermal analysis of F 2314 [34] have shown that the behavior of F 2314 binder with temperature falls into following steps: glassy state, the transition zone between glassy state and rubber state, and rubber state, followed by the change to the viscous flow state. In the viscous flow state, the binder loses more adhesion to explosive particles and the compression strength of PBX decreased significantly. Experimental results [32,[35][36][37] of PBX compression strength are shown in Figure 4 for comparison. Different experimental data come from either different techniques or different temperatures. The simulation results of this study are within the range of these data.  The aforementioned comparisons show that the results of MPM mesoscale simulation are reliable and can be applied to develop the constitutive model for the mechanical properties of TATB/F2314 PBX. Based on the stress-strain variations from MPM mesoscale simulation, the constitutive equation of the PBX can be established in the following steps [38].

Temperature (K) K (GPa) G (GPa) E (GPa
First, normalized stress n σ can be obtained as follows: = / n c σ σ σ (12) where c σ and σ are the compression strength and compression stress respectively. The stressstrain curves of PBX at different temperatures are normalized and plotted in Figure 5. The obtained stress-strain curves at different temperatures are basically coincident and therefore the stress-strain curve fitting can be performed at a fixed temperature.  The aforementioned comparisons show that the results of MPM mesoscale simulation are reliable and can be applied to develop the constitutive model for the mechanical properties of TATB/F 2314 PBX. Based on the stress-strain variations from MPM mesoscale simulation, the constitutive equation of the PBX can be established in the following steps [38].
First, normalized stress σ n can be obtained as follows: where σ c and σ are the compression strength and compression stress respectively. The stress-strain curves of PBX at different temperatures are normalized and plotted in Figure 5. The obtained stress-strain curves at different temperatures are basically coincident and therefore the stress-strain curve fitting can be performed at a fixed temperature.  The aforementioned comparisons show that the results of MPM mesoscale simulation are reliable and can be applied to develop the constitutive model for the mechanical properties of TATB/F2314 PBX. Based on the stress-strain variations from MPM mesoscale simulation, the constitutive equation of the PBX can be established in the following steps [38].
First, normalized stress n σ can be obtained as follows: = / n c σ σ σ (12) where c σ and σ are the compression strength and compression stress respectively. The stressstrain curves of PBX at different temperatures are normalized and plotted in Figure 5. The obtained stress-strain curves at different temperatures are basically coincident and therefore the stress-strain curve fitting can be performed at a fixed temperature.  The relationship between compression strength and loading temperature is expressed as: where T 0 is reference temperature and takes a value of 293 K, σ 0 is the compression strength at reference temperature, and α refers to a temperature effect factor. Each temperature has a corresponding compression strength that is listed in Table 3. σ 0 and α can be obtained by fitting the data in Table 3 into equation 13; the results are: σ 0 = 34.10 MPa and α = 2.357. The stress-strain relationship can be obtained by fitting the normalized stress-strain curve as follows: where the first term refers to the elastic stage, and the second term refers to visco-plastic stage. Parameters A, B, and n can be approximately regarded reflecting the essential properties of material and are independent of loading temperature. Normalized stress and strain data at 293 K in Figure 5 are used in the fitting process and the parameters were found to be: A = −0.1206, B = 1.084, and n = 0.2467. The constitutive equation eventually obtained by combining Equations (12)- (14) is: Figure 6 shows the variations of the stress of TATB/F 2314 PBX against strain calculated from the constitutive equation at the different temperature. The simulation results are also plotted for comparison. The stress-strain curve shows two stages: linear elastic deformation and a long plateau region after yielding. The plots calculated by the constitutive equation are in good agreement with simulation results. Experimental study [39] investigated the local damage initiation mechanism of polymer bonded sugar (PBS) and found that the stress-strain curve of PBS shows the similar behaviors. However, the experimentally found failure strain of PBS at about 6.5% was not observed in this study due to the absence of local failure strain for small strain deformation. The constitutive equations developed for TATB/F 2314 PBX for a mesoscale simulation can be a bridge for the numerical simulations at macroscale. Materials 2019, 12, x FOR PEER REVIEW 10 of 15  Figure 7 gives stress-strain curves of TATB-based PBX at different porosities. The stress-strain dependence on porosity follows the same trends. When the compression strain is smaller than 0.5%, the slope of the stress-strain curve increases with the increase of density, indicating that the elastic modulus of the PBX increases as the porosity of the sample decreases in the linear elastic stage. As the strain goes above 0.5%, the compression enters the plastic yield stage, where the stress varies slowly with the strain. At 99.0% theoretical density, all three porosity samples reach the stress value of 42.3 MPa. This result is supported by the experimental observation [40] that the pores between explosive particles almost disappeared when the loading force reached 40 MPa during the compaction of TATB granules.   Figure 7 gives stress-strain curves of TATB-based PBX at different porosities. The stress-strain dependence on porosity follows the same trends. When the compression strain is smaller than 0.5%, the slope of the stress-strain curve increases with the increase of density, indicating that the elastic modulus of the PBX increases as the porosity of the sample decreases in the linear elastic stage. As the strain goes above 0.5%, the compression enters the plastic yield stage, where the stress varies slowly with the strain. At 99.0% theoretical density, all three porosity samples reach the stress value of 42.3 MPa. This result is supported by the experimental observation [40] that the pores between explosive particles almost disappeared when the loading force reached 40 MPa during the compaction of TATB granules.  Figure 7 gives stress-strain curves of TATB-based PBX at different porosities. The stress-strain dependence on porosity follows the same trends. When the compression strain is smaller than 0.5%, the slope of the stress-strain curve increases with the increase of density, indicating that the elastic modulus of the PBX increases as the porosity of the sample decreases in the linear elastic stage. As the strain goes above 0.5%, the compression enters the plastic yield stage, where the stress varies slowly with the strain. At 99.0% theoretical density, all three porosity samples reach the stress value of 42.3 MPa. This result is supported by the experimental observation [40] that the pores between explosive particles almost disappeared when the loading force reached 40 MPa during the compaction of TATB granules.   Figure 8 presents the variations of bulk modulus, shear modulus, and Young's modulus for different porosities. All three moduli increase as the porosity decreases and this is consistent with experimental observations [41]. This study suggests that the smaller porosity of initial filling leads to higher compression strength and elastic modulus of PBX and, therefore, to better mechanical properties. Research found [42] that the initial damage is associated with pores in a compressed device and a poor density uniformity always appears in the parts where pores exist. These parts may become the failure zones owing to the existence of microcracks. Figure 8 shows experimental data for comparison. It is seen that the variation trend of PBX modulus with porosity obtained in this study is in good agreement with that of experimental results and that the simulated values are consistent with experimental values.  Figure 8 presents the variations of bulk modulus, shear modulus, and Young's modulus for different porosities. All three moduli increase as the porosity decreases and this is consistent with experimental observations [42]. This study suggests that the smaller porosity of initial filling leads to higher compression strength and elastic modulus of PBX and, therefore, to better mechanical properties. Research found [41] that the initial damage is associated with pores in a compressed device and a poor density uniformity always appears in the parts where pores exist. These parts may become the failure zones owing to the existence of microcracks. Figure 8 shows experimental data for comparison. It is seen that the variation trend of PBX modulus with porosity obtained in this study is in good agreement with that of experimental results and that the simulated values are consistent with experimental values. The stress-strain variations and modulus changes presented in Figures 7 and 8 reveal that the mechanical properties of PBX depend on the porosity systemically and that the compression strength and elastic modulus both significantly depend on the porosity of initial filling. Based on these factors, an improved Hashion-Shtrikman model [42] was used to fit the simulation data for the mechanical properties of different porosity.

Mechanical Properties with Different Porosity
For the improved Hashion-Shtrikman model, shear modulus and bulk modulus can be expressed as: where ( )  The stress-strain variations and modulus changes presented in Figures 7 and 8 reveal that the mechanical properties of PBX depend on the porosity systemically and that the compression strength and elastic modulus both significantly depend on the porosity of initial filling. Based on these factors, an improved Hashion-Shtrikman model [41] was used to fit the simulation data for the mechanical properties of different porosity.
For the improved Hashion-Shtrikman model, shear modulus and bulk modulus can be expressed as: where f (ρ/ρ 0 ) refers to interface bonding parameter, ρ 0 is the theoretical maximum density, ρ is the actual density of PBX, the theoretical calculation values, and G U and K U are upper bounds for effective shear modulus and effective bulk modulus in Hashion-Shtrikman model, respectively. These parameters take the values from reference [41]: G U = 10.84 GPa, and K U = 10.52 GPa.
The actual density is expressed as ρ = 1−φ v , where φ refers to porosity and v is the specific volume of material: v = σ p , of which σ p , ρ p , σ b , ρ b refer to mass percentage and the density of the particle and binder, respectively.
Considering the correlation between theoretical density and actual density, f (ρ/ρ 0 ) can be set as: where k and n are parameters related to each effective modulus. Three porosities of PBX (0.0932, 0.0798, and 0.0530) were used for this study. When porosity reaches 0, the theoretical maximum density of PBX is 1.942 g/cm 3 . Equations (19) and (20) can be obtained from the improved Hashion-Shtrikman Equations (16) and (17): ln k + n ln(ρ/ρ 0 ) = ln G e f f /G U (19) ln k + n ln(ρ/ρ 0 ) = ln K e f f /K U (20) where G e f f and K e f f are the values from this simulation study and k and n can be obtained from linear fitting of the simulation data. The resulting constitutive equations for shear modulus and bulk modulus with respect to porosity are: G/G U = 0.318(ρ/ρ 0 ) 12.11 (21) K/K U = 0.696(ρ/ρ 0 ) 12.01 .
Equations (21) and (22) can be used to predict the modulus of TATB/F 2314 PBX with different densities. The predicted variations of two moduli at different porosities are plotted in Figure 9a,b, respectively. The figures show that the shear modulus and bulk modulus of the PBX decreased exponentially with the density. A finite element simulation reported by Wei et al. [43] found that the effective elastic modulus of TATB-based PBX decreases exponentially with porosity, and these results were also shown in Figure 9a. It is worth noting that a discrepancy of predictive shear modulus occurs at lower porosity compared to finite element simulation results [43]. The discrepancy comes from the differences in the loading conditions and systems simulated: In the finite element simulation, TATB and binder have a mass ratio of 50:50. However, in the present calculation, the mass ratio between TATB and binder is 95:5. Increasing proportion of binder results in a significant decreasing of the effective elastic modulus. Moreover, the loading rates are completely different in two studies, which also cause a discrepancy of the shear modulus. However, the predictive shear modulus at low porosity is close to simulation result at 293 K [33], which indicates that the tendency of predictions is reliable.

Conclusions
Starting with the constitutive equations and the equations of state for both TATB grain and F2314 binder, the mechanical properties obtained from MPM mesoscale simulation can be applied to derive the constitutive relations of TATB/F2314 PBX. The constitutive equations of the PBX obtained from MPM simulation can be used to numerically determine mechanical properties of the PBX at a macroscale.
The stress-strain variations of TATB/F2314 PBX obtained from MPM mesoscale simulation for different temperatures were consistent with available experimental observations. The simulation results can be fitted numerically to the well-known constitutive law to determine the parameters

Conclusions
Starting with the constitutive equations and the equations of state for both TATB grain and F2314 binder, the mechanical properties obtained from MPM mesoscale simulation can be applied to derive the constitutive relations of TATB/F 2314 PBX. The constitutive equations of the PBX obtained from MPM simulation can be used to numerically determine mechanical properties of the PBX at a macroscale.
The stress-strain variations of TATB/F 2314 PBX obtained from MPM mesoscale simulation for different temperatures were consistent with available experimental observations. The simulation results can be fitted numerically to the well-known constitutive law to determine the parameters usually obtained from sophisticated experiments for the development of constitutive equations.
The stress-strain changes of the PBX for different porosities computed by MPM mesoscale simulation are in reasonable agreement with available experimental observations or other numerical simulations. The computed properties can be fitted to an improved Hashion-Shtrikman equation to determine the parameters that are important to the constitutive equations of the PBX which can be used to predict the shear modulus and bulk modulus of the PBX at different densities (and, therefore, different porosities).
The development of the constitutive relationship is an important research topic in the field and the method presented in this paper provides a numerical approach for the investigation of the constitutive equations. This would be useful for multiple scale simulation of PBX, where the computation is carried from microscale to macroscale.