A viscoelastic–viscoplastic constitutive model for polymer bonded explosives under low impact loading

Viscoplastic work is very important to explosive ignition under impact loading. At present, a large number of constitutive models only consider the viscoelastic and damage behavior of explosives, ignoring the plastic effect under low impact loading. A new viscoelastic–viscoplastic (VE–VP) model was developed and studied to describe the dynamic mechanical behaviors of polymer-bonded explosives (PBXs). The total strain was assumed to be the sum of the viscoelastic (VE) and viscoplastic (VP) components. A generalized Maxwell model was used to determine the VE responses. A VP model was developed by using the classical J2 rate-dependent model with isotropic hardening. Viscoplastic flow was considered in hyperbolic sinusoidal form. The explicit algorithms of VE model were proposed and assessed by using two different integration methods. The accuracy and efficiency of these two methods are similar at high strain rates. The coupled algorithms of VE–VP model were developed by referring to the classical elasto-viscoplasticity (EVP) provided and using the expression of incremental relaxation modulus. The proposed model was implemented in the ABAQUS using a user-subroutine (VUMAT) to predict the response behaviors of PBX 9501 under low impact loading. Several numerical simulations illustrated the computational efficiency and the accuracy of the proposed methods. The model predictions were compared with experimental data, and reasonable agreement was obtained.

Polymer-bonded explosives (PBXs) comprise energetic crystals that adhere to a polymer binder and can be used in numerous applications, including propellant grains and warheads 1 . Understanding the mechanical behavior and establishing a constitutive model of PBXs are of great interest for defense and commercial applications, and predictive constitutive model descriptions that can be used in finite element simulations of damage and deformation are in demand, with emphasis on understanding the dynamics of the localization phenomena and mechanical failure of PBXs 2 .
The violent reaction induced by low velocity impact (about tens to hundreds meter per second) is not the result of direct shock wave loading, since the pressure generated by the impact is relatively low  and the duration time of the pressure is relatively long (50 μs-10 ms). The response of solid explosives under low velocity impact includes several complex processes, such as the localization of plastic work, the conversion from plastic work to heat, the heat conduction, the ignition and the combustion to detonation transition. The split Hopkinson pressure bar (SHPB) test is the simplest and most commonly used method of assessing the sensitivity of low-velocity impacted explosives.
Up to now, many constitutive models have been proposed and applied to study the dynamic mechanical behavior of PBXs 1,3,4 . Dienes et al. 5 proposed a statistical crack mechanics model (SCRAM) considering anisotropic damage. Base on the SCRAM, Bennett et al. 6 and Hackett et al. 7 developed a viscoelastic cracking constitutive model (Visco-SCRAM) to simulate the combined VE and brittle damage behavior of PBX 9501. Liu et al. 8 developed a VE-VP damage model to describe the stress-strain responses and damage behaviors of PBX based on Bodner-Partom viscoplasticity. Liu et al. 9 developed the constitutive model which incorporates a porosity model and a VE-VP damage model to study the strain rate and pressure-related deformation and damage of PBXs. Le et al. 10  Viscoelasticity. The VE response is expressed in a hereditary integral form, which is Boltzmann's form if the VE model is linear and usually the Schapery form for nonlinear VE models [19][20][21] . Based on Boltzmann's superposition principle 22 , the integral form of the VE stress is as follows: For an isotropic material, the fourth-rank relaxation tensor is written as follows: where G(t) and K(t) represent the rate dependent shear and bulk relaxation functions, respectively. A Prony series representation was used to model the shear and bulk relaxation function. G(t) and K(t) can be expressed as a Prony series of shear elements and bulk elements: where G (n) and K (n) represent the shear and bulk stiffness, respectively, τ (n) is the shear and bulk relaxation times, G (∞) and K (∞) are the long-term elastic shear and bulk moduli, respectively, and N is the number of Maxwell elements.
(1)  (3), the Cauchy stress σ ij can be rewritten as follows: The deviatoric ( S ij (t) ) and hydrostatic ( σ V (t) ) parts of the stress tensor ( σ ij (t) ) may be expressed as the sum of the elastic and the Maxwell components: The elastic element is given as follows: Each Maxwell element is given as follows: Viscoplasticity. The viscous stress ( σ v ) considered in the subsequent formula derivation is a power-law function, which is defined as follows 23 : where k and m are material constants. The constant m is called the material strain rate sensitivity. The stress is given as follows 24 : where σ y is the initial yield stress, r p is the hardening stress, and p is the effective plastic strain. In the VP model, the uniaxial stress depends on the yield stress, the hardening of the yield stress, and the plastic strain rate.
According to the Eq. (11), the plastic strain rate is For a von Mises material, the uniaxial stress is identical to the effective stress and similar for the effective plastic strain rate. Equation (12) can therefore be written as follows: According of the normality hypothesis of plasticity, we also write the normality hypothesis for viscoplasticity as follows: For the isotropic hardening, the von Mises yield function is given as follows 25 : where σ eq is the von Mises equivalent stress.
Equation (14) can be rewritten as .
σ = σ y + r(p) + σ v = σ y + r(p) + kṗ m , (15) f (σ, R) = σ eq − σ y − R(p) and σ eq = 3 2 S : S, www.nature.com/scientificreports/ The effective plastic strain is an internal variable that keeps track of the past history of the material and is linked to the VP law as follows: The effective plastic strain is defined as The effective plastic strain rate is obtained by By substituting Eq. (19) into Eq. (16), the VP strain rate can be obtained as The VP multiplier ṗ is determined by the following conditions 23 : where a and β represent the VP modulus and exponent, respectively.
The hardening function r(p) is a linear law which is defined as 26 : where h is the hardening modulus.  V (t u+1 ) at time t u+1 in the VE model can be obtained as follows 27 : Two approaches are considered to compute the integrals over [t u , t u+1 ]. If �ε ij =ε ij (t)�t , integration of Eq. (24) leads to the following: Based on the mid-rectangle rule 28 , the numerical computation of the integrals in Eq. (24) leads to the following: Equations (25) and (26) can be rearranged as follows:  (7) can be written in the following incremental form: The incremental relaxation moduli G and K are defined as follows 26,29 : The influence of the choice between the two definitions of G and K is assessed in the following work.

Coupled VE-VP model computational algorithm. It is based on a combination of techniques used
separately for VE and elastic-viscoplastic (EVP) models. If the trial stress does not exceed the yield stress, the new stress is set equal to the trial stress. The trial stress can be written as follows: The new stress at t u+1 can be given by: If the yield stress is exceeded, viscoplasticity occurs in the increment. We then write the incremental analogs of the rate equations as 29 where L ijkl (t) =G(t) δ ij δ kl + δ il δ jk + K (t) − 2 3G (t) δ ij δ kl represents the incremental relaxation moduli. The relation between Cauchy and trial stresses remains form-identical to EVP, except that the predictor is VE and the relaxation modulus is a function of the time increment.
In the VP model, the effective plastic strain rate is This implies that Here, the VP strain tensor is written as follows: Equation (34) can be written as follows: Using the integrated flow rule, together with the Mises definition of the flow direction n ij , this becomes This implies that By Eqs. (37) and (41) if f is positive, time integration of the VE-VP law simply implies that an iterative solution of two scalar equations can be obtained: The two unknows, σ eq t n+1 and p , can be obtained by solving these equations. We solve these by the Newton-Raphson method. We have the following two scalar equations: At each iteration (r), the equations can be written as follows: We obtain: where ∂ψ α ∂�p = 1 − ∂φ ∂�p �t = 1 − φ �p �t , and ∂ψ α ∂σ eq = − ∂φ ∂σ eq �t = −φ σ �t. Substituting Eq. (45) the second into the first gives the following: We obtain: At iteration (r + 1) step, the p is obtained as Computation framework. The numerical algorithm was implemented in the finite element program ABAQUS by coding the user material subroutine VUMAT. The flow chart is shown in Fig. 1.
At the beginning of each time step, the strain increment ( �ε ij = �ε ve ij , �ε vp ij = 0 ) was imported. The trial stresses were calculated by using Eqs. (31) and (32). Then the yield function was calculated by Eq. (15). If f < 0, the stress state was within the yield stress space and yield criterion is not met; then the trial results were real (35) p = φ σ eq , �p = a sinh β 3 2 S : S − r(p) − σ y .

Model validation
Material parameters for PBX 9501. In the VE-VP model, the five Maxwell components have been used to describe the dynamic VE behavior. VE material parameters are obtained in two ways. By means of the timetemperature superposition principle, the master curve of E(t) versus t was constructed using relaxation tests at different temperatures. The material parameters of PBX 9501 are obtained by fitting the modulus master curve. Figure 2 shows the master Young's modulus curve of PBX 9501. The relaxation times for a GMM were obtained from the relationship that the relaxation times were approximately one-tenth the reciprocal strain rate of the test values. Figure 2 shows the experimental data for the relaxation time and the fit to the data for PBX 9501. The material parameters for PBX 9501 associated with the VE model are given in Table 2.

Obtain viscoelastic and viscoplastic parameters
Compute the trial stress tensor   Input: e ve ij , �ε ve V , t ; Params: G (n) , K (n) , g (n) , k (n) , σ y , h, β , and a Output: ε vp , ε ve , S ij (t u+1 ) , σ V (t u+1 ) 1: Read the strain tensor at the (u + 1)th step 2: Compute the trial stress tensor using Eqs. (25) and (26) 3: Compute the incremental relaxation modulus: 3: Compute the von Mises yield function f: 4: If f ≥ 0 , compute the plastic strain rate; otherwise, go to step (1) 5: Use New iteration to compute the effective plastic strain increment as follows: �p (r+1) = �p (r) + d�p 6: Compute the nominal VP strain tensor as follows: Compute the nominal stress tensor: Figure 2. Young's modulus versus relaxation time data for PBX 9501. www.nature.com/scientificreports/ An inverse method is proposed for parameter identification in VP deformation. The inverse solution procedure consists of an optimization method allowing the adjustment of the parameters so that the calculated response matched that measured in a uniaxial compressive test.
We consider a monotonic uniaxial compression test such that σ 11 is the only non-zero component. The von Mises equivalent stress is defined as follows: For the hyperbolic plastic flow, the VP function given by Eq. (50) is as follows: Finally, using kinematic hardening, the effective stress is written as follows: where σ exp is the experimentally measured stress under the quasi-static compressive tests (Fig. 3).
In the VP formulation, four parameters need to be obtained: σ y , h, a and β. The initial yield stress σ y is considered to be rate independent, and σ y was obtained under strain rate of 1/s. The optimization minimizes an objective function that is formed from the square of the difference between the measured and simulated stress responses for the considered uniaxial compressive test. Figure 3 shows the stress-strain curve under strain rate of 1/s. The material parameters for PBX 9501 associated with the VP model is given in Table 3.
(50) p = φ �p, r = a sinh β σ eq − r − σ y .  Figure 3. The stress-strain curves for PBX 9501 under the quasi-static compressive tests.  30 . To verify the simulation code and the effectiveness of the material constitutive law, we first performed simulations with the developed material laws and the existing linear VE model in ABAQUS to obtain the dynamic compressive responses of the PBX 9501. The SHPB system consists of the striker, an incident bar, a transmission bar and an absorption bar. All bars were made of aluminum, of which the Young's modulus and density were 73 GPa and 2.7 × 10 3 kg/m 3 , respectively, and the diameters of all the bars were 9.18 mm. The length of the impact bar was 75 mm, and the length of the incident and transmission bars were 1000 mm. The diameter of the specimen was 6.35 mm, and the thickness was 3.15 mm. The SHPB system is partitioned by hexahedron elements, and the mesh size of the specimen is 0.7 mm. The mesh size of the all bar is 1.2 mm, as shown in Fig. 4. Contact constraints are imposed to prevent interpenetration of the system components. The material parameters were obtained from Table 2, where the incremental time was 8.85 × 10 −8 s. Incident, reflected, and transmitted waves in the bars are shown in Fig. 5a. We compared the stress-strain curves in the compression direction calculated by the developed material models and the ABAQUS model. The stress-strain curves of the three methods are shown in Fig. 5b. The results obtained with the proposed VE constitutive law agreed fairly well with the results calculated using the linear VE material law in ABAQUS. Figure 6 shows a comparison of the numerical time integration schemes. The time of the incident wave is 50 µs. By decreasing the time increment, we note that the two definitions of the incremental relaxation modulus lead to very similar and rather accurate results. Due to small Δt in higher strain rate, both integration methods can be used to predict the dynamic mechanical behavior of materials. An important consequence is that if the total strain rate is constant and VP strains do not evolve, then the Cauchy stress computed with the first integration method is insensitive to the value of Δt.  Tables 2 and 3. As shown in Fig. 7, the stress-strain curves of the PBX 9501 at a strain rate of 2200/s were calculated for the VE, EVP, and VE-VP models. The comparison of the experimental results and the simulation results in Fig. 7 shows that the VP effect caused the VE and plastic stages to be smoothly connected. Compared with the simulation results of the three models, VE-VP can better describe the dynamic mechanical properties of PBX 9501. The VE-VP and EVP responses were very similar in the VP stage. As expected, the initial slopes of the VE-VP and VE curves coincided in the VE stage.

Comparison of the VE/EVP/VE
Uniaxial compression at various strain rates. The strain rate-related behavior of the presented model is examined by uniaxial dynamic compression. At about 194 μs, the stress wave reaches the left end face of the specimen, which was the same as the transmitted time calculated using the elastic wave velocity ( c 0 = E 0 ρ 0 = 5200 m/s, where E 0 = 73 GPa and ρ 0 = 2700 kg/m 3 are Young's modulus and mass density of the bar). Von Mises stress was nearly zero on the specimen, as shown in the Fig. 8a. The unit of von Mises stress cloud is MPa.
Part of incident wave is reflected back to the incident bar, and the remainder passes through the specimen to the transmission bar (the transmitted wave). To characterize the stress equilibrium, the relative stress difference was calculated as α = �σ/σ T , where σ T is the mean Von Mises stress. The assumption of stress equilibrium was nearly satisfied if α ≤ 5% . The specimen achieves stress equilibrium at about 203 μs, as shown in Fig. 8b. Von Mises stress difference at each point does not exceed 0.7 MPa. The relative stress difference was < 5%. The simulation model determined the stress equilibrium.
In fact, the stress-strain curve simulated by finite element model (FEM) is in good agreement with the experimental results under high strain rate. A point in the middle of the specimen surface was taken out, and

Incident bar T ransmitted bar Specimen Striker
Drawing of partial enlargement    Fig. 9. The uniaxial compressive stress-strain curves of PBX 9501 are compared by using the VE-VP model in Fig. 10. The strain rates in the loading direction were specified to be 700, 1500, and 2200/s for compression, respectively. The effect of rate on the material behavior are evident in the Figure. As shown in Fig. 10, the current model captures the main features of the experimental results, such as the initial slope, curvature and peak of the curve, quite well. As the plastic strain continued to increase, the plastic flow stress at different strain rates were parallel to each other, because the modulus of the hardening stage had the same value. Thus, this model could well describe the rate-dependent characteristics of PBX 9501.
Analysis and discussion. The algorithm established in this paper can accurately and effectively simulate the dynamic mechanical properties of PBX9501. Compared with the methods reported in the existing literature 14,[31][32][33][34] , there are two important differences: (I) For the discretization of the VE integral, an algorithm for the explicit update of a GMM was developed by using two schemes: constant VE strain rate, and mid-point rule. Due to small Δt under impact loading, the methods are exact that the VE strain rate is constant over the step. Viscoplastic flow was considered in hyperbolic sinusoidal form. (II) The coupled algorithms of VE-VP model were developed by using the expression of incremental relaxation modulus. The methods are remarkably simple and form-identical to EVP provided that the constant elastic stiffness operator is replaced with L ijkl (t) Eq. (31). (III) The display algorithm avoids the convergence in the calculation process, but this algorithm is satisfied by both integration methods for small Δt. This algorithm is very suitable for predicting the mechanical properties of materials under impact loading. Ellsiepen et al. 35 studied a class of diagonally implicit Runge-Kutta (DIRK) methods. This algorithm is much more advanced than the one presented in this paper. In a future work, it will be very interesting to develop a DIRK method instead, in order to achieve better accuracy and optimal time step control.

Conclusions
The present work aimed at the coupled VE-VP modeling of the dynamic compressive behavior of PBXs. The total strain was assumed to be the sum of VE and VP components. The stress was related to the history of the VE strains by a GMM. A VP model was developed for the history of the VP strains. This study presented numerical and analytical results. An algorithm for the explicit update of a GMM was developed by using two schemes: constant VE strain rate, and mid-point rule. Due to the time of loading was small, using the algorithm to predict the mechanical behavior was simple with high efficiency. It is worth to be recommended. If the Δt is not small, the algorithm can only be using under constant strain rate. The expression of incremental relaxation modulus was deduced by using the explicit algorithm proposed. The explicit algorithm of VE-VP model was developed by using the incremental relaxation modulus. The expression of VE-VP model was form-identical to classical EVP provided.
The VE-VP model was implemented in ABAQUS via a VUMAT. The VE, EVP, and VE-VP responses were also shown for comparison. The VE-VP model was more suitable for dynamic mechanical properties of PBX 9501. A SHPB numerical simulation was presented for the dynamic compressive responses of PBX 9501. The stress-strain curves of PBX 9501 were predicted by using the VE-VP model under different higher strain rates. Good convergence was achieved.

Data availability
All data generated or analyzed during this study are included in this published article.