Meso-scale failure simulation of polymer bonded explosive with initial defects by the numerical manifold method

Computational Materials Science xxx (xxxx) xxxx

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Meso-scale failure simulation of polymer bonded explosive with initial defects by the numerical manifold method ⁎

Ge Kanga, Youjun Ningb, , Pengwan Chena, a b

⁎

State Key Laboratory of Explosion Science and Technology, Beijing Institute of Technology, Beijing 100081, China School of Mechatronic Engineering, Southwest Petroleum University, Chengdu 610500, Sichuan, China

A R T I C LE I N FO

A B S T R A C T

Keywords: Polymer bonded explosives (PBXs) Initial defects Constitutive models Cohesive contact Numerical manifold method (NMM)

Polymer bonded explosive (PBX) is a composite consisting of the polymer binder and embedded explosive particles, along with a large number of the particle/binder interfaces as the third constituent. The particle volume fraction is often higher than 90%. In the present work, a visco-elastic constitutive model, an elastic viscoplastic constitutive model and a bilinear cohesive contact relationship model are implemented into the numerical manifold method (NMM) program, an open source code programmed with C language, to describe the deformations of the polymer binder, the explosive particles and the particle/binder interfaces, respectively. The fracturing of the polymer binder and explosive particles is described based on the maximum tensile stress and the Mohr-Coulomb criteria. Three categories of initial defects, including the initial interfacial debonding, initial voids in the polymer binder, and initial micro-cracks in the explosive particles, are considered in the PBX mesostructures under both uniaxial tensile and compressive conditions. The tension-compression asymmetry, the influence of the initial defects on the meso-structure failure modes/patterns and the macroscopic effective tensile/compressive strength of PBXs are investigated. The factors that cause the differences between the NMM results and other numerical or experimental results are analyzed and discussed. This work enables and proves the NMM to be an robust numerical tool for further simulation studies of the mechanical performances of PBXs, as well as other particle-filled composites, at the meso-scale.

1. Introduction Polymer bonded explosive (PBX) is a kind of energetic material mainly consisting of the polymer binder and filled explosive particles along with the particle/binder interfaces. Taking the HMX-based PBX9501 as an example, the volume fraction of HMX particles is about 92–93%, while the content of Estane (a kind of polymer binder) is only about 7–8% [1]. Due to the extremely lower elastic modulus of the binder [2], it allows PBXs to deform largely and absorb most of the deformation energy. Meanwhile, due to the high rigidity and fracture strength of the particle, failures/damages in PBXs are easier to evolve along the weaker particle/binder interfaces. Therefore, besides of the polymer binder and explosive particles, the particle/binder interfaces also have significant effects on the mechanical response of the PBXs [3]. It is a big challenge for the stability and safety of PBX because of its diverse and complex situated environments. In the past years, many researchers have been focusing on experimental characterization of PBXs to investigate the deformation, fracture behavior and fracture mechanism of the PBX at the meso-scale [4–8]. Some researches

⁎

devoted to the mechanical properties of the particle/binder interfaces and developed different micro-scale cohesive interfacial laws [3,9,10]. With the development of SEM and MCT Technologies, the formation reason, distribution pattern and evolution process of some main intrinsic defect forms (voids and micro-cracks) at the meso-scale could also be observed [11,12], especially the particle/binder delamination and flow of the particles [13]. In addition, in recent years, a lot of numerical simulation works were also conducted to investigate the meso-scale mechanical properties of composites. In the respect of PBX, different constitutive models were proposed for the polymer binder and explosive particles [14,15,16], respectively. The constitutive models of damaged interface, which is considered as a mixed binder-void interphase layer [17,18], were also proposed. Besides, some new numerical frameworks were proposed specially to investigate the deformation and damage sensing capabilities of nano-composite bonded explosives (NCBX) [19] and PBX [20]. Moreover, many researchers also had done a lot of numerical works to study effective modulus of PBX [21–26]. In the respect of other composites, a PERMIX software framework [27] and a methodology for stochastic modelling of the fracture [28] were

Corresponding authors. E-mail addresses: [email protected] (Y. Ning), [email protected] (P. Chen).

https://doi.org/10.1016/j.commatsci.2019.109425 Received 23 July 2019; Received in revised form 8 November 2019; Accepted 18 November 2019 0927-0256/ © 2019 Elsevier B.V. All rights reserved.

Please cite this article as: Ge Kang, Youjun Ning and Pengwan Chen, Computational Materials Science, https://doi.org/10.1016/j.commatsci.2019.109425

Computational Materials Science xxx (xxxx) xxxx

G. Kang, et al.

proposed for polymer/particle nanocomposites. At the meso-scale, the SEM results of PBX show that it always contains a large proportion of initial defects, such as interfacial debonding, voild in the binder and micro-cracks in the particles, and they will influence the mechanical performances of PBXs to a large extent [11,12]. The large number of displacement jumps and stress discontinuities across the debonded interfaces and micro-cracks in the PBX meso-structure bring challenges to continuum-based numerical methods, such as the FEM. In the present work, the NMM, a continuousdiscontinuous numerical method, is further developed on the basis of our previous work [29] to simulate the deformation and failure of PBXs at the meso-scale, by mainly considering the effects of initial defects. A visco-elastic constitutive model and an elastic visco-plastic constitutive model are implemented for the polymer binder and explosive particles, respectively, and a bilinear cohesive contact relationship (BCCR) model is implemented for the particle/binder interfaces. The fracturing of the polymer binder and explosive particles bases on the maximum tensile stress and the Mohr-Coulomb criterion [30,31]. The failure process of the PBX meso-structures is simulated under both uniaxial tensile and compressive conditions. Besides, the volume fraction of explosive particles in this present work is 90% more or less, much higher than that of other work [2].

2. Basic concepts of the NMM 2.1. Concepts of MC and PC The numerical manifold method (NMM) is a unified continuousdiscontinuous numerical method firstly proposed by Shi [32]. This method employs an implicit solving scheme and can be adopted to simulate the deformation and failure behaviors of solids in both static/ dynamic conditions. Especially, two cover systems, i.e., the mathematical cover (MC) and the physical cover (PC), are two important concepts in the NMM. For the example, in Fig. 1, the MC system is composed of regular triangles which completely cover the whole physical region (Ω). For each node of the MC system, the hexagonal piece, which is composed of six triangles, such as Mi, Mj, Mk and Ml at nodes “i”, “j”, “k”, and “l”, is called a mathematical cover (MC). The overlapping part of each MC and the physical region is called a physical cover (PC), such as Pi, Pj, and Pk (case (i), (ii), (iii) in Fig. 1). Besides, if a MC is completely divided into two or more totally separated parts by cracks, each part will form a PC, such as Pl1 and Pl2 derived from Ml (case (iv) in Fig. 1). Theoretically, in the NMM framework, arbitrary shapes of MCs are applicable [30]. However, for the convenience of the cover system generation, definition of the weight function, and integration of the

Fig. 2. Illustration of ME in the NMM (“□” denotes an MC completely split by a crack, while “○” denotes an MC isn't or not completely split by a crack).

weak form, at present, the triangular mesh topology, as shown in Fig. 1, is most adopted. Based on this kind of topology, the common area of three, and must only three PCs is defined as a manifold element (ME), which is also an important concept in the NMM.

2.2. Concept of ME The concept of ME in the NMM is illustrated in Fig. 2. For the triangle e1, its three nodes are i, j and k, respectively. Due to the existence of the crack, Mj and Mk are completely divided into Pj1 and Pj2, Pk1 and Pk2, respectively, while Mi is not completely divided, as shown in Fig. 2(b). Therefore, for ME e11, its three PCs are Pi, Pj1 and Pk1; while for ME e12, its three PCs are Pi, Pj2 and Pk2. Similarly, in Fig. 2(c), Ml, Mm and Mn are all completely divided by the crack into two separate parts. Thus, three PCs of the ME e51 are Pl1, Pm1 and Pn1, and for ME e52, its three PCs are Pl2, Pm2 and Pn2. In addition, in Fig. 2(d), Mr, Ms and Mt are all not completely divided by the crack. Thus, the triangle e8 forms only one ME with three PCs of Pr, Ps and Pt. It can be observed that an

Fig. 1. Illustration of MC and PC in the NMM. 2

Computational Materials Science xxx (xxxx) xxxx

G. Kang, et al.

clockwise, which satisfies

ME may be a triangle or an irregular polygon, however, always determined by three and only three PCs. The three PCs for an ME are also defined as the nodes of this manifold element. Furthermore, taking ME e8 for instance, its displacement field can be expressed by the displacements of its three nodes as follows

u (x) =

∑ wh (x) uh (x),

(x ∈ e8, h = r , s, t )

1 x1 + u1 y1 + v1 Δ = 1 x2 + u2 y2 + v2 1 x3 + u3 y3 + v3

where uh(x) is the displacements of the node h (h = r, s, t); wh(x) is termed as the weight function, which satisfies

e51

(2)

2

However, for MEs and e5 across the crack surface, the displacement jump across can be expressed as

uc (x) =

∑ wh1 u1h − ∑ wh2 u2h,

(h = l, m , n)

c

(3) uh1,

where u (x) represents the displacement jump across the crack; uh2 1 2 and wh , wh (h = l, m, n) represent the local independent unknowns and weight functions, respectively. 2.3. Concept of contact

K Δu = Δf The contact also is an important concept in NMM. As shown in Fig. 3(a), assume two MEs locate along two sides of a crack or interfacial surface, three categories of contact exist, in which category I and III will be transferred into category II in contact treatment in NMM. Category I can be transferred into two “vertex-edge” contacts (p1 to p2p3 and p4 to p2p3), while category III only can be transferred into one (p1 to p2p3 or p2 to p1p3), as marked by red dashed circles in the figure. For the particle/binder interface and micro-crack simulations in the present study, most contacts are category I and they are transferred into category II in contact treatment. Fig. 3(b) shows how the invasion of a vertex-edge contact is judged and treated. At time t, the sequence of vertexes p1, p2, and p3 is anti-clockwise (as marked by black arrows in Fig. 3(b)), which satisfies

1 x1 y1 Δ = 1 x2 y2 1 x3 y3

(5)

where (ui, vi) is the displacement increment of vertex pi (i = 1, 2, 3) in Δt. Therefore, Δ < 0 is a sign of invasion. If invasion happens, a normal penalty spring is added to constrain the normal movement of vertex p1. Meanwhile, if it does not slide (according to the Mohr-Coulomb criterion) in shear direction, a shear penalty is added to constrain its shear movement simultaneously; otherwise, no shear spring will be adopted. The normal and shear invasion displacements can be calculated by projecting vector p0p1′ into the local n-o-s coordinate system, in which p0 is the projection point of p1 on segment p2p3. If the contact is open (Δ > 0), a normal spring can also be added to represent the tensile strength on crack or interfacial surfaces. Besides the penalty method described above, the Lagrange multiplier method, or the augmented Lagrange multiplier method, also can be employed to prevent the invasion between the two sides of a discontinuity. Details can be found in relevant references [33]. The corresponding sub-matrices of contact are added into the NMM system equilibrium equation as follows

(1)

0 ≤ wh (x) ≤ 1, ∀ x ∈ Mh wh (x) = 0, ∀ x ∉ Mh, (h = r , s, t ) ∑ wh (x) = 1

<0

(6)

where K is the global stiffness matrix, Δu is unknowns of displacement increment, and Δf is global loading matrix. Then, by solving the system equilibrium equation, the displacement field of a problem can be obtained. Different from that in FEM, here, Δu can include the displacement jumps across discontinuity surfaces due to the two cover systems in NMM. 3. Contact and material models for PBX simulations In practice, the categories of energetic particles in PBXs mainly include HMX, RDX, TATB, and PETN, et al. Among these, HMX is widely used and has drawn intensive interest for several decades because of its high energy density. For HMX-based PBX9501, it mainly contains 92% HMX particles and 8% Estane [1] in volume. Its particle/binder interface can be treated as the third constituent with cohesive mechanical properties [3]. Studies [2,14] also indicated that the HMX particle is a hydrostatic strain rate dependent elastic visco-plastic material, and the polymer binder, such as Estane in the PBX9501, Nitrocellulose in the PBX9404, is a visco-elastic material which is extremely sensitive to strain rates and ambient temperatures.

>0 (4)

where (xi, yi) is the coordinates of vertex pi (i = 1, 2, 3). At moment t + Δt, assume p1 moves to p1′, then the sequences of p1′, p2, p3 turns

3.1. Bilinear cohesive contact relationship model for particle/binder interfaces For the particle/binder interfaces, the bilinear cohesive contact relationship (BCCR) model is applied both on the normal and shear direction. In Fig. 4, d is the open or sliding contact displacement; F is the corresponding contact force; k0 and k1 are the spring stiffness of the ascending and descending stages, respectively; d0 is the contact displacement corresponding to the interface tensile or shear strength; dc is the critical open or sliding displacement. A damage variable D is introduced to describe the damage degree of the contact when d > d0, and the damaged spring stiffness k′ can be calculated as k′ = (1-D)k0. When d > d0, if unloading occurs, k′ will be used in the unloading process and further used in the subsequent loading stage until D gets a new value. The area of the triangle below the two solid lines is defined as the fracture energy GF. Tan et al. [3] used the extended Mori-Tanaka method and the equivalence of cohesive energy on the macro-scale and micro-scale, respectively, to link the macro-scale compact tension experiment with the micro-scale cohesive law for particle/matrix interfaces. The

Fig. 3. Illustration of contact technique in NMM. 3

Computational Materials Science xxx (xxxx) xxxx

G. Kang, et al.

(i = 1, 2, …, 22) is fitted with a series of measured storage modulus in reference [34]. Then, the Young's relaxation modulus Ei and the relaxation time τi at the room temperature (25 °C) can be calculated, respectively, as shown in Table 2. According to generalized Maxwell elements, the stress–strain behavior of the binder is determined as [34]

∫0

σ=

t

N

E (t − τ ) ε ̇ (τ ) dτ =

∑ Ei ∫0

t

e−

t−τ τi ε ̇ (τ ) dτ

(13)

i=1

Extend equation (10) to the two dimensional (2D) state: 22

σ=

∑ Ei ∫0

t

i=1

[A]

∂ε − t − τ e τi dτ ∂τ

(14)

in which σ is the Kirchhoff stress tensor, ε is the Green strain tensor. Here, for the plane stress or plane strain, [A] can be expressed as below, respectively: Fig. 4. Illustration of the BCCR model.

1 ν 0 ⎞ 1 ⎛ ν 1 0 , 2 ⎜ 1 − ν 0 0 (1 − ν )/2 ⎟ ⎝ ⎠ ν ⎛ 1 1−ν ν 1−ν ⎜ 1 [A] = (1 + ν )(1 − 2ν ) ⎜ 1 − ν ⎜ 0 0 ⎝

[A] = parameters for particle/binder interfaces they measured are listed in Table 1. For the BCCR model, its implementation in NMM and the verification can be found in reference [29]. 3.2. Visco-elastic constitutive model for the polymer binder (Estane)

N

∑ Gi e−t / τi

(7)

i=1

where τi is the relaxation time of the ith Maxwell element, Gi is ith shear relaxation modulus, and Ge is the long term modulus when the binder is fully relaxed. Here, N = 22 and Ge = 0. The frequency-modulus relations for different temperatures are shifted and superposed by using the Williams-Landell-Ferry (WLF) shift function

log(αT ) = −6.5

(T − T0) 120 + T − T0

(8)

in which T0 = 19 °C is called the reduced temperature [34], and αT is the time-temperature shift factor. Then, the relaxation time τi can be calculated by αT, as follows

τi = 1.5·αT(7 − i)

(i = 1, 2, …, 22)

⎞ ⎟ ⎟ 1 − 2ν ⎟ 2(1 − ν ) ⎠ 0

(15)

where ν is the Poisson’s ratio of the polymer binder. The implementation of the above visco-elastic constitutive model in NMM was repoeted in detail in reference [29]. In the present work, the elastic term of the visco-elastic constitutive model is neglected in the description of Estane in PBX9501. Here, to verify the implemented visco-elastic constitutive model in NMM for the polymer binder of Estane, a single element mode l as shown in Fig. 5 is simulated, The length and width of the element model are 1 m and 0.1 m, respectively. The left boundary is fixed. On the right boundary, a displacement loading for different strain rates (0.004933 s−1, 0.04933 s−1, 0.4995 s−1, 1.568 s−1, 14.95 s−1) of the element is applied. The strain rates are identical to that in reference [34]. The larger length/width ratio (10:1) of the element model guarantees its one-dimensional deformation to some extent. The comparisons of NMM simulation results and theoretical/experimental results [34] are shown in Fig. 6. It can be found that quite satisfactory agreements are achieved.

A visco-elastic constitutive model, which contains N generalized Maxwell elements, is introduced into the NMM. As described in reference [34], the shear modulus is expressed as a Prony series form

G (t ) = Ge +

0

3.3. Elastic visco-plastic constitutive model for the explosive particles (HMX)

(9)

For the case of uniaxial stress, the relevant Young's modulus is given by

E (t ) = 2·G (t )·(1 + v )

In one-dimensional case, the elastic visco-plastic constitutive model can be expressed as

(10)

σ̇ = E (ε ̇ − εṗ )

in which ν = 0.495 is the Poisson’s ratio of the polymer binder. Therefore, the Young's modulus can be approximately expressed as:

Here, E is the elastic modulus, and ε i̇ s the total strain rate, which can be decomposed into an elastic term and a visco-plastic term as

(11)

E (t ) = 2·G (t )·(1 + v ) ≈ 3·G (t )

ε ̇ = εė + εṗ

Suppose there exists a series of Young’s relaxation modulus Ei (i = 1, 2, … , 22), according to Eq. (11), the Young’s relaxation modulus is given by

(17)

According to the flow rule of the Mises materials [35], under the complicated stress state, the visco-plastic strain rateε̇p satisfies

(12)

Ei ≈ 3·Gi

(16)

εṗ =

According to the WLF shift function, the shear relaxation modulus Gi

3¯εṗ 2σ¯

S , with σ¯ =

3 S: S 2

(18)

Table 1 Parameters for particle/binder interface [3]. Normal spring stiffness (GPa)

Shear spring stiffness (GPa)

Tensile strength (MPa)

Shear strength (MPa)

Fracture energy (J/m)

300

120

1.66

1.66

81

4

Computational Materials Science xxx (xxxx) xxxx

G. Kang, et al.

Table 2 The shear relaxation modulus Gi, the Young's relaxation modulus Ei, and the relaxation time τi [14]. Shear relaxation Prony terms (/MPa)

G∞ 0 G6 0.0891 G12 2.618 G18 346.737

G1 0.00417 G7 0.0056 G13 12.882 G19 251.188

G2 0.00741 G8 0.1622 G14 52.481 G20 177.83

G3 0.00159 G9 0.2218 G15 223.872 G21 117.489

G4 0.0380 G10 0.4753 G16 436.516 G22 75.8

G5 0.0676 G11 2.104 G17 457.088

Normal relaxation Prony terms (/MPa)

E1 0.0125 E7 0.0168 E13 38.647 E19 753.56

E2 0.0222 E8 0.4866 E14 157.44 E20 533.48

E3 0.0475 E9 0.6654 E15 671.616 E21 352.46

E4 0.114 E10 1.426 E16 1309.54 E22 227.573

E5 0.203 E11 6.312 E17 1371.26

E6 0.2673 E12 7.854 E18 1040.21

Relaxation time Prony terms (/s)

τ1 7.355e5 τ7 7.355e−1 τ13 7.355e−7 τ19 7.355e−13

τ2 7.355e4 τ8 7.355e−2 τ14 7.355e−8 τ20 7.355e−14

τ3 7.355e3 τ9 7.355e−3 τ15 7.355e−9 τ21 7.355e−15

τ4 7.355e2 τ10 7.355e−4 τ16 7.355e−10 τ22 7.355e−16

τ5 7.355e1 τ11 7.355e−5 τ17 7.355e−11

τ6 7.355 τ12 7.355e−6 τ18 7.355e−12

materials, respectively. Actually, the most critical term in Eq. (19) is the equivalent viscoplastic strain rate ε¯ṗ . It mainly reflects the plastic flow of the material. In the present work, according to references [36,37], ε¯ṗ is expressed as Fig. 5. Single element model.

ε¯ṗ = in which ε¯ṗ and σ¯ are the equivalent visco-plastic strain rate and the equivalent stress, respectively, which both are scalar quantity, and S is the Kirchhoff stress tensor. Thus, Eq. (16) can be further extended to two-dimensional as

σ̇ = M : (ε ̇ − εṗ ) = M : ε ̇ − M :

3¯εṗ 2σ¯

m

σ¯ ε¯1̇ = ε¯0̇ ⎡ g(¯ε , T ) ⎤ ⎣ p ⎦ ε¯2̇ = ε¯ṁ exp[−ag(¯εp, T )/ σ¯ ]

(

g(¯εp, T ) = σ0 1 +

S

ε¯p N ε0

)

{1 − β ⎡⎣ ( )

T κ T0

− 1⎤ ⎦

}

(21)

(19) where ε¯p is the equivalent plastic strain, ε¯0̇ and ε¯ṁ are reference strain rate for a low regime of strain rate and a high regime of strain rate, respectively, and m and a are relevant strain rate sensitivity parameters. σ0 is the quasi-static yield stress, ε0 is a reference strain, N is the strain hardening exponent, T0 is a reference temperature, and β and κ are thermal softening parameters. The function g(¯εp, T ) represents the quasi-static stress-strain response at ambient temperature. Obviously, Eq. (18) has considered the strain hardening and strain-rate dependence of plasticity. When the calculation proceeds from step n to step n + 1, suppose the time increment Δt is very small, then, Eq. (19) can be further expressed as

where M is the matrix of elastic modulus. At plane stress or strain state, it can be expressed as followings, respectively

1 ν 0 ⎞ E ⎛ ν 1 0 , 1 − ν 2 ⎜ 0 0 (1 − ν )/2 ⎟ ⎝ ⎠ ν ⎛ 1 1−ν E (1 − ν ) ⎜ ν 1 M= (1 + ν )(1 − 2ν ) ⎜ 1 − ν ⎜ 0 0 ⎝

ε¯1̇ ·¯ε2̇ ε¯1̇ + ε¯2̇

M=

0

⎞ ⎟ ⎟ 1 − 2ν ⎟ 2(1 − ν ) ⎠ 0

(20)

in which E and ν are the elastic modulus and the Poisson's ratio of

Fig. 6. Comparisons of NMM simulation results with theoretical/experimental results under different strain rates at T = 25°. 5

Computational Materials Science xxx (xxxx) xxxx

G. Kang, et al.

Δσ = M : Δε − M :

3¯εṗ 2σ¯

Δt: S

parameters (E and ν) are deduced through 13 elastic constants for the βHMX monoclinic single crystals at room temperature in [14], while the plastic parameters are derived by inversion of constitutive equation in [36,37]. Fig. 7 shows the comparisons of NMM simulation results with theoretical results under different strain rates and different temperatures. It can be found that, each stress-strain curve will undergo three different stages in the loading process. In Stage I, due to the small value of the visco-plastic strain rate, the deformation is dominated by the elasticity. The curve is close to a straight line with a slope of the elastic modulus E. In Stage II, as the increase of the visco-plastic strain, the slope of the stress-strain curve decreases gradually towards horizontal. In Stage III, the stress-strain curve finally approaches to horizontal and the final stress level is close to the yield stress of the material (260 MPa). Moreover, in the second and third stage, a higher loading strain rate, or a lower temperature, means a higher stress level of the material response. It can be found that good agreement between the theoretical and NMM simulation results is achieved, even under large strain (> 20%). The implemented elastic visco-plastic constitutive model is verified and it represents the mechanical character of nonlinear elasticity, deformation hysteresis, strain rate and temperature dependence of visco-plastic materials well.

(22)

in which Δσ is the increment of the stress tensor, while S is the accumulated total stress tensor. It is written in the matrix form as

σ ⎛ Δσ11 ⎞ ⎛ Δε11 ⎞ 3 ε¯ṗ ⎛ 11 ⎞ ⎜⎜ Δσ22 ⎟⎟ = M ⎜⎜ Δε22 ⎟⎟ − M · 2 σ¯ Δt·⎜ σ22 ⎟ ⎝ σ12 ⎠ ⎝ Δσ12 ⎠ ⎝ Δε12 ⎠

(23)

Assuming matrix Q as

Q = −M ·

σ 3 ε¯ṗ ⎛ 11 ⎞ Δt·⎜ σ22 ⎟ 2 σ¯ ⎝ σ12 ⎠

(24)

Then Eq. (23) can be written as

⎛ Δσ11 ⎞ ⎛ Δε11 ⎞ ⎜⎜ Δσ22 ⎟⎟ = M ⎜⎜ Δε22 ⎟⎟ + Q ⎝ Δσ12 ⎠ ⎝ Δε12 ⎠

(25)

Now, the strain energy of a manifold element can be deduced as:

⎛ Δσ12 ⎞ 1 = ∬A 2 (Δε11 Δε22 Δε12 ) ⎜ Δσ12 ⎟ dxdy ⎜ ⎟ ⎝ Δσ12 ⎠ ⎫ ⎧ ⎧ Δε11 ⎫ 1 = ∬A 2 (Δε11 Δε22 Δε12 ) M Δε22 + Q dxdy ⎬ ⎬ ⎨ ⎨ ⎭ ⎩ ⎩ Δε12 ⎭ ∐e

3.4. Fracturing criterion The maximum tensile stress criterion and Mohr-Coulomb criterion have been extensively used in rock mechanics due to the simplicity and the capability to predict peak strength and failure directions effectively. Assuming there is an inclined plane with an angle β to the direction of the minor principal stress σ3, according to the Mohr-Coulomb criterion, the shear strength on this inclined plane is

⎧ Δε11 ⎫ 1 1 = ∬A 2 (Δε11 Δε22 Δε12 ) M Δε22 dxdy + ∬A 2 (Δε11 Δε22 Δε12 ) Qdxdy ⎬ ⎨ Δ ε ⎩ 12 ⎭ 1

1

= 2 [De]T Ae [Be]T M [Be][De] + 2 [De]T Ae [Be]T Q (26)

s = c + σn tan φ

in which Ae is the area of the manifold element, [De]T=(De(1) De(2) De(3)) is the displacement unknowns of three PCs, [Be]=(Be(1) Be(2) Be(3)) is a deformation matrix, and e(i) is the ith PC, i = 1, 2, 3. Therefore, the following stiffness matrix for a manifold element is added to the global stiffness matrix K:

where c is the cohesion, φ is the internal friction angle, and σn is the normal stress. On this inclined plane, the normal stress is

σn =

1 1 (σ1 + σ3) + (σ1 − σ3) cos 2β 2 2

(32)

and the shear stress is

T

⎛ Be (1) ⎞ Ae [Be]T M [Be] = Ae ⎜ BeT(2) ⎟ M ( Be (1) Be (2) Be (3) ) ⎜⎜ T ⎟⎟ B ⎝ e (3) ⎠

1 (σ1 − σ3) sin 2β 2

τ=

Ae [Be (r) ]T M [Be (s) ] → [K e (r ) e (s)], r , s = 1, 2, 3

Substituting Eq. (32) to Eq. (31), and letting s = τ, the limit stress condition on the inclined plane is calculated as

(28)

σ1 =

and the following loading matrix for a manifold element is added to the global loading matrix Δf:

(29)

(34)

φ π + 4 2

β=

That is

1 e A [Be (r ) ]T Q → {Fe (r ) }, r = 1, 2, 3 2

2c + σ1 (sin 2β + tan φ (1 − cos 2β )) sin 2β − tan φ (1 + cos 2β )

Fig. 8 shows the Mohr-Coulomb criterion combined with the maximum tensile stress criterion. There must have a critical plane, on which the available shear strength will be firstly reached when σ1 increases. Fig. 8(a) also gives the orientation of this critical plane as

T

⎛ Be (1) ⎞ 1 1 − Ae [Be]T Q = − Ae ⎜ BeT(2) ⎟ Q 2 2 ⎜ ⎜ B T ⎟⎟ ⎝ e (3) ⎠

(33)

(27)

That is

−

(31)

(35)

which can also be obtained by solving d(s-τ)/dβ = 0. On the critical plane, sin2β = cosφ and cos2β = −sinφ, then Eq. (34) reduces to

(30)

σ1 =

The single element model in Fig. 5 is used again to verify the elastic visco-plastic constitutive model for the HMX particles. The parameters for the constitutive model are listed in Table 3, in which the elastic

2c cos φ + σ3 (1 + sin φ) 1 + sin φ

(36)

The linear relationship between σ3 and σ1 is plotted in Fig. 8(b).

Table 3 Parameters for the elastic visco-plastic constitutive model [14,36,37]. E (GPa)

ν

ε¯0̇ (s−1)

m

ε¯ṁ (s−1)

a (1/MPa)

σ0 (MPa)

ε0 (s−1)

T0 (K)

N

β

κ

13.75

0.32

1e-4

100

8e12

22.5

260

5.88e-4

293

0.05

2.4

0.15

6

Computational Materials Science xxx (xxxx) xxxx

G. Kang, et al.

Fig. 7. Comparisons of NMM simulation results with theoretical results under different strain rates and different temperatures.

If the Mohr-Coulomb envelop is extrapolated to σ1 = 0, it will intersect the σ3 axis (T0′ in Fig. 8(b)) at an apparent value of the uniaxial tension strength. However, the experimentally measured uniaxial tensile strength is generally lower. Therefore, a tension cutoff is applied at a selected value of the uniaxial tensile stress, which derives the maximum tensile stress criterion as

- σ3 = T0

simulation. 4.1. The influence of initial interface debonding Fig. 10 illustrates the initial interface debonding in PBX mesostructure, in which 4 randomly distributed debonded interfaces are indicated by arrows. The initial debonded interfaces are set with zero tensile and shear strengths. To describe the content of initial interface debonding quantitatively, an interface debonding coefficient α is defined as follows:

(37)

4. Simulations and analysis

n

α=

In reference [38], the construction process of meso-structures for PBX 9501 was described in great detail. In that work [38], 5 groups of meso-structures with different sizes (0.6 mm, 0.9 mm, 1.2 mm, 1.5 mm, 1.8 mm) had been simulated to investigate the influence of the mesostructure size on the effective modulus. It was found that the influence becomes negligible when size of the meso-structure reaches 1.2 mm, thus the model characteristic length was decided to be 1.2 mm, around 5 times of the size of the maximum particle. Therefore, in the present work, the meso-structure model size of 1.2 mm is used directly. In Fig. 9, a meso-structure model is displayed. The corresponding particle size distribution is shown in Fig. 9(c) and the particle volume fraction is about 90.06%. The boundary condition of the meso-structure model under uniaxial tension and compression are shown in Fig. 9(a) and (b), respectively. For the left boundary, the x-displacement is constrained and the ydisplacement is free. For the bottom boundary, the y-displacement is constrained and the x-displacement is free. For the right boundary, the kinetic-couple condition is applied to guarantee that the meso-structure always keeps a square shape during its deformation [24] as illustrated by the dashed lines. For the top boundary, a displacement loading with speed 0.01 mm/s is applied, thus the strain rate of the whole model in the y direction is about 0.008 s−1. The mechanical parameters of the PBX constituents in Table 1, Table 2, and Table 3 are used in the NMM

∑i = 1 li Lp

× 100%

(38)

In which, li is the length of the ith debonded interface, n is the total number of the debonded interfaces, and Lp is the total length of all the particle-binder interfaces in the whole meso-structure. As shown in Fig. 11, a group of PBX meso-structures with α varying from 0% to 50% is constructed. These models are simulated with the loading strain rate and environment temperature of 0.08 s−1 and 25 °C, respectively. The same loading strain rate and environment temperature are considered in all the following simulations in this paper. Fig. 12 shows the failure patterns of the PBX meso-structures with different α values under uniaxial tension. It can be found that, the micro-failure mainly occur on the particle/binder interfaces because of the obviously lower strength of the interface as compared with that of the other two bulk constituents. The fracture of the binder will take place thereafter due to the large deformation under tension. The interface debonding is a usual phenomenon observed in experiments [39]. The debonded interfaces gradually coalesce to form main cracks almost perpendicular to the loading direction and lead to the ultimate macroscopic tensile failure of the whole meso-structures. The location of the formed main cracks is quite different in the six models with different α values. This is because the different distributions of the initial debonded interfaces provide different opportunities for the

Fig. 8. Mohr-Coulomb criterion with a tension cutoff. 7

Computational Materials Science xxx (xxxx) xxxx

G. Kang, et al.

Fig. 9. PBX meso-structure model under uniaxial tension/compression.

bands are formed clearly, indicating an obvious shear failure in the macroscopic. However, when α ≥ 40%, the main shear crack bands turn vague with more scattered micro-cracks in the whole meso-structure. Moreover, in all the six cases, the new interface debonding also likes to take place around bigger explosive particles, along with a phenomenon of the flow of small particles together with the deformation of the polymer binder. Fig. 14 shows the macroscopic effective stress-strain curves of the meso-structures with different α values under uniaxial tension and compression. It can be found that, the effective tensile and compressive strength of the meso-structures (peak value of curves) has an obvious asymmetric character. This is because that the particle/binder interfacial mechanical properties play a dominant role under the tensile loading condition, while under the compressive loading condition, the mechanical properties of all the three constituents of PBXs play important roles. This asymmetric phenomenon is also widely observed in relevant experiments [39]. It is also found that, the increase of the initial debonding content decreases the macroscopic effective tensile and compressive strength dramatically. From Fig. 13, it can be found that the global large strain under tension and compression is mainly attributed to the large deformation of the polymer binder and failure of the particle/binder interfaces. Table 4 shows NMM simulated effective tensile and compressive strength along with some ABAQUS simulation and experimental results. For the NMM simulation results, when α = 0%, the effective tensile strength is close to 1.66 MPa (given interfacial tensile strength). It further verifies the dominant role of the interfacial mechanical property under tensile condition. When α = 0%, the ABAQUS and experimental

Fig. 10. Illustration of the initial interface debonding of PBX meso-structure (Only particle/binder interfaces are displayed, in which black lines represent debonded interfaces).

debonding to spread. However, it also can be observed that the new interface debonding generally likes to take place around bigger explosive particles. Fig. 13 shows the failure patterns of PBX meso-structures with different α values under uniaxial compression. It can be found that, the interfacial debonding is also the main failure mode along with large deformation of the binder, similar to that in tensile condition. Generally, the main cracks, or crack bands, formed by coalesced debonded interfaces eventually lead to the failure of the whole meso-structures in a compression-shear pattern. When α ≤ 30%, two main shear crack 8

Computational Materials Science xxx (xxxx) xxxx

G. Kang, et al.

Fig. 11. A group of PBX meso-structures with different α values (Red lines represent bonded interfaces; Black lines represent debonded interfaces).

For the effective compressive strength, the NMM results are obviously lower than other results. This can probably be attributed by the following factors: (1) The shear spring stiffness value on the particle/ binder interfaces may be unreasonable in the present BCCR; (2) The friction force on the particle/binder interfaces are neglected to achieve

tensile strength is lower than the NMM simulation result, but very close to that when α = 10%. The differences between NMM results and other results may be caused by the application of a complex PBX mesostructure in the NMM simulation, in which a large number of contacts exist on the particle/binder interfaces.

Fig. 12. Simulated failure patterns of PBX meso-structures with different α values under uniaxial tension. 9

Computational Materials Science xxx (xxxx) xxxx

G. Kang, et al.

Fig. 13. Simulated failure patterns of PBX meso-structures with different α values under uniaxial compression.

Fig. 15. Illustration of the initial void in the binder of PBX meso-structure (Only the polymer binder layer is displayed). Fig. 14. Effective stress-strain curves of the meso-structures with different α values under uniaxial tension and compression.

As compared with the case α = 0%, when α = 50%, the effective tensile and compressive strength decreases 74.26% and 60.53%, respectively. The initial interfacial debonding defect reduces the tensile strength more seriously.

better convergence speed of the computation; (3) Even thickness of the binder layer is assumed. Other factors, such as the size of the mesostructure, the particle size distribution and particle gradation, etc., may also have certain degree of influence on the results because of the complexity under the compressive condition. The specific influences of these factors could be investigated by simulations in the future.

4.2. The influence of initial voids Fig. 15 illustrates the initial void defect in the polymer binder of

Table 4 Comparison of NMM simulated effective tensile and compressive strength with other numerical or experimental results under different initial interface debonding coefficient. NMM results under different α values (MPa)

Tension Compression

Other results (MPa)

0%

10%

20%

30%

40%

50%

ABAQUS results

Experimental results

1.589 6.686

1.117 5.764

0.935 4.782

0.737 4.291

0.555 3.381

0.409 2.639

1.18 [14] 9.58 [14]

1.07 [40] 11.11 [41]

10

Computational Materials Science xxx (xxxx) xxxx

G. Kang, et al.

Fig. 16. A group of PBX meso-structures with different β values (Blank parts represent voids).

Fig. 17. Simulated failure patterns of PBX meso-structures with different β values under uniaxial tension.

11

Computational Materials Science xxx (xxxx) xxxx

G. Kang, et al.

Fig. 18. Simulated failure patterns of PBX meso-structures with different β values under uniaxial compression.

Fig. 20. Illustration of the initial cracks in the explosive particles (Only explosive particles are displayed). n

Fig. 19. Effective stress-strain curves of the meso-structures with different β values under uniaxial tension and compression.

β=

∑i = 1 Ai Ap

× 100%

(39)

In which Ai is the area of the ith void, n is the total number of the voids, and Ap is the total area of the whole meso-structure. The produce of voids will lead to the decrease of interface content obviously when the void content is high, which may influence the simulation results. However, the maximum β value considered below is only 2.0%. The study [42] pointed out that the β value for PBX9501 with the PVF of 92% is around 1%~3%. As shown in Fig. 16, a group of PBX

PBXs, in which 4 randomly distributed voids are indicated by arrows. The voids are generated by randomly deleting some manifold elements in the polymer binder domain. Similarly, to describe the content of voids in the meso-structure quantitatively, a void volume fraction β is defined as follows:

Table 5 CComparison of NMM simulated effective tensile and compressive strength with other numerical or experimental results under different initial void volume fraction. NMM results under different void volume fraction β

Tension (/MPa) Compression (/MPa)

Other results

0%

0.4%

0.8%

1.2%

1.6%

2.0%

ABAQUS results

Experimental results

1.589 6.686

1.354 5.961

1.173 5.344

0.953 5.288

0.899 3.946

0.851 3.692

1.18 [14] 9.58 [14]

1.07 [40] 11.11 [41]

12

Computational Materials Science xxx (xxxx) xxxx

G. Kang, et al.

Fig. 21. A group of PBX meso-structures with different d values (Black segments represent micro-cracks in the explosive particles).

Fig. 22. Simulated failure patterns of PBX meso-structures with different d values under uniaxial tension.

13

Computational Materials Science xxx (xxxx) xxxx

G. Kang, et al.

Fig. 23. Simulated failure patterns of PBX meso-structures with different d values under uniaxial compression.

there are more scattered micro-cracks generated simultaneously in the meso-structure. Fig. 18 shows the failure patterns of PBX meso-structures with different β values under uniaxial compression. It can be found that, interfacial debonding also is the main failure mode. When β ≤ 1.2%, obvious main shear cracks, or crack bands, are formed, which indicates compression-shear failure of the whole meso-structures. However, when β ≥ 1.6%, the main shear cracks become vague and there are more scattered micro-cracks generated in the whole model. Fig. 19 shows the macroscopic effective stress-strain curves of the meso-structures with different β values under uniaxial tension and compression. The tensile and compressive strength of the meso-structures also shows an obvious asymmetric character, and the increase of the void volume fraction β reduces the effect tensile and compressive strength of the meso-structures obviously. As shown in Table 5, when β = 0.8% or 1.2%, the NMM simulated effective tensile strength is close to the ABAQUS and experimental results, and the NMM simulated effective compressive strength is also obviously smaller than the ABAQUS and experimental results. As compared with the case with zero initial defects, when β = 2.0%, the effective tensile and compressive strength reduces 46.44% and 44.78%, respectively, a quite close proportion.

Fig. 24. Effective stress-strain curves of the meso-structures with different d values under uniaxial tension and compression.

meso-structures with β value varying from 0% to 2.0% is constructed. Fig. 17 shows the failure patterns of PBX meso-structures with different β values under uniaxial tension. It can be found that, interfacial debonding is still the main failure mode, and different content of the initial voids in the binder brings different locations of the approximately horizontal main cracks that lead to the final tensile failure of the whole meso-structures. Meanwhile, with the increase of the β value,

4.3. The influence of initial micro-cracks Based on the fact that random distribution of micro-cracks exists in real explosive particles, some manifold element nodes in the particles

Table 6 Comparison of NMM simulated effective tensile and compressive strength with other numerical or experimental results under different initial micro-crack density. NMM results under different crack density d

Tension (/MPa) Compression (/MPa)

Other results

0

1.0

2.0

3.0

4.0

5.0

ABAQUS results

Experimental results

1.589 6.686

1.432 6.183

1.406 5.945

1.392 5.879

1.275 5.755

1.054 5.748

1.18 [14] 9.58 [14]

1.07 [40] 11.11 [41]

14

Computational Materials Science xxx (xxxx) xxxx

G. Kang, et al.

and the initial micro-cracks in the explosive particles, are simulated with the extended NMM. Results indicate that:

domain are separated randomly in the initial PBX meso-structure model and zero tensile and shear strengths are set to their surfaces to simulate initial cracks in PBX. As shown in Figs. 20, 4 randomly distributed crack segments are indicated by arrows. To describe the density of crack in the meso-structure quantitatively, a parameter of crack density d (in mm/mm2) is defined as follows:

(1) Under the tension condition, the particle/binder interface property plays the dominant role in the determination of the macroscopic effective tension strength of the PBX meso-structure. Under the compression condition, along with that of the particle/binder interface, the properties of the explosive particles and polymer binder also play important roles to determine the effective compression strength. Therefore, the PBX performs an obvious tension-compression strength asymmetry. (2) The initial defects will lead to the decrease of the tension and compression strength of the PBX in different degree. For initial interfacial debonding and micro-cracks, the tension strength is more likely to be reduced as compared with the compression strength. For initial voids, its influence on the tensile and compressive strength is almost identical. In addition, the initial defects will lead to different microscopic failure modes and macroscopic failure patterns of the PBX meso-structures as well. Under both the tension and compression conditions, the interfacial debonding, along with large deformation and possible subsequent fracture of the binder, is a major microscopic failure mode. The initial micro-cracks in the explosive particles also bring an obvious failure mode of particle fracture. In the macroscopic, different types and content of the initial defects will lead to different locations of tensile main cracks under the tension condition, and different remarkableness of the shear main cracks (or crack bands) under the compression condition. However, an unchanged phenomenon is that failures are like to take place around bigger particles. (3) The microscopic deformation and failure behavior of PBX mesostructures under compression is of much complexity and brings challenges to precise reproduction of it through numerical simulation, which indicates a direction of our future work. Moreover, in the present work, to obtain high volume fraction of the explosive particles, an ideal PBX meso-structure with identical thickness of the polymer binder layer is used. This ideal model is different from the real PBX meso-structure and could influence the simulation results to some extent. How to obtain more realistic PBX mesostructures with high particle volume fraction is also a future work direction. Generally, the developed NMM in the present work should be a promising tool for the simulation of particle-filled composites at the meso-scale.

n

d=

∑i = 1 di Ap

(40)

In which di is the length of the ith crack segment, n is the total number of the cracks, and Ap is the total area of the whole mesostructure. As shown in Fig. 21, a group of PBX meso-structures with d varying from 0 to 5.0 is constructed. Fig. 22 shows the failure patterns of PBX meso-structures with different d values under uniaxial tension. It can be found that, along with the interfacial debonding, the fracture of explosive particles also becomes an obvious micro-failure mode, as signed with arrows. The whole meso-structures fail eventually due to the generation of tensile main cracks that are almost perpendicular to the loading direction. The main cracks like to evolve around bigger particles, and different d value brings different locations of the failure routes as well. Fig. 23 shows the failure patterns of PBX meso-structures with different d values under uniaxial compression. It can be found that, similar to that in the tensile condition, the interfacial debonding along with the fracture of explosive particles produce the main cracks, or crack bands, to lead to the macroscopic failure of the meso-structure models. Meanwhile, the macroscopic shear failure routes seem do not change apparently until d increases to 5.0. This is because the initial defect of micro-cracks is located in the explosive particles. Although particle fracture is involved under the compressive loading, interfacial debonding and binder large deformation are still the main failure factors. Moreover, similar to all the previous examples in the present study, the failures in the meso-structures are likely to take place around bigger particles and the small particles mainly perform as flow matter with the deformation of the polymer binder. Fig. 24 shows the macroscopic effective stress-strain curves of the meso-structures with different d values under uniaxial tension and compression. The tensile and compressive strength of the meso-structures also shows an obvious asymmetric character. The increases of the d value leads to the decease of the tensile and compressive strength simultaneously. As shown in Table 6, the NMM simulation result of the effective tensile strength with a d value of 5.0 is close to the ABAQUS and experimental results. However, similar to the previous examples in the present study, the NMM simulated effective compressive strength is obviously lower than the ABAQUS and experimental results. When d = 5.0%, the effective tensile and compressive strength of the mesostructure decreases 25.74% and 14.03%, respectively, as compared with case without initial defects. The influence of the initial microcracks in the explosive particles on the effective tensile strength is much more obvious than that on the effective compressive strength of the meso-structures. It indicates that the influences of the existence of closed micro-cracks on tensile failure are greater than that on compressive failure.

CRediT authorship contribution statement Ge Kang: Software, Formal analysis, Validation, Writing - original draft. Youjun Ning: . : Software, Methodology, Supervision, Writing review & editing. Pengwan Chen: Conceptualization, Resources, Funding acquisition, Project administration. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

5. Conclusions

Acknowledgements

In the present paper, a bilinear cohesive contact relationship model, a visco-elastic constitutive model, and an elastic visco-plastic constitutive model are implemented into NMM to simulate the three constituents of PBX, respectively. Based on the maximum tensile stress and the Mohr-Coulomb criteria, the microscopic fracturing of the particles and polymer binder is simulated. The deformation and failure behaviors of PBX meso-structures with three categories of initial defects including the initial interfacial debonding, the initial voids in the polymer binder,

This work was supported by the Science and Technology Program Project of Sichuan Province under Grant number 2017JY0128, and the NSAF project under Grant number U1330202. Data availability The data used to support the findings of this study are available from the corresponding author upon request. 15

Computational Materials Science xxx (xxxx) xxxx

G. Kang, et al.

References

[21] B. Banerjee, D.O. Adams, On predicting the effective elastic properties of polymer bonded explosives using the recursive cell method, Int. J. Solids Struct. 41 (2) (2004) 481–509. [22] H.L. Duan, J. Wang, Z.P. Huang, Z.Y. Luo, Stress concentration tensors of inhomogeneities with interface effects, Comput. Mech. 37 (7) (2005) 723–736. [23] H.K. Lee, S.H. Pyo, Micromechanics-based elastic damage modeling of particulate composites with weakened interfaces, Int. J. Solids. Struct. 44 (25) (2007) 8390–8406. [24] K.D. Dai, Y.L. Liu, P.W. Chen, Y. Tian, Finite element simulation on effective elastic modulus of PBX explosives, Trans. BIT 32 (11) (2012) 1154–1158. [25] J.C. Wang, J.R. Luo, Micromechanical study on the effective elastic moduli of polymer-bonded explosives with imperfect interfaces, J. Energ. Mater. 6 (2017) 1–13. [26] J.C. Wang, J.R. Luo, Predicting the effective elastic properties of polymer bonded explosives based on micromechanical methods, J. Energ. Mater. 11 (2017) 1–12. [27] H. Talebi, M. Silani, S.P.A. Bordas, P. Kerfriden, T. Rabczuk, A computational library for multiscale modeling of material failure, Comput. Mech. 53 (5) (2014) 1047–1071. [28] K.M. Hamdia, M. Silani, X.Y. Zhuang, P.F. He, T. Rabczuk, Stochastic analysis of the fracture toughness of polymeric nanoparticle composites using polynomial chaos expansions, Int. J. Fracture 206 (2) (2017) 215–227. [29] G. Kang, P.W. Chen, X. Guo, Y.J. Ning, G.W. Ma, Simulations of meso-scale deformation and damage of polymer bonded explosives by the numerical manifold method, Eng. Anal. Bound. Elem. 96 (2018) 123–137. [30] Y.J. Ning, X.M. An, G.W. Ma, Footwall slope stability analysis with the numerical manifold method, Int. J. Rock. Mech. Min. 48 (6) (2011) 964–975. [31] X.M. An, Y.J. Ning, G.W. Ma, Modeling progressive failures in rock slopes with nonpersistent joints using the numerical manifold method, Int. J. Numer. Anal. Met 38 (7) (2014) 679–701. [32] G.H. Shi, Discontinuous deformation analysis: a new numerical model for the statics and dynamics of block systems, Dept. Civ. Eng. 1988. [33] M. Hu, Y. Wang, J. Rutqvist, Development of a discontinuous approach for modeling fluid flow in heterogeneous media using the numerical manifold method, Int. J. Numer. Anal. Met. 39 (17) (2015) 1932–1952. [34] E.M. Mas, B.E. Clements, B. Blumenthal, C.M. Cady, A visco-elastic model for PBX binders, AIP. 2002:661-664. [35] S.R. Bodner, Y. Partom, Constitutive equations for elastic-viscoplastic strain-hardening materials, J. Appl. Mech. 42 (2) (1975) 385. [36] M. Zhou, A. Needleman, R.J. Clifton, Finite element simulations of shear localization in plate impact, J. Mech. Phys. Solids 42 (3) (1994) 423–458. [37] S. Kim, Y. Wei, Y. Horie, M. Zhou, Prediction of shock initiation thresholds and ignition probability of polymer bonded explosives using mesoscale simulations, J. Mech. Phys. Solids. (2018) S0022509617305513. [38] G. Kang, P.W. Chen, Y.J. Ning, K.S. Ni, Meso-structure construction and effective modulus simulation of PBXs, J. Energ. Mater. (2019) 1684594. [39] P.J. Rae, H.T. Goldrein, S.J.P. Palmer, J.E. Field, A.L. Lewis, Quasi-static studies of the deformation and failure of β-HMX based polymer bonded explosives, Phys. Eng. Sci. 2002 (458) (2019) 743–762. [40] P.J. Rae, S.J.P. Palmer, H.T. Goldrein, J.E. Field, A.L. Lewis, Quasi-static studies of the deformation and failure of PBX 9501, Proc. R. Soc., Phys. Eng. Sci. 458 (2025) (2002) 2227–2242. [41] G.T. Gray, D.J. Deanne, W.R. Blumethal, C.M. Cady, High-and low-strain rate compression properties of several energetic material composites as a function of strain rate and temperature, Off. Sci. Res. Inf. Tech. Rep. 1998. [42] T.M. Willey, L. Lauderbach, F. Gagliardi, T.V. Buuren, Mesoscale evolution of voids and microstructural changes in HMX-based explosives during heating through the β-δ phase transition, J. Appl. Phys. 118 (5) (2015) 055901.

[1] D.J. Benson, P. Conley, Eulerian finite element simulations of experimentally acquired HMX microstructures, Model. Simul. Mater. Sc. 7 (3) (1999) 333–354. [2] A. Barua, M. Zhou, A Lagrangian framework for analyzing microstructural level response of polymer bonded explosives, Model. Simul. Mater. Sc. 19 (5) (2011) 055001. [3] H. Tan, C. Liu, Y. Huang, P.H. Geubelle, The cohesive law for the particle/matrix interfaces in high explosives, J. Mech. Phys. Solids 53 (8) (2005) 1892–1917. [4] M. Li, J. Zhang, C.Y. Xiong, J. Fang, J.M. Li, Y. Hao, Damage and fracture prediction of plastic-bonded explosive by digital image correlation processing, Opt. Laser. Eng. 43 (8) (2005) 856–868. [5] P.W. Chen, F.L. Huang, Y.S. Ding, Microstructure, deformation and failure of polymer bonded explosives, J. Mater. Sci. 42 (13) (2007) 5272–5280. [6] Z.W. Liu, H.M. Xie, K.X. Li, P.W. Chen, F.L. Huang, Fracture behavior of PBX simulation subject to combined thermal and mechanical loads, Polym. Test. 28 (6) (2009) 627–635. [7] Z.B. Zhou, P.W. Chen, Z.P. Duan, F.L. Huang, Study on fracture behaviour of a polymer bonded explosive simulant subjected to uniaxial compression using digital image correlation method, Strain 48 (4) (2012) 326–332. [8] Z.N. Yuan, H. Chen, J.M. Li, B. Dai, W.B. Zhang, In-situ x-ray tomography observation of structure evolution in 1,3,5-Triamino-2,4,6-Trinitrobenzene based polymer bonded explosive (TATB-PBX) under thermo-mechanical loading, Materials 11 (5) (2018) 732–746. [9] L.J. Ling, F. Hua, T.D. Wang, L.F. Yun, Fracture behaviour investigation into a polymer bonded explosive, Strain (2012). [10] J.D. Yeager, K.J. Ramos, S. Singh, M.E. Rutherford, J. Majewski, D.E. Hooks, Nanoindentation of explosive polymer composites to simulate deformation and failure, Mater. Sci. Tech. 8 (9–10) (2012) 1147–1155. [11] L. Chen, D. Han, S.L. Bai, F. Zhao, J.K. Chen, Compressive behavior and damage Evaluation of a PBX substitute material, Mech. Adv. Mater. Struc. (2016) 737–744. [12] J.K. Chen, J.L. Li, L.M. Zhu, K.W. Li, F. Zhao, S.L. Bai, On the tension-induced microcracks’ nucleation in a PBX substitute material under impact compression loading, Int. J. Mech. Sci. 137 (2017) 263–272. [13] V.W. Manner, J.D. Yeager, B.M. Patterson, D.J. Walters, In situ imaging during compression of plastic bonded explosives for damage modelling, Materials 10 (6) (2017) 638–652. [14] Y.Q. Wu, F.L. Huang, A micromechanical model for predicting combined damage of particles and interface debonding in PBX explosives, Mech. Mater. 41 (1) (2009) 27–47. [15] H. Arora, E. Tarleton, J. Li-Mayer, M.N. Charalambides, D. Lewis, Modelling the damage and deformation process in a plastic bonded explosive microstructure under tension using the finite element method, Comput. Mater. Sci. 110 (2015) 91–101. [16] X.J. Wang, Y.Q. Wu, F.L. Huang, Numerical mesoscopic investigations of dynamic damage and failure mechanisms of polymer bonded explosives, Int. J. Solids. Struct. 129 (2017) 28–39. [17] K.C. Bennett, D.J. Luscher, M.A. Buechler, J.D. Yeager, A micromechanical framework and modified self-consistent homogenization scheme for the thermoelasticity of porous bonded-particle assemblies, Int. J. Solids. Struct. 139 (2018) 224–237. [18] K.C. Bennett, D.J. Luscher, Effective thermoelasticity of polymer-ponded particle composites with imperfect interfaces and thermally expansive interphases, J. Elasticity. 136 (1) (2019) 55–85. [19] N. Prakash, G.D. Seidel, Computational electromechanical peridynamics modeling of strain and damage sensing in nanocomposite bonded explosive materials (NCBX), Eng. Fract. Mech. 177 (2017) 180–202. [20] X. Zhang, C. Oskay, Material and morphology parameter sensitivity analysis in particulate composite materials, Comput. Mech. 62 (2018) 543–561.

16