Material Point Method Simulation of the Equation of State of Polymer-Bonded Explosive under Impact Loading at Mesoscale

: Mesoscale simulation using the material point method (MPM) was conducted to study the pressure–volume (PV) variations of Octahydro-1,3,5,7-Tetranitro-1,2,3,5-Tetrazocine (HMX) / Estane polymer-bonded explosive (PBX) under impact loading. The PV isotherms and Hugoniot data were calculated for the di ﬀ erent porosities and binder volume fractions. The PV isotherms were used to determine the parameters for the Birch– Murnaghan equation of state (EOS) for the PBX. From the EOS, the isothermal bulk modulus ( K 0 ) and its pressure derivative ( K (cid:48) 0 ) were calculated. Additionally, the pseudo particle velocity and pseudo shock velocity variations were used to obtain the bulk wave speed c and dimensionless coe ﬃ cient s for the Mie–Grüneisen EOS. The simulations provide an alternative approach for determining an EOS that is consistent with experimental observations. The ﬀ includes the e ﬀ ective factors of the microstructures of the PBX and the mass ratio of the energetic crystals and the binder contained in the PBX. The isothermal bulk modulus ( K 0 ) and its pressure derivative ( K (cid:48) 0 ) for the PBX were determined by analyzing the obtained isotherms using the Birch–Murnaghan equation of state. The results demonstrate that the established e ﬀ ective EOS meets the criteria based on the basic thermodynamic discussion. The bulk wave speed c and dimensionless parameter s in the Mie–Grüneisen EOS of the PBX obtained by ﬁtting u s - u p variations with the di ﬀ erent porosities and di ﬀ erent binder volume fractions are also e ﬀ ective.


Introduction
Polymer-bonded explosive (PBX), a polymer-matrix granular composite, is mainly composed of viscoelastic polymers and elastoplastic explosives particles. The huge difference between the properties of polymers and explosives makes it difficult to determine the thermodynamic and mechanical properties of the granular composite. An equation of state (EOS) is an equation that summarizes the relationship between several state functions (such as pressure P, volume V, temperature T, and internal energy E) of a system. An EOS can provide valuable microscopic insight into the bulk properties of energetic composite materials. Two types of EOS, reacted and unreacted, have been developed for energetic materials and have been the subject of extensive studies in the field. The commonly used model for reacted explosives is the Jones-Wilkins-Lee (JWL) EOS; for unreacted explosives, the Mie-Grüneisen EOS is common [1,2].
The determination of an EOS for a PBX is typically based on experimental studies, although the constitutive models of each component have also been studied by computational approaches. To obtain an EOS for an unreacted high explosive, various attempts have been reported involving static compression experiments [3][4][5][6], where the experimental data under static high pressure and room temperature conditions were derived from X-ray diffraction measurements. Analyses of isotherms from an empirical EOS allow for the determination of the isothermal bulk modulus (K 0 ) and its pressure derivative (K 0 ) for PBXs. The isotherms are transformed to a pseudo particle velocity (u s ) and pseudo shock velocity (u p ) plane using the Rankine-Hugoniot jump conditions [4].
Experimental determination of the isothermal high-pressure EOS for a PBX using conventional diffraction techniques is expensive and is not always feasible [7]. The experimental tests can be difficult to duplicate due to the complex environments and the sensitive behavior of explosives. Experimental studies also rarely include the effects of microstructures.
The high-pressure EOS of explosive crystals or binders has been intensively studied using molecular dynamics and other simulation methods [8][9][10][11][12]. Some fitting parameters of an isothermal EOS for different energetic explosives were also published [13]. However, the numerical simulations of an EOS for PBXs have rarely been reported.
The material point method (MPM) [14][15][16] has some advantages over other numerical techniques to study the thermodynamic and mechanical properties of composites with complicated morphologies and a wide range of particle sizes such as PBXs. As a meshless method, MPM is easier to use for discretizing the complex geometries of PBX composites compared to the periodical remeshing required for finite element method (FEM) calculations. Investigations reported by Banerjee et al. show that MPM is successful in large deformation and interface contact simulations of a PBX with a high-strain rate loading test [17][18][19]. Material points that are surrounded by background grids in an MPM avoid drawbacks for the costly searching for contact surfaces and/or remeshing by numerical approaches.
In this study, the simulation starts with the component properties of Octahydro-1,3,5,7-Tetranitro-1,2,3,5-Tetrazocine (HMX) grain and an Estane binder, such as the constitutive equations and equations of state. Interface interactions between the grain and binder are taken into consideration. Then, the MPM technique is used to compute the parameters for the EOS of HMX/Estane PBX composite at mesoscale. The parameters also account for the effects of porosity and binder volume fraction. The derived EOS can be used in a macroscale simulation for the thermodynamics study of the PBX.

Material Point Method
MPM, as described by Sulsky et al. [16], is a meshless-based method through the presence of a background grid and has a lot of advantages in computational mechanics. It is more suitable for solving problems such as large deformation, crack propagation, surface contact, and also shock loading simulation as it eliminates the need for periodic remeshing steps and the remapping of state variables [17].
In MPM, materials are discretized into a collection of particles, namely "material points", on a background mesh/grid, and each particle has its mass, coordinate, velocity, and other material properties. The discrete momentum equation can be written as: where m is the mass matrix, a is the acceleration vector, f ext is the external force vector, and f ext is the internal force vector generated by the divergence of the material stresses. In each step of an MPM simulation, all particle quantities, including particle position, velocity, and mass, are extrapolated to the nodes at the corners of the cell in the background grid, then using Equations (2) to (4), the mass matrix m, velocities v, and the external forces f ext are calculated at individual nodes using the material point values: where Σ summates values over all particles, i refers to individual grid nodes, m p is the particle mass, v p is the particle velocity, f ext p is the particle external force, and S ip is the nodal generalized shape function calculated at the each particle.
The nodal internal forces f int i are calculated as a volume integral of the divergence of particle stresses: where G ip is the nodal shape function gradient calculated at each particle, while σ p and V p are the particle stress and particle volume, respectively. The shape function S ip and shape function gradient G ip extrapolated from node i to particle p are as follows: where V p is the particle volume, χ p (x) is a particle basis function that equals to 1 within the particle's domain and zero elsewhere, N i (x) is a nodal shape function, and ∇ refers to gradient calculation. In each time step, all particle quantities and initial conditions are mapped to the background mesh nodes, and the momentum of each node is updated using Equation (1). Then, the information, such as the momentum carried on the node, is mapped back to the particle. When the next time step starts, the mapping process of the previous time is repeated until the specified time or the specified condition is reached, thereby completing the solution to the actual problem.

Interface Contact
An HMX-Estane PBX consists of HMX grains and an Estane polymeric binder matrix. Estane 5703 is a commonly used polymer in PBXs for various uses, such as being a component of binder in the typical highly filled energetic material PBX9501 [20]. HMX and Estane binders have different constitutive models and equations of state in microscale; therefore, we took into consideration the interface interaction between the two components in a mesoscale simulation. In the MPM, a multimaterial algorithm can realize the interaction calculation of material interfaces [21]. In a multimaterial MPM, each material extrapolates to its own velocity field on nodes. If a node in the background grid has more than one velocity field, contact detection should be determined. When nodes on surfaces are determined to be in contact, the contact mechanical laws are implemented to modify nodal momenta and forces; details of the multimaterial algorithm are as follows.
To perform friction contact law between HMX grain and the Estane binder, the tangential traction generated by frictional sliding is given as: where S slide is sliding traction, µ is the frictional coefficient, ∆v i is the relative sliding velocity, k refers to a linear coefficient of sliding velocity, N is the contact normal pressure, and S a is a shear adhesion strength. Some earlier studies [22][23][24] used the Coulomb friction law to investigate the frictional contact, intergranular interactions, and frictional heating during shock loading or shear deformation. The simple Coulomb friction is a special case [21], taking S a = 0, k = 0, and ∆v i = 0 for Equation (8). Then, the sliding force can be written as: where d n is the normal momentum changes required to stick, A c is the contact area, and ∆t is the time step. During frictional sliding, the momentum change ∆p becomes wheren andt are the normal and tangential vectors, respectively. The momentum change reflects the implemented sliding frictional contact mechanics. Otherwise, two materials would stick during the simulations. Such contact obeys a stick contact law, where the final momentum does not change, and two materials will move in a single velocity filed. The MPM automatically handles the stick contact law using a detecting contact approach, and the Coulomb friction contact conditions are applied via the grid analyses. The coefficient of friction can be taken from the literature where it is available. Some frictional contact parameters are appropriately inferred in previous work [25].

Simulations Performed
The interface model for the HMX/Estane PBX in this mesoscale analysis is a two-phase mesostructure that consists of HMX grains and the Estane binder. In the MPM code [26], two materials are described by different constitutive equations and equations of state, and they are discretized into material points over respective domains. These material points own the initial conditions of history-dependent state variables, such as stress and strain, and kinematic variables, such as position, velocity, and acceleration. Interfaces between grains and binders are taken into consideration, and the history-dependent state variables for the PBX are calculated using the updated kinematic variables on nodes of the two materials and their interfaces. HMX grains are described by an elastic-plastic model and the Mie-Grüneisen EOS. The elastic-plastic material model [27] for deviatoric stress is given as: where σ and σ 0 are the stress and initial yield stress of the material, respectively; . ε is the strain rate; C c and P c are Cowper-Symond strain rate parameters, respectively; ε P e f f is the effective plastic strain [28], which is the incremental plastic strain tensor and is defined as ε P e f f = 2 3 dε P ; and E P and δ are the plastic hardening modulus and hardening parameter, respectively. The Cowper-Symonds strain rate parameters were set to zero as a constant strain rate was applied and only one yielded stress for the PBX.
The Mie-Grüneisen equation [29] of the 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 atmospheric pressure, γ 0 is Grüneisen's gamma constant at the reference state, U is internal energy per reference volume, and s is the linear Hugoniot slope coefficient. The Estane binder is described by a viscoelastic constitutive model [30] and the Mie-Grüneisen equation of state. The Cauchy stress for the Estane binder is given as: where σ is 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 parts of the Eulerian strain tensor, respectively; t and τ refer to physical current and relaxation times, respectively; and I is the identity tensor. The parameters of HMX grain [31,32] and the Estane binder [13,32] used in this mesoscale simulation are listed in Table 1.

Construction of Simulation Systems
In this study, an idealized mesostructure was packed with circular particles and consisted of circular HMX grains and an Estane binder covering the outer layer of grains. All particles with a bimodal distribution of the diameters were chosen so as to be consistent with experimental observations [33] and were placed randomly in the packing box. The system designed for MPM mesoscopic calculation is shown in Figure 1. The system is composed of rigid mold, HMX grains, the Estane binder, and voids between particles. The void in the system was set as a vacuum. Moreover, interfacial properties described by frictional contact between grains and binders are considered in the model. The system model was initially set at a stress-free state, then impact loading was implemented by applying a constant vertical velocity at the top boundary. The other three boundaries were confined in fixed rigid walls to avoid lateral expansion. Liu [34] reported that the representative volume element (RVE) has a minimum size of 1.5 mm for a sample with an average HMX grain size of 125 µm. Specifically, we chose the HMX grain size in the range of 10-150 µm with a bimodal distribution, and the resulting RVE domain of the sample size is 0.5 × 1.5 mm. The distributions of the HMX grain size for all configurations are plotted in Figures S1 and S2 in the Supplementary Materials. The microstructural RVE has an aspect ratio of 3:1 (1.5 × 0.5 mm), allowing the stress wave transmission process that resulted from impact loading to be resolved. The Estane binder has a homogenized approximation in order for it to be considered an effective binder matrix. A homogenized binder with a different thickness covered the outer layer of the HMX grains, and the thickness of the binder was selected by the mass ratio between the HMX and the Estane binder. A frictional coefficient of 0.3 was used in the calculation. An earlier study used this value for the frictional contact between HMX grains and the binder [21], where the HMX grains and binder were simulated with a Mie-Grüneisen bulk response and elastic shear response. Its simulation results were close to the experimental observations. A parametric study was carried out in this study, focusing on the effects of (i) the binder volume fraction and (ii) initial porosity of the PBX. For all simulations, the initial temperature was set to 300K. The impact velocity at the top boundary of the specimen was chosen to be 250 m/s, resulting in an overall strain rate of 1.67 × 10 5 s −1 .
The four specimens, with the initial porosities of 0.113, 0.073, 0.033, and 0.003, which correspond to the densities of 1.654, 1.728, 1.803, and 1.859 g/cm 3 , respectively, as shown in Figure 2, were studied. The theoretical maximum density (TMD) of HMX/Estane PBX is 1.864 g/cm 3 . The thickness of the Estane binder ranges from 0.026 to 3.9 μm and was selected to cover the outer layer of HMX grains. The mass ratio between the HMX and the binder is 95:5, and the impact velocity for different porosities is 250 m/s. Simulations using different volume fractions of the binder were performed with an initial porosity of 0.003 and a constant velocity of 250 m/s. PBX9501 is theoretically made up of 95% HMX and 5% Estane binder and plasticizer by weight. The volume fraction of the HMX and polymeric binder was calculated to be 92.7% and 7.3%, respectively; however, the apparent binder volume A parametric study was carried out in this study, focusing on the effects of (i) the binder volume fraction and (ii) initial porosity of the PBX. For all simulations, the initial temperature was set to 300K. The impact velocity at the top boundary of the specimen was chosen to be 250 m/s, resulting in an overall strain rate of 1.67 × 10 5 s −1 .
The four specimens, with the initial porosities of 0.113, 0.073, 0.033, and 0.003, which correspond to the densities of 1.654, 1.728, 1.803, and 1.859 g/cm 3 , respectively, as shown in Figure 2, were studied. The theoretical maximum density (TMD) of HMX/Estane PBX is 1.864 g/cm 3 . The thickness of the Estane binder ranges from 0.026 to 3.9 µm and was selected to cover the outer layer of HMX grains. The mass ratio between the HMX and the binder is 95:5, and the impact velocity for different porosities is 250 m/s. A parametric study was carried out in this study, focusing on the effects of (i) the binder volume fraction and (ii) initial porosity of the PBX. For all simulations, the initial temperature was set to 300K. The impact velocity at the top boundary of the specimen was chosen to be 250 m/s, resulting in an overall strain rate of 1.67 × 10 5 s −1 .
The four specimens, with the initial porosities of 0.113, 0.073, 0.033, and 0.003, which correspond to the densities of 1.654, 1.728, 1.803, and 1.859 g/cm 3 , respectively, as shown in Figure 2, were studied. The theoretical maximum density (TMD) of HMX/Estane PBX is 1.864 g/cm 3 . The thickness of the Estane binder ranges from 0.026 to 3.9 μm and was selected to cover the outer layer of HMX grains. The mass ratio between the HMX and the binder is 95:5, and the impact velocity for different porosities is 250 m/s. Simulations using different volume fractions of the binder were performed with an initial porosity of 0.003 and a constant velocity of 250 m/s. PBX9501 is theoretically made up of 95% HMX and 5% Estane binder and plasticizer by weight. The volume fraction of the HMX and polymeric binder was calculated to be 92.7% and 7.3%, respectively; however, the apparent binder volume calculated to be 92.7% and 7.3%, respectively; however, the apparent binder volume fractions in the micrograph would actually be 23%-26% [35] because many of the fine particles were dislodged from the sample during its preparation. This is called the "dirty binder" effect by Xue et al. [36]; the "dirty binder" approximation is effective when it is required to reduce the computation for systems with very small particles. In our work, the HMX particle has a range of 10-150 µm; in such a case, the "dirty binder" approximation is not necessary, so we only considered the effect of the volume fraction of the pure binder on the PBX composite. Computationally generated microstructures with a binder volume fraction of 0%-20% are shown in Figure 3, and the corresponding volume ratios between the HMX and the Estane binder are 100:0, 95:5, 90:10, and 80:20. fractions in the micrograph would actually be 23%-26% [35] because many of the fine particles were dislodged from the sample during its preparation. This is called the "dirty binder" effect by Xue et al. [36]; the "dirty binder" approximation is effective when it is required to reduce the computation for systems with very small particles. In our work, the HMX particle has a range of 10-150 μm; in such a case, the "dirty binder" approximation is not necessary, so we only considered the effect of the volume fraction of the pure binder on the PBX composite. Computationally generated microstructures with a binder volume fraction of 0%-20% are shown in Figure 3, and the corresponding volume ratios between the HMX and the Estane binder are 100:0, 95:5, 90:10, and 80:20. Adopting the irregular circular particles for the HMX grains is for convenience. Barua et al. [37] showed that the grain morphology of explosives has significant effects on interface interaction, interface debonding, and energy distribution, and it may also affect the Hugoniot data compared to the real microstructure. However, Barua et al. also indicated that the little effects on the overall mechanical response were observed for grain morphology changing from a circular to diamond shape [37]. In this study, this effect is not considered but will be addressed in a future publication.

Justification of Model
For the justification of the model, the mesostructure effects on the compression behavior of HMX grains and the HMX/Estane PBX under impact loading were studied. At low pressure, weak shocks are predominated by compressional effects, and shock waves can be defined as a transmission wave front, across which a discontinuous adiabatic jump in state variables occurs [38]. Consequently, isothermal compression can be approximated by the adiabatic compression process, and the pressure-volume (PV) data can be transferred to pseudo particle velocity and pseudo shock velocity in terms of the method of Olinger and Cady et al. [39]. The Rankine-Hugoniot jump conditions are expressed as follows: where and are the initial pressure and volume of the system, respectively. P and V describe the pressure and volume at different compressions, respectively.
The calculated PV curve of HMX at room temperature is illustrated in Figure 4a and compared with available experimental observations [40][41][42][43] and numerical simulation results [44][45][46]. It is seen Adopting the irregular circular particles for the HMX grains is for convenience. Barua et al. [37] showed that the grain morphology of explosives has significant effects on interface interaction, interface debonding, and energy distribution, and it may also affect the Hugoniot data compared to the real microstructure. However, Barua et al. also indicated that the little effects on the overall mechanical response were observed for grain morphology changing from a circular to diamond shape [37]. In this study, this effect is not considered but will be addressed in a future publication.

Justification of Model
For the justification of the model, the mesostructure effects on the compression behavior of HMX grains and the HMX/Estane PBX under impact loading were studied. At low pressure, weak shocks are predominated by compressional effects, and shock waves can be defined as a transmission wave front, across which a discontinuous adiabatic jump in state variables occurs [38]. Consequently, isothermal compression can be approximated by the adiabatic compression process, and the pressure-volume (PV) data can be transferred to pseudo particle velocity u p and pseudo shock velocity u s in terms of the method of Olinger and Cady et al. [39]. The Rankine-Hugoniot jump conditions are expressed as follows: where P 0 and V 0 are the initial pressure and volume of the system, respectively. P and V describe the pressure and volume at different compressions, respectively. The calculated PV curve of HMX at room temperature is illustrated in Figure 4a and compared with available experimental observations [40][41][42][43] and numerical simulation results [44][45][46]. It is seen that the calculated results are in good agreement with the experimental data and match those from the molecular dynamics simulation studies by Sewell et al. [46]. that the calculated results are in good agreement with the experimental data and match those from the molecular dynamics simulation studies by Sewell et al. [46]. The us-up variations of HMX converted from the static PV data using Equations (14) and (15) are plotted in Figure 4b and also compared with experimental results [40,42,43] and numerical simulation data [44][45][46]. The calculated us-up data also show a good agreement with experimental observations. Pseudo particle velocity and pseudo shock velocity are related by a linear relation [47], which is a Hugoniot EOS: where the constant c and s are parameters of a Mie-Grüneisen EOS; c is related to the bulk wave speed, and s is a dimensionless coefficient. Curve fitting yields c = 2.759 and s = 2.545, which are in good agreement with the experimental values (c = 2.74 km/s and s = 2.60) [40]. Figure 5a shows a PV isotherm of the PBX (95% HMX: 5% Estane) compared with experimental and other previous simulation data [48][49][50][51][52]. The authors in [50] compared Hugoniot data for Estane and the PBX9501 binder, where 2.5% plasticizer is mixed, and the results showed similar changes for the two components. The pressure required for the PBX in this study is close to the experimental data of reference [48] at the same compression ratio and is slightly higher than the experimental data when the compression ratio is lower than 0.8. Divergence may be generated from different simulation techniques or different loading temperatures. The u s -u p variations of HMX converted from the static PV data using Equations (14) and (15) are plotted in Figure 4b and also compared with experimental results [40,42,43] and numerical simulation data [44][45][46]. The calculated u s -u p data also show a good agreement with experimental observations. Pseudo particle velocity u p and pseudo shock velocity u s are related by a linear relation [47], which is a Hugoniot EOS: where the constant c and s are parameters of a Mie-Grüneisen EOS; c is related to the bulk wave speed, and s is a dimensionless coefficient. Curve fitting yields c = 2.759 and s = 2.545, which are in good agreement with the experimental values (c = 2.74 km/s and s = 2.60) [40]. Figure 5a shows a PV isotherm of the PBX (95% HMX: 5% Estane) compared with experimental and other previous simulation data [48][49][50][51][52]. The authors in [50] compared Hugoniot data for Estane and the PBX9501 binder, where 2.5% plasticizer is mixed, and the results showed similar changes for the two components. The pressure required for the PBX in this study is close to the experimental data of reference [48] at the same compression ratio and is slightly higher than the experimental data when the compression ratio is lower than 0.8. Divergence may be generated from different simulation techniques or different loading temperatures.  Figure 5b is a linear fit of the simulation data in this study; the two dash dotted lines are us-up linear fits taken from [49,52]. The Hugoniot data can be fitted as us = 2.356up + 2.71 in this study, which is coincident with us = 2.3up + 2.65 given by Gustavsen et al. [49]. Moreover, simulation data in this study are close to the results from [50] at low up (0 < up < 0.6 km/s). The aforementioned comparisons indicate that results of the MPM mesoscale simulation are reliable and that this model is applicable to develop the equation of state for investigating the pressure-volume variation behaviors of HMX/Estane PBX, as well as microstructure effects on EOS of PBX in mesoscale.

PV Isotherms and Hugoniot Analysis of HMX/Estane PBX
The PV isotherms of HMX/Estane PBX at different porosities and different binder volume fractions are plotted in Figure 6. Figure 6a presents the PV data of the PBX at different porosities. The four porosities are 0.113, 0.073, 0.033, and 0.003, and their corresponding packing densities are 1.654, 1.728, 1.803, and 1.864 g/cm 3 . Each PV isotherm shows two trends: a linear variation at an early time (high V/V0), followed by rapid curved growth (low V/V0). The linear variation is associated with the rearrangement and deformation of particles squeezing out void space. At this stage, a small compression pressure results in a large volumetric change. After void space is removed, the PBX becomes highly compact, and the volume change slows down with the increase in compression pressure. It is seen from Figure 6a that at the same compression ratio the compression pressure decreases with the increase of porosity. A sample with smaller porosity has a higher packing density and therefore lower compressibility. It requires a large compression pressure for the same compression ratio, as shown in the experimental observation [3].   Figure 5b is a linear fit of the simulation data in this study; the two dash dotted lines are u s -u p linear fits taken from [49,52]. The Hugoniot data can be fitted as u s = 2.356u p + 2.71 in this study, which is coincident with u s = 2.3u p + 2.65 given by Gustavsen et al. [49]. Moreover, simulation data in this study are close to the results from [50] at low u p (0 < u p < 0.6 km/s). The aforementioned comparisons indicate that results of the MPM mesoscale simulation are reliable and that this model is applicable to develop the equation of state for investigating the pressure-volume variation behaviors of HMX/Estane PBX, as well as microstructure effects on EOS of PBX in mesoscale.

PV Isotherms and Hugoniot Analysis of HMX/Estane PBX
The PV isotherms of HMX/Estane PBX at different porosities and different binder volume fractions are plotted in Figure 6. Figure 6a presents the PV data of the PBX at different porosities. The four porosities are 0.113, 0.073, 0.033, and 0.003, and their corresponding packing densities are 1.654, 1.728, 1.803, and 1.864 g/cm 3 . Each PV isotherm shows two trends: a linear variation at an early time (high V/V 0 ), followed by rapid curved growth (low V/V 0 ). The linear variation is associated with the rearrangement and deformation of particles squeezing out void space. At this stage, a small compression pressure results in a large volumetric change. After void space is removed, the PBX becomes highly compact, and the volume change slows down with the increase in compression pressure. It is seen from Figure 6a that at the same compression ratio the compression pressure decreases with the increase of porosity. A sample with smaller porosity has a higher packing density and therefore lower compressibility. It requires a large compression pressure for the same compression ratio, as shown in the experimental observation [3].  Figure 6b shows the PV isotherms of the PBX at different binder volume fractions. The four volume fractions are 0%, 5%, 10%, and 20%, where 0% represents only the HMX component in the model. The PV isotherms show that at the same compression ratio, the compression pressure decreases with the increase in the binder volume fraction. Due to the addition of more binders, the elastic moduli of the PBX decrease, and the PBX becomes more elastoplastic. Bulk modulus, which is a form of elastic modulus and reflects the capacity of PBXs to resist external compression loading, decreases with the increasing binder volume fraction. Consequently, the compression pressure decreases at the same compression ratio. The simulations show that an increase in the porosity or the binder volume fraction in the PBX sample decreases in the compression pressure at the same compression ratio. Furthermore, comparisons between the curves of porosity and binder volume fraction suggest that porosity has a greater influence on compression pressure than binder volume fraction.
The Hugoniot data derived from the PV data of the PBX at different porosities and different binder volume fractions are plotted in Figure 7a,b, respectively. Figure 7a shows the Hugoniot data of the PBX at different porosities of 0.113, 0.073, 0.033, and 0.003, the corresponding packing densities of which are 1.654, 1.728, 1.803, and 1.864 g/cm 3 , respectively. The linear fitted parameters of the Hugoniot data are listed in Table 2. It is seen from Table 2 that c decreases gradually with increasing porosity. However, no significant changes are observed for s when porosity increases.
Experimental work by Gustavsen et al. [49] showed a substantial difference of 0.18 km/s in the parameter c of the calculated Hugoniots for PBX 9501 at typical densities (range from 1.80 to 1.837 g/cm 3 ) and at the theoretical maximum density (TMD = 1.86 g/cm 3 ).  Figure 6b shows the PV isotherms of the PBX at different binder volume fractions. The four volume fractions are 0%, 5%, 10%, and 20%, where 0% represents only the HMX component in the model. The PV isotherms show that at the same compression ratio, the compression pressure decreases with the increase in the binder volume fraction. Due to the addition of more binders, the elastic moduli of the PBX decrease, and the PBX becomes more elastoplastic. Bulk modulus, which is a form of elastic modulus and reflects the capacity of PBXs to resist external compression loading, decreases with the increasing binder volume fraction. Consequently, the compression pressure decreases at the same compression ratio. The simulations show that an increase in the porosity or the binder volume fraction in the PBX sample decreases in the compression pressure at the same compression ratio. Furthermore, comparisons between the curves of porosity and binder volume fraction suggest that porosity has a greater influence on compression pressure than binder volume fraction.
The Hugoniot data derived from the PV data of the PBX at different porosities and different binder volume fractions are plotted in Figure 7a,b, respectively. Figure 7a shows the Hugoniot data of the PBX at different porosities of 0.113, 0.073, 0.033, and 0.003, the corresponding packing densities of which are 1.654, 1.728, 1.803, and 1.864 g/cm 3 , respectively. The linear fitted parameters of the Hugoniot data are listed in Table 2. It is seen from Table 2 that c decreases gradually with increasing porosity. However, no significant changes are observed for s when porosity increases. Experimental work by Gustavsen et al. [49] showed a substantial difference of 0.18 km/s in the parameter c of the calculated Hugoniots for PBX 9501 at typical densities (range from 1.80 to 1.837 g/cm 3 ) and at the theoretical maximum density (TMD = 1.86 g/cm 3 ).    Figure 7b shows the Hugoniot data of the PBX at the binder volume fractions of 0%, 5%, 10%, and 20%. The parameters s and c obtained from the linear curve fitting using Equation (16) are presented in Table 3. It is seen that with the increase of binder volume fraction, the slope of the us-up curve, namely the dimensionless parameter s of the Mie-Grüneisen EOS, also decreases. On the other hand, the parameter c does not have a great decrease. Table 3 suggests that mixing a small content of the binder with HMX particles reduces the parameters c and s of the EOS obviously, and further addition of the binder has little effect on the parameters' values. Experimental shock loading analyses of the PBX [51] reveal that for the PBX with the amount of the binder beyond a threshold value, the binder's properties rather than the binder's volume fraction exert a measurable influence on the shock sensitivity of the PBX. The authors in [53] also reported that the EOS and constitutive relation of the PBX mainly depend on the properties of the binder (such as density, elastic modulus, etc.), but this has limited dependence on the binder volume fraction. This is consistent with the simulation results.
Simulated Hugoniot data for the PBX suggest that porosity and binder volume fraction both affect the equation of state under impact loading but in different ways. The porosity mainly affects parameter c, in that an increase in porosity leads to a decrease in c. However, the binder volume fraction mostly affects parameter s: An increase in binder volume fraction gives rise to a decrease in s. Therefore, the EOS of PBXs depends on both microstructures, such as porosity and mass ratio, of the HMX/binder.   Figure 7b shows the Hugoniot data of the PBX at the binder volume fractions of 0%, 5%, 10%, and 20%. The parameters s and c obtained from the linear curve fitting using Equation (16) are presented in Table 3. It is seen that with the increase of binder volume fraction, the slope of the u s -u p curve, namely the dimensionless parameter s of the Mie-Grüneisen EOS, also decreases. On the other hand, the parameter c does not have a great decrease. Table 3 suggests that mixing a small content of the binder with HMX particles reduces the parameters c and s of the EOS obviously, and further addition of the binder has little effect on the parameters' values. Experimental shock loading analyses of the PBX [51] reveal that for the PBX with the amount of the binder beyond a threshold value, the binder's properties rather than the binder's volume fraction exert a measurable influence on the shock sensitivity of the PBX. The authors in [53] also reported that the EOS and constitutive relation of the PBX mainly depend on the properties of the binder (such as density, elastic modulus, etc.), but this has limited dependence on the binder volume fraction. This is consistent with the simulation results. Simulated Hugoniot data for the PBX suggest that porosity and binder volume fraction both affect the equation of state under impact loading but in different ways. The porosity mainly affects parameter c, in that an increase in porosity leads to a decrease in c. However, the binder volume fraction mostly affects parameter s: An increase in binder volume fraction gives rise to a decrease in s. Therefore, the EOS of PBXs depends on both microstructures, such as porosity and mass ratio, of the HMX/binder.

Equation of State (EOS) Study
The isothermal compression simulations show that porosity and binder volume fraction have significant effects on the PV isotherms and Hugoniot data. Based on the PV variations obtained from the MPM mesoscale simulations, an effective EOS that can be used to study the PV isothermal properties of the PBX as well as porosity and binder volume fraction effective factors is obtained in the following procedures.
Normalized compression pressure P n can be obtained as follows: where P is the compression pressure of a system and P η is the pressure of a system when the compression ratio η = V/V 0 = 0.80. The PV isotherms of the PBX at different porosities and different binder volume fractions are normalized and plotted in Figure 8a,b, respectively. The normalized PV curves at different porosities are basically coincident, and the same can be said for different binder volume fractions. This allows the EOS fitting to be performed at any one of the four porosities and binder volume fractions.

Equation of State (EOS) Study
The isothermal compression simulations show that porosity and binder volume fraction have significant effects on the PV isotherms and Hugoniot data. Based on the PV variations obtained from the MPM mesoscale simulations, an effective EOS that can be used to study the PV isothermal properties of the PBX as well as porosity and binder volume fraction effective factors is obtained in the following procedures.
Normalized compression pressure can be obtained as follows: where P is the compression pressure of a system and is the pressure of a system when the compression ratio = / 0 = 0.80. The PV isotherms of the PBX at different porosities and different binder volume fractions are normalized and plotted in Figure 8a,b, respectively. The normalized PV curves at different porosities are basically coincident, and the same can be said for different binder volume fractions. This allows the EOS fitting to be performed at any one of the four porosities and binder volume fractions.
where 0 and 0 ′ respectively refer to the zero-pressure isothermal bulk modulus and its pressure , 0 is zero-pressure volume, and V is the volume of the PBX under compression loading. Based on Equation (18), the zero-pressure bulk modulus changes with the normalized pressure. One set of the normalized PV data is fitted to Equation (18)   The third-order Birch-Murnaghan EOS [54] is used to analyze the compression curves for the PBX. The equation is as follows: where K 0 and K 0 respectively refer to the zero-pressure isothermal bulk modulus and its pressure , V 0 is zero-pressure volume, and V is the volume of the PBX under compression loading.
Based on Equation (18), the zero-pressure bulk modulus changes with the normalized pressure. One set of the normalized PV data is fitted to Equation (18) as shown in Figure 8. The parameters are found to be K 0,n = 0.8765 and K 0,n = 23.27, where K 0,n is the normalized zero-pressure bulk modulus and K 0,n is the corresponding pressure derivative.
When the compression ratio η = 0.80, the relation between compression pressure and porosity is represented as Equation (19); also shown in Equation (20) is the relation between compression pressure and binder volume fraction.
where ϕ 0 and b 0 are the reference porosity and reference binder volume fraction, taking the values of 0.003and 5%, respectively. P 0,ϕ is the pressure at the reference porosity, and α refers to porosity effective factor. P 0,b is the pressure at the reference binder volume fraction, and β refers to the binder volume fraction effective factor. Each porosity or binder volume fraction has its corresponding reference pressure P η , as listed in Tables 4 and 5 for the four porosities and four binder volume fractions. P 0,ϕ = 9.668 GPa and α = 0.0168; these values are obtained by fitting the data in Table 4 into Equation (19). The values of P 0,b = 10.02 GPa and β = 0.0724 are obtained by fitting the data in Table 5 into Equation (20). Based on Equation (17), the relations among the parameters K 0 , K 0,n , and P η can be expressed as: Therefore, The parameters in the Equations (18) to (23) obtained by fitting the corresponding variables with porosity and binder volume fraction effective factors are given in Table 6. Using Table 6, the zero-pressure isothermal bulk modulus and its pressure derivative with porosity, K 0,ϕ and K 0,ϕ , are as follows: The zero-pressure isothermal bulk modulus and its pressure derivative with binder volume fraction, K 0,b and K 0,b , are as follows: where ϕ ϕ 0 and b b 0 are respectively the specific values of porosity and the binder volume fraction. Similar relationships between effective factors and the parameters c and s in the u s -u p plane of the PBX can be expressed as: where ϕ 0 and b 0 are the reference porosity and reference binder volume fraction, taking the values of 0.003 and 5% respectively; s 0,ϕ and c 0,ϕ are the parameters fitting to the u s -u p curve at the reference porosity; m ϕ and n ϕ refer to porosity effective factors. Here, s 0,b and c 0,b are the parameters fitting to the u s -u p curve at the reference binder volume fraction, and m b and n b refer to the porosity effective factors. The parameters obtained by fitting the data in Tables 2 and 3 to Equations (26) and (27) are given in Table 7. It is seen in Table 7 that n ϕ > n b , indicating that porosity has a greater influence on parameter c, i.e., the bulk wave speed of PBX, than the binder volume fraction. Here, c plays an important role in shock wave propagation. Comparisons with the parameters obtained from the effective factors of porosity and binder volume fraction suggest that porosity has a greater influence on Hugoniot data. Porosity is related to pores in a compressed system, and the authors in [55] found that density heterogeneity is always located at the parts where pores exist. These locations may become the failure zones for the existence of microcracks.
The pressure-volume data were used to determine certain thermodynamic parameters (bulk modulus, first pressure derivative of the bulk modulus) of the Birch-Murnaghan EOS, which can demonstrate the suitability of an EOS at different pressures, especially at high pressure. Hugoniot data derived from isothermal curves will provide quantitative measurements for shock parameters. Therefore, the effective EOS fitted in this research can be used to predict some thermodynamic and mechanical properties, such as the initiation process of the PBX under extreme conditions.

Bulk Modulus and Its First-Order Pressure Derivatives
The bulk modulus is given as K = −V(dP/dV) for different pressure, so the expression for bulk modulus from Equation (18) comes out as: The corresponding expression for the first-order pressure derivative of the bulk modulus K = dK/dP obtained from (28) comes out as: where K 0 and K 0 are respectively the zero-pressure bulk modulus and its pressure derivative. Variations of the isothermal bulk modulus K and its pressure derivative K with pressure for the PBX with 0.003 porosity and 5% binder volume fraction by using the Birch-Murnaghan EOS are shown in Figure 9a,b, respectively. Figure 9 shows that with the increase in pressure, K increases continuously, while K decreases progressively with the increase of pressure and gradually becomes asymptotic. Stacey et al. [56] suggested a few basic criteria of an EOS for its validity and applicability. These are as follows: i.
In the limit of infinite pressure, V/V 0 → 0. ii.
With the increase in pressure, the isothermal bulk modulus increases continuously, and in the limit of infinite pressure, K → ∞ iii.
K must decrease progressively with the increase in pressure, and K remains greater than 5/3 in the limit of infinite pressure.
The results of this study show that the effective EOS derived from mesoscale simulations satisfies the above three criteria based on thermodynamic relations. ii.
With the increase in pressure, the isothermal bulk modulus increases continuously, and in the limit of infinite pressure, → ∞ iii.
′ must decrease progressively with the increase in pressure, and ′ remains greater than 5/3 in the limit of infinite pressure.
The results of this study show that the effective EOS derived from mesoscale simulations satisfies the above three criteria based on thermodynamic relations.

Conclusions
Starting with the constitutive equations of each component of HMX/Estane PBX, an effective EOS of the PBX with microstructure effects was established by fitting the PV isotherms at different porosities and different binder volume fractions obtained by an MPM calculation. The equation of state is effective because it includes the effective factors of the microstructures of the PBX studied and the mass ratio of the energetic crystals and the binder contained in the PBX. The isothermal bulk modulus (K0) and its pressure derivative ( 0 ′ ) for the PBX were determined by analyzing the obtained isotherms using the Birch-Murnaghan equation of state. The results demonstrate that the

Conclusions
Starting with the constitutive equations of each component of HMX/Estane PBX, an effective EOS of the PBX with microstructure effects was established by fitting the PV isotherms at different porosities and different binder volume fractions obtained by an MPM calculation. The equation of state is effective because it includes the effective factors of the microstructures of the PBX studied and the mass ratio of the energetic crystals and the binder contained in the PBX. The isothermal bulk modulus (K 0 ) and its pressure derivative (K 0 ) for the PBX were determined by analyzing the obtained isotherms using the Birch-Murnaghan equation of state. The results demonstrate that the established effective EOS meets the criteria based on the basic thermodynamic discussion. The bulk wave speed c and dimensionless parameter s in the Mie-Grüneisen EOS of the PBX obtained by fitting u s -u p variations with the different porosities and different binder volume fractions are also effective.
The MPM mesoscale simulation technique presented in this study is an alternative approach for the development of an equation of state that can be used to predict the thermodynamic properties of PBXs for a wide range of pressures and temperatures. An effective EOS fitted at the different porosities and binder volume fractions in this study provides necessary input data for the description of particle dynamic behaviors at the continuum level, which is a key intermedium to connect microscopic and macroscopic properties in the study of energetic materials.
Author Contributions: Data curation, S.Y.; investigation, S.G., W.Z., and J.S.; methodology, S.Y. and Y.D.; writing-original draft, S.G.; writing-review and editing, G.V.L. All authors have read and agreed to the published version of the manuscript.