Thermal mechanical anisotropic constitutive model and numerical simulations for shocked β-HMX single crystals

European Journal of Mechanics A/Solids 36 (2012) 66e82

Contents lists available at SciVerse ScienceDirect

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

Thermal mechanical anisotropic constitutive model and numerical simulations for shocked b-HMX single crystals Wu Yan-Qing*, Huang Feng-Lei Department of Engineering Mechanics, State Key Laboratory of Explosion Science and Technology, Beijing Institute of Technology, 100081 Beijing, PR China

a r t i c l e i n f o

a b s t r a c t

Article history: Received 23 February 2011 Accepted 23 February 2012 Available online 5 March 2012

The present study develops a crystal plasticity model for low-symmetric b-HMX octahydro-l,3,5,7tetranitro-l,3,5,7-tetrazocine single crystals with only limited operative slip systems, accounting for nonlinear elasticity, volumetric coupled with deviatoric behavior, as well as thermo-dynamical consistency. Based on the decomposition of the stress tensor, the modified equation of state for anisotropic materials is adopted. Simulation results of the planar impact on b-HMX single crystals show good agreement with existing experimental data by Dick et al. (2004a). In addition to providing new perspective to a range of orientation-dependent shock behaviors of b-HMX single crystal, the present work also discusses dislocation density, shear stress, strain localization, and anisotropic temperature increase in shocked b-HMX single crystals under shock loading. The proposed formulation and algorithms can also be applied to other low-symmetric crystals under impact or shock loading which gives irrecoverable deformation by crystallographic slip. Temperature calculations with various characteristic features for different orientations based on numerical simulations are explained, but no comparison with available experimental data is possible to our knowledge. Future studies should also examine phase change and twining as they also often occur in b-HMX single crystals. Ó 2012 Elsevier Masson SAS. All rights reserved.

Keywords: b-HMX single crystal Crystal plasticity Orientation dependence Shock loading

1. Introduction Several studies have observed crystal anisotropy in the mechanical behavior of several energetic single crystals (Dick, 1982, 1997; Hooks et al., 2006). Moreover, some researchers argue that anisotropy in shock sensitivity may be explained using such purely mechanical mechanisms (Elban et al., 1984; Dick et al., 1991; Dick, 1993). Dick (1984) found that the time and run distance to detonation in pentaerythritol tetranitrate (PETN) depend strongly on the direction of shock propagation relative to the crystal axes under plane shock wave loading. In subsequent studies, a steric hindrance model was developed to explain the observed anisotropic deformation pathways and the sensitivity to shock initiation of PETN (Dick and Ritchie, 1994, 1997). Gruzdkov and Gupta, 2000 proposed a chemical mechanism to explain the observed anisotropy in the shock wave initiation of PETN single crystals based on semiempirical quantum chemical calculations. Jindal and Dlott (1998) showed that the orientation-dependence of shock-induced temperature increase may give rise to large anisotropies in thermo-chemical reactions.

* Corresponding author. Tel.: þ86 10 68918019; fax: þ86 10 68461702. E-mail address: [email protected] (W. Yan-Qing). 0997-7538/$ e see front matter Ó 2012 Elsevier Masson SAS. All rights reserved. doi:10.1016/j.euromechsol.2012.02.011

Using the information on the above mechanism explanations, more conventional, straightforward, and full-anisotropic mechanical model is needed to be developed for examining the directiondependent shocked behaviors of energetic single crystals. Several studies have identified the presence of dislocations in PETN (Sherwood et al., 1988) and hexahydro-1,3,5-trinitro-1,3,5-triazine (HMX) (Sheen et al., 1993), as well as their operative slip systems. Information on the actual dynamic dislocation- controlled process is required to expand the understanding of initiation mechanisms. The dominant role of shear-driven plastic deformation in an initiation scenario is evident in explosive single crystals at three global scales: microscopic, mesoscopic and macroscopic (Plaksin and Coffey, 2007). Several researchers have investigated the anisotropic behavior of energetic single crystals because of the plastic localization and grain fragment of energetic single crystals at subgrain scales, which may provide structural conditions for hot spots under shock loading (Field et al., 1982; Palmer and Field, 1982; Bardenhagen et al., 2006). The present study thus develops a thermo-mechanical model in which the single crystal accommodates plastic deformation through the dislocation motion on preferred slip systems. For a plastic-bonded explosive (PBX) high explosives, if shock wave propagation and damage or initiation occur at scales less than their grain size, anisotropy in HMX single crystals must be

W. Yan-Qing, H. Feng-Lei / European Journal of Mechanics A/Solids 36 (2012) 66e82

considered (Leiber, 2000). b-HMX is the most stable structure at room temperature and pressure. Dick et al. (2004a,b) performed a series of plate impact experiments for ”shock” waves propagating in a single-crystal b-HMX, and found that anisotropy of the strength and elastic precursor shock decay reflect its anisotropic elasticity as well as its plasticity. Menikoff et al. (2005) adopted the isotropic elasticeplastic model to analyze the data for b-HMX accounting for crystal anisotropy by letting the constitutive parameters vary with uniaxial flow of propagation. Nevertheless, constitutive equations that assume isotropic material responses are certainly not applicable for anisotropic b-HMX single crystals. To address limitations of assumption of isotropic constitutive equations for b-HMX single crystals, inelastic anisotropies can be mapped by crystal plasticity simulations which can provide a guide to anisotropic effects. The experimental observations by Dick et al. (2004a) suggested that the conventional mechanisms of plastic flow (i.e., dislocation slip) are also operative in this weak, brittle, molecular HMX single crystal. Sheen et al., 1993 provided the slip plane and the direction of slip for operative slip systems, (001) <100> and ð101Þ < 101 >, based on the monoclinic crystallographic axis (a, b, and c) in space group P21/n. Assuming Schmid’s law for small plastic deformation of b-HMX single crystal, a model developed by Zamiri and De (2011) can predict the anisotropic plastic yielding of b-HMX under uniaxial compression along different crystallographic directions. According to Coffey’s theory (1998, 2001), energy dissipation and localization by the moving dislocations are the means of generating the local hot spots within shocked crystals that will occur directly through molecular excitation. Coffey and Sharma (1999) demonstrated the dominant role of shock-induced shear at reaction localization and energy dissipation within b-HMX single crystals. A dislocation pile-up mechanism was proposed to explain the generation of hot spots in dropweight impact tests (Armstrong et al., 1982). The dislocation pileups in obstructed slip bands could easily span a crystal crosssection (Armstrong, 1995). He also indicated that adiabatic heating associated with dislocation pile-up avalanches play the important role for the initiation of rapid chemical decompositions (Armstrong, 2009). Armstrong and Elban., 1999 pointed out that the low symmetry of the crystal lattice structures leads to expectation of a large variety of dislocation Burgers vectors and line vectors within the total dislocation density. In low-symmetric b-HMX energetic single crystal, little studies have considered the three-dimensional anisotropic crystal plastic effects. Simulations based on crystal plastic theory may provide the conceptual framework for examining the extent to which dissipative energy can be deposited in shocked b-HMX single crystals. Consequently, the following question arises: To what extent does the anisotropic nature of b-HMX alter the technical expectations based on plane wave engineering approximations? Some anisotropic modeling developments utilized a linear elastic response for describing shock wave propagation in single crystal (Johnson et al., 1970, 1971; 1972). Barton et al. (2009) examined the crystal mechanics based model for interactions among wave motion, slip kinetics, defect generation kinetics at physical length scale. Pore collapse inside the HMX crystal and implications for experimental observations are predicted. Still, nonlinear elastic behavior has not been considered in their model. Such anisotropic models are adequate for weak shocks, however, for large compressions, nonlinear elasticity is needed. The question then arises regarding what happens when nonlinearly elastic anisotropic materials are shock-loaded in directions other than those specific ones which result in pure longitudinal motion. In anisotropic materials, transverse particle velocity can be produced by the coupling that exists between all components of stress and strain (Johnson, 1971). For some materials, the degree of coupling

67

may be small, either by virtue of a significant amount of symmetry possessed by the crystal or through the magnitude of the elastic constants being nearly representative of an isotropic material. For example, Winey and Gupta (2004, 2006a, 2006b) incorporated nonlinear elasticity and crystal plasticity in a thermodynamically consistent tensor formulation for problems involving shock waves in copper and LiF single crystals. Developing new models to extend crystal plastic anisotropy under quasi-static loading to simulations under shock loading of more widely-used materials is necessary (Gray et al., 2002; Vignjevic et al., 2002). Another challenge in developing a crystal plastic model under high pressure loading for single crystals is to couple their volumetric and deviatoric behaviors, and to obtain the required material data. Becker (2004) applied pressure-dependent crystal elastic moduli to capture large volume strains properly and to enable evolution of shocks from steep pressure gradients. Nevertheless, lattice rotation which is significant for low-symmetric shocked single crystals, and should thus not be neglected in such models. Hence, maintain material frame indifference of the constitutive equations, polar decomposition of elastic deformation gradient is employed in the stress and strain transformations for finite rotation in the model of the present study. To our knowledge, no model which incorporating the combinations of low-symmetric crystal plasticity and heat effect, have been developed. Existing experimental data available for b-HMX single crystals consist of the research conducted by Sewell et al. (2002) and Menikoff and Sewell (2001, 2003), as well as some preliminary data are also reported by Zaug (1998). The present study investigates the capability of current constitutive models to accurately predict the behavior of shocked b-HMX single crystals. Moreover, the present study examines the effect of orientation, and crystal anisotropy on the wave characteristics of bHMX single crystals, as well as explains their physical significance. Rate-dependent plastic flow in terms of dislocation approach in single crystals is described. Any phase changes in b-HMX that may be induced by the shock are neglected. Duvall and Fowles (1963) argued that, for convenience, the meaning of “shock” can be extended to include cases where stress and density change abruptly from one value to another, even though one or both of the connected states may be changing with time. Hence, to the present study simulates the wave propagation in b-HMX, and compares the simulation results with experimental data in the literature. After identifying the mechanical responses of b-HMX single crystals, the orientation-dependent temperature increases and the distributions are plotted, which allow us to evaluate the extent of heating induced by directional compression waves. Notation is based on the following conventions: Tensor scalar products of appropriate order are denoted by a raised dot. The superscripts “T”, “1”, “e” or “p” denote “transpose”, “inverse”, “elastic part”, and “plastic part”, respectively. The basis of the global orthonormal coordinate system is denoted by {e1 e2 e3}, and the crystallographic axes in group space P21/n are f a b c g for monoclinic b-HMX single crystals. 2. Model description The development of crystal plastic models can be traced to the studies of Hill and Rice (1966), Asaro and Rice (1977), Peirce et al. (1983), Asaro (1983a,b), and Asaro and Needleman (1985). Most of available crystal plasticity models consider materials with high symmetry, such as face-centered cubic, and body-centered cubic crystals. In the present study, assuming that the deformation induced by shearing along crystallographic slip planes is the only plastic mechanism for monoclinic b-HMX, the general framework

68

W. Yan-Qing, H. Feng-Lei / European Journal of Mechanics A/Solids 36 (2012) 66e82

described by Hill and Rice (1972) is adopted and integrated with nonlinear thermo-elasticity theory. 2.1. Kinematics and crystal plasticity Starting with the multiplicative decomposition and polar decomposition theorems to the deformation gradient, a deformed ~ is defined by F~ p ¼ U e $F p, and Ue is but unrotated configuration B the elastic right stretch tensor. All quantities with a bar “w” above ~ The rate of symmetric refer to the reference configuration B. ~ and the antisymmetric spin tensor W ~ is related stretching tensor D to the lattice rigid rotation tensor Re by,

~ ¼ ReT DRe D

(1)

~ eT W  U [ Re WR

(2)

_e

where U [ R $Re1 represents the rate of rigid-body rotation. To express the crystal constitutive equations through slip systems, the resolved shear stress based on Schmid law is given by, ðaÞ sðaÞ ¼ P~ : J sL

large shear stresses can exist momentarily; thus, the creation of dislocations at stress concentrations must be considered (John and Lothe, 1982). That is, at higher stresses, the multiplication of dislocations and nucleation will operate simultaneously. If multiple cross glide is dominant then the rate of mobile dislocation density rm is assumed to be proportional to the accumulated shear strain rate



ðaÞ r_ m ¼ M g_ pðaÞ exp  H g=sðaÞ

Pn



(9)

Rt

_ ðaÞ jdt denotes the total cumulative shear where g ¼ a ¼ 1 0 jg strain on each slip system. To incorporate the dislocation nucleation concept into the model, the rate of dislocation density is described by adding another stress-dependent nucleation term (Gupta et al., 1975; Winey and Gupta, 2006a),





ðaÞ r_ m ¼ M g_ pðaÞ exp  H g=sðaÞ þ Aðs  sc Þs_

(10)

where sc is a threshold stress for nucleation and A is a material constant.

(3)

where s is the Cauchy stress defined in the deformed but unro~ sL is related to Cauchy stress s in current tated configuration B. configuration by, L

s ¼ R sR L

eT

e

(4)

J is the third invariant of the elastic part of the deformation gradient. The slip direction of a particular system a in initial configðaÞ ðaÞ uration, is denoted by s0 and the slip plane with normal m0 . The ðaÞ ðaÞ ~ convect with the lattice to become, orthonormal pair ~s and m

~ ðaÞ ¼ ~sð0aÞ ,U e1 ~sðaÞ ¼ U e ,~s0ðaÞ ; m

(5)

~ ðaÞ and tensor It is convenient to introduce the Schmid factor P ðaÞ ~ as, W

i h ~ ðaÞ ¼ 1 ~sðaÞ m ~ ðaÞ~sðaÞ ~ ðaÞ þ m P 2 h i ~ ðaÞ ¼ 1 ~sðaÞ m ~ ðaÞ  m ~ ðaÞ~sðaÞ W 2

(6)

In the description of g_ ðaÞ induced by dislocation mobility, the Orowan equation (Gilman, 1969) is adopted, ðaÞ ðaÞ ðaÞ g_ ðaÞ ¼ rm b v

(7) ðaÞ

where b(a) is the Burgers vector, rm is the mobile dislocation density and vðaÞ is the average dislocation velocity for slip system a. Dislocation velocity is assumed to be related to resolved shear stress s(a) on a slip system by the relation (Gupta et al., 1975),

# " s ðaÞ vðaÞ ¼ v0 exp  ðaÞ d  s  s0 ðaÞ

(8)

where v0 is the sound speed for shear waves on the a -slip plane, and sd is the drag stress, and s0 is the threshold shear stress for dislocations motions. The dislocations are limited to move at velocities less than the shear wave velocity at the corresponding shear slip plane. Two distinct processes contribute to the strain rate: nucleation of dislocation loops and growth of existing loops (Hull and Bacon, 1984). At low strain rates, increases in dislocation density are usually attributed to multiplication of existing dislocations via a process called multiple cross glide (Gilman and Johnston, 1962). On the other hand, for solid under extremely high strain rates, very

2.2. Nonlinear thermo-elasticity For a monoclinic single crystal, when the elastic property is determined, an isochoric deformation may produce pressure, which is not consistent with what they were traditionally understood to be. To include finite volume changes under shock loading, a logarithmic strain measure is used to formulate the elasticity (Hoger, 1987; Becker, 2004),

~3e ¼ ln U e

(11)

The material time derivatives of elastic strain is written as, e ~e ~3_ ¼ D

(12)

In thermo-elasticity, the elastic strain and entropy specify the thermodynamic state. Internal energy per unit mass is expanded in terms of the elastic strain tensor and entropy (Wallace, 1980). At constant entropy, higher-order elastic constants may be used to approximate internal energy increments,

    S 1 S 1 S r0 e S;~3 eij ¼ r0 e S;0 þ C~ ij~3 eij þ C~ ijkl~3 eij~3 ekl þ C~ ijklmn~3 eij~3 ekl~3 emn 2!

3!

(13)

S

~ where C ijkl is the second-order isentropic elastic modulus with ~ The nonlinear elastic effects in Eq. (12) reference configuration B. may be used as a basis for the conditions that arise in shocked materials. Assuming small distortional elastic strains, Cauchy stress in an unrotated configuration times J ¼ det (Fe) is work conjugate to ~3e (Becker, 2004). When the entropy change is zero, the formula can be written as:,

J sLij ¼ r0



ve v~3 eij

 and J ¼

r0 r

(14)

S

An objective quantity is one which transforms in the same manner as the energy conjugate stress and strain rate pair under a superposed rigid-body motion. The fundamental advantage of the unrotated Cauchy stress sL is that the material derivative of sL is objective, whereas the material derivative of s is not. A similar stress rate, referred to as the Green-Naghdi rate (Johnson and Bammann, 1984) can be derived by transforming the rate of the unrotated Cauchy stress to the fixed global reference frame as follows:

W. Yan-Qing, H. Feng-Lei / European Journal of Mechanics A/Solids 36 (2012) 66e82



VL

s ¼ Re s_ L ReT ¼ s_  U$s þ s$U Green-Naghdi rate 

V

s ¼ s_  W$s þ s$W Jaumann rate





Although the Green-Naghdi stress rate

 VL 

s

(15)

avoids the finite V

strain problems that plague the Jaumann rate ðsÞ, the use of unrotated Cauchy stress rate directly is better when casting all constitutive models in the unrotated frame. In the case of positive entropy change, rate-form of the constitutive equation can be given based on Eqs. (13) and (14),

s_ Lij þ sLij tr



~ D

e



0

1  , ve e ~e ¼ r0 @v v~3 mn A D mn v~3 eij S 1 0  , ve @ vSA S_  sLik q_ kj þ q_ ik sLkj þ r0 v v~3 eij 3e   s 1 s e e ~e ~ D ~ ~ ¼ C 3 kl D mn ijkl kl þ C ijklmn~ 2   !!  ~  ~ vC vC ij  ijkl  e þ S_  sLik q_ kj þ q_ ik sLkj ð16Þ  þ  3 vS  e vS  e kl 3

where

3



q_ ¼ ReT  W  UÞ  Re

(17)

~ The three-order elastic tensor C ijklmn is assumed to be independent of strain and entropy in the present formulation. The Grüneisen tensor is defined as,

 ~  1 vC ij  ~ Gij ¼   rT vS  e

(18)

~3

C1111 C2222 C3333 C2323 C3131 C1212 C1122 C1133 C2233 C1131 C2231 C3331

3

(19)

S

~ is independent of strain and entropy for where we assume that G ij ~ s represents the isentropic components of the simplicity, and C ijkl elastic stiffness tensor of the material with orientation in the deformed but unrotated configuration. Hence, the complete expression of rate-form constitutive equation is given by,

s_ Lij ¼

  1~s e ~e ~ _ ~ ~ 3 e S_ ~s D ~e C 3 kl D mn  rT Gij S þ rT Gij Gkl~ ijkl kl þ C ijklmn~ kl 2   ~ e  sL q_ þ q_ sL  sLij tr D ik kj ik kj

(20)

A small amount of data is obtained on the higher-order elastic constants of single crystal b-HMX, which are herein defined as s ~s ~ 3 mn to a first approximation. Typically, the C ijklmn ¼ vC ijkl =v~ nonlinear elastic response of b-HMX single crystals is assumed to e ~ ~ 3e ¼ ðvC=vpÞðvp=v~ 3 V Þ. be dependent only on the pressure vC=v~ The pressure-dependent elastic moduli of b-HMX are a subject for future research. To obtain the pressure-dependent relation of ~ s =vp, molecular dynamic simulations are elastic constants vC ijkl performed using COMPASS field of force for b-HMX single crystals. According to the neutron diffraction data (Choi and Boutin, 1970), an b-HMX superlattice model with 4  2  3 molecular configuration is constructed and simulated using MS-DISCOVER NPT at a constant temperature of 295K, with pressure ranging from 1  104 GPa to 27 GPa.

¼ 9:02397 þ 11:0554 p  0:12183 p2 ðGPaÞ; ¼ 15:56815 þ 7:59921 p  0:10686 p2 ðGPaÞ; ¼ 18:39271 þ 5:90386 p  0:06176 p2 ðGPaÞ; ¼ 6:24105 þ 1:352751 p  0:01491 p2 ðGPaÞ; ¼ 1:59267 þ 5:16237 p  0:10296 p2 ðGPaÞ; ¼ 1:13088 þ 3:78434 p  0:09663 p2 ðGPaÞ; ¼ 1:32088 þ 6:52789 p  0:10792 p2 ðGPaÞ; ¼ 1:3792 þ 7:77819 p  0:1551 p2 ðGPaÞ; ¼ 0:22377 þ 5:7947 p  0:10229 p2 ðGPaÞ; ¼ 0:2347 þ 1:53836 p  0:02168 p2 ðGPaÞ; ¼ 3:09946  0:53953 p þ 5:02487 e  4 p2 ðGPaÞ; ¼ 0:7547ðGPaÞ;

C2312 ¼ 2:41649  0:21976 p þ 5:72921 e  4 p2 ðGPaÞ

> > > > > > > > > > > > > > > > > > > > > > > > ; (21)

~ s =vp can be obtained based on the above mentioned Thus, vC ijkl equations, where p is the hydrostatic pressure increment referred to in ambient conditions. The pressure-dependent elastic moduli for b-HMX is still a future experimental task. To complete the thermo-dynamic description, the temperature increasing rate is given by,

T_ ~ D ~e T_ ¼ T G ij ij þ S c3

(22)

The simple relationship between hydrostatic pressure and volumetric strain does not exist for anisotropic materials. To distinguish between thermodynamic response (equation of state or EOS) and the ability of the material to carry shear loads (strength), the stress increment derived from Eq. (20) is split as follows:

 1 j : C~ : I D~3eV  rT DSG0 3 1 _e0 þ rT DSj : G5G: 3 þ rT DSðj : G5G: IÞ~3eV 3   0  sL tr D~3e  sLik Dqkj þ Dqik sLkj

DsL ¼ j : C~ : D~3e þ 0

To maintain thermodynamic consistency, the derivatives ~ =vSÞ must satisfy, ðvC ijkl 3

 " # ~  vC vT ijkl  ~ ~ ~ G ¼ rT G  ¼ rGij ij kl vS  e v~3 kl

69

9 > > > > > > > > > > > > > > > > > > > > > > > > =

DsLm

0

1  e0 ~  D~3 : C : I þ ¼ 3

~ ~ 0 1 e0 dC 1 0 dC ~3 : : D~3e þ D~3e : :I 2 dp 3 dp   ~ m þ 1rT DS ~3e0 : G ~ 5G ~ : I  sL D~3 e  rT DSG V m 3 1 ~  e I : C : I D~3 V þ 9 |fflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl}

ð23Þ !

dp d3 eV

ð24Þ

replaced by DPEOS

where j [ J L (1/3) I 5 I, J is the fourth order identity and I is the second-order identity tensor. The last term in Eq. (24) deals with a portion of the mean stress that varies directly with the volumetric strain which should be replaced by an EOS. The high pressure thermodynamic EOS must be reduced to the correct relations at small volumetric strains. To account for this term with respect to volumetric strain, the Mie-Grüneisen EOS with principal Hugoniot as the reference curve is written as,

    G pjEOS ¼ PH m 1  m þ rGe 2

(25)

where G is the Grüneisen coefficient, PH (m) is the Hugoniot pressure which represents the locus of points that may be reached by a shock transition from an initial unshocked state; and m ¼ r/r0  1 is the relative change of volume. The Hugoniots are determined by a linear Us w up relationship. The following form for PH (m) is employed in the present study (Anderson et al., 1994),

  PH m ¼



A01 m þ A2 m2 þ A3 m3 ; A01 m;

m>0 m<0

(26)

70

W. Yan-Qing, H. Feng-Lei / European Journal of Mechanics A/Solids 36 (2012) 66e82

where A01 ; A2 ; A3 can be analytically determined through a Taylor’s series expansion of the Hugoniot curve (Alexander, 2008):

A01

¼

r0 c20 ;

A2 ¼

A01 ½1 þ 2ðs  1Þ;

A3 ¼

A01

where r0 is the reference density, c0 is the sound speed at zero pressure, and s is the constant in the approximated linear relationship between shock velocity Us and particle velocity up,

 1 ¼ I:C:I 9

(29)

Thus, c0 has the following definition,

c0 ¼

sffiffiffiffiffiffi A01

(30)

r0

The methodology described above can be applied for all anisotropic materials and represents a mathematically consistent generalization of the conventional isotropic case.

(35) _ ð aÞ

bðaÞ ¼ W

Constitutive equations presented in Section 2.2 is a ratedependent model as the slip rates are directly linked to the instantaneous resolved shear stresses. A variety of integration procedure for updating the stresses and different sets of governing equations have been used by different researchers (Peirce et al., 1983; Ling et al., 2005). Taking Ds(a) as the primary unknowns, the controlling nonlinear equation in the deformed but unrotated system is deduced, which is used to solve the slip rate iteratively. Following the exponential evolution of the plastic part of the ~ p proposed by Ortiz et al. (2001), deformation gradient F p ~ F~ tþDt ¼ exp Dt L

 p

p F~ t

p

" p F~ tþDt ¼

Iþ

n X

a¼1

# ~ ðaÞ F pt Dg~sðaÞ 5m

(32)

Based on Eq. (32), the elastic part of the deformation gradient at the end of the time step can be given calculated. The resolved shear stress increments Ds(a) is calculated by P _ðbÞ ðbÞ considering the relationship D3e ¼ D3  P Dg , b

Ds

ðaÞ

¼ r

ðaÞ



ðabÞ r1 DgðbÞ

ðaÞ

~ : D3 þ P

ðaÞ





h  .  ðaÞ ðaÞ ðaÞ ¼ bðaÞ v0 Dt rm t þ M Dg_ p exp  HgtþDt stþDt 2 3 i sd ðaÞ 4 5 exp  ðaÞ þAðs  sc ÞDs stþDt  s0

(37)

Substituting Eq. (38) into Eq. (34) leads to the controlling equation which can iteratively calculate the resolved shear stress ðaÞ Ds(a), and then the resolved shear stress at the end of each step stþ Dt is obtained. Temperature and entropy increments (DT, DS) are calculated by the following equations,

DT ¼ Tt G : D3 þ

Tt DS CV

(38)

The irreversible work due to shear slip strains is dissipated primarily as heat and results into entropy production,

DS ¼

1 X ðaÞ s DgðaÞ T a tþDt

(39)

After obtaining the shear strain Dg(a) from Eq. (38), the plastic strain increment with respect to local unrotated crystallographic axes D~3p is calculated as,

D~3ptþDt ¼

n X

a¼1

~ ðaÞ DgðaÞ P t tþDt

(40)

The spin tensor increment is given by, n X

p

~ DU tþDt ¼

a¼1

~ ðaÞ DgðaÞ W t tþDt

(41)

A method used to integrate Ue directly use,

 e  ~ Dt $U e U etþDt ¼ exp L t

(42)

where U et is the elastic stretch tensor at the beginning of the time ~e is the velocity gradient referred to the unroated step. L configuration. ðaÞ Average velocity vtþDt and dislocation density increment ramtþDt ðaÞ are updated by substituting stþDt into Eq. (8) and Eq. (9)(10) 0^ respectively. The increments of deviatoric stress DsL and hydrostatic stress Dsm are updated using Eq. (23) and (24) respectively,



ðaÞ sLtþDt  sLt ¼ DsL D3etþDt ; D3 eVtþDt ; Dgtþ Dt ; TtþDt 0

0

0



0

aÞ Dsm ¼ Dsm D3etþDt ; D3 eVtþDt ; Dgðtþ Dt ; TtþDt 0





¼ Dp

(43) (44)

(33)

where

r ðaÞ ¼ b

(36)

ðaÞ ðaÞ ðaÞ ðaÞ ðaÞ DgðaÞ ¼ g_ tþ Dt Dt ¼ rmjt þ Drm b vtþDt Dt

(31)

where F~ t is the plastic deformation gradient at the beginning of the time step. Kalidindi et al. (1992) approximated the exponential evolution equation as,

_ ðaÞ

$sL  sL $W

Express the shear strain increment based on Eqs. (7)e(10) as,

2.3. Numerical implementation



 X ðaÞ ðbÞ 1 ðaÞ vC ðbÞ ðaÞ ðbÞ b : P~ þP~ : C: P~ þ P~ : 3e : e : P~ DgðbÞ 2 v3 b

(28)

Given that the sound speed varies with different orientations of single crystals, an equivalent (bulk) sound used for c0 is introduced. For monoclinic b-HMX single crystal, the “effective” or “average” bulk modulus A01 is defined based on the uniform dilatation condition,

A01

¼

h  i  2 s  1 þ 3ðs  1Þ2 (27)

Us ¼ c0 þ sup

ðabÞ

r1

 ðaÞ  1 ~ ðaÞ e vC ~ : G rT DS : C: D3 þ P : 3 : e : D3  P 2 v3 (34)

The elastic constants in terms of the pressure increment are given by, s

~0 ~0 ~ C ijkl ¼ C ijkl þ C ijkl

~ vC ijkl vp

Dp i; j; k; l ¼ 1; 2; 3

(45)

W. Yan-Qing, H. Feng-Lei / European Journal of Mechanics A/Solids 36 (2012) 66e82

Cauchy stress increment with referred to the current configuration can be approximated by, e eT DstþDt ¼ RetþDt sLtþDt ReT tþDt  DR st DR

(46)

where st is the stress at t and defined with respect to configuration at time t, st þ Dt is defined with respect to configuration at time t þ Dt. The ABAQUS/explicit uses a forward Euler integration scheme and integrates through time by using small time increments. An approximation to the stability limit is often written as the smallest transit time of a dilatational wave across any of the elements in the mesh,

1 Lmin 3 cd

Dtzpffiffiffi

(47)

where Lmin isffi the smallest element dimension in the mesh and cd ¼ p ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðl þ 2mÞ=r is the dilatational wave speed. l, m are Lame elastic constants and r is material density. 3. Simulation of shocked b-HMX single crystal 3.1. Finite element model and parameters determination The constitutive equations presented in Section 2 for b-HMX single crystal are implemented in the user material subroutine (VUMAT) of ABAQUS/Explicit software code. In this approach, the algorithm for integration rate constitutive equations which is the objective for large rotation increments, is adopted (Hughes and Winget, 1980; Flanagan and Taylor, 1987), and the adaptive meshing technique, which is referred to as the Arbitrary LagrangianeEulerian (ALE) analysis is employed. The model geometry is a cube, 8 mm high and 1.2 w 4 mm thick which are derived from the experimental conditions as shown in Fig. 1. Periodic boundary conditions are applied to the four faces for the computational cell. Except for a few specific single crystal orientations, a shear accompanies the compression. The essential feature of such boundary conditions is that the center regions try to maintain uniaxial strain condition where measurements are made during the experimental recording time. Shock plane can be chosen to be consistent with arbitrary orientations. A time-varying

71

pressure-controlled loading is applied on the front surface normal to the z-axis, as shown in Fig. 1(b). The pressure is kept on its maximum value for a period of wave holding time. A frictionless contact condition is imposed on the surface between the tested sample and a support material with the same property as b-HMX single crystals. Quasishear reflection waves may exist at the interface between the target domain and the backed block. Material constants generally affect the capability of constitutive models to a large extent. The determination of constitutive parameters of b-HMX is challenging due to lack of sufficient experimental data. There are totally 14 parameters in the present ðaÞ model, including rm0, M, H, A, v0 , sd, s0, sc, A, b, cV, r0, c0, s, G .. Some can easily be obtained through searching references, e.g. cV, r0, c0, s for b-HMX. Among them, cV and r0 are the specific heat and density for b-HMX single crystal at ambient temperature and pressure. The parameters c0, s, G are sound speed, slope of UseUp Hugoniot, and Gruneisen parameter, which can be found in the literature (Menikoff et al., 2004; Menikoff and Kober, 2000). ðaÞ Some parameters such as rm0, v0 ,b, are chosen according to their physical interpretation. rm0 is initial dislocation density in single crystal, and b is the Burgers vector which can be defined as the shortest lattice translation vectors. Following the work of Armstrong (1982, 2009), their value are taken to be 1.4  10 m2 and 0.713 nm, respectively. Gallagher et al. (1992) have provided values of Burgers vector for Pentaerythritol Tetranitrate (PETN) and Cyclotrimethylene Trinitramine (RDX) single crystals, which are in the range 1.08e1.48 nm that are comparable to that of b-HMX. Large Burgers magnitude of elastic stiffness gives relatively large dislocation self-energies, consonant with low dislocation densities ðaÞ being observed in such crystals (Armstrong and Elban, 2004). v0 refers to shear wave speed limit on two slip planes {001} and {101}. ðaÞ For convenience, v0 has been taken to be the same value of 3280.0 m/s. Others consist of M, H, A, sd, s0, sc, are tuned using an optimization method with initial values chosen to be physically reasonable. Among them, sd, s0, sc are the characteristic drag stress, threshold shear stress for dislocation motion and dislocation nucleation, which should satisfy sd > s0 > sc. A series of plane impact simulations are performed to calibrate the parameters. M, H are parameters influencing dislocation growth mechanisms, and A, sc are incorporated into the dislocation nucleation mechanisms. sd,

Fig. 1. a) Finite element mesh for the single crystal cube with uniform pressure suddenly applied to the shock plane; b) time-varying pressure ramped from zero to a certain level, keeping constant and then decreases to zero.

72

W. Yan-Qing, H. Feng-Lei / European Journal of Mechanics A/Solids 36 (2012) 66e82

s0 are related to dislocation velocity. Lower M, A, or larger sd, H will bring an increase in the elastic precursor peak value. s0, sc is taken to be in the order of 106 according to the works of Zamiri and De (2010). It’s found that a slight variation in dislocation density will lead to an intense change in shear stress. All the calibrated simulations have shown that dislocation density evolution and shear stress changes are strongly interdependent. We calibrate the adjustable parameters by minimizing differences with the experimental measured velocity histories at the back surface of the crystals. ~ s =vp is calcuIn addition, dependence of elastic constants vC ijkl lated by molecular dynamics which have been provided in Section 2.2. The method to determine Grüneisen tensor G in Eq. (22) has been provided in Appendix. As in the P21/n space group, the second molecule in the b-HMX lattice is cell centered. For lower symmetry crystals, each of the slip system of (001) <100> and ð101Þ < 101 > are unique cases that are not like simple cubic. Elastic constants and slip systems data should be transformed from the crystallographic system to the Cartesian system, which are given in the Appendix. All parameters listed in Table 1 and Table 2 are related to the plastic responses and EOS properties of b-HMX single crystal. 3.2. Mesh convergence To check the mesh convergence, several simulations are performed for the computational cell with the same sample size 8  8  3 mm (i.e. the thickness of 3 mm) with shocked direction along (110) orientation. Several cases are investigated: the samples consisting of 8, 12, 17, 22 and 26-elements along the shocked direction, which have the mesh sizes of 375.0, 250.0, 176.0, 136.0, 115.0 mm, respectively. A shocked pressure was ramped from zero to 1.5 GPa in 0.01 ms, which is held constant for 0.1 ms and then decrease to zero as shown in Fig. 2(a). A short period pressure wave is chosen to examine the effects of an evolving shock spanning the thickness. The simulated elastic wave speed and plastic wave speed observed at the rear face for the five mesh-sized sample are plotted against the mesh numbers along the sample thickness, as shown in Fig. 2(b). It’s noted both the elastic wave speed and the plastic wave speed are influenced by mesh size. The elastic wave speed converges to 4.83 mm s1, and the plastic wave speed to 3.85 mm s1 respectively with mesh size diminishes. From Fig. 2(b), when the mesh size is finer than 136.0 mm (corresponding to 22elements along thickness), the results differences are minor. From a physical standpoint, mesh convergence of a model should be reflected through distribution of some critical parameters results. Fig. 2(c) and (d) show the distribution of the z-direction stress component szz at 0.3 ms and 0.5 ms along one specific path, which is from the central position at the loading face to the center of the rear face. The stress wave moves from left to right at a bit faster speed and the stress gradient has sharpened in the case of finer mesh size. The curve patterns become finer when the mesh resolution is improved. Fig. 2(c) and (d) show that the curve discrepancies on stress component distributions for 8-element and 12-element are significant. There exist 15%e20% variations between the simulated stress at the same position using 17-element model and that by 22element and 26-element models. If one mesh size makes the relative errors between the magnitude of stress values are less than a certain maximum error (e.g. 1%),

Table 2 Parameters related to b-HMX physical properties. CV (J kg1 K1)

r0 (kg m3)

EOS parameters c0 (m/s)

s

G

1.559e3

1900

2565.0

2.38

1.1

we can say that such a mesh size is an available mesh size. Results in Fig. 2(bed) have implications for mesh generation that a relatively low number of elements in which there are 22-elements along 3 mm thickness sample can provide convergent description of the stress distribution and wave motion. 3.3. Particle velocity histories Stress wave propagation in an isotropic linear elastic continuum has been well understood (Kolsky, 1953). In contrast, the stress wave propagation in their anisotropic nonlinear elastic media is more complex. b-HMX single crystals exhibit monoclinic anisotropy in its elasticity and preferential slip shearing in plasticity. These anisotropies lead to a more complex shock-induced deformation as well as wave propagation behavior. The plate-impact experimental configuration produces conditions of roughly uniaxial strain at the center of b-HMX single crystal samples. As a check on the constitutive model, a series of plane shock simulations is performed on the b-HMX crystals with input stresses of 1.5 GPa and 2.4 GPa respectively. Dick et al. (2004a,b) performed a series of plate impact experiments and measured Lagrangian particle velocity time histories, using a velocity interferometry system for reflector diagnostic, of “shock” waves propagating in a single-crystal HMX. In Fig. 3, the particle velocity time histories of the center point on the back surface of the b-HMX samples are compared between simulated results and experimental data. The time origin corresponds to the start of the impact on the b-HMX samples. Fig. 3(aed) correspond to the cases wherein the shocked strength is 1.5 GPa and samples with (011), (110) and (010) orientations of different thicknesses. When the shock strength increases to 2.4 GPa shown in Fig. 3(c), the plastic shock amplitude for each sample with different thicknesses is increased by the higher input stress. Upon impact, an elastic stress jump occurs followed by a stress relaxation due to rate-dependent plasticity. Obviously, the anisotropic wave structures exhibited by b-HMX can be grasped using the present constitutive model during the distance of run. Simulated elastic precursor strengths in Fig. 3 also show orientation-dependent elastic precursor limits which are closely related to the plastic slip mechanisms. In Fig. 3(aec), for (110) and (001) orientations, the simulated results for the elastic wave precursor and the plastic wave regime are in good agreement with the experimental data. A maximum difference in the arrival time of the elastic wave precursor is no more than 0.04 ms and at most 14.2 m/s difference between simulated particle velocity and experimental one. Nevertheless, it should be mentioned that the peak values in the elastic precursor wave measured by experiments are much higher than the simulated ones in the case of (010) orientation. The maximum differences are 53.3 and 67.9 m/s between experimental data and the simulated results for (010) orientated sample with 3.95 mm and 4.4 mm, respectively. This inconsistency may imply that a different deformation mechanism

Table 1 Parameters for dislocation motion model which are related to b-HMX plastic properties. ðaÞ

rm0 (m2)

M (m2)

H (Pa)

v0

1.4e10

1.23e12

1.0e10

3280.0

(m/s)

sd (Pa)

s0 (Pa)

sc (Pa)

A (m2 Pa2)

b (m)

0.17e9

1.2e6

1.56e6

1.0e-9

7.1332e-10

W. Yan-Qing, H. Feng-Lei / European Journal of Mechanics A/Solids 36 (2012) 66e82

73

b 5.5

a

-1

Wave speed (mm.μs )

5.0

Pressure (GPa)

1.5

4.5 loading direction along (110) orientation plane Elastic wave speed Plastic wave speed

4.0 3.5 3.0

0.0

0.01

0.11 0.101 μ s 2.5

8

c

0.3

d

at 0.3μs

-0.3

at 0.5 μs

-0.3

The number of elements along the sample thickness

-0.6

8-elements 12-elements 17-elements 22-elements 26-elements

-0.9 -1.2 -1.5 0.0

28

0.0

σzz (GPa)

σzz (GPa)

0.0

0.3

12 16 20 24 Mesh numbers along the thickness

The number of elements along the -0.6 sample thickness 8-elements 12-elements 17-elements 22-elements 26-elements

-0.9 -1.2 -1.5

0.5

1.0 1.5 2.0 2.5 Distance along shocked direction (mm)

0.0

3.0

0.5

1.0 1.5 2.0 2.5 Distance along loading direction (mm)

3.0

Fig. 2. Mesh convergence analysis for (110) b-HMX single crystal model a) a short time maintaining pressure wave is applied to examine the effects of an evolving shock, b) through plotting elastic and plastic wave speeds against mesh number along the shocked direction, c) and d) z-direction stress component distributions along the wave propagation path at 0.3 ms and 0.5 ms, respectively.

Shock strength: 1.5GPa

b

300

Shock plane (110) with thickness of 3.18mm

200

Experimental data Simulated results Shock plane (110) with thickness of 1.23mm Experimental data Simulated results

100

0 0.0

0.4

0.8

1.2

Particle Velocity (m/s)

Particle Velocity (m/s)

a

300

Shock strength: 1.5GPa

Shock plane (011) with thickness of 3.00mm

200

Experimental data Simulated results

100

Shock plane (011) with thickness of 1.39mm

Experimental data Simulated results

0 0.0

1.6

0.4

Time (μs)

500

Shock strength : 2.4GPa

d 300

1.2

1.6

Shock strength: 1.5GPa

400 Shock plane (011) with thickness of 3.11mm Experimental data Simulated results

300

200 Shock plane (110) with thickness of 3.57mm Experimental data Simulated results

100

0 0.6

0.8

1.0 1.2 Time (μs)

1.4

1.6

Particle Velocity (m/s)

Particle Velocity (m/s)

c

0.8 Time (μs)

200

100

0 0.0

Shock plane (010) with thickness of 3.95mm Experimental data Simulated results Shock plane (010) with thickness of 4.40mm Experimental data Simulated results

0.4

0.8

1.2 Time (μs)

1.6

2.0

Fig. 3. Comparisons between the simulated and experimental particle velocity histories for shocked planes parallel to, a) {110} orientation with 1.23 mm and 3.18 mm thickness, subjected to shock strength of 1.5 GPa ramped from zero in 0.01 ms; b) {011} orientation with 1.39 mm and 3.00 mm thickness, subjected to shock strength of 1.5 GPa ramped from zero in 0.01 ms; c) {110} orientation with 3.57 mm thickness and {011} orientation with 3.11 mm thickness subjected to shock strength of 2.4 GPa ramped from zero in 0.01 ms; d) {010} orientation with 3.95 mm and 4.40 mm thickness, subjected to shock strength of 1.5 GPa ramped from zero in 0.01 ms.

74

W. Yan-Qing, H. Feng-Lei / European Journal of Mechanics A/Solids 36 (2012) 66e82

from that for (110) and (011) orientations may be operative. Barton et al. (2009) also pointed out that there may not be any easily activated slip systems for shots along (010) so that there is need to consider twinning and cleavage mechanisms in future work. In Fig. 3, decreases in the elastic wave velocity with the increase in sample thickness are evident, which can also be obtained by the present simulations. Some differences in particle velocity also exist when the plastic wave arrives with the same strength of 1.5 GPa and different orientations. In Fig. 3(a), the (110) orientation shows a slightly higher particle velocity (about 262e280 m/s) than that of the (011) orientation (about 256e262 m/s). when subjected to the shock strength of 2.4 GPa, the particle velocity increases from 392 to 417 m/s. The simulated mean elastic wave velocity for each case can be estimated by the sample thickness and the transit time, i.e., a) subjected to 1.5 GPa shocked loading, Uel (110) ¼ 3.18 mm/0.723 ms ¼ 4.40 mm/ms; and Uel (110) ¼ 1.23 mm/0.23 ms ¼ 5.35 mm/ms; b) subjected to 1.5 GPa shocked loading, Uel (011) ¼ 3.00 mm/0.695 ms ¼ 4.32 mm/ms; Uel (011) ¼ 1.39 mm/0.271 ms ¼ 5.13 mm/ms; and c) subjected to 2.4 GPa shocked loading, Uel (011) ¼ 3.11 mm/0.726 ms ¼ 4.28 mm/ ms; Uel (110) ¼ 3.57 mm/0.85 ms ¼ 4.18 mm/ms, and d) subjected to 1.5 GPa shocked loading, Uel (010) ¼ 3.95 mm/0.93 ms ¼ 4.25 mm/ ms; Uel (010) ¼ 4.40 mm/1.196 ms ¼ 3.68 mm/ms. To probe further into the effects of anisotropy on wave propagation, numerical simulated results are compared for samples with (001), (010), (011), (100), (110), and (111) orientations. Subsequently, the variations in elastic wave velocities with different orientations are examined. Taking the 3.18 mm-thick samples under 1.5 GPa shock loading as example, the elastic wave speeds are determined for each multiple locations along propagation distance. The elastic wave velocities as a function of the propagation distance are shown in Fig. 3(a). The elastic wave velocity for the (010) orientation is the fastest, followed by the (110) and (011) orientations when arriving at the rear surface. Among these curves, the elastic precursor velocity for the (011) orientation decreases most rapidly, which varies from 5.55 to 4.20 mm/ms along the propagation distance. Elastic wave variations with different thicknesses may confirm the effects of the nonlinear elastic constitutive model for b-HMX single crystals. 3.4. Anisotropic wave behavior and effects of pressure pulse shape In the present study, the elastic precursor velocity for the (010) orientation does not seem to correspond very well with the experimental data of Dick et al. (2004a,b) literature. This inconsistency may imply that other dominant deformation mechanisms are operative for the (010) orientation, and that the dislocation slip shearing mechanism is not enough. Although the present constitutive model accounts for nonlinear elasticity and pressure-

5.6

b

(001) (010) (011) (100) (110) (111)

5.4 5.2 5.0

Plastic wave velocity(mm/ms)

Elastic wave velocity (mm/μs)

a

dependent elastic constants, the simulated results do not ensure that all orientations are completely consistent with the experimental results. To our knowledge, Previous experiments have not been performed for many orientations, such as {111}, {100}, and {001}, as a possible demonstration for anisotropy in elastic precursor velocity. For the {111} orientation, the simulated mean elastic wave velocity decays from 4369.60 to 4020.70 mm/ms going from the distance of 0.65e3. 02 mm. The {001} orientation decays from about 4586.2 to 3875.88 mm/ms, and the {100} orientation decays from 4844.22 to 4138.40 mm/ms. Likewise shown in Fig. 4(b) are the plastic wave velocity inferred from the wave profile rise times, the elastic wave velocity, and the distance between front surface and multiple target locations. Fig. 4(b) shows that plastic wave velocities are anisotropic, which increase with propagation distance and become slower, indicating that the steady plastic wave can be achieved. No direct correlativity seems to appear between elastic and plastic wave velocities for each orientation. The {111} orientation impact plane exhibits the highest plastic wave speed while the {110} orientation is the lowest. As elastic velocity increase with higher stresses, the simulations are performed under different shock strengths: 1.0, 1.5, 2.0, 2.5, and 3.5 GPa. Fig. 5(a) shows the mean elastic wave speed vs. particle velocity and their variations with orientations. Data for the elastic wave speed and the particle velocity are all deduced from the rear surface of the samples. The simulated elastic precursor speeds, which are deduced directly from this full-anisotropic constitutive model are orientation dependent. For each orientation, a linear fit to the data for elastic velocity versus. particle velocity can be written by U ¼ c þ s up. Fig. 5(a) presents the simulated results for the (010), (011), and (110) orientations, in which the linear fit relationships are U(010) ¼ 3286.45 þ 2.94 up,U(011) ¼ 3652.87 þ 5.50 up, and U(110) ¼ 3865.94 þ 8.07 up, respectively. Also shown in Fig. 5(a) are the experimental measured data which are in a very limited range; hence, demonstrating the simulated results using these experimental data is difficult. Nevertheless, the fit coefficients for the simulated results are very close to those of the experimental data as has been pointed out by Dick et al. (2004a,b). The large coefficient of the particle velocity (up) indicates stronger nonlinearity for the elasticity. Thus, the (011) and (110) orientations show stronger nonlinearities than the other orientations. Fig. 6 shows the contour plots of the accumulated slip ðgÞ at 0.8 ms. The plots show that no reflection waves arrived at the rear surfaces for all orientations. In all cases, the pressure is ramped from0 to 1.5 GPa in 0.01 ms. The contour plots shows that the slip deformation gradient sharpens at the edge zones for each orientation. Apparently, non-uniformity of the plastic deformation exists around the boundary zones because of the effects of anisotropy. The contour plots in Fig. 6 also illustrated nearly uniform deformation at the center, indicating that each element in it experiences

4.8 4.6 4.4 4.2

3.6 3.4 3.2

(001) (010) (011) (100) (110) (111)

3.0 2.8 2.6

4.0 3.8 0.4

0.8

1.2

1.6

2.0

2.4

Wave propagation distance (mm)

2.8

3.2

2.4 0.4

0.8

1.2

1.6

2.0

2.4

2.8

3.2

Wave propagation distance (mm)

Fig. 4. a) Elastic wave velocity and b) plastic wave velocity for (001), (010), (011), (100), (110), and (111)-orientations reflecting the anisotropic wave structure for HMX single crystal.

W. Yan-Qing, H. Feng-Lei / European Journal of Mechanics A/Solids 36 (2012) 66e82

6000 5500 5000

b 4800

Simulated results: {010} {011} {110} Experimental data: {010} {011} {110}

{100} {001} {111}

4600 Elastic wave velocity (m/s)

Elastic wave velocity (m/s)

a

4500 4000 3500

75

U = 3485.384 + 2.8236up

U = 3842.27 + 3.40796up

4400 4200 4000

U = 3382.615+3.0244up

3800 3600 3400

3000 0

50

100

150

200

250

300

Particle velocity (m/s)

0

50

100

150

200

250

300

350

Particle velocity (m/s)

Fig. 5. Elastic wave speed vs. particle velocity for a) (011), (110) and (010) orientations; b) (100), (001) and (111) orientations which reflect the anisotropic elastic wave structure for HMX single crystal.

identical loading. Even if the elements want to shear in the plane of impact, all elements near the center shear identically. Variations in the slip deformation field due to the different orientations can also reflect the effects of anisotropy on wave structures. The observations on six different orientations in Fig. 6 clearly showed that the impact plane of {100} orientation exhibits the largest slip deformation. When the pressure was ramped from zero to a same peak value within different time and then held constant, the deformation patterning may be influenced by pressure pulse shapes to some extent. Herein, three cases are investigated: the slowest loading rate with 0.05 ms ramp time, the medial with 0.01 ms ramp time, and the fastest loading rate with 0.002 ms ramp time. The peak value of the loading pressure is chosen to be 2.0 GPa applied to the crystal plane with (110) orientation. With the stress wave propagates from left to right, the leading stress wave front comes into being from zero to the peak value, thus creating a plastic process zone in which the plastic deformation occurs. This is illustrated by the equivalent plastic strain distribution along the centered path line at various times (i.e., 0.2 ms, 0.5 ms, and 0.8 ms) in Fig. 7. At the earlier time (0.2 ms), the plastic strain is lower and the plastic deformation process covers a zone of about 0.7 mm. The loading face has higher plastic strain since it is subjected to loading for longer time. At 0.2 ms and 0.5 ms, with the slowest loading rate, the plastic strain through the deformation zone is lower than that in the case of faster loading rate with 0.01 and 0.002 ms ramp times. With the stress wave propagation proceeds the plastic profiles tend to show multiple fluctuations due to nonlinear elasticity and anisotropic effect. Such characteristics are evident at the time 0.8 ms. A wide strain variation and higher value plastic strain at loading face can be obtained for the slowest loading rate at 0.8 ms. The wave for the slowest loading rate which originally lags behind the other two may catch up with them and even reach higher strain rate after a period of time. The plastic strain distribution curves for the two faster loading rates tend to become close to each other, which means that choosing faster loading rate in experiments can avoid much uncertainty due to different pulse shapes. 3.5. Dislocation density and resolved shear stress When the applied pressure exceeds the Hugoniot elastic limit, the material deforms plastically as a result of the dislocation motion and multiplication. As dislocation nucleation and growth mechanisms have been incorporated into the present constitutive model, the localized deformation and resolved shear stress continue to

influence the evolution of dislocations. The time histories of the dislocation density for the (001), (010), (011), (100), (110), and (111) orientations corresponding to the two activated slip systems under the same loading conditions are plotted in Fig. 8. The slip system {101}< 101 > is denoted by ‘s1’ and {001} <100> refers to ‘s2’. In the present study, the slip systems are updated for each increment, and then transformed to the global Cartesian coordination systems. These plots show that the dislocation density is very sensitive to crystal orientation with the highest value in the (001) orientation, followed by (100), and then (110) followed respectively. While the wave has not yet reached the target back surface, the dislocation density maintains its initial value. For the (100) orientation, the values of rs1 increase rapidly from their initial value to 1.29 times its initial value (at Point A) and then remain almost unchanged in the subsequent 0.22 ms; this reflects a steady state. Subsequently, rs1 increases rapidly again, and soon reaches a second saturation stage, i.e., after Point B, rs1 becomes 1.33 times its initial value. In the (001) orientation, no obvious transition stage exists while the rs1 also reaches about 1.33 times its initial value at the end. In the (011) orientation, rs1 increases initially, and then enters a gradually growth stage. Finally, the {011} orientation reaches a saturation value at Point ‘D’. In both the (110) and (011) orientations, t two steady stages for evolutional curves of dislocation density rs1 exist, as shown in Fig. 8(a). The former reaches a saturation value of 1.11 times while the latter reaches 1.08 times the initial value. At 2.0 ms (i.e. the maximum value of the x-coordinate), no longitudinal reflection wave has yet reached the target rear surface. The rs1 for the (010) orientation develops gradually. In the (111) orientation, the resolved shear stress for the slip system ‘s1’ is not higher than the threshold shear stress; thus, the dislocation density maintains its initial value. In the slip system {001} <100> (denoted by ‘s2’) shown in Fig. 8(b), the highest value of dislocation density increases to 2.0 times its initial value, which is higher than that of the slip system ‘s1’. In Fig. 8(b), the curves for the {100} orientation represents an obvious steady stage, while the (001) orientation corresponding to the slip system ‘s2’ also presents the highest dislocation density at 2.0 ms. As the analysis above demonstrates, the dislocation density is strongly influenced by the orientations, which may be attributed to the resolved shear stress along the slip systems. The same result is inferred from the constitutive equations, including the observation that the dislocation density is strongly influenced by the activated slip systems and the shocked plane. To understand the interactions

76

W. Yan-Qing, H. Feng-Lei / European Journal of Mechanics A/Solids 36 (2012) 66e82

Fig. 6. Contour plots of accumulated slip ðgÞ for different orientated single crystal samples at 0.8 ms.

between the dislocation density evolution and stress relaxation, the time histories of the resolved shear stress for different orientations are shown in Fig. 8(c) and (d). ss1 ; ss2 refer to the resolved shear stress along the two specified slip directions. As shown in Fig. 8(c) and (d), the resolved shear stresses increased to some peak value once the elastic wave reaches the rear surface. At the same time, the dislocation densities for all the orientations increase once the resolved shear stress exceeds the threshold value for the dislocation motion. in the (111) orientation, the resolved shear stress ss1 fluctuates below the threshold value; hence, the increment of the dislocation density rs1 remains 0. The

dislocation density in the {100} reaches saturation stage at point ‘A’ in Fig. 8(a); this state may have resulted from a relaxation of shear stress during the period from A0^to B0^in Fig. 8(b). This reduction in shear stress leads to a full loss of shear strength (Point ‘B0^’), so that the dislocation density begins to increase gradually. A slight variation in dislocation density leads to an intense change in shear stress. In the {110} and {011} orientations, the fluctuations in shear stress shown in Fig. 8(c) and (d) coincide with the gradually increasing dislocation density evolutions in Fig. 8(a) and (b). Two turning points (denoted by C, D and E, F) exist for each evolution

W. Yan-Qing, H. Feng-Lei / European Journal of Mechanics A/Solids 36 (2012) 66e82

77

3.6. Shock-induced temperature increase anisotropy

Fig. 7. Effects of loading pressure pulse shapes corresponding to 0.002, 0.01 and 0.05 ms ramp times on the equivalent plastic strain distribution along the wave propagation distance.

curve of dislocation density in the {011} and {110} orientations. In the second increasing stage, i.e., after D and F, the growth rates become slower and smoother. In the (001) and (010) orientations, the resolved shear stresses in both Fig. 8(c) and (d) increase to a peak value, followed by a steady process. Compared with the dislocation density evolution in Fig. 8(a) and (d), the {001} and {010} orientations have only one smoothly increasing stage, respectively, throughout the entire process. Such changing laws in shear stress and dislocation density for all orientations are reasonable.

In addition to typically measured mechanical variables (stress or pressure, shock and particle velocities), the ability to determine temperature under shock loading can also provide insight into the physical conditions that govern shock-induced structural changes and chemical reactions. As is now well known, a shock heating associated with some localized effects in an energetic material may be a significant source of hot spots in the regime of weak initiation. Note that the present model includes temperature increase because of reversible adiabatic compression; hence, the model can be used in predicting bulk temperature rise in b-HMX single crystal under impact. For adiabatic analyses, the temperature of each interior interface increases as the wave fronts pass, and the thermal interactions between wave fronts are excluded in the model. The temperature histories with time and distributions along the wave propagation direction for the b-HMX single crystal samples in the (001), (010), (011), (100), (110), and (111) orientations are shown in Fig. 9. For given similar loading conditions, the temperature histories of the center at the front and rear surfaces in Fig. 9(a) and (b) (with initial temperature of 300K). The strong orientation sensitivity of the increase in temperature is apparent. The greatest temperature increases occur in the {010} orientation and the smallest in the {100} orientation. Jindal and Dlott (1998) proposed that temperature differences may arise due to differences in compressibility, strength and Gruneisen parameters along each orientation. They used thermo-chemical data of energetic materials to show that temperature increase anisotropy is large enough to lead to reaction rates which differ by several orders of magnitude. For the (010) orientation, the temperature initially increases by 14.42K, and then remains nearly unchanged. In the (100) orientation, the temperature increases rapidly from the first step to the

a

b

c

d

Fig. 8. (a)(b) Dislocation evolutions and (c)(d) resolved shear stress evolutions for the center point at rear surface corresponding to shocked planes (011), (110) and (010), (100), (001) and (111) respectively. Among these figures, a) c) corresponds to slip system s1: {101}< 101 >; b) slip system s2: {001} <100>.

78

W. Yan-Qing, H. Feng-Lei / European Journal of Mechanics A/Solids 36 (2012) 66e82

b

a 320

c

d

Fig. 9. Temperature histories for a) center Lagrangian particle at the front surfaces; b) center Lagrangian particle at the rear surfaces and temperature profiles along wave propagation distances at c) 0.6 ms; d) 1.0 ms.

end of wave propagation. Notice that the trends of the curves in the (011) and (110) orientations are very close, in which the temperature increases by 10K in a short time, and then increases gradually along with wave propagation. In Fig. 9(a), the (001) and (100) orientations have the same final temperature (308.63K) at the end of time (2.0 ms) though their time histories are different. In Fig. 9(a) and (b), the temperature appears to increase in the center at the rear surface, similar with that at the front surface in the same orientation. To demonstrate how the model for b-HMX single crystals developed in the present study may be helpful in predicting temperature wave profiles, temperature distributions along the sample thickness in different orientations at 0.6 and 1.0 ms are shown in Fig. 9(c) and (d), respectively. The temperature distributions show that the degree and pattern of temperature uniformity in different orientations have great variations. Given that explosive grains are crystalline, the anisotropy of a crystal may be one of the sources of nonuniform temperature distribution. The compressive wave does not reach the back surface at 0.6 ms, and thus the predicted variations in temperature are not identical to those given at 1.0 ms through the compaction zone. This difference in the predicted temperature is explained by the anisotropy in both the uniaxial irreversible compressibility and the bulk reversible thermo-elastic coupling. This trend coincides with the flexibility to transfer energy and plastic deformation during shock or impact loading along different orientations. In the (100) and (001) orientations, the compressive waves have significant fluctuations in the temperature between front and rear surface. In the present study, given that the temperature-dependences of parameters are difficult to measure and are neglected during

the simulation, the model may not be suitable for very strong shock loading because higher temperature may lead to substantial changes in the physical property of b-HMX single crystals. 4. Summary and conclusions The present study develops a model of full elasticeplastic anisotropy. In comparison, previous researchers developed their models by varying the parameters with orientations. The foundation of shock initiation is based on the observations and analyses of mechanical responses of solid samples to shock loading. In the regimes of mechanical responses, a nonlinear elastic compression, crystal plasticity model and an EOS for a lowsymmetric single crystal with limited operative slip systems have been developed. Only dislocation slip mechanisms are incorporated into the constitutive model despite possible twinning and fracture mechanisms. The present study considers the volumetric couple with deviatoric behaviors, as well as thermo-elastic heating with no temperature-dependent flow stress because of lack of experimental data. The present work provides insight into the elastic and inelastic wave propagation in anisotropic b-HMX single crystal. Specially, the formulas and algorithms developed herein can be applied to study other low-symmetric crystals under high pressure loading that deform primarily by crystallographic slip.. Plasticity is important for granular HMX as well as for b-HMX single crystal. A range of complex, elasticeplastic behaviors of HMX crystalline under shock loading are encountered. Elastic anisotropy and plastic anisotropy can be instructive if they are described in the high-pressure state. Thus, the present study

W. Yan-Qing, H. Feng-Lei / European Journal of Mechanics A/Solids 36 (2012) 66e82

develops a more realistic description of the mechanical responses of b-HMX single crystal using a crystal plasticity model combined with nonlinear elasticity for shock loading simulations. The use of rate-dependent crystal plasticity based on a dislocation approach has typically led to descriptions of single crystal response containing more adjustable parameters that can be independently verified. In the simulation of the present study, a b-HMX plate sample is subjected to a shock loading of 1.5 GPa and 2.4 GPa, which coincides with the experimental conditions in literature. The simulated wave structures are then compared with the measured data by Dick et al. (2004a). An important feature of the model is its capability to capture the evolution of elasticeplastic wave profiles. Simulation results of the planar impacts on b-HMX single crystal are in good agreement with experimental data on particle velocity histories in the (110) and (011) orientations. Consequently, wave properties and wave number of the single crystals also become complex. The particle velocity vector consisting of nonzero transverse components with respect to the longitudinal axis in the for (110) and (011) orientations are examined. Results clearly manifest the effects of crystal anisotropy on the deviations from pure mode wave propagation. Controlled by the monoclinic crystal symmetry, longitudinal loading causes the propagation of quasi-longitudinal and quasi-shear waves whose amplitudes are limited by their shear strength. To describe further its anisotropic property, simulations are performed for arbitrary orientations (001), (010), (011), (100), (110), and (111) using the same loading conditions. The contour plots on the orientation-dependent accumulated slip deformation are compared. The plastic deformation appears to exhibit non-uniformity around the boundary zones because of the anisotropy. Variations in slip deformation field due to different orientations can also reflect the effects of anisotropy on wave structures. The variation in elastic wave speed with particle velocity and orientation is worth observing. The fit coefficients of the simulation results are very close to those of the experimental data by Dick et al. (2004a,b). The (011) and (110) orientations show stronger nonlinearities than the other orientations. Although the existing experiments have not been performed for many orientations, such as {111}, {100}, and {001}, simulated anisotropy in the elastic precursor speed versus particle velocity for these orientations may be demonstrated by future experiments. In contrast to elastic wave velocity, plastic wave velocity increases with propagation distance and becomes slower, indicating that a steady plastic wave can be achieved. Seemingly, no direct correlation between elastic and plastic wave velocity for each orientation exists; that is, the {111} orientation impact plane exhibits the highest plastic wave speed while the {110} orientation exhibits the lowest. Given that dislocation nucleation and growth mechanisms are incorporated into the present constitutive model, slip deformation and resolved shear stress continue to influence the evolution of dislocations. The dislocation density is very sensitive to crystal orientation which may be attributed to the resolved shear stress along the slip systems. The present study indicates that a slight variation in dislocation density leads to an intense change in shear stress. All these simulation results demonstrate that dislocation density evolution and changes in shear stress are strongly interdependent. Hot spots are subgrain to an extent, consequently, simulations of hot spot initiation require that the temperature inside an individual grains is determined. The model developed in the present study can indirectly account for the role contributed by bulk thermal heating to hot spot mechanisms. Although the mechanism for the initiation of decomposition chemistry has not been determined, the behavior

79

observed may be responsible for the possible anisotropy of shock initiation sensitivity in HMX crystals. The present model presents a nonlinear elastic and anisotropic plastic model, which can be used to establish correlations between orientation-dependent wave behavior and dislocation slip mechanisms. More complete experimental data, including the shocked mechanical and chemical properties of b-HMX along other orientations, should be the main task of future research. Fortunately, theoretical work toward describing the decomposition and deformation mechanisms is progressing (Strachan et al., 2003; Margetis et al., 2002). What is more, in addition to slip system, the twin system and the fracture should also be involved in the model. The present model may be easily extended to include twinning and cleavage fracture and then integrated in detonation initiation of single crystal explosives. Acknowledgment The authors would like to thank the Chinese National Nature Science Foundation (Grant No. 10902016 and 10832003), and the project of State Key Laboratory of Science and Technology (Beijing Institute of Technology) with Project No. ZDKT10-03A for supporting this project. Appendix The details on the transformation between monoclinic crystallographic slip systems and orthogonal coordination system are outlined here. The lattice parameters are specified in space groups of P21/n. For the monoclinic b-HMX single crystals with symmetrical axis, the 13 elastic constants are defined with reference three orthogonal axes e1, e2, e3 as shown in Fig. 1. Zaug (1998) determined the elastic constants for the b-HMX monoclinic crystals at 24  C.

Fig. A1. Orthogonal coordination system for defining initial elastic coefficient tensor of b-HMX single crystal with three axes of (e1, e2, e3), where e2 corresponds to the b axis and c to e3 respectively, a-axis is perpendicular to b and c.

Elastic relationship described by a 6  6 elastic stiffness matrix form is given by,

2

3

2

s11 C1111 6 s22 7 6 C1122 6 7 6 6 s33 7 6 C1133 6 7 6 6 s12 7 ¼ 6 0 6 7 6 4 s13 5 4 C1113 s23 0

C1122 C2222 C2233 0 C2213 0

C1133 0 C2233 0 C3333 0 0 C1212 C3313 0 0 C1223

32 3 C1113 0 3 11 7 6 C2213 0 7 76 3 22 7 7 6 C3313 0 7 76 3 33 7 6 3 12 7 0 C1223 7 76 7 C1313 0 54 3 13 5 0 C2323 3 23

(A1)

80

W. Yan-Qing, H. Feng-Lei / European Journal of Mechanics A/Solids 36 (2012) 66e82

C1111 C2222 C3333 C1212 C1313 C2323

where

¼ ¼ ¼ ¼ ¼ ¼

20:8; C1122 ¼ 4:0; C1133 ¼ 13:0; C1113 ¼ 0:6 26:9; C2233 ¼ 6:6; C2313 ¼ 1:5; 17:6; C3313 ¼ 0:1; 2:9; C1223 ¼ 3:0; 6:6; 3:8 ðunits : GpaÞ

Sheen et al., 1993 observed the slip systems of (001) <100> andð101Þ < 101 >. Dick et al. (2004a) suggested that the conventional mechanisms of plastic flow, i.e., dislocation slip are operative in this weak, brittle, molecular crystal. Given that the index (001) <100> and ð101Þ < 101 > are for Bravais primitive cell based on monoclinic crystallographic axis (a, b, and c) in the space group P21/ n. the transformation from crystallogrophic system to rectangular coordination system (e1, e2, e3) should be performed through a rotation matrix Rs,

0

Rs ¼ @

sin l 0 cos l

0 1 0

1

0 0A 1

(A3)

after Zaugð1998Þ

dða cosðl  p=2ÞÞ da dl ¼ sin l þ a cos l dT dT dT

ax ¼

0

sin l @ 0 cos l 0

sin l @ 0 cos l

0 1 0

2 12 3 3 1 0 sin l ð1Þ 0 A4 0 5/s0 ¼ 4 0 5 0 cos l 1

(A4)

0 1 0

2 3 12 3 0 0 0 ð1Þ 0 A4 0 5/m0 ¼ 4 0 5 1 1 1

(A5)

Thus, three expansion coefficients in global rectangular coordination system e1, e2, e3 can be calculated by the means described by Wang et al. (2005),

2

2

3

2 12 3 3sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 sin l sin l 0 0 1 5 @ 0 1 0 A4 0 5/sð2Þ ¼ 4 0  0 2 2 sin l þ ð1  cos lÞ cos l 0 1 ðcos l  1Þ 1 (A6) 2 12 3 3sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 sin l sin l 0 0 1 ð2Þ 4 A 4 5 5 @ 0 0 /m0 ¼ 1 0 0 2 2 l sin þ ð1 þ cos lÞ 1 cos l 0 1 1 þ cos l 0

(A7) The thermal expansion properties of b-HMX are measured using x-ray diffraction. They were found to be highly anisotropic. The linear expansion coefficients are (0.29, 12, and 2.3)  105 K1 along the (a, b, and c) axes, and d (ln l)/dT ¼ 2.6  105 K1 (Herrmann et al., 1993, 1992). For convenience in numerical calculation, a second-rank thermal expansion tensor a for b-HMX in orthogonal coordination system (e1, e2, e3) is written as,

a11

0

a13

3

a ¼ 4 0 a22 0 5 a13 0 a33

3

2

3

a11 ac 0:283  105 4 a13 5 ¼ ½T r 4 ax 5 ¼ 4 0:307  105 5 a33 aa 2:3  105 2

0 6 6 1 ½Tr ¼ ½x ¼ 6 0:0483 4 0:4289 1

1 0:9517 0:4289 0

0

(A10)

3

7 7 1 7  0:4289 5 0

(A11)

2

0

2

(A9)

¼ 0:2829  105

where, The primitive monoclinic crystallographic constants for the unit cell of b-HMX are: a ¼ 0.654 nm, b ¼ 1.105 nm, c ¼ 0.73 nm, l ¼ 102.06 , respectively (Palmer and Field, 1982). The initial slip ðaÞ ðaÞ direction denoted by s0 and the slip plane with normal m0 are then rotated to global orthogonal coordination system (e1, e2, e3) as,

(A2)

(A8) 

In the (010) plane, three orientations are chosen: c-axis (xc ¼ 0 ), a  axis (xa ¼ 102.7 ), and the direction xx ¼ 90 perpendicular to the c axis. Linear thermal expansion coefficient along the xx ¼ 90 direction is calculated as follows:

3 sin2 xc sin 2xc cos2 xc x ¼ 4 sin2 xx sin 2xx cos2 xx 5 sin2 xa sin 2xa cos2 xa 2 3 0 0 1 5 ¼ 4 1 0 0 0:9517 0:4289 0:0483

(A12)

In summary, the second-rank thermal expansion tensor a can be written as

2

0:283 a¼ 4 0 0:307

0 12:0 0

3 0:307 0 5  105 K 1 2:3

(A13)

For anisotropic materials, the Grüneisen tensor will change with varying properties in space which means that the assumed EOS becomes a function of the orientation of the crystallites. (Leiber, 2000),

Gij ¼

1

rCV

Cijkl akl

(A14)

References Alexander, A. Lukyanov, 2008. Constitutive behaviour of anisotropic materials under shock loading. International Journal of Plasticity 24 (1), 140e167. Anderson, C.E., Cox, P.A., Johnson, G.R., Maudlin, P.J., 1994. A constitutive formulation for anisotropic materials suitable for wave propagation computer programII. Computational Mechanics 15, 201e223. Armstrong, R.W., Elban, W.L., 2004. Dislocation in energetic crystals. In: Dislocation in Solids, vol. 12, p. 406 (Chapter 69). Armstrong, R.W., Elban., W.L., 1999. Dislocation characteristics in energetic crystals. In: Furnish, M.D., Chhabildas, L.C., Hixson, R.S. (Eds.), Shocked Compression of Condensed Matter. 2000. American Institute of Physics. Armstrong, R.W., Coffey, C.S., Elban, W.L., 1982. Adiabatic heating at a dislocation pile-up avalanche. Acta Metallurgica 30, 2111e2116. Armstrong, R.W., 1995. Dislocation mechanisms for shock-induced hot spots. Journal of Physics IV. France 05 C4-89-C4-102.

W. Yan-Qing, H. Feng-Lei / European Journal of Mechanics A/Solids 36 (2012) 66e82 Armstrong, R.W., 2009. Dislocation mechanics aspects of energetic material composites. Reviews on Advanced Materials Science 19, 13e40. Asaro, R.J., Needleman, A., 1985. Overview .42. Texture development and strainhardening in rate dependent polycrystals. Acta Metallurgica 33 (6), 923e953. Asaro, R.J., Rice, J.R., 1977. Strain localization in ductile single crystals. Journal of the Mechanics and Physics of Solids 25 (5), 309e338. Asaro, R.J., 1983a. Micromechanics of crystals and polycrystals. Advances in Applied Mechanics 23, 1e115. Asaro, R.J., 1983b. Crystal plasticity. Journal of Applied Mechanics Transactions of the ASME 50, 921e934. Bardenhagen, S.G., Brydon, A.D., Williams, T.O., Collet, C., 2006. Coupling grain scale and bulk mechanical response for PBXs using numerical simulations of real microstructures. In: Furnish, M.D., Elert, M., Russell, T.P., White, C.T. (Eds.), Shock Compression of Condensed Matter-2005, pp. 479e482. Barton, N.R., Winter, N.W., Reaugh, J.E., 2009. Defect evolution and pore collapse in crystalline energetic materials. Modelling and Simulation in Materials Science and Engineering 17, 1e19. Becker, R., 2004. Effects of crystal plasticity on materials loaded at high pressures and strain rates. International Journal of Plasticity 20 (11), 1983e2006. Choi, C.S., Boutin, H.P., 1970. A study of the crystal structure of b-Cyelotetramethylene tetranitramine by neutron diffraction. Acta Crystallogr B 26, 1235e1240. Coffey, C.S., Sharma, J., 1998. Initiation of crystalline explosives due to energy dissipated during plastic flow. In: Proc of 11th Symp (Int) on Detonation, Snowmass, CO, Aug 4eSept 4, pp. 751e757. Coffey, C.S., Sharma, J., 1999. Plastic deformation, energy dissipation, and initiation of crystalline explosives. Physcial Review B 60 (13), 9365e9371. Coffey, C.S., Sharma, J., 2001. Lattice softening and failure in severely deformed molecular crystals. Journal of Applied Physics 89, 4797e4802. Dick, J.J., Ritchie, J.P., 1994. Molecular mechanics modeling of shear and the crystal orientation dependence of the elastic precursor shock strength in pentaerythritol tetranitrate. Journal of Applied Physics 76 (5), 2726e2737. Dick, J.J., Mulford, R.N., Spencer, W.J., et al., 1991. Shock response of pentaerythritol tetranitrate single crystals. Journal of Applied Physics 70 (7), 3572e3587. Dick, J.J., Hooks, D.E., Menikoff, R., Martinez, A.R., 2004a. Elastic-plastic Wave Profiles in Cyclotetramethylene Tetranitramine Crystals. LA-UR-04e1229. March 5. Dick, J.J., Hooks, D.E., Menikoff, R., Martinez, A.R., 2004b. Elastic-plastic wave profiles in cyclotetramethylene tetranitramine crystals. Journal of Applied Physics 96 (1), 374e379. Dick, J.J., 1982. Plane shock initiation of detonation in g-irradiated pentaerythritol tetranitrate. Journal of Applied Physics 53 (9), 6161e6167. Dick, J.J., 1984. Effect of crystal orientation on shock initiation sensitivity of pentaerythritol tetranitrate explosive. Applied Physical Letters 44, 859e861. Dick, J.J., 1993. Orientation-dependent explosion sensitivity of solid nitromethane. The Journal of Physical Chemistry 97 (23), 6193e6196. Dick, J.J., 1997. Anomalous shock initiation of detonation in pentaerythritol tetranitrate crystals. Journal of Applied Physics 81 (2), 601e612. Duvall, G.E., Fowles, G.R., 1963. In: Bradley, R.S. (Ed.), 1963. High Pressure Physics and Chemistry, vol. 2. Academic, New York, pp. 209e291. Elban, W.L., Hoffsommer, J.C., Armstrong, R.W., 1984. Journal of Materials Science 19, 552. Field, J.E., Swallowe, G.W., Heavens, S.N., 1982. Ignition mechanisms of explosives during mechanical deformation. Proceedings of the Royal Society of London. A, Mathematical and Physical Science 382 (1782), 231e244. Flanagan, D.P., Taylor, L.M., 1987. An accurate numerical algorithm for stress integration with finite rotations. Computer Methods in Applied Mechanics and Engineering 62 (3), 305e320. Gallagher, H.G., Halfpenny, P.J., Miller, J.C., Sherwood, J.N., Tabor, D., 1992. Dislocation slip systems in Pentaerythritol Tetranitrate (PETN) and Cyclotrimethylene Trinitramine (RDX) [and Discussion] 339, 293e303. Gilman, J.J., Johnston, W.G., 1962. In: Seitz, F., Turnbull, D. (Eds.), 1962. Solid State Physics, vol. 13. Academic, New York. Gilman, J.J., 1969. Micromechanics of Flow in Solids. McGraw-Hill, New York. Gray III, G.T., Lopez, M.F., Bourne, N.K., Millett, J.C.F., Vecchio, K.S., 2002. Influence of microstructural anisotropy on the spallation of 1080 eutectoid steel. In: Furnish, M.D., Thadhani, N.N., Horie, Y. (Eds.), Shock Compression of Condensed Matter-2001. AIP Press, Melville, NY, pp. 479e483. Gruzdkov, Y.A., Gupta, Y.M., 2000. Shock wave initiation of pentaeythritol tetranitrate single crystals: mechanism of anisotropic sensitivity. The Journal of Physical Chemistry A 104 (47), 11169e11176. Gupta, Y.M., Duvall, G.E., Fowles, G.R., 1975. Dislocation mechanisms for stress relaxation in shocked LiF. Journal of Applied Physics 46 (2), 532e546. Herrmann, M., Engel, W., Eisenreich, N., 1992. Thermal expansion, transitions, sensitivities and burning rates of HMX. Propellants, Explosives, Pyrotechnics 17, 190e195. Herrmann, M., Engel, W., Eisenreich, N., 1993. Thermal analysis of the phases of HMX using x-ray diffraction. Zeitschriff fur Kristallographie 204, 121e128. Hill, R., Rice, J.R., 1966. Generalized constitutive relations for incremental deformation of metal crystals by multislip. Journal of the Mechanics and Physics of Solids 14 (2), 95e102. Hill, R., Rice, J.R., 1972. Constitutive analysis of elasticeplastic crystal at arbitrary strain. Journal of the Mechanics and Physics of Solids 20 (6), 401e413. Hoger, A., 1987. The stress conjugate to logarithmic strain. International Journal of Solids and Structures 23 (12), 1645e1656.

81

Hooks, D.E., Ramos, K.J., Martinez, A.R., 2006. Elasticeplastic shock wave profiles in oriented single crystals of cyclotrimethylene trinitramine (RDX) at 2.25 GPa. Journal of Applied Physics 100. 024908-1w5. Hughes, T.J.R., Winget, J., 1980. Finite rotation effects in numerical integration of rate constitutive equation arising in large-deformation analysis. International Journal of Numerical Methods in Engineering 15 (12), 1862e1867. Hull, D., Bacon, D.J., 1984. Introduction to Dislocations, third ed. Pergamon Press Ltd, Heading ton Hill Hall. Oxford OX 3 0BW. England. Jindal, V.K., Dlott, D.D., 1998. Orientation dependence of shock-induced heating in an harmonic molecular crystals. Journal of Applied Physics 83 (10), 5203e5211. John, H., Lothe, J., 1982. Theory of Dislocations, second ed. McGraw-Hill, New York. Johnson, G.C., Bammann, D.J., 1984. A discussion of stress rates in finite deformation problems. International Journal of Solids and Structures 20 (8), 725e737. Johnson, J.N., Jones, O.E., Michaels, T.E., 1970. Dislocation dynamics and singlecrystal constitutive relations: shock-wave propagation and precursor decay. Journal of Applied Physics 41 (6), 2330e2339. Johnson, J.N., 1971. Shock propagation produced by planar impact in linearly elastic anisotropic media. Journal of Applied Physics 42 (13), 5522e5530. Johnson, J.N., 1972. Calculation of plane wave propagation in anisotropic elasticeplastic solids. Journal of Applied Physics 43 (5), 2074e2082. Kalidindi, S.R., Bronkhorst, C.A., Anand, L., 1992. Crystallographic texture evolution in bulk deformation process of fcc metals. Journal of the Mechanics and Physics of Solids 40 (3), 537e569. Kolsky, H., 1953. Stress Wave in Solids. Clarendous Press, Oxford. Leiber, Carl-Otto, 2000. Aspects of the mesoscale of crystalline explosives. Propellants, Explosive, Pyrotechnics 25 (6), 288e301. Ling, X., Horstemeyer, M.F., Potirniche, G.P., 2005. On the numerical implementation of 3D rate-dependent single crystal plasticity formulations. International Journal for Numerical Methods in Engineering 63, 548e568. Margetis, D., Kaxiras, E., Eisner, M., Frauenheim, Th., Manaa, M.R., 2002. Electronic structure of solid nitromethane: effects of high pressure and molecular vacancies. Journal of Chemistry Physics 117, 788. Menikoff, R., Kober, E., 2000. Compaction waves in granular HMX. In: Furnish, M.D., Chhabildas, L.C., Hixson, R.S. (Eds.), Shock Compression of Condensed matter1999. American Institute of Physics, New York, pp. 397e400. Menikoff, R., Sewell, T.D., 2001. Equation of State Data and Fitting Forms. In: Detonation Theory & Application, T-14. Theoretical Division. A portion of LAUR-01-1847, May. Menikoff, R., Sewell, T.D., 2003. Complete Equation of State for b-HMX and Implications for Initiation. Shock Compression of Condensed Matter-2003. AIP Conference Proceedings, Melville, NY. 157e160. Menikoff, R., Dick, J.J., Hooks, D.E., 2004. Analysis of Wave Profiles for Single Crystal HMX. LA-UR-04-3928. Menikoff, R., Dick, J.J., Hooks, D.E., 2005. Analysis of wave profiles for single-crystal cyclotetramethylene tetranitramine. Journal of Applied Physics 97. 023529-1w6. Ortiz, M., Radovitzky, R.A., Repetto, E.A., 2001. The Computation Of The Exponential And Logarithmic Mappings And Their First And Second Linearizations. Cal. Tech ASCI Technical Report 80. Palmer, S.J.P., Field, J.E.,1982. The deformation and fracture of b-HMX. Proceedings of the Royal Society of London. A, Mathematical and Physical Science 383 (1785), 399e407. Peirce, D., Asaro, R.J., Needleman, 1983. Material rate dependence and localized deformation in crystalline solids. Acta Metallurgica 31, 1951e1976. Plaksin, I., Coffey, C.S., 2007. Reaction localization in crystalline explosive subjected to strong shock. In: Proceedings of the International Autumn Seminar on Propellants, Explosives and Pyrotechnics 259e272. Xi’an, Shaanxi, China. Sewell, T.D., Bedrov, D., Menikoff, R., Smith, G.D., 2002. Elastic properties of HMX. In: Furnish, M.D., Thadhani, N.N., Horie, Y. (Eds.), Shock Compression of Condensed Matter 2001, pp. 399e402. Sheen, D.B., Sherwood, J.N., Gallagher, H.G., Littlejohn, A.H., Pearson, A., 1993. An Investigation Of Mechanically Induced Lattice Defects In Energetic Materials, Office of Naval Research Final Report N00014-87-J-1173. 53 pp. Sherwood, J.N., Sheen, D.B., ONR, 1988. Workshop on Energetic Material Initiation Fundamentals. CPIA Publication 516. 638 p. Strachan, A., van Duin, A.C.T., Chakraborty, D., Dasgupta, S., Goddard III, W.A., 2003. Shock waves in high-energy materials: the initial chemical events in nitramine RDX. Physical Review Letters 91, 098301e098304. Vignjevic, R., Bourne, N.K., Millett, J.C.F., De Vuyst, T., 2002. Effects of orientation on the strength of the aluminum alloy 7010-T6 during shock loading: experiment and simulation. Journal of Applied Physics 92 (8), 4342e4348. Wallace, D.C., 1980. Irreversible thermodynamics of flow in solids. Physical Review B 22 (4), 1477e1486. Wang, Kunpeng, Zhang, Jianxiu, Wang, Jiyang, Zhang, Huaijin, Wang, Zhengping, Yu, Wentao, Wang, Xuping, Lu, Qingming, Ba, Mingfang, 2005. Anisotropic thermal expansion of monoclinic potassium lutetium tungstate single crystal. Journal of Applied Physics 98. 046101-1w3. Winey, J.M., Gupta, Y.M., 2004. Nonlinear anisotropic description for shocked single crystals: thermoelastic response and pure mode wave propagation. Journal of Applied Physics 96 (4), 1993e1999. Winey, J.M., Gupta, Y.M., 2006a. Nonlinear anisotropic description for the thermomechanical response of shocked single crystals: inelastic deformation. Journal of Applied Physics 99. 023510-11w11. Winey, J.M., Gupta, Y.M., 2006b. Anisotropic modeling for shocked single crystals. In: Furnish, M.D., Elert, M., Russell, T.P., White, C.T. (Eds.), Shock Compression of Condensed Matter-2005, pp. 367e372.

82

W. Yan-Qing, H. Feng-Lei / European Journal of Mechanics A/Solids 36 (2012) 66e82

Zamiri, A.R., De, S., 2010. Deformation distribution maps of b-HMX molecular crystals. Journal of Physics D: Applied Physics 43, 035404e035410. Zamiri, A.R., De, S., 2011. Modeling the anisotropic deformation response of b-HMX molecular crystals. Propellants Explosives Pyrotechnics 36, 247e251.

Zaug, J.M., 1998. Elastic constants of b-HMX and tantalum, equation of state of supercritical fluids and fluid mixtures and thermal transport determinations. In: 11th International Detonation Symposium. Snowmass Conference Center, Snowmass Village, Colorado, pp. 498e509.