Microcrack and microvoid dominated damage behaviors for polymer bonded explosives under different dynamic loading conditions

Mechanics of Materials 137 (2019) 103130

Contents lists available at ScienceDirect

Mechanics of Materials journal homepage: www.elsevier.com/locate/mechmat

Research paper

Microcrack and microvoid dominated damage behaviors for polymer bonded explosives under different dynamic loading conditions

T

Kun Yang, Yanqing Wu , Fenglei Huang ⁎

State Key Laboratory of Explosion Science and Technology, Beijing Institute of Technology, Beijing 100081, PR China

ARTICLE INFO

ABSTRACT

Keywords: PBX explosives Statistical microcracks and microvoids Void collapse Void distortion Dynamic failure modes

A damaged viscoelastic-plastic model incorporating with both microcrack and microvoid evolution laws has been developed and implemented in the hydrodynamic code DREXH, to predict overall mechanical behavior and gain insights into possible failure mode of polymer-bonded explosives (PBXs) under complicated dynamic loading conditions. Multiple stress-state loadings will motivate evolution modes of microcracks, i.e., frictionlocked, shear with friction, pure shear, mixed shear and open, and normal open, or microvoids, i.e., void collapse, void distortion. The simulated results connected the asymmetry of PBXs under uniaxial tension and compression with the critical stresses for open-crack and shear-crack initiation. Under high-confinement compression, crack growth is inhibited due to friction-locked (stable) state while more voids become collapsed accompanying with a strong plastic flow, corresponding to the occurrence of brittle-to-ductile transition for PBXs. Under simple shear loading, growth of pure-shear crack contributes more on failure rather than void distortion for pressed PBXs. Effects of void distortion on failure become more important for PBXs containing higher-fraction ductile binder. The simulated results agree with experimental observations and contribute to deepening our understanding of underlying microcrack and microvoid mechanisms for dynamic mechanical failure behavior of PBXs.

1. Introduction Polymer bonded explosives (PBXs) are kind of particulate-reinforced composite materials that consist of high volume-fraction (∼90%) energetic crystals held together with a viscoelastic polymer binder as well as additives such as nitro-plasticizer antioxidant and radical inhibitor (Wang et al., 2016). They are widely used in civil and defense related fields, such as mining, rocket propellants, and warheads (Asay, 2010). The heterogeneous microstructural defects in PBXs (i.e., cracks and voids) not only contribute to the macroscopic failure of material but also exacerbate energy localization process, thereby affecting hotspots formation and ignition (Barua et al., 2012). In general, several possible micromechanisms could jointly result in the material failure or hotspots formation; and the dominant mechanism (if any) will change with stress states and constituent contents (Malcher et al., 2012; Brünig et al., 2018). Therefore, understanding the overall damage response of PBXs and its non-unique potential micromechanisms is fundamental to the safe use of PBXs. Both brittle and ductile failure could occur in typical pressed PBXs (e.g., PBX9501, 95wt% HMX, 2.5wt% estane, and 2.5wt% BDNPA/F) since that such materials not only have high-fraction of brittle crystals

⁎

containing large quantities of inherent microcracks but also small quantities of high-ductile binder (Trumel et al., 2010). Generally, they often exhibit brittle-failure behavior especially without confinement (Parab et al., 2016; Ravindran et al., 2016; Ravindran et al., 2017). Parab et al. (2016) pointed out that, under dynamic uniaxial compression, most of the cracking in the crystals were formed due to the tensile stress generated by the diametral compression resulted from the contacts between the crystals. Ravindran et al. (2017) studied the dynamic multi-scale deformation of a PBX simulant with relatively highbinder content (85wt% sugar crystals, 15wt% binder). Their results suggested that delamination of the binder from crystals and binder cracking are the main local failure modes. The damage mechanisms of pressed PBXs exhibit a strong sensitivity to loading conditions. Under high confinement, a brittle-to-ductile transition will occur in PBXs (Wiegand et al., 2005 and 2011; Le et al., 2010; Bailly et al., 2011), which is similar to many other quasi-brittle materials, i.e. concretes, rocks, and ceramics. Both the experimental results of Wiegand et al. (2005) and Bailly et al. (2011) indicated that high-confined PBXs exhibit a strong pressure-dependent plasticity. The decrease of brittle behavior could be the result of partly-inhibited microcracking. The mesoscale factor responsible for the macro-plastic

Corresponding author. E-mail address: [email protected] (Y. Wu).

https://doi.org/10.1016/j.mechmat.2019.103130 Received 4 April 2019; Received in revised form 22 June 2019; Accepted 22 July 2019 Available online 23 July 2019 0167-6636/ © 2019 Elsevier Ltd. All rights reserved.

Mechanics of Materials 137 (2019) 103130

K. Yang, et al.

deformation is not sole, including plastic flow of binder, collapse of voids, contact flattening between grains, and crystal twinning. Further experimental or numerical evidences are required to determine the main reason for this type of brittle-to-ductile transition (Wiegand et al., 2011). In contrast to axisymmetric stress states (uniaxial or triaxial compression with lateral confinement), which is characterized by a high stress-triaxiality, PBXs also often undergo a severe shear-dominated deformation with a low triaxiality, i.e., penetration and punch (Peterson et al., 2001). The experimental results of Liu et al. (2016) indicated the initiation and development of shear bands in a low-rate punched PBX. Skidmore et al. (2000) observed the formation of a wedge structure and the evidence of chemical reaction in open cracks generating from the high shear region for an impacted PBX9501 sample. The fundamental micromechanisms controlling damage of shearing PBXs remain poorly understood due to few available experimental evidences (Hanina et al., 2018). In the current study, respective contributions by shear crack growth and void distortion to shear failure of PBXs were studied to produce the abundant damage modes. Besides loading conditions, whether PBXs mainly exhibit brittle or ductile damage behavior is affected by two other factors, namely, manufacture techniques and constituent contents (Boddy et al., 2016; Drodge et al., 2016). For example, unlike pressed explosives, some casted explosives (e.g., PBXN-110, 88wt% HMX and 12wt% binder) contain few incipient microcracks and mainly exhibit ductile-damage (Çolak et al., 2004). Moreover, solid propellants usually consist of lowfraction (75–80 wt%) crystal and a very soft rubber binder so that they are enough tough to prevent the formation of macro-cracks and withstand the extreme vibration during launch (Asay, 2010). The low modulus of the binder allows PBXs to suffer a larger deformation and to absorb most of mechanical energy resulting from shock loading, thus contributing to the increase of structural toughness and the decrease of impact sensitivity of PBXs. Several micromechanical approaches have been proposed to model the dependence of PBXs’ damage on the constituent fraction and interfacial strength (Tan et al., 2005; Xu et al., 2008; Wang and Luo, 2018). Tan et al. (2005) used the microscale cohesive law and the extended Mori-Tanaka method to study the debonding of particle-matrix interface for PBX9501. Xu et al. (2008) proposed a general 3D nonlinear macroscopic constitutive law for propellants to model microstructural damage evolution upon straining via continuous void formation and growth. Wang and Luo (2018) treated PBXs as pseudo polycrystals consisting of equivalent composite particles (crystals with binder coating) and developed a multistep micromechanical approach. Clements and Mas (2001) derived a constitutive model for filledpolymer composites via Schapery's theory based on the effective medium method. Since that (i) the high contrast in stiffness between crystal and binder (at least 104:1) (Wu et al., 2009); (ii) the volume fraction of crystals (∼90%) in PBXs is much higher than the typical filler content (30–40%) limited in micromechanics approach (Wang et al., 2016), it's still a challenge to accurately clarify the contributions by each constituent to damage response by these micromechanical approaches. For studying the effects micro-defects of on damage, lots of explicit mesoscale simulations have been performed (Barua et al., 2012, 2013; Arora et al., 2015; LaBarbera and Zikry, 2015; Wang et al., 2017a and Wang et al., 2017b; Kang et al., 2018; Grilli et al., 2018). Wu et al. (2009) applied a nonlinear cohesive law for grain-binder interface to study the effect of interfacial debonding on failure of quasistatic loaded PBXs. Using the cohesive finite element method (CFEM), Barua et al. (2013) tracked the contributions of individual constituents, fracture and frictional contact along failed crack surfaces to the heating of PBXs. LaBarbera and Zikry (2015) utilized a dislocation densitybased crystalline plastic model and a finite viscoelastic model to study dynamic nucleation and propagation of fracture of RDX/Estane aggregates. Wang et al. (2017b) employed a full field method to consider

the local damage behavior by modeling the cleavage damage of HMX crystals and the damage due to binder rupture and debonding. Although mesoscale modeling is able to explicitly describe micro-defects evolution, it's still impractical to extend these models to study the continuum-scale damage responses of PBXs under dynamic complicated loading conditions (Yang et al., 2018). Thus, a robust physically-based continuum damage mechanics (CDM) model is still desirable to implicitly account for the effects of factors (i.e., loading conditions and constituent contents) on the change of damage; and to help understand the relation between the macroscale and mesoscale damage under complicated loadings. In recent years, plentiful micromechanics motivated CDM models have been developed to predict the stress-state dependent damage behavior for a wide class of materials, i.e., metals, concretes, and polymers (Zhu et al., 2010; Lecarme et al., 2011; Malcher et al., 2012, 2014; Guo et al., 2014). For predicting ductile fracture of materials, many researchers introduced more effects in the constitutive formulation or in the damage evolution law of original Gurson model (Gurson, 1977), such as the Lode angle dependence, viscoplastic effects, and shear mechanisms (Hutchinson and Pardoen, 2000; Zuo and Rice, 2008; Gao et al., 2009). Malcher et al. (2014) coupled a new shear mechanism and a new micro-defects nucleation mechanism into the Gurson-TvergaardNeedleman (GTN) model to predict ductile fracture under a low stress triaxiality. To study the damage mechanisms of fiber-reinforced composites, Chen and Ghosh (2012) employed the GTN model with a Johnson-Cook type hardening law for ductile matrix and an isotropic CDM model for brittle fiber. Gao et al. (2009) presented a detailed analysis of the effects of hydrostatic stress and Lode angle on the fracture of aluminum alloy. Lu et al. (2013) presented a coupled microcrack-based damage model and flow modeling approach for progressive fracturing in permeable rocks. Brünig and Michalski (2017) proposed a damage rule accounting for different stress-state-dependent damage mechanisms based on the laws of irreversible thermodynamics for concrete materials. Specific to describe the deformation and damage responses of PBXs under dynamic loadings, few physical models have been proposed and most of them were derived based on the statistical crack mechanical model (SCRAM) (Bennett et al., 1998; Dienes et al., 2006; Zuo et al., 2006; Liu et al., 2018). Dienes et al. (1978, 2006) originally proposed the SCRAM model which considered the anisotropic damage by tracking the growth and coalescence of cracks with several typical orientations. Zuo et al. (2005, 2006) developed the Dominant Crack Algorithm (DCA) where damage evolution is described by the most unstable orientation for cracks in the material. Superior to the phenomenological models applied to PBXs (Gruau et al., 2009; Reaugh, 2010), the SCRAM-like models could specifically describe the microcrack-related damage and introduced the interfacial friction-heating of shear-crack as a single active hotspot mechanism. Whereas, the evolution of microcrack does not always play a dominant role on the damage and ignition of PBXs. For the type of PBXs with high-ductility, the nucleation, growth, and coalescence of microvoids are usually treated as the dominant ductile failure mechanism. Furthermore, when PBXs undergo a high-confined compression, crack growth is inhibited and void collapse will become a significant hotspot mechanism (Wiegand et al., 2011; Akiki et al., 2015). Therefore, a constitutive model incorporating with both microcracks and microvoids dominated damage or hotspot mechanisms is required to exactly predict the failure and ignition of PBXs under complex mechanical insults. The present study seeks to develop a physical model that describes both the micro-crack (brittle) and micro-void (ductile) related damage for multiple types of PBXs or same type of PBXs under various dynamic loading conditions, which may present different damage mechanisms. Several cracking-damage modes including friction-locked, shear with friction, pure shear, mixed shear and open, and normal open were incorporated in the model by utilizing the DCA model (Zuo et al., 2005, 2006). In addition, two void-damage modes including void collapse and 2

Mechanics of Materials 137 (2019) 103130

K. Yang, et al.

void distortion were also considered using the modified Gurson model developed by Nahshon and Hutchinson (2008). The features of the current model are illustrated by a study of mechanical and damage behavior of PBXs over a series of loading conditions (uniaxial stress/strain, different lateral confinements, hydrostatic compression, and simple shear loading) in Section 4. The dependence of damage-mode on defects-related parameters is presented in Section 5. Two main goals of the current work are as follows: (i) to extend our understanding of the primary micromechanisms for the overall damage response of PBXs under dynamic loadings; (ii) to lay the first stone for establishing a coupled mechanical-thermal-chemical model which could predict multiple damage-related hotspots formation and initiation of PBXs.

binders [Fig. 1(b)]. Several defects-related parameters (i.e., initial crack density, crack size, and void fraction) are introduced to characterize initial defect conditions and damage evolution of PBXs. Five crack-damage modes (friction locked, shear with friction, pure shear, mixed shear and open, and normal open) and two void-damage modes (void collapse and distortion) are considered. In addition, viscoelastic effects of PBXs, which are resulting from binders, is described by a generalized Maxwell model. The plastic flow of porous PBXs is described by the Gurson yield surface; and plastic deformation of solid matrix (void-free) is described by an isotropic-hardening rate-dependent yield model. A conceptual diagram of all kinds of considered mechanisms in the model is plotted in Fig. 2. Note that different types of PBXs could exhibit a variety of characteristic initial damage morphologies (i.e., high porosity, low crackdensity) due to the variations of manufacturing process, constitute contents and material properties. In the following discussion, we mainly focus on the mechanical behavior of typical pressed PBXs mentioned above. To evaluate the model's potential in describing other PBXs with different initial defects, dependences of damage-mode on defects-related parameters were investigated in Section 5.

2. Theory and constitutive models 2.1. Material microstructural morphology and modelling Fig. 1(a) shows a typical micrograph of a kind of pressed PBXs (PBX9501) obtained from (Tan et al., 2005). From the figure, the individual explosive crystals exhibit defects including large quantities of randomly distributed microcracks which are formed in pressed process; and a small amount of intragranular voids which are generated from crystal growth (Skidmore et al., 1997). Moreover, extragranular voids are always located at the interface between crystals and binder or between prills. The formation of these extragranular voids is due to several reasons: (i) machining process, such as pouring and non-uniform solidification of melt explosives; (ii) storage process, loss of volatiles in aging of explosive; (iii) service process, mechanical or thermal cycling loadings. From the above, in the current model, we assume that the cracks are mainly distributed in crystals and voids are mainly distributed in

2.2. Decompositions of stress and strain tensor The total strain rate ( ), which is calculated by differentiating the velocity field from the momentum equation in numerical analysis, can be decomposed into deviatoric and hydrostatic components,

=e+

vi

(1)

where e is the deviatoric strain, ɛv is the hydrostatic strain and calculated by v = kk /3, the dot denotes time derivative, i is the second-order identity. The evolution of stress rate tensor ( ) is divided into corresponding

Fig. 1. (a) Typical micrograph of PBX 9501, which is obtained from Tan et al.’s publication (Tan et al., 2005); (b) a schematic of PBXs microstructures in the current model.

Fig. 2. A conceptual diagram of all kinds of considered mechanisms in the current model. 3

Mechanics of Materials 137 (2019) 103130

K. Yang, et al.

deviatoric and hydrostatic components,

=s+

2.4. Hydrostatic response (2)

mi

To describe the irreversible compaction behavior of porous PBXs at low pressures compression and reduce to the correct Hugoniot description at high pressures, a porosity-dependent equation of state for porous PBXs is proposed following the method of Herrmann (1969). The pressure in the porous PBXs can be related to that in the void-free solid as follows,

where s is the deviatoric stress, σm is the mean stress whose value equals to the negative pressure, namely, m = P . 2.3. Deviatoric stress-strain formulations Based on the idea of Dienes et al. (2006), material deformation is often the consequence of numerous independent physical processes, i.e., elastic deformation, plastic flow, and opening/sliding of crack, which can be superimposed to obtain an overall deformation. In the model, the total deviatoric strain rate (e ) is decomposed into the viscoelastic deviatoric strain rate (e ve ), the cracking deviatoric strain rate (e cr ), and the plastic deviatoric strain rate (e p ), as given by Eq. (3),

P ( , e, f ) = (1

where f is the current void fraction, ρ and e are the density and specific internal energy of the porous material, respectively. The corresponding density and specific internal energy for the solid, ρs and es, are expressed as, s

(3)

e = e ve + e cr + e p

N

s=

(n )

(4)

e

c¯ a

2

e c¯ 3 c¯ s+ s a 3 a

e

P=

64 (1 = 15(2

e p)

Bs e

N

s (n) = 2G (n) A (e

e p)

G (n ) Bs + (A G

N

1) n=1

s (n) (n )

A

(14)

f ) Ps

(15)

1

:

=

1

( P

v

+ sij eij )

(16)

The plastic deformation of porous PBXs was described by the original Gurson's theory (Gurson, 1977). The yield surface of the Gurson Model F is given in terms of the effective stress σe and mean stress P by,

F ( e, P, f ) =

e

YM

2

+ 2f cosh

3P 2YM

f2

1

(17)

The yield strength of the matrix solid material (f = 0) is defined by Eq. (18), which considers the hardening response and rate-dependence of the material under dynamic loading.

YM = (

0

+ h ¯Mp )[1.0 + C ·ln(1 + *)]

(18)

where σ0 is the initial yield stress at 0 =10 s−1; ¯Mp is the effective plastic strain of matrix; * = / 0 is the normalized equivalent strain rate; h is the hardening modulus; C is the coefficient of strain rate. Note that the Gurson surface was first developed for void growth under a tensile stress state (P<0). In the current model, it is assumed that a reflection of the Gurson surface about the porosity-effective stress plane (P = 0) can be used for void compaction when the material is under a compressive stress state (P<0). When f→0, the Gurson surface

s (n )

s (n ) (n )

(13)

s s)2

(1

2.5. Gurson yield surface

c a)3) 1, B = A e ( a )2 a , and C = A n = 1 (n) . where A = (1 + 3 (¯/ The expression for the deviatoric stress rate for the nth Maxwell element is, c¯

(12)

v

es = e =

(8) c¯

s s es

The rate of the specific internal energy is given by the first law of thermodynamics (neglecting heat flux and heat sources),

(7)

C

+

2 s c0 s

f Ps + (1

=

(6)

) N0 )G

2

The rate of porous density can be calculated according to the law of mass conservation,

Combining Eqs. (4) and (6) provides an expression for the deviatoric stress rates as,

s = 2GA (e

s s

where c0 and s are the constants of solid material. The rate of pressure is obtained by taking time derivative of Eq. (10),

where c¯ is the mean crack radius, which is isotropic during loading (c¯ (n, t ) = c¯ (t ) ); a is the initial flaw size and expressed as a 3 = 6G ; β is a parameter relating to shear moduli, Poisson's ratio ν, and initial crack distribution N0; αe is a parameter relating to the state of open crack.

3, if P 0 = , 5 , if P < 0

(11a) (11b)

PH ( s ) =

(5)

f ) Gs

1 2G

,

where ηs=1-ρs0/ρs=-εv with ρs0 being the initial density of the solid, Γs is the Grüneisen coefficient, and the Hugoniot pressure PH is expressed,

Consider an ensemble of penny-shaped microcracks (with size c and normal n) randomly distributed within a statistically homogeneous volume of the material. The applied stress may result in the opposite faces of cracks to slide and/or open, thus causing an additional strain (crack strain) and an increase of macroscopic compliance. Following the works of Addessio and Johnson (1990), Bennett et al. (1998), and Yang et al. (2018), the total crack strain due to an ensemble of opencrack and shear-crack of all sizes (0 ≤ c≤∞) and orientations can be written as Eq. (6) with several assumptions: (i) the initial crack number density is isotropic (N0(n)=N0); (ii) the crack number density is an exponential function of c; (iii) neglecting the interactions between cracks.

e cr =

f

Ps ( s , es ) = PH ( s ) 1

where N is the number of elements in the generalized Maxwell model;s(n)and τ(n) are the deviatoric stress component and relaxation time N for the nth Maxwell element; Gs = n = 1 (Gs(n) ) is the shear modulus for the intact matrix (void-free) and the effects of void damage on shear modulus are described by,

G = (1

1

It should be noted that Eq. (11b) is the result of neglecting the surface energy of the voids. The pressure-volume relation in the solid follows the Mie-Grüneisen EOS as,

s (n )

n=1

=

es = e

The viscoelastic response of PBXs is described by a generalized Maxwell model, in which the strain is common for all elements of the model and the stress for an individual element is additive. The relation of deviatoric stress rate to deviatoric strain rate and stress can be expressed as, (Bennett et al., 1998)

2Ge ve

(10)

f ) Ps ( s , es )

(9) 4

Mechanics of Materials 137 (2019) 103130

K. Yang, et al.

reduces to the classical isotropic plasticity theory based on the VonMises yield surface. The plastic strain rate is calculated by, p ij

where J3 and σe denote a third invariant of stress and von-Mises stress, respectively. Two variables are written as,

J3 = det(s ) =

F ( e, P , f )

=

(19)

ij

(1

p

f ) YM ¯M =

p ij ij

(20)

Whilst PBXs undergo a high-pressure confined compression, crack growth is inhibited and collapse of void becomes vital for both damage and hotspots formation of explosive materials (Wiegand et al., 2011; Akiki et al., 2015). In addition, PBXs are often encountered in a severe shear-dominated damage and failure with a low triaxiality (e.g., punching). The voids in material may exhibit several features under shear loading, such as shape deformation, reorientation, and inter-void interaction (void sheeting), thereby contributing to the failure of material (Morin et al., 2016; Guo et al., 2018). Hence, we adopted the modified Gurson model developed by Nahshon and Hutchinson (2008) to characterize two competing mechanisms for void evolution behavior, namely, compression causing a decrease in void fraction [Fig. 3(a)] and shear causing an effective increase [Fig. 3(b)]. Nahshon and Hutchinson (2008) modified the original Gurson model to include a contribution to the void growth rate for pure-shear stress states in a manner which leaves the relation unaltered for axisymmetric stress states. The evolution of void fraction is expressed as Eq. (21). The modification is motivated by that the volume of voids undergoing shear may not increase, but void deformation and reorientation contribute to softening and constitute an effective increase in damage.

f)

p v

+ fk w ( )

( )=1

m )( 3

m)

(23a) (23b)

f

p v

v

s

1

f

kw ( )

sij eijp (24)

e

Substituting Eqs. (16), (21) and (24) into Eq. (14) and the rate of pressure is obtained,

P= =

(21)

f Ps + (1 K

v

+

f)

s sij eij

+

s

Ps s p v

+

Ps es es

+ fk w ( )

sij eijp e

(25)

where the bulk modulus of porous material K is related to that of the solid Ks by K = (1 f ) Ks , and = K (1 + s) P , = Ks (1 + s) Ps . The value of Ks increases with the pressure, corresponding to the development of shock wave in the material under high-velocity impact loading,

p /3 is the plastic volumetric strain rate, kw is a material where vp = kk constant related to the shear states, ω is the parameter related to the current stress state (Eq. (22)). To ensure the consistency of subsequent derivation in EOS, the second term sij ijp in the Nahshon et al.’s model is replaced by sij eijp in the current model with minor effect.

27J3 2 e3

=

s

sij eijp e

m )( 2

The value of ω lies in the range, 0≤ω≤1, with ω=0 for all axisymmetric stress states and ω=1 for all states comprised of a pure shear stress plus a hydrostatic contribution. Note that the method relating the void distortion to the stress-state parameter ω is relatively phenomenological. However, this method has been widely-used in describing shear failure of metals under low-triaxiality loading and obtained reasonable results (Nahshon et al., 2009; Xu et al., 2014). This method is predictably reasonable to be used to describe the effects of void distortion in binder on PBXs damage. In general, the nucleation, growth, and coalescence of microvoids are the dominant ductile failure mechanism for metallic materials under tension (Malcher et al., 2014). Whereas, the study of Ellis et al. (2005) for pressed PBX indicated that PBXs under tension exhibit brittle failure (with little plastic deformation). The growth of Mode-I microcracks rather than void coalescence plays a more important role on the failure. Thus, increase of void evolution due to void growth and coalescence in PBXs under tension is not considered in the current model. Applying Eqs. (11a) and (21) into Eq. (15) gives the rate of the solid density,

2.6. Evolution of void fraction

f = (1

1

3 sij sij 2

=

e

where the plastic multiplier determined by the consistency condition. As in the Gurson model, plastic work in the matrix is equated to macroscopic plastic work according to Eq. (20) and the effective plastic p strain rate of the matrix ( ¯M ) can be obtained.

1 sij sik sjk = ( 3

Ks

2

(22)

s

Ps

+

s Ps

s

(26)

The Eq. (25) relates the pressure rate to the applied volumetric

Fig. 3. Two types of void evolution considered in the current model (a) void collapse; (b) void distortion. 5

Mechanics of Materials 137 (2019) 103130

K. Yang, et al.

strain rate ( v ), the strain rate arising from void collapse ( vp ), and an effective strain rate due to void distortion [ fk w (sij / e) eijp ]. Moreover, the work done by shear deformation influences the volumetric response of material as a result of using a Mie-Grüneisen EOS. That has an explicit physical meaning, namely, shear deformation of the material increases the internal energy, which in turn increases the pressure in material.

Appendix A.where ¯ is the effective surface energy of the material,

L ( , n) c¯ G

(27)

where L(σ, n) is given by,

L ( , n) =

(1

v 2

(sn + µ

n

)

)2

2 n

+ sn2,

sn + µ

n

>0

,

n

n

(28)

0

3. Trial-correction numerical algorithm

where μ is the friction coefficient; the normal and shear components of remote traction are n = n · n and sn = [n · 2 n (n · n)2]0.5 , respectively. Given the crack size and stress state, the critical crack orientation (nc) maximizes the function L(σ, n) and the energy-release rate g ( , n , c¯) . The crack with that orientation becomes unstable at the lowest applied stress and is defined as the dominant crack. Energy-release rate is defined by applying Eq. (27) to the critical orientation g ( , c¯) = g ( , nc , c¯) . The crack growth modes can be divided into five-types based on applied stress states (divided in the maximum σ1 and minimum σ3 principle stress plane by stress biaxiality r=σ3/σ1), namely, [I] ) r 1); [II] Mixed tensile and shear Pure tension ( (1 (1 ) ); [III] Pure shear ( Hµ2 < r 1); [IV] Shear ( 1
A trial-correction time-integration algorithm was established for solving the current modeling equations. At the beginning of time step (t), the state variables of material in the last cycle, such as the stress, mean crack size, and porosity, will be treated as input variables. The objective is to update the corresponding values at the end of the step (t +∆t), for a given total strain rate Δɛ over the step (∆t). The trial state is first calculated assuming that the step is purely elastic, namely, ijp = 0 . In view of Eqs. (8), (25), and (21), the trial deviatoric and hydrostatic stress and porosity are expressed as Eq. (30a–c), respectively.

s tr = s t + 2GAt e

P tr = P t

withHµ = µ2 + 1 + µ . The critical orientation (nc) and energy-release rate g ( , c¯) in each crack-growth mode are listed in Table 1. The complete details of the method calculating the critical orientation are introduced in

Remarks

Critical orientation (nc)

Energy-release rate g ( , c¯)

[I]

Pure tension Mixed tensile and shear Pure shear Shear with friction Friction locked

nc = e1

2¯ 1 Scr (c¯) / 1

[II]

[III] [IV] [V]

1 2

1 {( 2

o

nc , e1 =

s

nc , e1 = tan

1 (1/

c

nc , e3 = tan

1 (H ) µ

−

cos

1)

1+r } 1 r

(

2

+

t s sij

(30b)

eij

(30c)

s (n) tr = s (n) t + 2G (n) At e G (n ) t t B s + (At G

N

1) n=1

s (n) t (n )

+ At

s (n ) t (n )

t

(31)

The trial state is then checked to determine whether or not it lies outside the yield surface F ( ijtr ) . If F ≤ 0, then the stress state is indeed elastic and the calculation for the step is complete. Otherwise (F > 0), the step involves plasticity and the trial stress state should be corrected. Plastic strain for the step can be calculated using the trial stress. The obtained plastic strain is then used to correct the trial stress rate, and the corrected result is named as the updated stress rate at the end of time step. From Eqs. (8), (25), (21), and (30a–c), we can obtain,

/2

+ 3 2 3 )2 1)( 1 ) +( 1 2 2 2¯ 2 Scr (c¯)

Hµ ( 1

v

where superscript tr denotes the trial results. Using Eq. (8), the trial deviatoric stress for the nth Maxwell element is expressed as,

2¯ 1 3 2 (c¯) Scr

r)

K

(30a)

(B t s t + C t ) t

f tr = f t

Table 1 The critical crack orientation (nc) and energy-release rate g ( , c¯) in each type of crack growth. (Zuo et al., 2005). Type

(29)

where cmax is the terminal speed for crack growth. The damage of bulk material due to microcrack growth is defined as dcr = c¯3/(a3 + c¯3) to connect the evolution of microcracks with damage response of PBXs at the macroscale. It can be expected that interactions between micro-defects, i.e., cracks and voids, cracks and cracks, are important to the evolution of damage. However, to the authors’ knowledge, few continuum damage models could be used to detailedly describe these interactions due to the difficulties in theoretical derivation. And at least until now, it's still a challenge for mesoscopic simulations which could explicitly model these defects. In the model, evolution of microcracks in the crystal is assumed to be the main reason for brittle damage of PBXs and microvoid in the binder is for ductile damage. These two micro-defects damage parameters are mutually independent to some extent. However, the evolution of microcrack and microvoid will interact each other in an indirect way through stress states. Namely, crack growth will affect the calculation of deviatoric stress and von-Mises stress σe, thereby influencing the calculation of plastic part and void evolution by the Gurson surface. Conversely, void growth affects the calculation of hydrostatic stress and further affects the calculation of current stress states and crack growth.

If the applied stress is sufficiently large, some cracks in the material will become unstable and grow in size. In the current model, we adopted the Dominant Crack Algorithm (DCA) developed by Zuo et al. (2006) to describe the evolution of cracks. The basic assumption of the DCA model is that damage is related to the growth of crack with the most unstable (dominant) orientation rather than the mean crack size over all crack orientations. In this section, a brief outline of the DCA model is presented and the complete details of the model should be referenced as (Zuo et al., 2005, 2006, and 2010). The energy-release rate for a single crack with size c¯ and normal n is,

41 2

2¯ g ( , c¯)

c¯ = cmax 1

2.7. Dominant microcrack growth model

g ( , n , c¯) =

G¯

2

Scr (c¯) . 21 c¯ The crack growth law is calculated by Griffith instability criterion,

3) + 2µ 3 ¯ Scr (¯) c

−

s t+ 6

t

= s tr

2GAt (e p)

(32a)

Mechanics of Materials 137 (2019) 103130

K. Yang, et al.

Pt+

f t+

t

= P tr +

t

p kk )

(

= f tr + (1

+ kw f t+ p kk )

f t+ t ) (

t

t

s t+

t

(e p )

F( (32b)

t+ t e

+ k w f t+

t+ ts

t

t

(32c)

The increment in the plastic strain is assumed to be normal to the Gurson surface in the stress space, t+ t,

F(

( p) =

f t+ t )

sijt+

(33)

t+ t

The increment in the deviatoric and hydrostatic plastic strain can be written as,

(e p ) =

(

p kk )

t+ t ,

F(

st +

=

f t+ t ) t

t+ t ,

F(

=

f t+ t )

P t+ t

3 t+ s 2 YM

=

t

3f t + YM

sinh

3P t + t 2YM

(34b)

Substitution of Eq. (34a) into Eq. (32a) yields,

1+

6At G 2 YM

s t+

t

6At G 2 YM

t+ t e

= s tr

(36a)

Using the associative flow equation Eq. (34a,b) into Eq. (32b) yields,

Pt+

t

+

3f t + YM

t

sinh

3P t + t 2YM

kw f t+

t

t2

2 YM

t+ t e

= P tr

(36b)

Similarly, substituting Eq. (34a,b) into Eq. (32c) for the porosity gives,

f t+

t

+ (1

f t+ t )

3f t + YM

t

sinh

3P t + t 2YM

k w f t+

t+ t e tr tr sij e

=

sijtr 1+

6At G 2 YM

(37)

To illustrate the main features of the current model, several typical loadings (i.e., uniaxial stress, different lateral confinements, uniaxial strain compression, hydrostatic compression, and simple shear) have been applied to PBX material and corresponding deformation and damage behaviors were studied. PBX9501 was selected as a numerical sample because it's a typical explosive which could exhibit both brittle and ductile behaviors under different loadings. Parameters related to the cracking model and Mie-Grüneisen EOS for PBX9501 were verified by Dienes et al. (2006) using the multiple shock experiment data. Parameters for the viscoelastic and plastic model were cited from (Bennett et al., 1998) and (Yang et al., 2018). PBX9501 generally has an initial porosity (f0) of 1∼3% (Millett et al., 2017) and in this section, f0 is determined as 1.0%. The effects of the change of parameter on simulated results will be discussed in Section 5. The void distortion related parameter kw was adopted from (Nahshon et al., 2009) and approximated as 2.0. In the future, more experimental stress-strain curves of PBX9501 under shear loading are required to validate the value of kw. All parameters employed in the model, and the sources for these values, are listed in Table 2.

(35)

tr e

=

=

4. Numerical examples for different loading conditions

Then the effective stress can be obtained by,

1+

t

As such, the other state variables (e.g., mean crack size) at the end of time step can be calculated. The procedure discussed above was implemented by a set of subroutines and integrated into a Lagrangian finite element code Drexh (Yang et al., 2017). The numerical results for typical pressed PBXs under different loading conditions will be presented in Section 4.

(34a) t

(36d)

P t+ t , f t + t ) = 0

A Newton iteration algorithm was employed to solve four nonlinear coupled equations Eq. (36a–d) with the variables ( et + t , P t + t , f t + t , ). Given that the calculation of Von-Mises stress (Eq. (23b)), the deviatoric stress sijt+ t (orsij(n) t + t ) in the final state can be updated once t+ t are obtained, e

(e p )

t+ t e

t+ t , e

t

t2

2 YM

t+ t e

= f tr (36c)

The consistency equation is given by,

Table 2 Summary of material properties for PBX9501. (Dienes et al., 2006; Bennett et al., 1998; Yang et al., 2018; Millett et al., 2017; Nahshon et al., 2009). Conditions Matrix (uncracked solid) Viscoelasticity

Property Initial mass density Shear modulus (solid) Poisson's ratio Shear modulus and relaxation time for the nth Maxwell element

Plasticity

Initial yield strength

Microcracks

Microvoids Equation of state

Symbol ρs0 Gs ν G(1)∼G(5) 1/τ(1) ∼1/τ(5) 0 y

Strain rate coefficient Hardening modulus Initial crack size Crack number density Specific surface energy Rayleigh wave speed Friction coefficient Initial porosity Void distortion parameter Sound speed Us-Up slope Grüneisen constant

C h c¯0 N0 ¯ cmax μ f0 kw c0 s Гs

7

Value Units 1860 kg∙m−3 3.2 GPa 0.3 944, 173.8, 521.2, 908.5, 687.5 MPa 0, 7.32×103, 7.32×104, 7.32×105, 2.00×106 s 12.0 MPa 0.76 500 30 30 0.5 300 0.5 0.01 2.0 2500 2.26 1.50

MPa μm cm−3 J∙m−2 m/s

m/s

Mechanics of Materials 137 (2019) 103130

K. Yang, et al.

4.1. Asymmetry in tension and compression from microcrack initiation

hardening due to plastic deformation. The crack is initiated at point A’ in tension and at A* in compression at 2000 s−1 strain-rate. The critical stress for crack growth under compression (A*, 42.12 MPa) is more than twice that for tension (A’, -14.13 MPa), indicating the asymmetry in tension and compression. According to the dominate crack growth model described in Section 2.7, under uniaxial tension loading (σ1≥0, σ2=σ3=0), the crack growth mode is Mode [I]. In contrast, under uniaxial compression loading (σ1=σ2=0, 0≥σ3), the crack growth mode is Mode [IV]. The differences between the local stress states near the crack surface result in the different crack growth modes, thereby causing the overall asymmetry shown in the uniaxial tension and compression loading. The critical stresses for Mode [I] crack growth is determined by the Eq. (29) and given by,

The uniaxial tensile and compressive stress-strain curves reflecting the damage characteristics of PBX9501 are compared in Fig. 4(a). The loading strain rates in loading direction are specified to be −1 for tension and compression, re11 = ±500 , ± 1000, and ± 2000 s spectively. The rate-effects on material behavior can be seen from the figure. Several characteristic parameters of the curves, including the elastic modulus, the yield point, and the peak stress, increase as the strain rate increases from 500 to 2000 s−1. The calculated compressive curve at 2000 s−1 strain-rate was compared with the experimental result obtained from the study of Gray et al. (1998). From Fig. 4(a), the current model captures the main features of the experimental result quite well, such features include the initial slope, bend over, and peak of the curve. From Fig. 4(a), in the stage of A’-B’, the tension curve shows a decrease of stress level due to the degradation of elasticity by crack growth. Then a post-peak softening (B’-C’) occurs when the mean crack size is sufficient large. Unlike the tension curve, in the stage of A*-B*, the compression curve first arrives at the yield point D* (e.g., 2000 s−1: ε11≈1.1%, σ11≈51.45 MPa). Following D*, the material behavior is the result of combined effects of softening due to crack growth and

t

=

t 11

= Scr (c¯0)/ 1

/2 (Tension)

(38a)

The critical stresses for Mode [IV] crack growth is given by, c

=

c 11

= 2Scr (c¯0)( µ2 + 1 + µ) (Compression)

(38b)

Then the ratio of compressive critical stress to tensile critical stress is obtained by,

Fig. 4. (a) Comparisons of uniaxial stress-strain responses for PBXs under tensile and compression loading at 500, 1000, and 2000 s−1 strain-rates; (b) uniaxial tensile and compression stress-strain curves for PBXs under quasi-static loading (∼10−3 s−1) from Ellis et al. (2005) experiment; (c) crack growth rate versus uniaxial strain under tensile and compression loading at 2000 s−1 strain-rate; (d) void fraction (left) and equivalent plastic strain (right) versus uniaxial strain under compression loading at 2000 s−1 strain-rate.

8

Mechanics of Materials 137 (2019) 103130

K. Yang, et al. c/ t

=2 1

/2 ( µ2 + 1 + µ)

As seen from Fig. 5(a) and (b), the stress-strain curve with 500 s−1 lateral strain rate begins with an elastic response (O-A1) and then shows a elasticity degradation resulted from crack-growth (A1-C1) followed by crack-softening (C1-D1). During O-D1 stage, both the lateral stress (σ22) and mean stress (σm=(σ11+2σ22)/3) are positive, corresponding to a tension-dominated stress state. The yield behavior of material is not considered in such tension stress state. As the lateral strain rate decreases to 400 s−1, the confined stress (-σ22) is approximate to zero (unconfined) in OeC2 stage and increases in C2-E2 stage. The stress state is compression-dominated (σm<0). The strain hardening occurs in B2eC2 stage. In C2-E2 stage, the crack-softening is interrupted at point D2 and a secondary ‘strain-hardening’ (D2-E2) response appears. For the highest level of confinement (300 s−1), the curve only shows an elastic stage (O-B3) and a strain-hardening (B3-E3); the crack-related softening doesn't occur. With increasing confining pressure, the peak stress disappears and the ductility of material is improved, namely, a brittle to ductile transition occurs. Fig. 5(c) plots the experimental stress-strain curves for PBS9501 (94 wt% sucrose and 6 wt% binder) under three different confining pressures (0.1, 34, and 138 MPa) at 10−3 s−1 quasi-static loading strain-rate (Wiegand et al., 2005). From the figure, the simulated results are consistent with the experimental curves in two aspects, (i) increasing confining pressure improves the stress level of PBXs; (ii) a change from strain softening at the lower confining pressure (0.1 MPa) to work hardening at the higher confining pressure (138 MPa), namely, brittle to ductile transition. The simulated results help to reveal the underlying micro-mechanism of this brittle to ductile transition, which is still unclear in experiments. As marked in Fig. 5(d), at 500 s−1 lateral strain rate, cracks grow in a mode of pure-shear without friction in OeC1 stage. The maximum (e.g., point C1, σ22=24.36 MPa) and minimum principle stress (e.g., point C1, σ11=-11.56 MPa) are mixed signs; and its ratio (r = σ22/σ11=2.11) is located in the range of [III] region (defined in section 2.4). Following point C1, the value of ratio (r) decreases to the [IV] region so that the growth mode of crack transfers into a mixed shear and friction mode. At 400 s−1 lateral strain rate, the confined stress (-σ22) starts to increase at point C2, thereby causing the decrease of ratio r. The crack growth mode correspondingly transfers from [IV] unstable growth state into [V] friction-locked stable state at point D2. Consequently, the axial response changes from softening (C2-D2) to a secondary hardening (D2E2). At 300 s−1 lateral strain rate, the stress state applied on crack is always in [V] friction-locked state so that non-cracks grow. In other words, higher confinement inhibits the growth of crack. In Fig. 5(e), at 500 s−1 lateral strain rate, void fraction remains unchanged because void growth is not considered in tension-dominated stress state. At 400 s−1 lateral strain rate, void collapse doesn't occur until the curve reaches a point (E2) located in the secondary hardening stage (D2-E2), corresponding that the stress triaxiality increases to a critical high-value (q=|σm|/σe=0.95). At 300 s−1 lateral strain rate, the stress triaxiality is sufficiently high that void collapse occurs immediately once the stress reaches the yield point (B3). Thus, unlike the inhibited effects of confinements on crack evolution, voids are easier to be compressed under higher confinement.

(39)

For PBX9501 material (μ=0.5, ν =0.3), the ratio of compressive to tensile critical stresses is 2.98. Eq. (39) shows the value of σc/σt increases with the friction coefficient so that asymmetry in compression and tension become more obvious for the materials with more surface cohesion. Note that the asymmetry under a simple uniaxial tension-compression loading is just a special case that describes influences of crack growth on mechanical behavior. The current model is capable of describing the difference of mechanical behavior resulted from crack growth, which is presented in arbitrary complicated loadings, for example, a triaxial tension-dominated and compression-dominated loading conditions. The prediction of the tension-compression asymmetry could be verified by the uniaxial tensile and compression stress-strain curve obtained from Ellis et al., (2005) quasi-static experiment (∼10−3 s−1), as plotted in Fig. 4(b). In their experiment, a kind of UK pressed PBXs (named as EDC 37, 91wt% HMX and 9wt% binder), whose composition and mechanical behavior are both very similar to PBX9501 (Govier et al., 2008), was studied. From the curve, the asymmetry in tension and compression was clearly observed; and the ratio of peak compressive to tensile stresses is 3.10, which is consistent with calculated results. Fig. 4(c) shows the evolution of crack growth rate under tension and compression at 2000 s−1 strain-rate. The uniaxial strain at the starting point of crack growth rate under compression (A*, 0.51%) is larger than the strain under tension (A’, 0.19%). Moreover, the rate of shear-crack growth under compression (e.g., maximum B*, 0.84cR) is slower than the rate of open-crack growth under tension (B’, 0.90cR). These two phenomena indicate the crack growth-related reasons for asymmetry in tension and compression. Fig. 4(d) shows the void fraction (left) and equivalent plastic strain (right) versus uniaxial strain under compression at 2000 s−1 strain-rate. Void collapse starts from the yield point D* and void faction only shows a slight decrease (0.008%, approximately) due to the relatively low stress triaxiality (|σm|/σe=0.33). At softening point B*, void collapse is ended because the cracks grow into a sufficient large size so that brittle failure occurs and plastic deformation is inhibited (plastic parts of residual strain, 0.53%). Under uniaxial loading, void distortion doesn't occur since that ω(σ) equals to zero under axisymmetric stress states. 4.2. Microcrack inhibited and microvoid promoted by higher confinement Both the experimental data obtained by Wiegand et al. (2005) and Bailly et al. (2011) indicated that lateral confinement strength has a strong influence on deformation and damage behaviors of PBX9501. In this section, a constant compressive strain rate is applied along axial direction ( 11 = 1000s 1). Three different positive strain rates ( 22 = 33 = 300, 400 , and 500 s−1) were specified to apply three strengths of lateral confinement on the compressive PBXs. Fig. 5(a) and (b) present axial (σ11) and lateral (σ22) stress-strain responses for PBXs under different confinements. Correspondingly, Fig. 5(d) and (e) show the evolution of crack growth rate and void fraction as a function of axial strain, respectively.

9

Mechanics of Materials 137 (2019) 103130

K. Yang, et al.

Fig. 5. (a) Axial stress and (b) lateral stress-strain responses for PBXs under triaxial loading ( 11=-1000 s−1, 22 = 33 = 300, 400, and 500 s−1 strain rates); (c) Wiegand et al.’s experimental stress-strain curves of PBS9501 under different confining pressure (0.1, 34, and 138 MPa) at 10−3 s−1 quasi-static loading strain-rate; (d) crack growth rate and (e) void fraction versus axial strain (ε11) under triaxial loading ( 11=-1000 s−1, 22 = 33 = 300, 400, and 500 s−1 strain rates).

4.3. Microvoid dominated yield under uniaxial strain compression

state. Based on the linear elasticity theory, the difference between the axial (σ11) and lateral stress (σ22=σ33) under uniaxial strain loading is expressed as σ11-σ22 = (1–2ν)/(1-ν)σ11 with Poisson's ratio ν, and shows a linear relation with the axial strain (ε11). Thus, the nonlinear variation of σ11-σ22 with ε11 [O-A path in Fig. 6(b)] results from the viscoelastic effects. Following the yield point A, the material exhibits strain-hardening response and voids are compressed [Fig. 6(b)]. According to the von-

In the early stages of PBXs under high-velocity plate impact, the stress states of materials are similar to the states under uniaxial strain loading. Fig. 6(a) presents the axial and lateral stress-strain response under uniaxial strain compression ( 11 = 1000s 1, 22 = 33 = 0.0 ). Both axial and lateral stresses show a monotonously increase and non-strain softening occurs, corresponding to the stable crack in friction-locked 10

Mechanics of Materials 137 (2019) 103130

K. Yang, et al.

Fig. 6. (a) The stress-strain responses for PBXs and (b) axial/lateral stress (left) and void fraction (right) versus axial strain (ε11) under uniaxial strain compression loading ( 11=-1000 s−1 strain rate, 22 = 33 = 0.0 ).

Mises yield theory, the yield surface under uniaxial strain loading can be simplified as σ11-σ22 = ± YM; and then the Hugoniot elastic limit (YH) is defined as YH = (1-ν)/(1–2ν) YM. In the current model, the Gurson yield surface, which considers the dependence of the yield process on pressure and porosity, can be described by Eq. (40) under uniaxial strain loading. The nonlinear relation between σ11-σ22 and ε11 is shown in the stage of A-B in Fig. 6(b). The Hugoniot elastic limit (YH) is not only dependent on yield stress of matrix (YM) but also on void fraction (f) and pressure [P=-(σ11+2σ22)/3]. 11

22

= ±YM f 2 + 1

2f cosh

3P 2YM

along with hydrostatic strain in all three stages (OeC) is plotted in Fig. 7(c). It's easy to note that in the second stage, the bulk modulus first shows an abrupt decrease and then increases to the maximum. Then the bulk modulus decreases again until it reaches to a value nearby the initial value. This phenomenon is related to the change of plastic volume strain rate and collapsed rate of void. The increase of plastic volume strain rate ( vp ) leads to the abrupt decrease of K¯ around the yield point A, which is defined as a singular point mathematically. After that, voids are compressed in a rapid rate and the capacity of material for resistance to volume change increases. When void fraction decreases to a critical value (e.g., εv =0.93%, f = 0.2% for 1000 s−1), the rate of void collapse correspondingly decreases to a slow value. Whereas, plastic volume strain rate still continues to increase so that the value of K¯ gradually decreases. In Fig. 7(c), the fluctuation of K¯ for the curve with 1000 s−1 loading hydrostatic strain rate (maximum K¯ , 10.45 GPa) is more drastic than the curve with 2000 s−1 loading rate (maximum K¯ , 8.98 GPa). Furthermore, the duration of void collapse with 2000 s−1 (13.25μs) is longer than the duration with 1000 s−1 (10.10μs). That can be attributed to the rate-dependence of the material (void inertia effects). Namely, higher loading rate leads to a higher compressive strength of the material, thereby prolonging the collapsed process of voids. Fig. 7(d) shows the hydrostatic response of PBXs under hydrostatic compression with two initial void fractions (f0 =1 and 5%). The compressive strength for 5% porosity is smaller than that for 1% porosity since that higher initial void fraction reduces the initial Gurson surface. Moreover, the concavity of curve in void collapse stage (A-B) increases with increasing initial porosity. That indicates that the material containing higher initial porosity has lower ability to resist permanent hydrostatic deformation, thereby collapsed process of voids being accelerated. The nonlinear effects of void collapse on mechanical behavior are exhibited in many other quasi-brittle materials, i.e. concretes and ceramics (Burlion et al., 2000). As for PBXs, Wiegand et al. (2011) also proposed that collapse of voids in binder is a possible mechanism for the macro-plastic deformation. While, few experimental volumetric stress-strain data of PBXs under triaxial dynamic loadings (high strainrate and high-pressure) are available for quantitively validating this trend. More experimental evidence are required to validate the calculated results under hydrostatic compression loading.

(40)

When the voids are completely collapsed (point B), porosity decreases towards zero, thereby causing the Gurson yield surface turning into the von-Mises yield surface. From path B-C in Fig. 6(b), the absolute value of σ11-σ22 shows a proportional increase with the increase of axial strain (ε11), corresponding to the strain hardening responses. 4.4. Effects of void collapse on nonlinear hydrostatic compression Porous material often exhibits a nonlinear hydrostatic response due to the effects of void evolution. Fig. 7(a) presents the effects of void collapse on dynamic hydrostatic response of porous PBXs under lowpressure (0.05∼0.1 GPa) stage, in which the cracks are stable. Furthermore, the nonlinear hydrostatic response of PBXs under highpressure (0.1∼6 GPa) is presented in Fig. 7(b). In Fig. 7(a), according to the evolution of voids, the pressure-volume response can be divided into three stages. The first stage (O-A) is linear elastic hydrostatic response with an initial slightly damaged bulk modulus. The second stage is defined as the transition stage (A-B), in which the voids are gradually compressed out of the material and plastic volumetric strain increases. The third stage (B-C) is defined as a locked stage, in which all voids have been removed from the PBXs (fully dense). In this stage, the bulk modulus is increased than the initial value in the first stage. As the pressure increases into the high-pressure stage, the bulk modulus continually increases and the hydrostatic response described by nonlinear EOS is shown in Fig. 7 (b). The evolution of averaged bulk modulus [K¯ = P / v = P /( ve + vp) ]

11

Mechanics of Materials 137 (2019) 103130

K. Yang, et al.

Fig. 7. (a) Hydrostatic stress-strain responses (left) and void fraction (right) for PBXs under low-pressure hydrostatic compression ( v =-1000 s−1); (b) hydrostatic stress-strain responses under high-pressure hydrostatic compression ( v =-1000 s−1); (c) volumetric modulus versus volumetric strain (εv) with v =-1000, 2000 s−1 strain rates; (d) pressure (left) and void fraction (right) versus volumetric strain (εv) for PBXs with initial porosity (f0 =1 and 5%) at v =-1000 s−1.

4.5. Joint contribution by microcrack and microvoid under simple shear

brittle crystals and ductile binders and debonding of crystal-binder interface. Thus in this section, joint contributions by shear-crack and void distortion on shear failure of PBXs were investigated. Fig. 8(a) and (b) show the predicted deformation response and damage-mode for PBXs under simple shear loading, respectively. From Fig. 8(a), the curve with 2000s−1 strain-rate shows the characters of elasticity degradation and crack-softening. The yield behavior of

In general, void distortion and inter-void linking are the main reasons for shear failure of ductile metallic material (Malcher et al., 2012); and shearing crack is the main reason for shear failure of quasi-brittle material (Zuo et al., 2006). The overall shear damage and failure of composite explosives is definitely the combined results of failure of

Fig. 8. (a) The shear stress-strain responses for PBXs ( 23 =2000, 3000, and 4000 s−1 strain rates); (b) crack (left) and void fraction (right) growth rate versus shear strain (ε23) under simple shear loading ( 23 =4000 s−1 strain rate). 12

Mechanics of Materials 137 (2019) 103130

K. Yang, et al.

material is not observed due to the effects of rapid-growth of cracks. As we mentioned in Introduction, the mechanical behavior of PBXs is controlled by the individual behavior of the crystal and binder components and also by their bond strength. With the loading strain rate increasing to 4000s−1, the matrix (binder) of material exhibits high plastic flow so that PBXs show a high ductility (with relatively high fracture strain). Meanwhile, the rate-dependence of energetic crystal strength results in the increase of PBXs strength. As marked in Fig. 8(b), under simple shear loading (σ23 =σ32 =τ), the stress state applied on the crack (r =-1.0) is located in [III] pureshear region. When voids start to be distorted (point B), crack growth rate has reached a high value (0.9cR); and void distortion is interrupted by crack-softening at point C. From the calculated results, distortion of voids only causes a small increase of effective void volume fraction (0.05%). The above analysis indicates rapid-growth of pure-crack plays a more significant role on failure of PBXs rather than void distortion (a potential mechanism for ductile failure). Note that shear dilatancy of the material, which results from crack opening under a pure-shear state, is not described in the current model. In simple shear, the parameter ω(σ) is unity and plastic volumetric strain rate ( vp ) is zero. The evolution law of void fraction described in the section 2.6 reduces to Eq. (41),

f = fk w

sij eijp e

= fk w

p e

(41)

is defined as where the effective plastic strain rate Then the void fraction can be integrated as, p e

f = f0 e kw

(1) microcracks are assumed to be randomly distributed within a statistically volume of material; (2) the initial number density N0 keeps unchanged over time; (3) the number density of cracks with radius larger than c decreases exponentially with the crack radius, namely, N (c, t ) = N0 exp( c / c¯ (t )) . Both the initial number density N0 and mean crack size c¯0 of PBXs strongly depend on the pressing pressure in hotpressing process. As analyzed in section 4.1, the crack-softening strain (εc) can be utilized to characterize the brittleness for PBXs under uniaxial stress compression. However, changing the values of N0 and c¯0 may affect the softening response of material, thus changing its failure mode. To investigate the combined effects N0 and c¯0 on predicted brittle behavior of PBXs, Fig. 9 shows the surface contour plotted by the log function of N0 (x-axis, N0: 0.1∼1000 cm−3), c¯0 (y-axis, 0∼40 μm), and εc (z-axis, 1%∼4%). From the figure, as the values of N0 and c¯0 decrease, the softening strain increases (e.g., from N0 = 1000 cm−3, c¯0 =40 μm, εc = 1.05%; to N0 = 0.1 cm−3, c¯0 =0.1 μm, εc = 3.89%), indicating that the material exhibits a brittle-to-ductile transition. According to Eq. (38b) mentioned in section 4.1, the PBXs material with smaller-size of initial microcracks will show a higher critical stress for crack growth, thereby causing a harder crack initiation. Moreover, according to Eqs. (6) and (7), the portion of cracking-strain in total strain will decrease for PBXs with lower initial density of cracks. Thus, with decreasing N0 and c¯0 , the contribution by crack evolution on material damage is reduced and PBXs material increasingly shows a high ductility.

p e

p e

=

2 ep. 3 23

(42)

The Gurson yield surface can be correspondingly reduced to Eq. (43), which shows that the yield of material under simple shear is a function of initial porosity, current effective plastic strain, and the yield strength of matrix. e

= (1

f ) YM = (1

f0 e kw

p e

) YM

(43)

Since that the rate of void distortion depends on the effective plastic strain rate, voids are easier to be distorted in the matrix with higherductility. In other words, as for the PBXs with higher-fraction or softer binder, increasing volume fraction of voids will be observed in damaged material. Note that the evolution of shear-crack and shear-void are coupled and interacted with each other through the current stress states. The crack and void behaviors simulated by the current model are useful for understanding the underlying damage mechanism of sheared PBXs. Further experimental observations for sheared PBXs will be useful to provide direct evidence for above predictions.

Fig. 9. The 3D surface contour plotted by the log function of initial crack number density N0 (x-axis), initial mean crack size c¯0 (y-axis), and softening strain εc (z-axis) for PBXs under uniaxial compression at 1000 s−1 strain rate.

5. Dependence of damage-mode on microstructural parameters Many factors including material composition, fabrication process, and loading conditions, could cause the variations of PBXs microstructure, i.e., grain shape and size distribution, grain/binder interface strength, and defects (Xu et al., 2008). These microstructural variations have significant influences on mechanical property, ignition mechanism, and detonation performance of PBXs. In the present model, several microstructures-related parameters (initial crack number density, initial crack size, initial porosity, and void distortion parameter) could be utilized to characterize the initial microstructures of PBXs statistically. In this section, the dependence of predicted damage response on these microstructural parameters was studied.

5.2. Effects of f0 and kw on ductile behavior As mentioned in section 4.5, compared to the rapid-growth of cracks, void distortion only has a small effect on failure for pressed PBXs under simple shear. Unlike pressed PBXs, some types of PBXs containing few microcracks generally exhibit ductile behavior due to the effects of soft polymer matrix (Asay, 2010). Undoubtedly, crack evolution plays little roles on damage for such type of PBXs. In this section, we idealized c¯0 =0.0 μm and decoupled the crack and void growth. Using a series of mathematically parameterization analysis, two basic features of the current model are discussed, (i) the effects of initial porosity on PBXs damage; (ii) the capability for describing the mechanical response of PBXs with low initial density microcracks, which mainly exhibits ductile damage. Fig. 10 displays the shear stress-strain response for PBXs with three initial porosities (f0 = 1, 2, and 5%). In the figure, the curves don't

5.1. Combined effects of N0 and c¯0 on failure mode In the SCRAM-type model (Dienes et al., 2006; Bennett et al., 1998), 13

Mechanics of Materials 137 (2019) 103130

K. Yang, et al.

6. Discussion The failure of PBXs is usually due to the combined effects of several possible micro-mechanisms; and the dominant mechanism (if any) will change with lots of factors, such as loading conditions, constituent contents, and manufacturing techniques. This paper describes a viscoelastic-plastic damage model which accounts for multiple stress-state motivated evolution-modes of microcracks (friction-locked, shear with friction, pure shear, mixed shear and open, and normal open) and microvoids (collapse and distortion). In the experimental research of Ellis et al., (2005), the damage evolution for a kind of UK pressed PBX (EDC37) was studied by Acoustic Emission (AE) method. They suggested that failure of PBX is less sensitive to the accumulation of observable microcracks under compression than tension. Our simulated results presented the influence of crack initiation on the asymmetry in tension and compression. The critical stress for initial shear-crack growth under compression σc is much higher than the critical stress for open-crack growth under tension σt. And σc/σt is a function of friction coefficient μ. In the model, contribution by the debonding at binder-crystal interface to tensile failure is modelled by cracking damage (crack open and growth) rather than void damage (void growth and coalescence). This is reasonable for pressed PBXs which fail in a nearly brittle way under tension. The void nucleation and growth should be considered if the model is applied to types of PBXs with high ductility. As noticed by Wiegand et al. (2011) and Bailly et al. (2011), a brittle-to-ductile transition occurs in PBXs under high confined pressure; and plastic behavior of material depends on confined pressure. The results revealed the effects of evolution of cracks and voids on brittle-to-ductile transition. As the lateral confinement level increases, crack growth becomes less important to failure while more voids are collapsed accompanying with a strong plastic flow. The dependence of plastic deformation on pressure is described by the decrease of initial porosity due to the collapse of voids in binder matrix. If all voids are compressed out of binder matrix, some factors related to explosive particles, such as dislocation slip or crystal twinning, maybe responsible for the further pressure-dependence of plastic flow. Consideration for the pressure dependence of yield strength and hardening modulus is the next task in the model. In the current work, the nonlinear hydrostatic response of porous PBXs under dynamic hydrostatic compression was studied. The collapse of voids in binder matrix is responsible for the nonlinearity of hydrostatic response under low-pressure compression (0.05∼0.1 GPa). As the pressure continuously increases, all voids are collapsed and the material is compacted. In the high-pressure stage (0.1∼6 GPa), the nonlinear

Fig. 10. The shear stress-strain responses for PBXs under simple shear loading ( 23 =4000 s−1 strain rate) with initial porosity (f0 = 1, 2, and 5%).

exhibit crack-softening behavior (brittle failure) anymore (referenced curve with c¯0 =30.0 μm also plotted in the figure). Instead, shape deformation and reorientation of voids contribute to material softening; and the curves show features of ductile failure (defined the critical void fraction at facture fv as 20%). As the initial void fraction increases, the fractured strain (εf) decreases (e.g., from f0 = 1.0%, εf =7.16%; to f0 = 5.0%, εf =5.05%). That indicates the material with higher initial void volume corresponds to a higher speed of void distortion (or intervoid interaction), thus speeding up the ductile fracture of material. Fig. 11(a) and (b) show the shear stress and void fraction versus shear strain with three values of kw (given f0=1.0%; c¯0 =0.0 μm) under simple shear. From Fig. 11(a), increasing kw causes a decrease of the strain at which shear localization or fracture occur (e.g., from kw = 1.0, εf =12.0%; to kw = 3.0, εf = 6.0%). From Fig. 11(b), the distorted rate of voids shows a notable increase with increasing kw, corresponding to the increasing velocity of ductile failure. In the case of kw = 0, the yield surface reduces to the original Gurson model; void-related damage never accumulates and no failure is observed in shear-dominated states. From the above, kw is a constitutive parameter related to localization of plastic flow in shear. To sum up, predicted deformation response and damage-mode of PBXs changes with the values of microstructural parameters, demonstrating the potential of the model in describing mechanical behavior for types of PBXs with different initial defects.

Fig. 11. (a) Shear stress-strain responses; (b) void fraction versus shear strain for PBXs under simple shear loading ( 23 =4000 s−1 strain rate) with void distortion parameter (kw =1.0, 2.0, and 3.0). 14

Mechanics of Materials 137 (2019) 103130

K. Yang, et al.

hydrostatic response of PBXs is related to the volume change of energetic crystal lattice and described by a porosity-independent EOS. Moreover, the inertia effects of void evolution resulted from the ratedependence of matrix was also presented, namely, voids are harder to be collapsed under higher loading strain-rate. Under a combined highpressure and high strain-rate condition (e.g., shock pulse), the localized heat will be deposited around the collapsed void boundary so that a hotspot could be induced, and thus a hotspot model related to void collapse should be developed in the next. The underlying micromechanisms of the softening, localization and fracture for PBXs in shear are still not well understood (Reaugh, 2010). Joint contributions by two potential mechanisms, namely, shear crack growth and void distortion, to shear failure of PBXs were investigated in this work. The simulated results show that growth of shear-crack plays a dominant role in failure as for PBXs with large quantities of inherent microcracks. Void distortion becomes important to failure for PBXs with higher-fraction or high-ductile binder. This prediction can be indirectly proved by the experimental observations by Ravindran et al. (2017). By applying a dynamic compressive load to a type of PBX with high-binder content (15 wt%), they proposed that microvoid evolution and void coalesce in the polymer binder caused the formation of large cracks and eventually failed the specimen. Further experimental observations for shearing PBXs will be useful to provide more direct evidence for above predictions. Moreover, mesoscale simulations should be performed to obtain more straightforward information on the evolution of micro-structures. For reactive explosives, shear localization not only contributes to the damage and failure of the material but also has significant effects on localized heating deposition, thereby influencing ignition and initiation response especially when PBXs undergo low-velocity impact (Wang et al., 2017a; Timms et al., 2018). Many hotspot mechanisms related to shear localization were proposed including shear-crack friction, shear bands, and friction between grains (Field et al., 1992; Dienes et al., 2006); however, the leading mechanism is still to be determined (Barua et al., 2013; Kim et al., 2018). In the present model, it's reasonable to suppose that shearing-voids in binder will increase their effective collective cross-sectional area by linking each other; thus promoting to the formation of micro shear-bands and leading to a high heat deposition rate inside the bands (a possible hotspot). This supposition is partly consistent with the viewpoint of Ravindran et al. (2017) concluded from experiments. Namely, shear band formation in PBX is due to the combined mechanism of thermal softening due to high strain localization in the binder and mechanical softening due to crack growth in the material. In the future, more shear localized parameters (not just kw) should be introduced to develop an ignition model; and more experimental evidences are required to verify this hotspot mechanism. In practical application, PBXs often undergo a series of dynamic complicated loadings. For example, high mechanical confinement, including charged in warheads, plate impact, and pressurized by hot gases during burning; and high shear loading, including projectile penetration, cropping, needling, and punching. In view of the limitations of mesoscale simulation, a physically-based constitutive model is still

required to predict overall damage response of PBXs charge under complicated stress states; and meanwhile to gain an insight into the micromechanisms for damage. The current study has presented evolutions of two primary microdefects (cracks and voids) in PBXs under various typical simple loadings. In the future, a numerical example for explosive charge under a more complicated stress state should be designed to verify the applicability of the model. Note that despite the model is developed for composite explosives, the model seems to be applicable for a wide class of materials with these two potential damage mechanisms, i.e. concretes and rocks, on the basis of reasonable material constants. 7. Conclusion A constitutive model incorporated with multiple stress-state motivated evolution-modes of microcracks (friction-locked, shear with friction, pure shear, mixed shear and open, and normal open) and microvoids (collapse and distortion) has been developed to study overall damage behavior of PBXs under different dynamic loading conditions. Under uniaxial stress compression, rapid-growth of crack has a dramatic effect on the failure; PBXs exhibit low ductility and little voids are compressed. With the increase of confinement level, effects of crack growth on failure decrease while more voids become collapsed accompanying with a strong plastic flow. Growths of pure-shear crack and void distortion jointly cause the failure of PBXs in shear. The simulated results show that both decreasing the values of initial crack number density or crack size will cause the material exhibiting lower brittleness and higher ductility. Under hydrostatic compression, increasing initial void fraction will reduce the initial Gurson yield surface and decrease the ability of the material to resist volumetric deformation, thereby enhancing the rate of void collapse. In simple shear, PBXs containing high initial void volume fraction and few inherent cracks exhibit a high rate of void distortion, thus speeding up the ductile fracture of material. Similarly, increasing the value of kw also causes a notable increase of void distortion rate and contributes to shear fracture. In the future work, a series of hotspot and reactive sub-models should be integrated into the model to describe damage-related possible hotspot mechanisms (e.g., shear-crack friction, void collapse, or microshear bands due to void sheeting); and to predict ignition and initiation responses of PBXs. In addition, more experimental works and mesoscale or molecule simulations are required to determine the accurate statistical parameters and confirm the merits of the model. Acknowledgements The authors would like to thank the Science Challenging Program (TZ2016001) and China National Nature Science Foundation (Grant nos. 11872119 and 11572045) and Pre-research Project of Armament (No. 1421002020101–01) for supporting this project. The authors would like to thank the PhD Nursery Foundation of Beijing Institute of Technology for supporting this project.

Appendix A. The calculation of the critical crack orientation According to the theory of the DCA model (Zuo et al., 2005, 2006), incipient failure of a brittle material, which contains a statistical, isotropic distribution of penny-shaped cracks, is assumed to occur when the crack with the critical orientation becomes unstable with two basic assumptions (i) neglecting the interactions between the cracks; (ii) the local stress states applied near a single crack is equal to the far-field Cauchy stress states (σ). Based on the Griffith energy-balance criterion, as described by Eq. (29)) in section 2.7, the crack instability surface is defined as Eq.(A.1) and is plotted in the σn–sn plane, as shown in Fig. 12.

F c ( , n , c¯) = g ( , n , c¯)

2 ¯ = L ( , n)

2 0

(A.1)

=0

where

15

Mechanics of Materials 137 (2019) 103130

K. Yang, et al.

Fig. 12. The failure surface for combined shear (sn) and normal (σn) stresses on a penny-shaped crack. Mohr circles are also shown.

0

G¯ , c¯

2 21

=

0

=

1

0

/2

(A.2) c

Given a loading state (σ), the critical crack orientation n is the one that maximizes the function L(σ, n). The critical orientation can be determined in the principal basis of the stress tensor σ.

=

i ei

ei

(A.3)

i = 1,3

where σ1, σ2, σ3 are the principal stresses (σ1≥σ2≥σ3), and e1, e2, e3 are the principal stress directions. The relation between the local stress state near a single crack (σn and sn) and the principal stresses and the crack orientation could be expressed by the Mohr circles. The point closest to the crack instability surface refers to the critical crack orientation, as shown in Fig. 12. Since the Mohr circle for the σ1–σ3 plane encloses the other two Mohr circles, the closest point to the instability surface must be on this circle, namely, the critical crack normal lies in the σ1–σ3 (the maximum and minimum principal stress) plane. The normal and shear components of the traction are expressed as, n

=

sn =

1 (( 2

1 ( 2

1

+

3)

+(

3)cos 2

1

)

(A.4)

3)sin 2

1

(A.5)

where θ denotes the angle between n and the e1 axis (0≤θ≤π/2). According the relation σ3≤σn≤σ1, cracks can be divided into three types, (i) open when σ1 > 0, (ii) closed when σ3<0, and (iii) pure-shear when σ1> 0 and σ3< 0. Accordingly, the target function L(σ, n) (mentioned in Section 2.7, Eq. (28)) to be maximized with respect to the crack orientation also have different forms related to the sign of σn,

L( ) =

1 4

{ (1 ) [ v 2

1 {(( 1 4

1

+

3

+(

3)sin 2

1

+ µ[

3)cos 2 1

+

3

]2 + (

+(

1

1

3)

2 sin2 2

3)cos 2

]) 2 },

}, n

n

>0

<0

(A.6)

By setting L ( ) = 0 , the candidates for the critical crack orientation could be solved. Moreover, the candidate for the critical pure-shear crack orientation crack is obtained from the equation, n(

s)

(A.7)

=0

The solutions of the candidate orientations, varying with the status of the critical crack (open, closed, pure-shear), are a function of biaxiality (σ1/ σ3) of the loading stress states (σ). Then, the energy-release rate g ( , c¯) in each crack growth mode is obtained using the critical orientation. Both the solutions for the critical orientation and energy-release rate are listed in Table 1 in Section 2.7.

of HMX-Estane PBXs. Comput. Mech. 52, 151–159. Bennett, J.G., Haberman, K.S., et al., 1998. A constitutive model for the non-shock ignition and mechanical response of high explosives. J. Mech. Phys. Solids 46, 2303–2322. Boddy, R.L., Gould, P.J., Jardine, A.P., et al., 2016. Damage in polymer bonded energetic composites: effect of loading rate. J. Dyn. Behav. Mater. 2 (1), 157–165. Brünig, M., Michalski, A., 2017. A stress-state-dependent continuum damage model for concrete based on irreversible thermodynamics. Int. J. Plast. 90, 31–43. Brünig, M., Michalski, A., 2018. Damage and failure at negative stress triaxialities: experiments, modeling and numerical simulations. Int. J. Plast. 102, 70–82. Burlion, N., Gatuingt, F., Pijaudier-Cabot, G., et al., 2000. Compaction and tensile damage in concrete: constitutive modelling and application to dynamics. Comput. Methods Appl. Mech. Eng. 183, 291–308. Chen, Y.L., Ghosh, S., 2012. Micromechanical analysis of strain rate-dependent deformation and failure in composite microstructures under dynamic loading conditions. Int. J. Plast. 32-33, 218–247. Clements, B.E., Mas, E.M., 2001. Dynamic mechanical behavior of filled polymers. I. Theoretical developments. J. Appl. Phys. 90, 5522–5534.

References Addessio, F.L., Johnson, J.N., 1990. A constitutive model for the dynamic response of brittle materials. J. Appl. Phys. 67, 3275–3286. Akiki, M., Menon, S., 2015. A model for hot spot formation in shocked energetic materials. Combust. Flame 162, 1759–1771. Arora, H., Tarleton, E., Li-Mayer, J., et al., 2015. . Modelling the damage and deformation process in a plastic bonded explosive microstructure under tension using the finite element method. Comput. Mater. Sci. 110, 91–101. Asay, B.W., 2010. Shock Wave Science and Technology Reference Library (Volume5): Non-Shock Initiation of Explosives. Springer, Berlin. Bailly, P., Delvare, F., Vial, J., et al., 2011. Dynamic behavior of an aggregate material at simultaneous high pressure and strain rate: SHPB triaxial tests. Int. J. Impact Eng. 38, 73–84. Barua, A., Horie, Y., Zhou, M., 2012. Energy localization in HMX-Estane polymer-bonded explosives during impact loading. J. Appl. Phys. 111, 054902. Barua, A., Zhou, M., 2013. Computational analysis of temperature rises in microstructures

16

Mechanics of Materials 137 (2019) 103130

K. Yang, et al.

Eur. J. Mech. A/Solids 27 (1), 1–17. Nahshon, K., Xue, Z., 2009. A modified Gurson model and its application to punch-out experiments. Eng. Frac. Mech. 76 (8), 997–1009. Parab, N.D., Roberts, Z.A., Harr, M.H., et al., 2016. High speed X-ray phase contrast imaging of energetic composites under dynamic compression. Appl. Phys. Lett. 109 (13), 3725–3744. Peterson, P.D., Mortensen, K.S., Idar, D.J., et al., 2001. Strain field formation in plastic bonded explosives under compressional punch loading. J. Mater. Sci. 36, 1395–1400. Ravindran, S., Tessema, A., Kidane, A., 2016. Local deformation and failure mechanisms of polymer bonded energetic materials subjected to high strain rate loading. J. Dyn. Behav. Mater. 2, 146–156. Ravindran, S., Tessema, A., Kidane, A., 2017. Multiscale damage evolution in polymer bonded sugar under dynamic loading. Mech. Mater. 114, 97–106. Reaugh, J.E., Jones, A.G., 2010. Mechanical damage, ignition, and burn: experiment, model development, and computer simulations to study high-explosive violent response (HEVR). In: Proceedings of the 14th Symposium (International) on Detonation. Coeurd’ Alene, ID. Skidmore, C.B., Phillips, D., Crane, N.B., 1997. Microscopical examination of plasticbonded explosives. Microscope 45, 127–136. Skidmore, C.B., Phillips, D.S., Asay, B.W., et al., 2000. Microstructural effects in PBX 9501 damaged by shear impact. Shock Compression of Condensed Matter-1999. AIP Publishing, pp. 659–662. Tan, H., Huang, Y., Liu, C., Geubelle, P.H., 2005. The Mori–Tanaka method for composite materials with nonlinear interface debonding. Int. J. Plast. 21 (10), 1890–1918. Timms, R., Purvis, R., Curtis, J.P., 2018. A mechanism for hot-spot generation in a reactive two-dimensional sheared viscous layer. Combust. Theo. Model. 22 (3). https://doi.org/10.1080/13647830.2018.1438670. Trumel, H., Lambert, P., Belmas, R., 2010. Mesoscopic investigations of the deformation and initiation mechanisms of a HMX-based pressed composition. In: Proceedings of the 14th Symposium (International) on Detonation. Coeurd’ Alene, ID. Wang, J., Luo, J., 2018. Predicting the effective elastic properties of polymer bonded explosives based on micromechanical methods. J. Energ. Mater. 36, 211–222. Wang, X.J., Wu, Y.Q., Huang, F.L., et al., 2016. Mesoscale thermal-mechanical analysis of impacted granular and polymer-bonded explosives. Mech. Mater. 99, 68–78. Wang, X.J., Wu, Y.Q., Huang, F.L., 2017a. Thermal-mechanical-chemical responses of polymer-bonded explosives using a mesoscopic reactive model under impact loading. J. Hazard. Mater. 321, 256–267. Wang, X.J., Wu, Y.Q., Huang, F.L., 2017b. Numerical mesoscopic investigations of dynamic damage and failure mechanisms of polymer bonded explosives. Int. J. Solids Struct. 129, 28–39. Wiegand, D.A., Reddingius, B., 2005. Mechanical properties of confined explosives. J. Energ. Mater. 23, 75. Wiegand, D.A., Reddingius, B., Ellis, K., Leppard, C., 2011. Pressure and friction dependent mechanical strength-cracks and plastic flow. Int. J. Solids Struct. 48, 1617–1629. Wu, Y.Q., Huang, F.L., 2009. A micromechanical model for predicting combined damage of particles and interface debonding in PBX explosives. Mech. Mater. 41, 27–47. Xu, F., Aravas, N., Sofronis, P., 2008. Constitutive modeling of solid propellant materials with evolving microstructural damage. J. Mech. Phys. Solids 56, 2050–2073. Xu, F., Zhao, S.D., Han, X.L., 2014. Use of a modified Gurson model for the failure behaviour of the clinched joint on Al6061 sheet. Fatig. Fract. Eng. Mater. Struct. 37, 335–348. Yang, K., Wu, Y.Q., Huang, F.L., Li, M., 2017. Numerical simulations of mechanical and ignition-deflagration responses for PBXs under low-to-medium-level velocity impact loading. J. Hazard. Mater. 337, 148–162. Yang, K., Wu, Y.Q., Huang, F.L., 2018. Numerical simulations of microcrack-related damage and ignition behavior of mild-impacted polymer bonded explosives. J. Hazard. Mater. 356, 34–52. Zhu, Q.Z., Shao, J.F., Mainguy, M., 2010. A micromechanics-based elastoplastic damage model for granular materials at low confining pressure. Int. J. Plast. 26, 586–602. Zuo, Q.H., Dienes, J.K., 2005. On the stability of penny-shaped cracks with friction: the five types of brittle behavior. Int. J. Solids Struct. 42, 1309–1326. Zuo, Q.H., Addessio, F.L., Dienes, J.K., et al., 2006. A rate-dependent damage model for brittle materials based on the dominant crack. Int. J. Solids Struct. 43, 3350–3380. Zuo, Q.H., Rice, J.R., 2008. An implicit algorithm for a rate-dependent ductile failure model. J. Appl. Phys. 104, 083526. Zuo, Q.H., Disilvestro, D., Richter, J.D., 2010. A crack-mechanics based model for damage and plasticity of brittle materials under dynamic loading. Int. J. Solids Struct. 47, 2790–2798.

Çolak, O.U., 2004. Mechanical behavior of polymers PBXW-128 and PBXN-110 under uniaxial and multiaxial compression at different strain rates and temperatures. J. Test. Eval. 32 (5), 390–395. Dienes, J.K., 1978. A statistical theory of fragmentation. In: Proceedings of the 19th US Symposium on Rock Mechanics. University of Nevada, pp. 51–55. Dienes, J.K., Zuo, Q.H., Kershner, J.D., 2006. Impact initiation of explosives and propellants via statistical crack mechanics. J. Mech. Phys. Solids 54, 1237–1275. Drodge, D.R., Williamson, D.M., 2016. Understanding damage in polymer-bonded explosive composites. J. Mater. Sci. 51, 668–679. Ellis, K., Leppard, C., Radesk, H., 2005. Mechanical properties and damage evaluation of a UK PBX. J. Mater. Sci. 40 (23), 6241–6248. Field, J.E., 1992. Hot spot ignition mechanisms for explosives. Acc. Chem. Res 25, 489–496. Gao, X., Zhang, T., Hayden, M., et al., 2009. Effects of the stress state on plasticity and ductile failure of an aluminum 5083 alloy. Int. J. Plast. 25 (12), 2366–2382. Govier, R.K., Gray, G.T., Blumenthal, W.R., 2008. Comparison of the influence of temperature on the high-strain-rate mechanical responses of PBX 9501 and EDC37. Metallurg. Mater. Trans. A (Phys. Metallurg. Mater. Sci.) 39, 535–538. Gray III, G.T., Idar, D.J., Blumenthal, W.R., Cady, C.M., Peterson, P.D., 1998. High- and low strain rate compression properties of several energetic material composites as a function of strain rate and temperature. In: Proceedings of the 11th Symposium (International) on Detonation. Snowmass, CO. Grilli, N., Duarte, C.A., Koslowski, M., 2018. Dynamic fracture and hot-spot modeling in energetic composites. J. Appl. Phys. 123, 065101. Gruau, C., Picart, D., Belmas, R., et al., 2009. Ignition of a confined high explosive under low velocity impact. Int. J. Impact Eng. 36, 537–550. Guo, T.F., Wong, W.H., 2018. Void-sheet analysis on macroscopic strain localization and void coalescence. J. Mech. Phys. Solids 118, 172–203. Guo, X., Ji, R., Weng, G.J., et al., 2014. Micromechanical simulation of fracture behavior of bimodal nanostructured metals. Mater. Sci. Eng. A 618, 479–489. Gurson, A.L., 1977. Continuum theory of ductile rupture by void nucleation and growth: path I-yield criteria and flow rules for porous ductile media. J. Eng. Mater. Tech. 99, 297–300. Hanina, E., Partom, Y., Havazelet, D., et al., 2018. Prediction of low-velocity-impact ignition threshold of energetic materials by shear-band mesoscale simulations. J. Energ. Mater 1–14. https://doi.org/10.1080/07370652.2017.1421726. Herrmann, W., 1969. Constitutive equation for the dynamic compaction of ductile porous materials. J. Appl. Phys. 40 (6), 2490–2499. Hutchinson, J.W., Pardoen, T., 2000. An extended model for void growth and coalescence. J. Mech. Phys. Solids 48, 2467–2512. Kang, G., Chen, P.W., Guo, X., et al., 2018. Simulations of meso-scale deformation and damage of polymer bonded explosives by the numerical manifold method. Eng. Anal. Bound. Elem. 96, 123–137. Kim, S., Wei, Y., Horie, Y., Zhou, M., 2018. Prediction of shock initiation thresholds and ignition probability of polymer-bonded explosives using mesoscale simulations. J. Mech. Phys. Solids 114, 97–116. LaBarbera, D.A., Zikry, M.A., 2015. Dynamic fracture and local failure mechanisms in heterogeneous RDX-Estane energetic aggregates. J. Mater. Sci. 50, 5549–5561. Le, V.D., Gratton, M., Caliez, M., et al., 2010. Experimental mechanical characterization of plastic-bonded explosives. J. Mater. Sci. 45, 5802–5813. Lecarme, L., Pardoen, T., 2011. Void growth and coalescence in ductile solids with stage III and stage IV strain hardening. Int. J. Plast. 27, 1203–1223. Liu, Z.W., Zhang, H.Y., Xie, H.M., Li, K.X., 2016. Shear band evolution in polymer bonded explosives subjected to punch loading. Strain 52, 459–466. Liu, R., Chen, P.W., 2018. Modeling ignition prediction of HMX-based polymer bonded explosives under low velocity impact. Mech. Mater. 124, 106–117. Lu, Y.L., Elsworth, D., Wang, L.G., 2013. Microcrack-based coupled damage and flow modeling of fracturing evolution in permeable brittle rocks. Comput. Geotech. 49, 226–244. Malcher, L., Pires, F.M.A., César de Sá, J.M.A., 2012. An assessment of isotropic constitutive models for ductile fracture under high and low stress triaxiality. Int. J. Plast. 30-31, 81–115. Malcher, L., Pires, F.M.A., César de Sá, J.M.A., 2014. An extended GTN model for ductile fracture under high and low stress triaxiality. Int. J. Plast. 54, 193–228. Millett, J.C.F., Taylor, P., et al., 2017. The strength of two HMX based plastic bonded explosives during one dimensional shock loading. J. Dyn. Behav. Mater. 3, 100–109. Morin, L., Leblond, J.B., Tvergaard, V., 2016. Application of a model of plastic porous materials including void shape effects to the prediction of ductile failure under sheardominated loadings. J. Mech. Phys. Solids 94, 148–166. Nahshon, K., Hutchinson, J.W., 2008. Modification of the Gurson model for shear failure.

17