Analysis of thermomechanical response of polycrystalline HMX under impact loading through mesoscale simulations

We investigate the response of polycrystalline HMX (Octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine) under impact loading through a 3-dimensional mesoscale model that explicitly accounts for anisotropic elasticity, crystalline plasticity, and heat conduction. This model is used to quantify the variability in temperature and stress fields due to random distributions of the orientations of crystalline grains in HMX under the loading scenarios considered. The simulations carried out concern the response of fully dense HMX polycrystalline ensembles under impact loading at imposed boundary velocities from 50 to 400 m/s. The polycrystalline ensemble studied consists of a geometrically arranged distribution of bi-modally sized and shaped grains. To quantify the effect of crystalline slip, two models with different numbers of available slip systems are used, reflecting differing characterizations of the slip systems of the HMX molecular crystal in the literature. The effects of microstructure and anisotropy on the distribution of heating and stress evolution are investigated. The results obtained indicate that crystalline response anisotropy at the microstructure level plays an important role in influencing both the overall response and the localization of stress and temperature. The overall longitudinal stress is up to 16% higher and the average temperature rise is only half in the material with fewer potential slip systems compared to those in the material with more available slip systems. Local stresses can be as high as twice the average stresses. The results show that crystalline anisotropy induces significant heterogeneities in both mechanical and thermal fields that previously have been neglected in the analyses of the behavior of HMX-based energetic materials.


I. INTRODUCTION
Under dynamic loading, energetic crystals may ignite due to the localization of dissipated energy into small, intensely heated regions termed "hotspots".2][3][4][5] Heating due to each of these mechanisms depends on the histories of the states of local stress and rate of deformation in the material.8][19][20][21] Analyses utilizing discrete methods have yielded understanding of many of the thermo-mechanical properties of the polymorphs of HMX including crystal structure and lattice parameters, 11 isotherms, 10 elastic constants, 7 viscosity, 6 thermal conductivity, and behavior in compression. 8On the other hand, the continuum analyses have considered anisotropic yield surfaces with 17 or without [18][19][20][21][22] circular (spherical) voids.The studies that consider PBXs as collections of heterogeneously distributed particles with isotropic properties are able to capture many of the complex interactions that arise from the mesostructure of the granular composites.2D and 3D Eulerian simulations like those of Benson 23 and Baer 24 use isotropic material descriptions to analyze the stress and temperature fields generated by shock loading of porous granular HMX.Lagrangian simulations have been carried out to analyze the effect on ignition sensitivity of inter-particle friction (Soulard et al. 28 and Panchadhara and Gonthier) 29 and the effects of frictional dissipation along both inter-and intra-granular cracks on hotspot formation (Barua et al.). 25Specifically, the study of Panchadhara and Gonthier concerns ensembles of initially porous granular HMX.The results showed the importance of both plasticity and friction along the surfaces of HMX grains in heating and hotspot formation.The studies of Barua et al. have concerned several issues, including the development of a framework to account for fracture, contact and frictional heating in microstructures of PBXs, 25 contributions to heating of dissipation through fracture, friction in both HMX grains and the polymer binder (as well as along the grain-binder interfaces) 30 and the effects of the transient nature of loading on the contributions of the various heating mechanisms. 31The study of Mas et al. 32 accounts for grain scale features in a coarse grained simulation by using the method of cells to model split Hopkinson pressure bar tests on PBXs.The results show that fracture plays a significant role at this scale on the deformation of the PBXs, due to large stress gradients between different phases.Even with this level of grain detail, it is still not possible to isolate specific localization sites due to the interaction of mesostructural organization, friction, and plasticity.Bardenhagen and Brackbill 33 investigated stress bridging and stress fingering in lattices of circular (cylindrical) grains.Using the particle and cell computational technique, they found the properties of the binder material to be the dominant factor affecting the stress distribution and localization in the grains.On the other hand, Bardenhagen et al. 34 embedded microstructures of a PBX obtained by X-ray microtomography in General Interpolation Material Point method simualtions of response under compressive loading.The method offers an explicit link between micromechanical information and macroscale models.The technique uses the stochastic transformation field analysis to capture multi-scale effects.The result is a statistical macroscale model capable of modeling the bulk material effects and strain distributions in the material.Taken together, these studies have clearly shown that material heterogeneities at different scales, elastic and inelastic deformation mechanisms, and dissipation through a range of bulk and interfacial avenues combine to affect the heating and development of hotspots in the materials.Despite the wide range of analyses reported, thus far no analysis has been carried out to account for the crystalline nature of HMX grains, the primary constituent in many PBX composites and granular explosives (GX).This lack of study for molecular crystals is in contrast to the extensive work that has been carried out for polycrystalline metals, including Cu, [35][36][37][38] Ta, 39,40 titanium alloys, 41,42 and steels 43 among many others.These studies have shown that the polycrystalline nature of materials significantly affect local mechanical and thermal behaviors.
The primary objective of this paper is to analyze and quantify the effect of crystalline anisotropy and polycrystalline microstructure on the response of HMX-based PBXs and GXs.To this effect, the 3D Lagrangian framework developed by Rimoli et al. 44 for polycrystalline PETN is adopted.The specific form of HMX analyzed is β-HMX which has a monoclinic structure. 27The analyses focus on the effects of polycrystalline slip and microstructure attributes on the early-time histories and heterogeneous variations of stress and temperature fields for imposed boundary velocities from 50 to 400 m/s.This paper does not deal with initiation.The model used does not account for the binder phase in PBXs, voids, cracks or other defects.In essence, what we analyze are polycrystalline aggregates of HMX grains with perfect bonding.This choice of framework allows us to fully quantify of the effects of (1) the structure of crystalline slip systems in HMX and (2) the microstructural level crystalline anisotropy on the conditions in energetic materials under prescribed loading and boundary conditions.In doing so, we provide a measure for the contribution to heating of crystalline plasticity as one of the many sources that hitherto has not been fully quantified for granular or polymer-bonded energetic materials.

II. ENERGETIC MOLECULAR CRYSTALS
The range of possible bonding strength between grains in PBXs and GXs can vary greatly, from no bonding (packed granular particles) to perfect bonding.Here, we only consider the case of perfect bonding, therefore, the systems can be regarded as polycrystalline ensembles in which the grain boundaries do not fail under conditions analyzed.
The polycrystalline model here can be regarded as an idealization of PBXs in which the thickness of the binder between adjacent grains is much smaller than the size of the grains so that it is taken to be zero.This approach simplifies the material condition and allows us to delineate the effect of intrinsic material anisotropy on the behavior of the overall material.The ensemble of crystals is generated using the relaxed dual complex (RDC) method proposed by Rimoli and Ortiz. 45The RDC method takes as input an initial triangulation of the domain.The triangulation is then refined appropriately to provide suitably sized grains.Barycentric subdivision is performed on the entire refined triangulation and grains are defined as the barycentric dual of the nodes in the refined triangulation.The grain boundaries are then relaxed using a grain boundary energy minimization scheme to achieve the final grain morphologies.At this point, the grain ensemble is geometrically arranged into a regular packing of two distinct grain morphologies and sizes, as seen in Figure 2. The resulting microstructure is a polycrystalline HMX with each grain having a random crystallographic orientation.It should be pointed out that, although this method may not yield grain morphologies that perfectly match those observed in experiments, it provides an efficient way to generate microstructures with desired attributes such as grain sizes and grain size distributions.These microstructures can allow trends in material behavior to be explored.More realistic generation of 3D microstructures (see, e.g.,) 46,47 is not the focus of this study.Regardless, the simulation framework adopted here lends itself to studies which could use microstructures generated using other methods.
HMX (octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine) is the primary energetic constituent in many PBXs, such as PBX 9501.This material is a polymorphic molecular crystal which has α, β, δ, and γ phases. 26For determining the critical behaviors of HMX, the β and δ phases are of the greatest importance.The β phase has a monoclinic crystal structure and is stable at room temperature and pressure.It is commonly used in the initial formulations for PBXs 27,48 and is the polymorph considered in this study.At 438 K, the β phase begins to transform to the δ phase, which has a hexagonal crystal structure, according to the Arrhenius kinetics. 49,50 t 521 K, the melting temperature, the δ phase is more stable than the β phase at any pressure. 51Under certain conditions, cracking can play a significant role in the response of β-HMX 2,17,52,53 and extensive analyses have been carried out to quantify the dissipation associated with fracture and post-fracture crack-face friction.The consideration of fracture and friction requires a more complicated model and is beyond the scope of this paper.Such a model also has a disadvantage in that the cross influences between fracture and other dissipation mechanisms can make the delineation of the effects of specific mechanisms (crystalline slip in this paper) challenging.Here, we opt for a framework that focuses on the anisotropic crystalline response by accounting for both monoclinic elasticity and crystalline plasticity and by ignoring phenomena associated with fracture.Acknowledging that this is an idealization, we note that the benefit of this approach is the opportunity to systematically quantify the scale and extent of oscillations in the stress, deformation, and temperature fields in the materials resulting from the anisotropy of material response and microstructure heterogeneity.

A. Elasticity -Monoclinic Crystal Symmetry
The elastic properties of β-HMX have been thoroughly studied by way of nondestructive measurement techniques such as impulsive stimulated light scattering (ISLS), 54 Brillouin scattering, 55 and MD simulations. 7Table I contains sets of elastic constants obtained using each of these three methods.
Due to the monoclinic crystal symmetry in β-HMX, thirteen elastic constants fully define the elastic stiffness tensor.The elastic strain energy density can be written in the form of  where ε e is the recoverable (elastic) portion of the strain tensor and C e is the elastic stiffness tensor which can be expressed as At high loading rates, pressure gradients cause gradients in the elastic moduli in the loaded and unloaded regions of the material.Therefore, a full representation of the anisotropic material behavior would include the evolution of each of the elastic moduli with respect to pressure.Such an analysis is beyond the scope of this investigation, but should be pursued when reliable, pressure dependent, anisotropic moduli are identified.
It should be noted that obtaining a complete set of temperature dependent constants would be advantageous due to the fact that the material undergoes significant heating during impact loading.Zaug 54 determined five of the necessary constants for two different temperatures using ISLS, but a lack of certainty on the remaining constants leads us to consider the other two data sets.An analysis of potential elastic constants has been done by Vial et al., 56 leading to the identification of the constants of Stevens and Eckhardt 55 as better than those of Zaug 54 because of the sample sizes used in the analyses.Zamiri 18 compares the responses obtained using these three sets of constants to the response experimentally measured by Rae et al. 52 in the [110] and [001] crystalline directions.It is found that the constants of Sewell et al. 7 provide the best fit to the experimental data.Based on this, we adopt the constants of Sewell et al. 7 in this study.The incorporation of fully pressure-and temperature-dependent elastic moduli is likely a future task when more complete data is available.

B. Crystalline Plasticity and Thermal Response
Relative to elasticity, considerably less information concerning the inelastic deformation of β-HMX exists in the literature.Analysis of impact experiments on specially grown single crystals has previously led to a belief that slip along the (001)[100] and ( 102) [201] systems (in P2 1 /c notation) as well as twinning and cleavage 57 can occur in β-HMX.The slip along these two systems was subsequently considered in a continuum model to analyze the yield surface of HMX single crystals under quasi-static loading. 18,19 ore recently, MD simulations show that there are five additional potential slip systems: ( 11 17 Barton et al. have implemented a crystalline plasticity model with all seven slip systems to study the development of hotspots in the vicinity of an intragranular pore in a single crystal under shock loading with pressures of the order of 10 GPa. Figure 1 illustrates all seven slip systems relative to the unit cell in the P2 1 /c space.The slip systems labeled in green also correspond to the two systems used by Zamiri.Due to the uncertainty concerning the slip behavior of β-HMX, we will consider cases with both two-and seven-slip systems and compare the differences in behavior for both cases.For clarity, we will refer to the case with two slip systems as the Zamiri case and the case with seven slip systems the Barton case.All seven slip systems and relevant material parameters are summarized in Table II.Note that the two slip systems in the Zamiri case are among the seven slip systems in the Barton case.
Configuration for analysis and microstructure model.
The crystal plasticity model outlined below follows the formulation of Rimoli et al.. 44 The local deformation condition is described by the displacement gradient where u is the displacement field and ∇ is the gradient operator.The plastic portion of the displacement gradient is In the above relation, N is the total number of slip systems, α is the number of a slip system, γ α is the magnitude of the plastic shear strain on the α th system, s α and m α are the unit vector in the slip direction and unit vector normal to a slip plane, respectively, corresponding to the α system.The irrecoverable (plastic) part of the linearized Green strain tensor can then be represented as and the elastic part of the strain tensor is The work conjugate stress tensor for the linearized Green strain is the Cauchy stress, σ .As a result of the linearization of the strain tensor, crystal orientations are constant and do not rotate under strains.The critical resolved shear stress (CRSS), τ α c , on each slip system can be described by the relation where τ α 0 is the initial CRSS listed in Table II and h is the hardening matrix for the crystal.Due to limited experimental characterization, the hardening matrix, h, is taken to be 0, so hardening is not explicitly considered.The collection of functions given by defines what is known as the yield surface in classical plasticity theory.The rate of plastic deformation is constrained by Equation ( 4) such that The admissible stresses are defined by the Kuhn-Tucker conditions, which specify that the slip rate is non-negative and that plastic flow can only occur on the yield surface.The elastic predictor -plastic corrector method proposed by Cuitino and Ortiz 58 for solving the constitutive equations above is used.This method is similar to the projected Newton method when line search is avoided.The plastic work density is The heat equation accounting for both thermal generation and conduction is where θ is temperature, ρ is mass density, c v is specific heat, k is thermal conductivity, and η is the fraction of plastic work that is converted into heat.For the durations of the simulations carried out here (∼1 μs), the characteristic distances of thermal conduction is smaller than the size of a single finite element, 59 therefore the diffusive term in (12) does not have a significant role in the results reported here.

C. Statistical Variation of Critical Resolved Shear Stress (CRSS)
In addition to the effect of the number and structure of slip systems, analyses are also carried out to study the effect of the variation of the critical resolved shear stress (CRSS) on the impact response of HMX.The average CRSS for each slip system, τ α 0 , is shown in Table II.In the analyses, two scenarios are considered.The first involves the CRSSs being constant spatially and equal to the averages shown in Table II.The second involves randomly varying the CRSS spatially throughout the microstructure in an element-by-element manner.The variation of each CRSS is relative to the average in Table II and follows a normal distribution, as illustrated in Figure 3.The standard deviation of the variations is 15% of the mean for each slip system to approximate the stochastic spatial variation of local strength measures.Comparing the uniform CRSS case and the case with random statistical variations of the CRSS allows the effect of stochastic variations of the material property to be delineated.The numerical implementation of the CRSS variations is achieved by assigning randomly varying CRSS values for each finite element in the microstructure.

III. NUMERICAL TECHNIQUE
Calculations carried out concern ensembles of β-HMX crystal grain subject to loading under planar impact.The calculations are performed using a Lagrangian finite element framework with 4-noded, linear strain tetrahedral elements.As shown in Figure 2, the specimen analyzed here has overall dimensions of 800 × 400 × 400 μm.The loading is effected by imposing an initially ramped (ramp time ∼5 ns) velocity boundary condition in the longitudinal direction at the left face of the specimen.This loading condition involves the specification of the particle velocity history along the left face of the specimen.The profile imposed includes a steeply ramped region of about 5 ns followed by a constant velocity regime for the remainder of the calculation.The lateral faces of the specimen are prevented from moving transversely, therefore, the conditions in the specimen can be regarded as those of nominally uniaxial strain at the scale of the specimen size.The imposed boundary velocities range from 50 m/s to 400 m/s.Each specimen consists of 817 approximately equiaxed grains of either 100 μm or 64 μm in size (as labeled in Figure 2).In the single crystal baseline cases, all grains are crystallographically aligned in the same orientation.In the polycrystal cases, each grain has a randomly aligned crystalline orientation (creating a polycrystal with a random texture).Due to the small-strain assumption implicit in the constitutive equations, the grain orientations assigned initially do not change during the deformation.This assumption is reasonable since the magnitude of plastic strain observed is small.

IV. RESULTS AND DISCUSSIONS
In order to establish a baseline thermal response due to crystal plasticity, we have analyzed the response of single crystals of β-HMX to imposed boundary velocities in the range of 50 to 400 m/s.Each specimen is loaded in the {110}, {0 11}, or {010} crystal direction.For each orientation, the temperature rise is analyzed as a function of imposed velocity for both the Barton case and Zamiri case.The characteristic temperature rise for each case shown in Figure 4 is the average temperature rise behind the wave front in the three crystalline orientations.The solid line represents the average characteristic temperature of the three orientations.The error bars denote the highest and lowest characteristic temperatures among the three orientations.This figure shows that the characteristic temperature rise increases from 0.1 K at 50 m/s to 11.2 K at 400 m/s for the Barton case and from 0.14 to 4.0 K for the Zamiri case.Overall, the temperature increases are relatively low at several Kelvin, due to the homogenous nature of the samples.The low temperature increase in single crystal specimens echoes the idea that localizations or "hotspots" are necessary for high temperature rises in energetic materials.In the following sections, we will use these characteristic temperatures as baselines for comparison to evaluate the effects polycrystalline microstructure, numbers of slip systems, and spatial variations of the CRSS on the fluctuations in stress and temperature fields.

A. Effect of Number of Slip Systems
The effect of the number of slip systems on the behavior of polycrystalline HMX is analyzed.Specifically, the two different material cases as described in Section II B, i.e., the Barton case with seven slip systems and the Zamiri case with two potential slip systems are used in the analysis.The calculations involve polycrystalline ensembles and imposed velocities from 50 to 400 m/s.
Figure 5 shows the evolutions of the longitudinal stress (a-c) and temperature (d-f) as the stress wave traverses a specimen of the Barton material.The imposed velocity is 200 m/s.The stress is non-uniform in the material, varying between 510 MPa and 1.65 GPa, due to the heterogeneous material microstructure.Although the stress level fluctuates spatially, the average longitudinal stress along the direction of wave propogation is approximately constant behind the stress wave front.In Figure 5(d)-5(f), the temperature field is characterized by distinct regions of thermal fluctuations, reflecting the influences of both the non-uniform stress field and the anisotropic and heterogeneous material response.The localized nature of the stress and temperature fields is the focus of this study.
The mean, standard deviation, maximum, and minimum of the axial stress on planes parallel to the impact face are calculated at different distances from the impact face.Figure 6(a) shows the resulting values for the Barton case at an imposed velocity of 200 m/s.The black curve represents the mean stress.The blue and green curves are the mean stress ± one standard deviation and the mean stress ± two standard deviations, respectively.The red asterisks denote the minimum and maximum stresses.Figure 6(b) shows the stress variation in the material in the form of a probability distribution function (PDF) as a function of distance from the impact face.At any given distance, a value of 1 (red) in the PDF represents the statistically most probable stress at that position in the specimen.The current analysis focuses on the behavior within the quasi-steady portion of the stress wave behind the wave front.In this plateau region, the mean and standard deviation of the axial stress have constant values of 1.1 GPa and 108 MPa, respectively, for this microstructure and the loading conditions.Figure 7 shows the mean longitudinal stress behind the stress wave front.This average stress is nearly linearly related to the imposed velocity for both material descriptions.The figure also shows that the mean stress is consistently higher in the Zamiri case than in the Barton case.This difference increases with imposed velocity from 4.5% at 50 m/s to 16.2% at 400 m/s.Physically, this is due to the Zamiri case having limited avenues for energy dissipation via plastic slip, resulting in higher stresses due to higher fractions of elasticity in the material.The Barton case, however, has more avenues for accommodating plastic deformation through slip in more orientations, resulting in lower mean stresses.Figure 8 shows a comparison between the numerically calculated and experimentally measured 60 hydrostatic stress as a function of imposed piston velocity.The calculated value shown (red squares) is the average behind the stress wave front.The experimental data is for solvent-pressed 99.6% dense HMX. 60The calculated stress levels are in good agreement with the experimental data.
In order to put the scatter of the stress states in perspective and quantify how the dispersion of the data changes with loading and material, the coefficient of variation (standard deviation normalized by the mean) and the normalized maximum stress (maximum stress normalized by the FIG. 7. Plateau stress level behind the stress wave front as a function of imposed boundary velocity.The two slip system case (left) has a lower plateau stress level than the seven slip system case (right) at all velocities.Both increase with the impact velocity essentially linearly.FIG. 8. Comparison of calculated hydrostatic stress (red squares) for polycrystalline HMX (Barton case) and experimentally measured hydrostatic stress (blue diamonds) for solvent-pressed 99.6% dense HMX. 60an) are calculated.The coefficient of variation allows trends affecting a large percentage of the total volume of sample to be identified, and the normalized maximum provides insight into the extreme conditions in the domain.Figure 9 quantifies the dispersion of stress states in the materials by looking at the coefficient of variation of stress as a function of imposed velocity for both material cases.For all velocities, the Zamiri case displays greater variations in stress state than the Barton case.The coefficient of variation for the Zamiri case varies from a minimum of 0.118 at 50 m/s to a maximum of 0.145 at 400 m/s.The Barton case exhibited a maximum coefficient of variation of 0.113 at 100 m/s that decreases to a minimum of 0.064 at 400 m/s.This indicates a relationship between the heterogeneity in the material and the dispersion in the stress field: a higher level of heterogeneity corresponds to higher levels of expected variations in stress states.The lower bound of this phenomenon is simply the case of a homogeneous single-crystal under uniaxial loading.In such a limiting case, the stress state is uniform on cross-sections perpendicular to the loading axis.On the other hand, for polycrystals anisotropy and mismatch in orientations between the individual grains.In the Zamiri case, this effect is stronger.In the Barton case, however, the difference in response between different directions in each grain is less pronounced, leading to stress fluctuation levels that are lower than those in the Zamiri case.Another observation from Figure 9 is that the dispersion increases monotonically with imposed velocity in the Zamiri case but shows no simple trend in the AIP Advances 4, 097136 (2014) FIG. 9. Dispersion of the longitudinal stress field in the quasi-steady region behind the stress wave front.As heterogeneity in the material increases due to material anisotropy, wider ranges of stress states are seen.FIG. 10.Maximum longitudinal stress normalized by averaged stress in the plateau region behind the stress front as a function of impact velocity for both Barton and Zamiri materials.This plot demonstrates the spread of the extreme values observed in the materials.
Barton case.Specifically, the decrease in stress field variations at higher imposed velocities (above 100 m/s) in the Barton case is in agreement with the results by Trott et al. 61 for a mock PBX under shock loading.In contrast, the increasing levels of stress fluctuations in the Zamiri case suggest that this material representation may be missing essential dissipation mechanisms represented by the additional slip systems in the Barton case or fracture planes which are not modeled here.
Figure 10 shows the normalized maximum stress for both cases as a function of imposed velocity.The maximum stress is between 1.6 and 1.9 times the bulk average stress for the Zamiri case and between 1.3 and 1.6 times the average stress for the Barton cases.The Zamiri case has greater variability in stress states for all impact velocities, as well as a wider range of stress states as indicated by the higher normalized maximum stress.The dispersion of the stress field indicates that greater anisotropy of the underlying crystal structure renders a material less capable of homogenized response.This quantitative analysis of the stress field reveals the influence of the number of slip systems on the variability of material states.Next, we analyze the effect of anisotropy on the variability of the temperature field.
Figure 11 shows the average (red) and peak (blue) temperature as functions of distance from the impact face for the Barton material case at an imposed velocity of 200 m/s.The average temperature is maximum at the impact face and decreases steadily away from the impact face.The peak temperature profile shows significant oscillations and does not follow a monotonic trend, suggesting a degree of stochasticity in the distribution of heating in the material.This observation echoes those of Refs.25, 29-31, 62, and 63 An additional feature of Figure 11 is the plateau-like regions in the peak temperature profile.The are not exact flat regions of temperature.Rather, they have very similar temperatures.These regions correspond to grains (size of the regions on the order of grain size).Each grain has a specific grain orientation (dominant slip system) and the slips at different locations of a grain are similar, resulting in similar temperature levels.This is not to say temperature is uniform over any grain.It is only to say that it is possible that, for grains located and oriented in certain ways, temperature change may not be extremely non-uniform.
Figure 12 summarizes the temperatures for both slip system cases with constant and variable CRSS. Figure 12(a) shows that the average temperature increases monotonically for both cases as imposed velocity increases and the rate of temperature increase is slightly higher at higher imposed velocities.The Barton case exhibits higher average temperatures at all imposed boundary speeds, ranging from 293.2 K at a speed of 50 m/s to 308.0 K at a speed of 400 m/s.The average temperature for the Zamiri material increases from 293.0 K at to a 50 m/s to 300.6 K at 400 m/s.The higher rate of average heating in the Barton case is due to the enhanced availability of slip systems relative to the Zamiri case. Figure 12(b) indicates that the peak temperatures are not strongly dependent on the number of slip systems.The peak temperatures increase from 293.4 K at 50 m/s to 325.4 K at 400 m/s for the Zamiri case and 293.6 K to 320.0 K for the Barton material.This similarity is due to both material descriptions sharing the (001)[100] system, which has the highest CRSS.Therefore, the peak temperature rise in both material cases is associated with the strongest slip system which yields the highest temperature change for the same amount of slip.
In order to quantify the effect of the polycrystalline microstructure on the temperature rise in the specimens, we have compared the average and maximum temperature rises in the polycrystalline materials with the characteristic temperature rises for single crystals discussed earlier in Section IV.The average temperature rise in the polycrystalline specimens is 90.8% and 41.1% higher than that in the corresponding single crystal specimens for the Zamiri and Barton materials, respectively.The polycrystalline microstructure causes the maximum temperature rise to be 7.5 times higher relative to the single crystal in the Zamiri case and 1.65 times higher in the Barton case.From these results, it is clear that the resulting temperature field in the Zamiri material, which has only two available slip systems, is much more strongly affected by the polycrystalline environment than in the Barton material.For each material case, however, the effect of the polycrystalline microstructure on heating is fairly constant.Specifically, this is seen over the imposed velocity range of 100 to 400 m/s for Barton case and over 200 to 400 m/s for the Zamiri case.Data from this regime is used to characterize the effect of polycrystalline microstructure in this section and the effect of varying the CRSS in the next section.
As a measure of variability in the temperature field, the ratio between the peak temperature increase and the average temperature increase (normalized peak temperature rise) in the form of is used.Here, θ max (x) is the maximum temperature change and θ avg (x) is the average temperature change at distance x from the impact surface.The normalized peak temperature change is shown in Figure 13.Over the range of imposed velocity from 50 to 200 m/s, the normalized peak temperature decreases from 8.01 to 4.04 for the Zamiri case and from 4.49 to 2.11 for the Barton case.So, although the Barton case exhibits overall more significant bulk heating under these loading conditions, the Zamiri case exhibits a consistently larger deviation from the uniformly heated case.This result reinforces the idea that when extreme material states are of interest and dominate material response, as is the situation in studies on the sensitivity of energetic materials, isotropic and homogeneous material properties lead to more uniform thermal and mechanical fields such as stress and temperature.The results presented here clearly indicate the effect of microstructure is more pronounced as the degree of anisotropy in the underlying crystal structure increases.In order to quantify this effect, it is advantageous to develop a functional form for the normalized peak temperature rise, θ n , that captures the evolution of θ n as a function of the parameters studied: material anisotropy (A) and imposed velocity (v).A value of θ n = 1 corresponds to the case where the peak temperature rise is equal to the average temperature rise, indicating uniform heating along planes perpendicular to the loading direction.To choose an appropriate form for this relationship, two cases are considered: 1) in the limit of an isotropic material (defined asA = 0), θ n = 1; and 2) as the load intensity increases, microstructural details become less significant, i.e., as v → ∞, θ n → 1.These conditions motivate the use of a relationship in the form of where c and k are constants, v is imposed velocity, and A is a scalar measure of the anisotropy of the plastic behavior of the underlying crystal structure.In order to use this relationship, a measure of material anisotropy as a function of available slip systems for both the Barton and Zamiri materials must be adopted.We have chosen to consider the anisotropy as a function of the number of available slip systems (n ss ) and the coefficient of variation of the available CRSS c C RSS v .The value of c C RSS v is calculated by using the magnitudes of the available slip systems listed in Table II.For the Zamiri case, the coefficient of variation of its two available slip systems is 0.63, and the coefficient of variation for the seven available systems in the Barton case is 0.39.The form of the relationship for A is determined by the limiting case of an isotropic material for which A = 0.An isotropic material can be thought of as one that has an infinite number of slip systems which are uniformly distributed in all directions and have the same constant CRSS, consequently, we take A = 0 if and only if n ss = ∞ and c C RSS v = 0.One form that meets the considerations is For the analyses carried out, A Barton = 0.53 and A Zamiri = 1.13.For reference, the value for a bcc material having 48 slip systems with the same CRSS is 0.02.The normalized peak temperature rise in Figure 13 leads to c = 6.962 and k = 0.004.The dash lines in this figure show the predicted value of θ n for both materials over the entire range of imposed velocity considered.This relationship for θ n may be used to estimate the expected peak temperature in polycrystalline HMX.

B. Effect of Variation in CRSS
In this section, we will discuss the effect of spatially varying critical resolved shear stress (CRSS) on the material response.Figure 11 and Figure 13 compare the thermal behaviors of microstructures with constant CRSS and those microstructures with CRSS that changes from location to location in the material.The spatial variations of the CRSS follow the normal distribution, as shown in Figure 3.The two figures show that variation in CRSS has minimal effect on the average temperature, peak temperature, or the normalized peak temperature.Also, the average stress for the variable CRSS case differs from that for the constant CRSS case by less than 1.5% at all imposed velocities, with the coefficient of variation and normalized maximum stress both falling within 8.0% of the constant CRSS values as well.
To further quantify the contribution of the variation in CRSS, we have also compared the thermal response from the cases to the characteristic single crystalline cases in Figure 4.The average temperature rise is 94.4% higher than the characteristic temperature rise in the single crystals for the Zamiri material and 42.2% higher for the Barton material.These constitute a 3.8% and 2.8% change from the heating in polycrystalline microstructures with constant CRSS for the Zamiri and Barton cases, respectively.The peak temperature rise exhibits a similar trend, with the Barton material showing a peak temperature rise that is 2.65 times the rise in the single crystals and the Zamiri material showing a rise 8.4 times that in the single crystals.The contribution of variation in the CRSS to peak temperature is 1.9% for the Zamiri case and 0.3% for the Barton case.
The fact that variations in the CRSS has very little effect on the thermo-mechanical response is a strong indication that the material response in the loading regime considered is dominated by microstructural features other than the spatial variation of the CRSS.Factors include the distribution of crystalline orientations, grain size and size distribution.Because the strongest slip system (highest CRSS) has a CRSS that is ∼4.5 times that for the weakest system, varying any CRSS by as much as 50% has less influence on the mechanical response of a polycrystalline specimen than rotating grains to more/less favorable orientations for slip.In this paper, the results certainly show that the orientation of individual grains has a more significant effect on the response of the specimen than the variation of the CRSS.However, if the misorientation between the grains was limited, as in a textured specimen or a single crystal, the spatial variation of the CRSS may have a larger effect, this aspect is not analyzed here but may be considered in the future.

V. CONCLUDING REMARKS
The behavior of a regularly structured polycrystalline HMX is investigated by considering perfectly bonded ensembles of single crystallites which have distinct individual crystalline orientations.The response characterization has focused on the evolution and variability of stress and temperature fields arising from impact loading.The quantification has primarily concerned the inhomogeneity in stress and temperature.It is found that crystalline anisotropy has a large effect on both the homogenized (averaged) material response and the dispersion in stress and temperature states in each specimen.The analysis is based on an implementation of the crystal plasticity models due to Barton et al. 17 and Zamiri and De, 18 which contain seven and two potential slip systems, respectively.The results presented here indicate that the results based on the description of Barton et al. more closely matches previously published data. 61Therefore, the Barton model is likely a more realistic representation of the inelastic behavior of HMX crystals in the absence of fracture or other dissipative mechanisms.The calculations carried out also show that spatial variations of the critical resolved shear stresses (CRSS) of up to 15% in a material have minimal effect on the bulk or local material response for the polycrystalline specimens analyzed.
The two key results of this work are (1) the determination of the magnitude of heating due to crystal plasticity and (2) the quantification the effect of crystal anisotropy on the variability of stress and temperature states.With regard to the first point, a polycrystalline sample of HMX is shown to undergo average heating of 0.04 to 15.0 K and maximum local heating of 0.4 to 32.4 degrees K for imposed velocities up to 400 m/s.Even though the temperature rise due solely to crystal plasticity is significantly lower than that required to cause initiation, others have shown that plastic heating may be significantly exacerbated by the presence of other stress concentrators 17 or friction along crack surfaces 64 not explicitly modeled here.
Finally, the anisotropy inherent in polycrystalline microstructures has been shown to contribute significantly to the dispersion in the stress and temperature fields.This dispersion means that localized regions within microstructures experience significantly higher temperatures and stresses than those predicted by a homogenized model.These regions provide favorable locations for hotspots to form.In order to accurately account for the hotspot development, it is essential to model the fluctuations in the local stress state and temperature field due to, not only microstructural heterogeneities analyzed here, but also other processes such as friction, viscous dissipation, and pore collapse.66] 2158-3226/2014/4(9)/097136/17 C Author(s) 2014 4, 097136-1 All article content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 Unported license.See: http://creativecommons.org/licenses/by/3.0/Downloaded to IP: 130.207.153.228On: Thu, 25 Sep 2014 14:25:01 1)[ 10 1], (011)[0 11], ( 102)[010], (011)[100], and (010)[100].

FIG. 1 .
FIG. 1.The seven potential slip systems of HMX crystal in the P2 1 /c space.**The Zamiri case involves only the two slip systems labeled in green.The Barton case involves all seven potential slip systems shown here.

FIG. 3 .
FIG. 3. Statistical distribution of the spatial variations of the critical resolved shear stress (CRSS) considered.The variations follow a normal distribution about the reported mean values in TableII.

FIG. 4 .
FIG. 4. Characteristic temperature behind the stress wave front for single crystalline specimens of HMX with different orientations.The behaviors considered are described by the Barton model (blue) and the Zamiri model (red).Each datapoint represents the average temperature in a single crystal specimen under impact along the {110}, {011}, and {010} directions.

FIG. 5 .
FIG. 5. Distributions of longitudinal stress (a-c) and temperature (d-f) at different times in a specimen with seven slip systems.The domain has been cut along a longitudinal midplane in order to show the stress and temperature states on the interior of the specimen.The imposed boundary velocity is 200 m/s.

FIG. 11 .
FIG. 11.Peak and average temperature profiles for the Barton material at an imposed boundary velocity of 200 m/s.The average is over cross-sections normal to the loading direction.

FIG. 12 .
FIG. 12. (a) Average temperature in the plateau region and (b) Peak specimen temperature as functions of impact velocity and CRSS.

TABLE I .
Elastic Properties of B-HMX.

TABLE II .
17ip systems and critical resolved shear stresses for β-HMX.17Allseven slip systems are active in the Barton case, but only the two systems labeled with an asterisk (#1 and 7) are used in the Zamiri case.