Damage and hotspot formation simulation for impact–shear loaded PBXs using combined microcrack and microvoid model

European Journal of Mechanics / A Solids 80 (2020) 103924

Contents lists available at ScienceDirect

European Journal of Mechanics / A Solids journal homepage: http://www.elsevier.com/locate/ejmsol

Damage and hotspot formation simulation for impact–shear loaded PBXs using combined microcrack and microvoid model Kun Yang, Yanqing Wu *, Fenglei Huang State Key Laboratory of Explosion Science and Technology, Beijing Institute of Technology, Beijing, 100081, PR China

A R T I C L E I N F O

A B S T R A C T

Keywords: PBXs Impact–shear loading Void distortion Shear crack hotspot Void collapse hotspot

A physically–based damage–hotspot formation framework incorporating multiple stress–state motivated evolu­ tion–modes of microcracks and microvoids is developed to study the overall damage behavior of polymer­ –bonded explosives (PBXs). Localized heating sub-models of shear–crack friction and void collapse hotspot mechanisms are described to predict impact–shear ignition of PBXs. Several features of microdefect evolution under a combined shear and compression loading are predicted as follows. (i) Crack growth causes an elasticity deterioration and softening response; (ii) void distortion begins at the yield point and ends at the softening point; (iii) the softening stage will be interrupted if the increasing lateral pressure is sufficiently large to inhibit crack growth; and (iv) void collapse occurs if the lateral pressure continually increases to a critical high value. Simulated results of a punched PBX charge show that shear–crack friction heating plays a critical role in ignition under low–velocity impact (<400 m/s). Under high–velocity impact (>400 m/s), the heating due to void collapse dominates ignition because the timescale to void hotspot formation (~1 μs) is considerably shorter than that of crack hotspots (~10 μs).

1. Introduction

microdefects on the impact ignition of PBXs (Idar and Straight, 2000; Borne and Beaucamp, 2002; Parab et al., 2016; Ravindran et al., 2017). Ravindran et al. (2017) investigated the failure mechanisms for a PBX simulant under dynamic loadings and proposed that crystal frictional heating and binder debonding are two possible hotspot mechanisms. Borne and Beaucamp (2002) pointed out that the intragranular voids between explosive crystals exhibit a strong influence on the sensitivity for a type of casted PBXs. Idar and Straight (2000) performed the modified Steven test for PBX9501 and found that the impact threshold for a damaged sample shows a significant reduction relative to the pristine threshold. The defect evolution and hotspot formation in PBXs remain to be inaccurately observed in in–situ experiments due to the limitation in obtaining images at high– resolution and speed. Adequate numerical simulations are required to provide additional insights into the defect–related localized thermomechanical response in PBXs. Strong shock and mild impact (nonshock) loadings may cause initi­ ation events of PBXs (Yang et al., 2017). That is, the ignition response of PBXs and its underlying hotspot mechanisms are greatly influenced by loading conditions. If the strength of initial shock waves is sufficient (�2 GPa), then the pressure and temperature in explosives will rapidly in­ crease (~0.1 μs). Accordingly, shock wave grows into a detonation

Polymer bonded explosives (PBXs) consist of energetic crystals bonded with a polymer binder and exhibit a wide range of applications including propellants, fuels, and warheads. Microstructural defects in PBXs, such as micro–cracks and micro–voids, are unavoidable sources of microstructural heterogeneities; these defects play a crucial role in the formation of localized temperature spikes (known as hotspots) of ma­ terials (Wei et al., 2018), thereby significantly affecting the safety and survivability of PBXs in their transportation, storage, and application (Tarver and Chidester, 2005; Asay, 2010; Wang et al., 2017a). In the past decades, several defect–related hotspot mechanisms, such as, void collapse, intragranular or interfacial crack friction heating, and heating at crack tips, have been proposed (Field, 1992; Akiki and Menon, 2015; Bassett et al., 2017). However, the dominant mechanism is yet to be generally agreed upon (Picart et al., 2013). These mechanisms likely simultaneously or sequentially operate under certain loading conditions. Thus, analysis of combined microcrack and microvoid–related hotspot mechanisms is conducive to deeply understand the roles of defects on damage and ignition of PBXs. Considerable experimental studies have focused on the effects of

* Corresponding author. E-mail address: [email protected] (Y. Wu). https://doi.org/10.1016/j.euromechsol.2019.103924 Received 23 July 2019; Received in revised form 5 November 2019; Accepted 6 December 2019 Available online 9 December 2019 0997-7538/© 2019 Elsevier Masson SAS. All rights reserved.

K. Yang et al.

European Journal of Mechanics / A Solids 80 (2020) 103924

wave, namely, shock to detonation transition (SDT). The evolution of microvoids is generally regarded as the dominant hotspot mechanism for shock initiation. Bowden and Yoffe (1952) identified adiabatic compression of gas entrapped within voids as a cause of hotspot. The experimental study of Bourne and Field (1991) indicated that asym­ metrical void collapse plays an important role in ignition. During void collapse, a high–speed jet flow is initially formed at the upper wall. The jet then crosses the void and penetrates the downstream wall to induce a high–transient temperature rise. In contrast to the high–velocity shock initiation, the low–velocity impact, such as, dropping, punching, and cropping, produces a rela­ tively long–duration (~1 ms) and low–amplitude (~100 MPa) load. In this case, shear heating becomes a critical aspect in hotspot formation (Reaugh and Jones, 2010). Skidmore et al. (2000) studied the effect of shear impact on the microstructure of PBX9501 and found that chemical reaction occurs in the high–shear region. Hanina et al. (2018) consid­ ered that localized adiabatic shear banding (ASB) leads to unwanted ignition of explosives and incorporated the simulation of shear band formation at the mesoscale in an ignition model. Browning and Scam­ mon (2002) assumed that the frictional work resulting from grain–to–grain slip contributes to the spike of localized temperature for PBXs in the Steven test. The fundamental mechanisms controlling the ignition processes of mild impacted PBXs remain to be determined (Bennett et al., 1998; Zhu et al., 2009; Liu et al., 2016). In recent years, several continuum mechanics–based numerical models have been developed to predict the ignition response of PBXs (Ma et al., 2013; Liu and Chen, 2018; Gruau et al., 2009). Gruau et al. (2009) developed an ignition criterion which is a function of pressure and plastic shear rate. The specific hotspot mechanism for ignition at the mesoscale level is ignored for this type of model. Hence, a large number of explicit mesoscale simulations have been performed (Barua and Zhou, 2013; LaBarbera and Zikry, 2015; Wang et al., 2017a, 2017b; Kim et al., 2018). Barua and Zhou (2013) used a cohesive finite element method (CFEM) to quantify the energy dissipated from binder viscoelastic deformation, grain–binder interfacial debonding, grain fracture, and friction between crack faces for shocked PBXs. Wang et al. (2017b) developed a mesoscopic framework to study the thermal­ –mechanical–chemical responses of PBXs under impact loading. Kim et al. (2018) used mesoscale simulations to predict the shock initiation thresholds and ignition probability of PBXs. These types of method exhibit difficulty in predicting the ignition response of PBXs under complicated loading conditions. A physical model that presupposes one or more hotspot mechanisms is still required for elaborately describing damage and ignition behaviors of PBXs to bridge the gap between mesoscale and continuum–scale simulations. A typical physical model applied to PBXs is the statistical crack mechanical (SCRAM) model, which was originally proposed by Dienes (1978, 1983, and 1985) to model the dynamic behavior of oil shale. In the SCRAM model, an ensemble of penny–shaped microcracks is assumed to be randomly distributed within a statistically homogeneous volume of the material. The anisotropic damage is obtained by tracking the growth and coalescence of cracks with several typical orientations. Dienes et al. (2006) extended the original SCRAM model to predict the ignition of PBXs by introducing a set of submodels that account for frictional heating, melting, ignition, and burning for shear–crack. Ben­ nett et al. (1998) simplified the SCRAM model and developed the Vis­ co–SCRAM model by introducing a generalized Maxwell model to describe the binder viscoelastic effect. Liu and Chen., (2018) used the Visco–SCRAM model to predict the ignition threshold of PBXs under the Steven test. Zuo et al., (2006) developed the dominant crack algorithm, wherein damage evolution is described by the most unstable orientation for cracks in the material. Since that most of SCRAM-type models just consider an individual presupposed micro-mechanism, a physically-based damage model which incorporates both microcrack and microvoid laws is required to efficiently describe multiple potential failure and hotspot modes of PBXs.

In the present work, a physical thermomechanical model that ac­ counts for multiple microcrack and microvoid growth modes and hot­ spot mechanisms is developed to predict the defect–related damage and ignition of PBXs under complicated dynamic loadings. In Section 3, the damage modes of PBX9501 under a typical combined shear and compression loading are studied. In Section 4, the damage and hotspot formation behaviors of punched PBX samples are investigated to illus­ trate the main model features. The current work aims to provide a better understanding of the respective contributions of the two defect–related hotspot mechanisms on impact–shear ignition of PBXs. 2. Model description 2.1. Constitutive formulations A damage viscoelastic–plastic model that considers several possible evolution–modes of microcracks and microvoids is described in this section. A conceptual diagram of all kinds of considered mechanisms in the model is plotted in Fig. 1. 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 consid­ ered. In addition, the 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 is described by an isotropic-hardening ratedependent yield model. 2.1.1. Deviatoric stress-strain response _ and hydrostatic The total strain rate (ε_ ) is divided into deviatoric (e) components (_εv ). The corresponding total, deviatoric, and hydrostatic _ respectively. _ and P, stress rates are denoted as σ_ s, _ is decomposed into In the model, the total deviatoric strain rate (e) viscoelastic (e_ve ), cracking (e_cr ), and plastic (e_p ) components. The viscoelastic response of PBXs is described by a generalized Maxwell model, and the total crack strain due to an ensemble of open–crack and shear–crack is adopted from Bennett et al. (1998). The deviatoric stress–strain formulation is given by s_ ¼ 2GAðe_

e_p Þ

Bs

C

(1)

� � 1 � �2 _ P e ðnÞ Where A ¼ 1 þ α3 ðc=aÞ3 , B ¼ Aαe ac ac, and C ¼ A Nn¼1 τsðnÞ ; c is the average crack radius; a is the initial flaw size (a 3 ¼ 6Gβ); G is the shear modulus,the sliding and opening of cracks will cause an additional crack strain and a decrease of macroscopic stiffness, namely, the effective shear modulus is G ¼ AG ¼ ð1 dcr ÞG with dcr ¼ c3 =ða3 þ c3 Þ; and sðnÞ and τðnÞ are the deviatoric stress component and relaxation time for the nth Maxwell element, respectively. � 64πð1 νÞN0 3; ​ ​ ​ ​ ​ ​ if ​ P � 0 ; β¼ αe ¼ (2) 5 ν; if ​ P < 0 15ð2 νÞG where P is the pressure, N0 is the initial crack number density, and ν is the Poisson’s ratio. 2.1.2. Hydrostatic response The pressure in the porous material is a function of porosity and described as Eq. (3). Herein, the pressure–volume relation in the solid is described by the Mie–Grüneisen equation of state. � � ρs c20 ηs � Γs ηs � Pðρ; e; f Þ ¼ ð1 f Þ 1 þ Γ s ρs e s (3) 2 2 ð1 sηs Þ where ρ and e are the density and specific internal energy of the porous material, corresponding variables for the solid are expressed as ρs ¼ ρ/(1 f) with void fraction f and es ¼ e; ηs ¼ 1 ρs0/ρs; Γs is the Grüneisen coefficient; and c0 and s are material constants. 2

K. Yang et al.

European Journal of Mechanics / A Solids 80 (2020) 103924

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

The rate of pressure is then obtained by P_ ¼

_ pv

K ε_ v þ Γ s sij e_ij þ αε þ χ fkw ωðσ Þ

sij e_pij

σe

Table 1 The energy-release rate gðσ ; cÞ in each type of crack growth (Zuo and Dienes (2005)).

(4)

where the bulk modulus of the porous material K is related to that of the solid Ks by K ¼ ð1 fÞKs ; α ¼ K ð1 þ Γs ÞP; χ ¼ Ks ð1 þΓ s ÞPs and Ks � ρs ð∂Ps =∂ρs Þ þ Γs Ps . 2.1.3. Gurson yield surface The plastic deformation of porous PBXs is described by the Gurson’s theory (Gurson, 1977); and the yield function is expressed as follows: � � �2 � σe 3P Fðσe ; P; f Þ ¼ f2 1 þ 2f cosh (5) 2YM YM

YM ¼ ðσ0 þ hεp Þ½1:0 þ C1 lnð1 þ ε_ Þ� þ C2 P

(6)

[I]

Pure tension

[II]

Mixed tensile and shear

[III]

Pure shear

[IV]

Shear with friction

[V]

Friction locked

2γσ1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Scr ðcÞ= � 1 ν=2 � �σ þ σ �2 �σ 2 σ 3 �2 1 3 1 1 þ ν 2 2 2γ S2cr ðcÞ 2γσ1 σ3 S2cr ðcÞ Hμ ðσ1 σ3 Þ þ 2μσ3 γ Scr ðcÞ –

f Þ_εpv þ fkw ωðσ Þ

sij e_pij

(8)

σe

where ε_ pv is the plastic volumetric strain rate; kw is a shear state related 2

constant; ω is defined as ωðσ Þ ¼ 1 ð27J3 =2σ3e Þ with a third invariant pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi of stress J3 ¼ detðsÞ and von–Mises stress σ e ¼ 1:5sij sij . The crack–related damage fraction is defined as dcr ¼ c3 =ða3 þc3 Þ to characterize the degree of brittle fracture due to crack growth. The damage fraction of bulk material is defined as D ¼ dcr f=f0 to connect the evolution of microcrack and microvoid with damage response of PBXs at the macroscale level.

2.1.4. Evolution law of microcrack size The dominant crack algorithm (DCA) developed by Zuo et al. (2006) is utilized to describe the evolution of cracks. The basic assumption of the DCA model is that damage is related to the crack growth with the most unstable (dominant) orientation. The law of crack growth is a function of the applied energy–release rate gðσ ; cÞ and expressed as follows: 2γ 〉 gðσ ; cÞ

Energy-release rate gðσ; cÞ

f_ ¼ f_vc þ f_vd ¼ ð1

where σ0 is the initial yield stress at ε_ 0 ¼ 10 s 1; εp is the effective plastic strain; ε_ * ¼ ε_ =_ε0 is the normalized equivalent strain rate; h is the hard­ ening modulus; and C1 and C2 are the coefficients of strain rate and pressure, respectively.

c_ ¼ c_max 〈1

Remarks

collapse (f_vc ) and void distortion (f_vd ).

The yield strength of the solid material is dependent on the loading strain rate and pressure, which is defined by *

Type

2.1.6. Evolution of bulk temperature To complete the thermomechanical description, the rate of bulk temperature (Tbulk) generating from the inelastic work is given by � (9) ρcv T_ bulk ¼ kr2 Tbulk þ ϑ w_ ve þ w_ p þ w_ cr þ w_ pr

(7)

where r2 is the Laplace operator; k and cv are thermal conductivity and heat capacity, respectively; and ϑ is the heat converted fraction of the _ due to viscous amount of inelastic work. The inelastic work rates (w) effects, plastic effects, crack damage, and adiabatic compression are denoted by subscripts ve, p, cr, and pr, respectively (Yang et al., 2018).

where < > is the Macauley notation, c_max is the maximum crack growth speed; and γ is the effective surface energy. 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] Pure tension ( ð1 νÞ � r � 1); [II] Mixed tensile and shear ( 1 < r � ð1 νÞ); [III] Pure shear ( H2μ < r � 1); [IV] Shear with friction pffiffiffiffiffiffiffiffiffiffiffiffiffi 4 (Hμ � r2 ); [V] Friction locked (1 � r < H2μ ) with Hμ ¼ μ2 þ 1 þ μ. The calculation of energy-release rates gðσ ; cÞ in each crack-growth mode are listed in Table 1 (Zuo and Dienes (2005)),where γ is the effective surface qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi energy of the material, Scr ðcÞ � 2π 21 νν Gγ c

2.2. Frictional shear-crack hotspots The shear–crack friction hotspot model used in this study follows that developed by Dienes et al. (2006) and is summarized in this section. Full details of the model and its parameterization should be referenced in the literature (Dienes et al., 2006; Yang et al., 2018). The temperature rise due to the frictional heating at shear–crack faces is expressed as follows: � � ∂ ∂Tcr ρcv T_ cr ¼ (10) k þ ρQr Ze E=RTcr þ fm μv ε_ 2 ∂x ∂x

2.1.5. Evolution law of void fraction Two competing mechanisms, namely, compression causing a decrease in void fraction and shear causing an effective increase, for void evolution behavior are described by the modified Gurson model developed by Nahshon and Hutchinson (2008). The evolution of void _ is the additive of the void fraction change due to void fraction (f)

where Tcr and Qr are the shear–crack hotspot temperature and chemical heat release, respectively; Z and E are the Arrhenius reaction rate 3

K. Yang et al.

European Journal of Mechanics / A Solids 80 (2020) 103924

parameters; fm is the melting fraction; μv is the melt viscosity; and ε_ is the shear strain rate. The heating per unit area due to interfacial friction is taken as the thermal boundary condition on the crack surface, � ∂Tcr �� k ¼ μ〈 σn 〉ð1 fm Þvc (11) ∂x �x¼0

where r, θ, and ϕ are spherical coordinates. The von–Mises yield criterion determines the yield limit of the ma­ trix, as shown as follows: pffiffiffiffiffiffi σ *e ¼ 3J2 ¼ σ*r σ *θ ¼ s*r s*θ ¼ YM þ ηε_ *r (14) where s*r and s*θ are the microscopic deviatoric stress components, and η is the material viscosity that describes the rate–sensitivity of the local­ ized material (matrix) around the void (Whitworth, 2008). The microscopic deviatoric stresses can be determined by � � � 2 2b3 s*r ¼ YM þ η 3 ε_ pv ; and s*θ ¼ s*ϕ ¼ s*r 2 (15) 3 r

where vc is the crack–sliding velocity, which is calculated by the central velocity of shear crack. If the temperature at the crack interface reaches the melting point (Tm), then a solid–liquid phase change is initiated. The melting rate is _ where Q_ is the rate of heating generated from the defined as f_ ¼ Q=L, m

chemical reaction and mechanical dissipation, and L is the latent heat of melting.

In view of material incompressibility, we have b30

2.3. Void-collapse related hotspots

b3 ¼ 2 ε_ pv r

2b3 b3 ε ¼ 3 ε_ pv ; and ε_ *θ ¼ ε_ *ϕ ¼ 3 ε_ pv r r

1 1

f0 3 b fvc 0

(16)

where the two terms on the right of the equation denote the work rates due to the plastic and visco effects in the matrix around the voids, respectively. The hotspot equation from the material around the void (Tvo) is, � 2 � ∂ Tvo 2 ∂Tvo ρcv T_ vo ¼ k þ (18) þ w_ *vp þ ρQr Ze E=RTvo ∂r2 r ∂r The rate of local reaction fraction related to the void hotspot (λvo) is calculated by Z b dΛ ðr; tÞ4πr2 dλvo dt �dr (19) ¼ 4 dt b3 a3 π a 3 where the reaction rate around the void is defined by

(12)

dΛ ¼ ð1 dt

The microscopic strain rates are then calculated by _ *r

a3 ; b3 ¼

where a and b are the radii of void and matrix, respectively; and the subscript “0” is the initial state. The rate of viscoplastic work per unit volume around the void at the microscale is expressed as follows: �2 p �2 � 2YM ε_ pv 1 f0 3 4η ε_ v 1 f0 * * * _ w_ vp ðrÞ ¼ sij εij ¼ b0 þ b60 (17) 3 6 r 1 fvc r 1 fvc

Theory of Gurson’s model states that void–matrix aggregate is idealized as a unit–cell, namely, a finite, spherical, incompressible ma­ trix containing a concentric, spherical void, as shown in Fig. 2(a). The stress and strain in the matrix are non–uniform and defined as the microscopic stress (σ*) and strain (ε*), respectively. By contrast, the uniform, effective, compressible medium of the homogenized void–matrix is referred to the macroscale, as shown in Fig. 2(b). The average stress and strain over the entire cell are defined as the macro­ scopic stress (σ) and strain (ε), respectively. The microscopic velocity field (v*) is divided into two parts, namely, shape change at a constant volume (vs)* and volume change at a con­ stant shape (vv)*. In the present model, we mainly focused on the heat generated from the collapse of the spherical symmetric void, in which the volume of spherical voids will decrease and its shape will remain unchanged. Only the volume change part of the microscopic velocity field (vv)* is considered. The microscopic velocity field (v*) must satisfy the external boundary condition in terms of the macroscopic volumetric plastic strain � rate (_εpv ) by v*r �S ¼ b_εpv . The non–vanishing component of v* is expressed as follows: v*r

a30 ¼ b3

ΛÞZe

E=RTvo ðr;tÞ

(20)

The proposed models are implemented as subroutines in the Drexh finite element code. The Drexh code uses a forward Euler integration scheme and that is integrated through time by using small time in­ crements. A typical pressed explosive named PBX9501 (95 wt% HMX,

(13)

Fig. 2. (a) Unit–cell of a finite, spherical matrix containing a concentric void; (b) uniform, effective medium of the homogenized matrix–voids at the macroscale. 4

K. Yang et al.

European Journal of Mechanics / A Solids 80 (2020) 103924

2.5 wt% estane, and 2.5 wt% BDNPA/F) is used in the current work. The parameters related to physical properties (i.e., density and modulus) and crack and void damage model for PBX9501 are adopted from (Yang et al., 2018, 2019; Bennett et al., 1998; Dienes et al., 2006) and listed in Tables 2(a) and 2(b). These parameters are verified using uniaxial compression stress–strain curves (Yang et al., 2018) and multiple shock experimental data (Dienes et al., 2006). The parameters related to shear–crack and void–collapse hotspots (η) are selected from the work of Dienes et al. (2006) and Whitworth (2008), respectively.

dramatic decrease in stage C1–D1 because the increasing lateral pressure inhibits the crack growth. At point D1, c_ arrives at a low value (0.08cmax); thus, the stress level begins to increase. At point E1, the increasing lateral pressure is sufficiently large to inhibit the crack growth entirely. The mode of crack growth transfers from the mixed shear and friction into the friction locked state. The strain at point D1 for the curve with a low compressive strain rate (500 s 1) is increased because the low lateral pressure exhibits a weak effect on the crack growth inhibition. As shown in Fig. 3(c), the void fraction under simple shear begins to increase at yield point B0 and ends at softening point C0. This result implies that void distortion is accompanied with the yield of material and is interrupted by crack–related softening. Under a combined shear and compression loading (1000 s 1), void distortion occurs in stage B1–C1. Void collapse initiates at point E1, which is the result of high lateral–pressure compression. Stress states are shear–dominated at a lower compressive strain rate (500 s 1). That is, the lateral pressure is insufficiently high, thereby preventing void collapse. Under uniaxial strain compression, stress states are compression–dominated; thus, voids are non–distorted. Voids are compressed once the material yields at point B2. At point G2, the voids are totally collapsed and porosity de­ creases toward zero. Fig. 3(d) shows the relation of macroscale damage evolution and its underlying micro–defects. The damage factor defined byD ¼ c3 f=ða3 þc3 Þf0 continually increases under simple shear. Thus, shear–crack growth and void distortion contribute to material damage. The appearance of a platform in D1–E1 for the curve of combined shear and compression (1000 s 1) corresponds to that where neither crack growth nor void distortion occurs. In point E1, the decrease in damage factor is attributed to the occurrence of void collapse. Under uniaxial strain compression, crack growth is inhibited, and voids are compressed. Accordingly, the damage factor show no increase and stays around the initial value (D0 ¼ 2.69 � 10 5).

3. Damage responses of PBXs under combined shear and compression In non–shock loading scenarios, such as, dropping, punching, and needling, PBXs often undergo a series of complicated shear and compression stress states (Chen et al., 2018). In this section, deformation and damage behaviors of PBXs under a combined shear and compression loading were investigated. A constant shear strain rate was specified in plane (_ε23 ¼4000 s 1) and two compressive strain rates (_ε11 ¼ 500 and 1000 s 1) were applied on the shearing PBXs. Fig. 3(a) presents the calculated effective stress–strain curves for PBXs under combined shear and compression loading. The stress–strain curves under simple shear (non–compression) and uniaxial strain compression (com­ pression–dominated) loadings were also plotted in Fig. 3(a). Fig. 3(b)– (d) show the evolution of crack growth rate, void fraction growth rate, and damage factor, respectively, as a function of effective strain under different loading conditions. As shown in Fig. 3(a), the stress–strain curve under simple shear begins with an elastic response (O–A0) and then shows an elasticity deterioration caused by crack–growth (A0–B0). The material begins to yield at point B0 and then enters a stage (B0–C0) in which the material behavior depends on the combined effects of softening due to crack growth and hardening due to plastic deformation. At point C0, the mean crack size is sufficiently large that a post–peak softening (C0–D0) occurs. Two significant differences between the curves under simple shear and combined shear and compression (1000 s 1) are observed after soft­ ening point (C1); that is, (i) the softening stage is interrupted at point D1; and (ii) the stress level increases in stage D1–E1–F1. The strain at inter­ rupted point (D1) increases with the decrease in compressive strain rate to 500 s 1. The stress level of the curve under uniaxial strain compres­ sion (O–B2–G2–F2) continually increases. The non–softening stage ap­ pears because the crack growth is inhibited by high lateral confinement. As shown in Fig. 3(b), the crack initiates at point A0 (εeff ¼ 0.14%) under simple shear and then rapidly grows in stage A0–B0. At point C0, _ reaches the maximum, which corresponds to the crack growth rate (c) the occurrence of softening. In stage C0–D0, the crack still grows rapidly (>0.8cmax) although the value of c_ shows a slight decrease. Under combined shear and compression loading (1000 s 1), c_ exhibits a

4. Damage and ignition responses for rod-punched PBXs A typical impact–shear (punch) loading was selected for modeling the typical stress states and wave propagation conditions experienced in explosive charge during accidental ignition scenarios. The spatial and temporal distribution of microcrack– and microvoid– related damage and ignition for PBXs under such complicated impact and shear loading was investigated in this section. A flat–fronted steel rod (Φ20�25 mm) normally impacts an uncov­ ered cylindrical explosive charge (Φ60�50 mm; Fig. 4(a)). Due to the axial symmetry of the charge and rod, only half of the geometry model was constructed by using a 2D finite element model. The displacements in the x-direction of the nodes on the symmetrical axis were constrained. The rod is discretized with 2000 quad elements. The explosive charge is

Table 2 (a) Parameters for physical properties of PBX9501 (Yang et al., 2018; Bennett et al., 1998; Dienes et al., 2006). ρs0 1860 kg∙m

3

G

ν

μ

cv

k

3.2 GPa

0.36

0.5

1000 J∙kg 1∙K

1

0.406 J∙m 1∙s

1

∙K

1

σ0

C1

C2

h

c0

N0

12.0 MPa

0.76

0.01

500 MPa

30 μm

30 cm

γ 3

50 J∙cm

2

f0

kw

0.01

0.5

(b) Parameters for viscoelastic model (units: MPa, s) (Bennett et al., 1998) G(1)

G(2)

G(3)

G(4)

G(5)

944 1/τ(1) 0

173.8 1/τ(2) 7.32 � 103

521.2 1/τ(3) 7.32 � 104

908.5 1/τ(4) 7.32 � 105

687.5 1/τ(5) 2.00 � 106

(c) Parameters for Mie–Grüneisen EOS and hotspot model (Yang et al., 2018; Dienes et al., 2006; Whitworth, 2008). c0 2500 m/s

s 2.26

Г 1.50

Tm 519 K

μv 0.2 Pa∙s

L

Qr 5

2.08 � 10 J∙kg 1

2.09 � 10 J∙kg 1

5

E 6

2.25 � 10 J∙mol 1

Z 5

19

5.0 � 10 s 1

a0

η

ϑ

10 μm

100 Pa∙s

0.9

K. Yang et al.

European Journal of Mechanics / A Solids 80 (2020) 103924

Fig. 3. (a) Effective stress–strain responses; (b) crack growth rate; (c) void fraction growth rate; (d) damage factor for PBXs under combined shear and compression loading (_ε11 ¼4000 s 1, ε_ 11 ¼ 500 and 1000 s 1 strain rates) versus effective strain.

Fig. 4. (a) Finite element model of rod–punched loading; (b) the histories of pressure at region A (left axis) and crack hotspot temperature at region B (right axis) with mesh sizes of 0.25, 0.5 and 0.75 mm with 200 and 600 m/s impact velocities.

discretized with three different mesh sizes of 0.25, 0.5, and 0.75 mm (triangular elements number: 24000, 12000, and 8000) to check the mesh independence. Since that crack hotspot temperature is calculated related to the crack growth rate and heat conduction algorithm, this variable is assumed to be more sensitive to mesh size. Fig. 4(b) displays the his­ tories of pressure at region A (200 m/s) and crack hotspot temperature at region B with different mesh sizes (200 and 600 m/s). Note that three

pressure curves exhibit insignificant difference. Compared to 200 m/s impact velocity, the curves between for 0.75mm and 0.5mm mesh size with 600 m/s shows a more distinct difference. The difference between the time when crack hotspot runs away with 0.5mm and 0.25mm mesh size for two impact velocities are not very significant (3.7% for 200 m/s and 4.8% for 600 m/s). That indicates that 0.5 mm mesh size is suffi­ ciently small to get reasonable results. Moreover, since that heat con­ duction equations for hotspots (Eqs.(10) and (18)) are solved using the 6

K. Yang et al.

European Journal of Mechanics / A Solids 80 (2020) 103924

forward difference algorithm with sub time-steps, the actual CPU time needed for each time-step (0.5mm) was relatively long (~1.0s). Increasing the mesh size from 0.25 to 0.5mm will significantly reduce the total needed CPU time. Thus, the explosive charge is discretized with 0.5 mm mesh size in the following calculations. The Johnson–Cook strength model was utilized for steel rod, and its parameters were adopted from (Yang et al., 2017) and listed in Table 3. The limited integration time increment is estimated based on the char­ acteristic length of the element (l) and sound velocity of the material (C), where Δtlim¼θl/C with the courant sound velocity speed fraction θ (~0.5). The real time-step for the calculation process is selected as Δt¼0.1Δtlim¼1.0 ns for the accuracy of the calculation. Two regions, labeled as region A and B, within the sample were selected to efficiently describe the following analysis (Fig. 4(a)). Regions A and B are located at the center of the sample and the junction between the side of the rod and the sample, respectively.

severe crack damage regions propagate from region B to approximately the entire sample, which corresponds to the travel of side rarefaction. Under 600 m/s, the maximum damage fraction and the occupied area of severe damage at the rod periphery show an increase due to the driving effect of high shear stress (Fig. 7). The high–strength incident waves (2.0 GPa) causes a larger displacement beneath the rod front. Thus, the hemispherical friction–locked region nearly disappears at 15 μs. Fig. 8(a) and (b) show the shear strain field of a PBX9501 inert simulation material (90 wt% KCL, 10 wt% binder) under low-rate punch loading (2 mm/min) obtained from Liu et al.‘s test (Liu et al., 2016). The Liu et al.’s test results provide several qualitative evidence for simulated results. (i) A hemispherical region with low shear-strain is formed beneath the impact surface (Fig. 8(b)). (ii) Two concentration regions with high shear strain initiate simultaneously following the boundary of the hemispherical region (Fig. 8(a)). (iii) The test results show that the material will finally crack along these two high shear regions. The simulated results capture the main features of tests including the for­ mation of hemispherical region and high-shear crack region. The open-crack damage predicted by the simulated results is not observed in Liu et al.‘s test because the side rarefaction effects are insignificant during quasi-static loading.

4.1. Comparisons of microcrack and microvoid related damage Microcracks and microvoids are two main microstructural defects influencing the damage responses in PBXs. The evolution modes of crack or void will change with the varying of impact strengths. In addition, the modes of these two microdefects will exhibit a significant difference under the same impact condition. Thus in this section, we made a comparison between the evolutions of cracks and voids under different impact velocities.

4.1.2. Microvoid related damage under low and high impact The combined shear and compression stress state at region B will create two competing mechanisms namely, compression causing a decrease in void fraction (collapse) and shear causing an increase (distortion), for void evolution. Fig. 9(a) and (b) show the void growth rate histories under 200 and 600 m/s impact velocities, respectively. In Fig. 8(a), the void growth rate at region B initially shows a positive value during 0.05–0.12 μs, indicating that distortion is the dominant behavior of voids. Void collapse is then dominated during 0.12–0.65 μs. Mean­ while, only void collapse occurs at region A. The maximum of void collapse rate at region A ( 8.45 � 104 s 1) is nearly twice the rate at region B ( 4.3 � 104 s 1) due to the effect of a stress state being close to hydrostatic compression. As shown in Fig. 9(b), the maximum of void fraction growth rate at region B under 600 m/s (0.4 � 104 s 1) is more than twice the rate under 200 m/s (0.18 � 104 s 1), because the sample at region B undergoes a higher shear stress under 600 m/s. The durations of void collapse at region A and B under 600 m/s are shorter than that under 200 m/s because a high–strength compression will result in a high rate of void collapse. In Figs. 6 and 9, an important feature observed from the evo­ lution curves of crack and void growth rate needs to be given attention. That is, under the same impact velocity and region, the timescale for crack growth (~10 μs) is two orders higher than that for void growth (~0.1 μs). Figs. 10 and 11 present the contours of microvoid fraction change (Δf) under 200 and 600 m/s impact velocities, respectively, at 5, 10, 15, and 20 μs. In Fig. 10, void collapse initially occurs beneath the impact surface. At 10 μs, the area undergoing a void collapse is mainly distributed at an elliptically shaped region near the impact surface and propagates after the travel of incident wave front. The decrement of void fraction at the periphery of the rod front is smaller than that at the elliptically shaped region. This condition is attributed to the high–shear stress state at the rod periphery, which contributes to void distortion and blocks the decrease of void fraction. Under 600 m/s (Fig. 11), the impact strength of incident shock wave is sufficiently high that a large area of material at the impact surface undergoes void collapse. The shear stress plays a less important role in the change of void fraction. After the propagation of incident wave front, a teapot–shaped area undergoing void collapse is formed. At 20 μs, nearly the total sample undergoes void collapse.

4.1.1. Microcrack related damage under low and high impact Crack growth can be divided into five modes on the basis of the biaxiality of applied stress states (σ3/σ1); these modes are [I] pure ten­ sion, [II] mixed tensile and shear, [III] pure shear, [IV] mixed shear and friction, and [V] friction locked (Zuo et al., 2006). Fig. 5(a)–(b) present the crack growth rate histories under two typical impact velocities of 200 and 600 m/s. As shown in Fig. 5(a), the applied shear stress on the crack surface (sn) can sufficiently overcome the interfacial friction (μσn) given the high–level stress concentration at region B. Accordingly, the cracks become unstable and rapidly grow in Mode [IV]. At 21 μs, the crack no longer grows because damage fraction of material arrives at a high value (dcr >0.995), indicating that the material is entirely fractured. In contrast to region B, the cracks at region A are stable and fric­ tion–locked during 0–24 μs. In this duration, a stress state close to hy­ drostatic compression forms at region A. Accordingly, the applied shear stress is always less than the interfacial friction (sn � μσ n). After 24 μs, the rarefaction waves from the rear of the rod arrive at region A, thereby contributing to the growth of open–crack. Fig. 5(b) illustrates that the crack growth mode at region A under 600 m/s (Mode [IV]) is different from that under 200 m/s (Mode [I]). The time to crack initiation under 600 m/s is 19 μs, which is 5 μs earlier than that under 200 m/s. The material at region A under 600 m/s un­ dergoes a larger vertical displacement compared with that under 200 m/ s. Thus, in comparison with the rarefaction from the rear of the rod, the side rarefaction generated at the periphery of the rod front plays a more important role in the transfer of stress state from Mode [V] to Mode [IV]. Figs. 6 and 7 show the contours of the microcrack–related damage fraction under 200 and 600 m/s at 5, 10, 15, and 20 μs. Severe crack damage (dcr >0.75) is mainly distributed at region B at 5 μs (Fig. 6). Then, a hemispherical region beneath the impact surface forms at 15 μs. The cracks in such a hemispherical region are friction–locked. At 20 μs, Table 3. Parameters of Johnson–Cook strength model for steel rod. Density (g/cm3)

Shear modulus G (GPa)

Yield stress A (MPa)

B

C

n

m

7.8

80

410

539

0.05

0.6

1.0

4.2. Thermal response related to microcrack and microvoid The microstructural defects in PBXs including cracks and voids not 7

K. Yang et al.

European Journal of Mechanics / A Solids 80 (2020) 103924

Fig. 5. Evolution of microcrack growth rate with (a) 200 and (b) 600 m/s impact velocities.

Fig. 6. Microcrack damage contours with 200 m/s impact velocity at (a) 5.0, (b) 10.0, (c) 15.0, and (d) 20.0 μs.

Fig. 7. Microcrack damage contours with 600 m/s impact velocity at (a) 5.0, (b) 10.0, (c) 15.0, and (d) 20.0 μs.

only lead to the macroscopic failure of material but also exacerbate energy localization, thereby significantly affecting hotspots formation and ignition of PBXs. It is required to study the localized thermal re­ sponses related to crack and void evolution, and to determine which one plays a dominant role under different loading conditions.

PBXs but also exhibits significant influences on the localized heating deposition, thereby influencing their ignition and initiation responses. In this section, evolution and distribution of localized temperature generated by surface friction of shear–crack (Tcr) is delineated. Fig. 12 (a) shows the history of Tcr at the riskiest location (region B) under a large range of impact velocities. Figs. 13 and 14 show the contours of Tcr under 200 and 600 m/s, respectively, at 5, 10, 15, and 20 μs. In Fig. 12(a), the value of Tcr rapidly rises at the initial impact

4.2.1. Microcrack friction hotspots under low and high impact The evolution of crack and void not only contributes to the failure of 8

K. Yang et al.

European Journal of Mechanics / A Solids 80 (2020) 103924

Fig. 8. The shear strain field obtained from the PBX9501 inert simulation low-rate punch test of Liu et al. (2016) using DIC method after (a) 6s and (b) 9s loading.

Fig. 9. Evolution of microvoid growth rate with (a) 200 and (b) 600 m/s impact velocities.

Fig. 10. Change of microvoid fraction (Δf) contours with 200 m/s impact velocity at (a) 5, (b) 10, (c) 15, and (d) 20 μs.

because high–shear stress at region B contributes to the unstable growth and frictional heating of shear–crack. The curve of Tcr for a relatively low–impact (e.g., 125 m/s) shows a peak point, indicating that heat generation reaches a maximum, and heat conduction begins to control the following behavior. With impact velocity increasing to 150 m/s, the peak point disappears and Tcr reaches a plateau corresponding to the melting process. After melting, the continual rise of temperature is caused by chemical reaction and intense viscous heating inside the molten crack surface. The runaway feature of the temperature curve demonstrates that the critical velocity of crack hotspot formation (vcr) is determined as 125–150 m/s. Fig. 12(b) plots the variations of time to crack hotspot formation (tcr) with impact velocity at region B. The value of tcr initially presents a

dramatic decrease as the impact velocity increases from the critical value to 200 m/s. Subsequently, tcr decreases toward a plateau when the impact velocity increases to a relatively high value (e.g., 400 m/s). The high intensity impact contributes to the deposition of viscous heating inside molten crack surface; hence, less time is required to form a crack hotspot. In Fig. 13, the area with a relatively high value of Tcr (�400 K) is mainly gathered around the periphery of the rod front at 5 μs under 200 m/s. Crack hotspot (Tcr � 750 K) forms at this area at 10 μs. Simulta­ neously, the area undergoing melting (Tcr� 519 K) forms a semi–ring green zone beneath the rod front. A wedge–shaped zone located be­ tween the rod front and melting zone exhibits few temperature increases (ΔTcr � 0 K) caused by a friction–locked crack. At 15 μs, the area with a 9

K. Yang et al.

European Journal of Mechanics / A Solids 80 (2020) 103924

Fig. 11. Change of microvoid fraction (Δf) contours with 600 m/s impact velocity at (a) 5, (b) 10, (c) 15, and (d) 20 μs.

Fig. 12. (a) Evolution of microcrack–related hotspot temperature with different impact velocities; (b) variations of time to crack hotspot formation (tcr) with impact velocity at region B.

Fig. 13. Microcrack–related hotspot temperature contours with 200 m/s impact velocity at (a) 5.0, (b) 10.0, (c) 15.0, and (d) 20.0 μs.

crack hotspot formation is spread along the boundary of the wedge zone and forms a small–radius semi–ring red zone, which is defined as a transition zone. The cracks in the transition zone experience severe frictional heating caused by high shear stress. The occupied area of melting and hotspot formation zone expands as impact velocity increases to 600 m/s (Fig. 14) because the high shear stress accelerates the temperature increase of shear–crack hotspot. As mentioned in Sec. 4.1.1, the area of wedge–shaped friction–locked

region is squeezed. Accordingly, the hotspot formation zone converges at the rod front. 4.2.2. Microvoid collapse hotspots under low and high impact Fig. 15(a) and (b) show the time–curves of the void collapse hotspot temperature (Tvo) at regions A and B under different impact velocities, respectively. In Fig. 15(a), Tvo at region A initially shows a distinct rise at ~0.1 μs, which corresponds to the accumulation of viscoplastic heating 10

K. Yang et al.

European Journal of Mechanics / A Solids 80 (2020) 103924

Fig. 14. Microcrack–related hotspot temperature contours with 600 m/s impact velocity at (a) 5.0, (b) 10.0, (c) 15.0, and (d) 20.0 μs.

Fig. 15. Evolution of microvoid–related hotspot temperature at regions (a) A and (b) B with different impact velocities.

generated from void collapse. Heat conduction is crucial to the change of Tvo after the peak point when the impact velocity is relatively low (<500 m/s). With the increase of the impact velocity to 550 m/s, the heat generated from viscoplastic deformation of voids and chemical reaction is larger than the heat loss due to conduction, thereby resulting in the continual rise and final runaway of Tvo. In Fig. 15(b), the increase of Tvo at region B is considerably higher than that at region A under the same impact. Under 400 m/s, the maximum of Tvo at high–shear region B (462 K) is 30 K higher than the value at high–compression region A. That is attributed to the high shear

stress at region B, which causes void distortion, thereby reducing the yield strength and increasing the rate of heat generation due to void collapse. This outcome points to a new shear–related hotspot mecha­ nism, namely, shearing–voids. At low stress triaxiality, shearing–voids will increase the effective collective cross–sectional area by linking each other and promote to the formation of micro shear–bands, which results in a high heat deposition rate inside the bands. In the future, experi­ mental evidence is required to verify this hotspot mechanism. Similar to the determination of vcr, the critical velocities of void hotspot formation at regions A (vvo–A) and B (vvo–B) are 500–550 m/s and 300–400 m/s,

Fig. 16. Microvoid–related hotspot temperature contours with 200 m/s impact velocity at (a) 5.0, (b) 10.0, (c) 15.0, and (d) 20.0 μs. 11

K. Yang et al.

European Journal of Mechanics / A Solids 80 (2020) 103924

respectively. Figs. 16 and 17 show the contours of Tvo under 200 and 600 m/s at 5, 10, 15, and 20 μs. In Fig. 16, the area with a relatively high value of Tvo (�400 K) is beneath the rod front at 5 μs under 200 m/s. Void hotspot (Tvo � 750 K) doesn’t form in the entire time due to insufficient impact strength. The area undergoing a distinct rise (ΔTvo � 50 K) forms an elliptical green zone and propagates after the travel of incident wave front. The area with a high rise of Tvo is concentrated on region B at a high impact velocity of 600 m/s (Fig. 17). As previously mentioned, void hotspot initially forms around region B due to the indirect effect of void distortion on heat generation. At 20 μs, the area with ΔTvo � 50 K spread to nearly the total sample, which is consistent with the void fraction change distribution.

Figs. 19 and 20 present the contours of Ttot under 200 and 600 m/s at 5, 10, 15, and 20 μs? Under 200 m/s, the distribution of maximum Ttot is consistent with the distribution of Tcr (Fig. 13), indicating the dominant role of Tcr on total temperature rise. The area undergoing a distinct temperature rise (Ttot � 400 K) becomes larger compared with the contours of Tcr due to the effect of Tvo. Under 600 m/s, the increases of Tcr, Tvo, and Tbulk show a maximum at region B. Thus, region B is the riskiest location where ignition occurs. At the top center of fric­ tion–locked wedge zone (region A, ΔTcr � 0 K), the total temperature shows a high value (Ttot � 700 K) due to the rapid collapsed rate of voids. Hence, the void collapse hotspot dominates ignition of the sample at this region. The simulated temperature distribution can be supported by the experimental results of Skidmore et al. (2000) (Fig. 21(a)) and Clements (2011) (Fig. 21(b)). In their tests, a plunger impacted into the PBX9501 sample from the left side of the figure. The micrograph shown in Fig. 21 (a) provides evidence that shear crack initiates and melting or reaction ever occurs around the friction–locked zone under a relatively low impact velocity (100 m/s). Void hotspot will form in the friction–locked zone under a critical impact velocity, as shown in Fig. 21(b). In this case, both void and crack hotspots exhibit strong influences on ignition of PBX, which demonstrates the predictions of simulated results.

4.2.3. Contributions of two hotspot mechanisms to ignition To delineate the respective contributions of localized heating due to crack friction and void collapse and bulk homogenous heating to igni­ tion, a variable named total temperature of sample (Ttot) is defined as functions of Tcr, Tvo, and ΔTbulk, Ttot ¼ MaxðTcr þ ΔTbulk ; Tvo þ ΔTbulk Þ

(21)

Fig. 18(a) shows the evolution of Ttot in a relatively long duration (0–40 μs) at region B under 125, 150, and 600 m/s. The behavior of Ttot is dominated by Tcr in a relatively large timescale (~10 μs) under the low–velocity impact (e.g., 125 and 150 m/s). The difference between the curves of Ttot and Tcr characterizes the effects of bulk temperature rise. A detailed view of Fig. 18(a) at the early stage (0–2 μs) is presented in Fig. 18(b). From the figure, the increase rate of Tvo is larger than that of Tcr under 150 m/s and during 0–1.0 μs. Thus, Tvo dominates the increase of Ttot. After 1.0 μs, heat conduction inhibits the continual rise of Tvo. Accordingly, Tcr begins to dominate. Under 600 m/s, Tvo runs away at 0.48 μs; hence, the sample will be immediately ignited, although Tcr also runs away at 7.0 μs (Fig. 18(a)). Thus, the crack hotspot dominates ignition under low–velocity impact (<400 m/s), whereas the void hot­ spot becomes dominant under high–velocity impact (>400 m/s). The evolution of bulk temperature increment (ΔTbulk) is shown in Fig. 18(c). The maximum of ΔTbulk is 86 K under a relatively low impact velocity of 200 m/s, which is insufficiently high to ignite explosives. The maximum of ΔTbulk shows a rapid increase when the impact velocity increases to a high value (i.e., 249 K and 525 K for 400 m/s and 600 m/s, respectively). This result corresponds to an accelerated accumulation of inelastic work rates. Bulk homogenous heating becomes increasingly important to ignition as impact strength increases. However, the local­ ized temperature increase due to heterogeneous defects (voids) is still a dominant factor on ignition behavior given that the timescale of Tvo increase under 600 m/s (~1 μs) is considerably smaller than that of bulk temperature rise (~10 μs).

5. Conclusions In the present work, a physical model that can describe multiple growth modes of cracks, collapse and distortion of voids, crack frictional heating, and viscoplastic heating during void collapse, is developed for PBXs under complicated dynamic mechanical insults. The study aims to delineate the respective contribution of cracks and voids to the damage and ignition responses of PBXs. The simulated results show that crack growth causes elasticity deterioration and softening response of PBXs under a combined shear and compression loading with 1000 s 1 compressive strain rate. Void distortion begins at the yield point and ends at the crack–related soft­ ening point. The softening stage will be interrupted when the increasing lateral pressure is sufficiently large to inhibit crack growth entirely. Void collapse occurs as lateral pressure continually increases to a critically high value. The lower lateral pressure at a low compressive strain rate (500 s 1) presents a weak effect on the inhibition of crack growth. Accordingly, the strain at the interruption point of the softening stage is increased, and no void collapse occurs. The effect of microcrack and microvoid evolution on the damage and ignition behaviors of punched PBXs are quantified. Crack is fric­ tion–locked at a wedge–shaped region beneath the rod front (region A). By contrast, void collapse mainly occurs at this region due to high­ –confinement. At the sides of the wedge zone, high shear stress results in

Fig. 17. Microvoid–related hotspot temperature contours with 600 m/s impact velocity at (a) 5.0, (b) 10.0, (c) 15.0, and (d) 20.0 μs. 12

K. Yang et al.

European Journal of Mechanics / A Solids 80 (2020) 103924

Fig. 18. At the riskiest region B with different impact velocities to predict the (a) evolution of total temperature during 0–40μs; (b) detailed view of the total temperature history at the early stage (0–2μs); (c) evolution of bulk temperature increment.

Fig. 19. Total temperature contours with 200 m/s impact velocity at (a) 5.0, (b) 10.0, (c) 15.0, and (d) 20.0 μs.

unstable crack growth and void distortion. Shear–crack friction heating plays a critical role in ignition under low–velocity impact (<400 m/s). Under high–velocity impact (>400 m/s), the heating due to void collapse dominates ignition because the timescale of void hotspot for­ mation (~1 μs) is shorter than the timescale of crack hotspots (~10 μs). The present study only separately examines the effect of shear–crack and void collapse heating on ignition. In the future, the interactions between the two types of defect–related hotspot should be considered. The other possible hotspot mechanisms, i.e., adiabatic shear bands, should be incorporated in the model to better understand the role of

each hotspot mechanism on ignition. Additional experimental evidence of punched PBXs is needed to fully validate the proposed model.More­ over, how the simulation trends would change if the composition of the material changes should be studied to reveal the effects of material type on its mechanical behavior and micro-defect evolution. Declaration of competing interest We have no conflict of interest to declare.

13

K. Yang et al.

European Journal of Mechanics / A Solids 80 (2020) 103924

Fig. 20. Total temperature contours with 600 m/s impact velocity at (a) 5.0, (b) 10.0, (c) 15.0, and (d) 20.0 μs.

Fig. 21. (a) Micrograph showing the evidence of shear crack hotspot formation along the high shear region for 100 m/s plunger-impacted PBX 9501, which is obtained from Skidmore et al. (2000); (b) micrograph showing the evidence of void and crack hotspot formation for plunger-impacted PBX, which is obtained from Clements (2011).

Acknowledgements

Browning, R.V., Scammon, R.J., 2002. Influence of mechanical properties on non-shock ignition. In: Proceedings of the 12th Symposium (International) on Detonation, San Diego, CA. Chen, P., Yuan, B., Chen, R., Qu, K., 2018. Compression and shear experimental study of PBX explosive. Propellants, Explos. Pyrotech. 43, 1–7. Clements, B., 2011. Merging the Mechanical Constitutive Nodel viscoSCRAM with the Henson-Smilowitz Thermal Ignition Model to Address Weak Shock Initiation. 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., 1983. On the stability of shear cracks and the calculation of compressive strength. J. Geophys. Res. 88, 1173–1179. Dienes, J.K., 1985. A statistical theory of fragmentation processes. Mech. Mater. 4, 325–335. 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. Field, J.E., 1992. Hot spot ignition mechanisms for explosives. Acc. Chem. Res. 25, 489–496. 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. 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. Technol. 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. Energetic Mater. 1–14. Idar, D.J., Straight, J.W., et al., 2000. Low amplitude impact of damaged PBX 9501. In: Shock Compression of Condensed Matter-1999. AIP Publishing, pp. 655–658. 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. Heterogeneous thermo-mechanical behavior and hot spot formation in RDX-estane energetic aggregates. Int. J. Solids Struct. 62, 91–103. Liu, R., Chen, P.W., 2018. Modeling ignition prediction of HMX-based polymer bonded explosives under low velocity impact. Mech. Mater. 124, 106–117.

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. References Akiki, M., Menon, S., 2015. A model for hot spot formation in shocked energetic materials. Combust. Flame 162, 1759–1771. Asay, B.W., 2010. Shock Wave Science and Technology Reference Library (Volume5): Non-shock Initiation of Explosives. Springer, Berlin. Barua, A., Zhou, M., 2013. Computational analysis of temperature rises in microstructures of HMX-Estane PBXs. Comput. Mech. 52, 151–159. Bassett, W.P., Johnson, B.P., Neelakantan, N.K., et al., 2017. Shock initiation of explosives: high temperature hot spots explained. Appl. Phys. Lett. 111, 061902. 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. Borne, L., Beaucamp, A., 2002. Effects of explosive crystal internal defects on projectile impact initiation. In: Proceedings of the 12th Symposium (International) on Detonation, San Diego, CA. Bowden, F.P., Yoffe, Y.D., 1952. Initiation and Growth of Explosion in Liquids and Solids. Cambridge University Press. Bourne, N.K., Field, J.E., 1991. Bubble collapse and the initiation of explosion. Proc. R. Soc. Lond. A 435, 423–435.

14

K. Yang et al.

European Journal of Mechanics / A Solids 80 (2020) 103924 Wang, X.J., Wu, Y.Q., Huang, F.L., 2017. Thermal-mechanical-chemical responses of polymer-bonded explosives using a mesoscopic reactive model under impact loading. J. Hazard Mater. 321, 256–267. Wei, Y., Kim, S., Horie, Y., et al., 2018. Quantification of probabilistic ignition thresholds of polymer-bonded explosives with microstructure defects. J. Appl. Phys. 124, 165110. Whitworth, N., 2008. Mathematical and Numerical Modelling of Shock Initiation in Heterogeneous Solid Explosives. Cranfield University. 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. Yang, K., Wu, Y.Q., Huang, F.L., 2019. Microcrack and microvoid dominated damage behaviors for polymer bonded explosives under different dynamic loading conditions. Mech. Mater. 137, 103130. Zhu, W., Xiao, J., Zhu, W., Xiao, H., 2009. Molecular dynamics simulations of RDX and RDX based plastic-bonded explosives. J. Hazard Mater. 164, 1082–1088. 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.

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. Ma, D.Z., Chen, P.W., Zhou, Q., Dai, K.D., 2013. Ignition criterion and safety prediction of explosives under low velocity impact. J. Appl. Phys. 114, 1–5. Nahshon, K., Hutchinson, J.W., 2008. Modification of the Gurson model for shear failure. Eur. J. Mech. A Solid. 27 (1), 1–17. 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, 3725–3744. Picart, D., Ermisse, J., Biessy, M., et al., 2013. Modeling and simulation of plastic-bonded explosive mechanical initiation. Int. J. Energetic Mater. Chem. Propuls. 12, 487–509. 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.S., Asay, B.W., et al., 2000. Microstructural effects in PBX 9501 damaged by shear impact. In: Shock Compression of Condensed Matter-1999, vol. 505. AIP Publishing, pp. 659–662, 1. Tarver, C.M., Chidester, S.K., 2005. On the violence of high explosive reactions. J. Press. Vessel Technol. 127, 39–45. Wang, X.J., Wu, Y.Q., Huang, F.L., 2017. Numerical mesoscopic investigations of dynamic damage and failure mechanisms of polymer bonded explosives. Int. J. Solids Struct. 129, 28–39.

15