Experimental study of anisotropic constitutive behavior of β-HMX crystals via nanoindentation and small-scale dynamic impact

For energetic crystals such as HMX, the sensitivity of the material to shock, the possibility of initiation, and the subsequent reaction is known to be controlled by processes occurring at the crystal level. The anisotropic nature of β-HMX can be critical in determining the performance of HMX based polymer bonded explosives, which are widely used across multiple industries as propellant or explosives. In this work, we experimentally obtain constitutive parameters for characterizing the response of multiple crystalline planes of β-HMX crystals to external loading. Nanoindentation and small-scale dynamic impact experiments were performed on multiple planes of β-HMX crystals to comparatively measure the indentation moduli in multiple orientation directions. Anisotropic material behavior, involving constitutive elastic and non-elastic parameters, was measured and studied. Findings regarding material properties for the (100), (010), (001), {110}, and {011} planes of β-HMX are presented and compared with literature data.


Introduction
Accurate characterization and prediction of the mechanical behavior of energetic materials (EMs) requires indepth understanding of their physical properties and of their response to external stimuli. There are many material responses whereby it is important to understand non-initiating behaviors, including deformation and failure, as damage may alter the safety parameters of these explosive material and set a new initial condition for future performance [1]. For an energetic material such as HMX, sensitivity to shock, the possibility of initiation, and the subsequent reaction is known to be controlled by processes occurring at the crystal level. HMX occurs in four polymorphic forms, α, β, δ and ε. Of these, the most thermodynamically stable phase is β-HMX [2]. Thermal stability of the β-HMX crystalline phase makes it a good candidate to characterize effect of anisotropy on HMX material behavior. Previous experiments have pointed to the anisotropic nature of β-HMX being critical in determining the performance of HMX based polymer bonded explosives (PBX) [3]. The main failure mechanisms in such materials have been identified as particle fracture, interfacial failure between particle and binder, and interfacial failure between particles [4]. With respect to particle fracture, due to the anisotropic nature of HMX, it is known that deformation in single crystals of β-HMX is orientation dependent and fracture tends to occur along cleavage planes. For instance, indentation experiments performed by Palmer and Field [5] on cleaved crystals determined the cleavage in β-HMX to take place along the {011} crystallographic planes.
The ability to understand and, in turn, accurately predict the behavior of EM critically relies, in part, on understanding the elasticity of the material, as the structure and compressibility of the material are relevant regardless of initiation shock mechanics [1]. Description of the anisotropic elastic behavior of β-HMX has been widely investigated. However, a significant portion of the description has been obtained via molecular dynamics related numerical methods [6]. Experimentally, it has proven difficult to measure elastic properties of energetic crystals such as β-HMX due to their anisotropic nature requiring complex geometric considerations and their brittle nature yielding relatively minute strains before failure [7]. Previous experimental methods involving elastic constant determination for β-HMX have involved the use of techniques such as impulsive stimulated light scattering (ISLS) [8], impulsive stimulated thermal scattering (ISTS) [9] and Brillouin scattering [10]. However, these methods, when compared, produce a variance in values that is yet to be fully understood or explained [1]. Nanoindentation as an experimental procedure has also been used in attempts to characterize the hardness and modulus of β-HMX in previous studies [5,[11][12][13]. The results from such studies provide unexplained variance which may only be fully understood with more work in the area.
With respect to EMs, hot spots are regions in the microstructure whereas the initiated volume is sufficiently large, and temperature is sufficiently high to start the chemical reactions leading to detonation. Friction, or interfacial interaction between particles as mentioned above, among other mechanisms, is known to be a cause of hot spot initiation. Susceptibility of an energetic material to such form of failure, most adequately described by the frictional sensitivity of the material, has been shown to be dependent on multiple factors, of which hardness has been speculated to be the most important [14]. In addition to elastic properties, this work also presents a study of the anisotropic hardness values and viscoplasticity parameters of β-HMX.
Currently, with few exceptions, simulation of the mechanical behavior of HMX-based EMs at the mesoscale generally overlooks the anisotropic properties. However, an argument can be made that incorporation of anisotropic properties of HMX crystals may improve the current ability to predict the energetic behavior of such materials. Examples of such exceptions can be seen in work done by Wang et al. [15], in which a full field method was employed to model the local damage behavior of a HMX-based EM. Other works simulate the anisotropic material behavior HMX to determine anisotropic material responses. Work by Conroy et al. [16] modelled the behavior of β-HMX crystals using first principles density functional theory (DFT). They found that the anisotropic material model predicted experimental results that adequately matched existing data. Zamiri et al. [17] also worked on numerically predicting the anisotropic plastic response of modelled single crystals of HMX based on experimental data. This work focuses on experimentally characterizing anisotropic quasistatic as well as dynamic local constitutive material properties of β-HMX crystals using indentation based static and dynamic schemes.

Methods
This section covers brief explanation of the experimental methods and theories to be used in this work. These methods have been covered in detail in previous works [18][19][20][21].

Sample preparation
Untwinned single crystals of β-HMX were formed by controlled solvent evaporation, at room temperature, from a saturated solution of raw HMX in HPLC grade 99.8% pure acetone. The solution was mixed at ratio of 0.024:1 by weight of HMX to acetone respectively and measured into vials with perforated caps for slow evaporation. The descriptions here are of monoclinic single crystal β-HMX in space group P2 1 /n. with lattice parameters a = 0.654 nm, b = 0.736 nm, c = 1.104 nm, and B = 102.66°. Two main morphologies of single crystal were observed, and the respective orientations were identified through x-ray diffraction shown in Fig. 1 below.
The first morphology observed, denoted 'morphology 1' , exposed facets corresponding to the {010}, {011}, and {110} crystallographic planes, while the second morphology displayed exposed facets corresponding to {011} and {110} planes only. In both morphologies, the facets corresponding to the {011} planes appeared to be of preferred growth, displaying larger surface areas in most crystals observed. The morphologies observed also grew longitudinally along the 'a' crystallographic axis as shown in Fig. 2. The length distribution (along the preferred a-crystallographic axis) of fully formed crystals varied from 200 μm to 4 mm with an average length of 2 mm and width (approximately perpendicular to the a-axis) of 1 mm. Typically crystals equal to or larger than 1 mm were chosen for experiment purposes to provide sufficient surface area for multiple indentations.
Experiments for further characterization of elastic behavior required access to indentation axes perpendicular to planes that were not exposed in the single crystal morphologies that presented, such as the (100) and (001) planes corresponding to the a and c crystallographic axes. In such cases, crystals were embedded in a hardened epoxy and cleaved to expose the necessary planes and sanded down with decreasing particle sizes. In particular, the hardened epoxy used here was the 'Aluminite Amazing Casting Resin' due to its lower volumetric contraction on curing. The surfaces were then polished using silicon carbide sandpaper with a < 8 μm average particle diameter.

Methods to obtain anisotropic elastic constitutive properties
This work presents experimentally obtained indentation material properties to describe the anisotropic elastic behavior of β-HMX as a response to external loading applied to different planes of the crystal. This section describes the experimental and theoretical methods employed.

Nanoindentation
Nanoindentation has been used in multiple studies of material elastic behavior and has proven to be a reliable technique to measure local elastic material properties [21][22][23][24], extending to elastic properties of Ems [5]. Properties such as elastic modulus and hardness can be easily obtained via nanoindentation due to the underlying fact that displacements recovered during unloading are known to be mainly elastic. In this work, a NanoTest platform (Micro Materials Ltd., UK) [25] is used to perform the indentation experiments. Multiple indentations were performed, at various indentation depths, on each plane at room temperature with a Berkovich type indenter having a diamond tip radius of 20 nm. Figure 3a shows the indentation test setup. Sample crystals were oriented with the loading plane of interest perpendicular to the indenter and indentations were spaced with a minimum of 100 μm between each indent to limit the possibility of coalescence between indents. The machine compliance is corrected as detailed in previous work [21] (Fig. 3b).

Hardness and elastic modulus evaluation
An indentation experiment is characterized by loading applied by the indenter until the desired maximum load is achieved, whereby the maximum load is held for 5 s and then unloaded. Each indentation experiment produces a load-depth plot, an example of which can be observed in Fig. 4. The Oliver-Pharr's method [26] is used to calculate the indentation modulus and hardness. In the case of isotropic materials, the Oliver-Pharr method has been shown to estimate material properties such as Young's modulus to a high accuracy, about 4% of literature values. However, this is not the case for anisotropic materials as the indentation modulus represents more of a weighted average value as deformation involves multiple directions [27]. Therefore, in this work, the indentation moduli for each direction will be obtained individually, as weighted averages, and compared to actual elastic constants for those directions. To obtain the necessary parameters, the values obtained from the indentation experiment are fit the load-depth curve obtained from the experiment.
where A m and m are parameters obtained empirically from the power law equation, h is the indentation depth corresponding to a specific load P, and h f is the plastic indentation depth remaining after the load has been completely removed. The indentation hardness can be calculated as,  For a Berkovich indenter, the contact area can be estimated as 24.5h 2 c , where h c is the contact depth. The reduced indentation modulus (M r ) is calculated as, The indentation modulus of the material is then obtained from the relation, Constants υ and M are the Poisson's ratio and indentation modulus of the sample, while υ i and E i represents Poisson's ratio and Young's modulus of the indenter, respectively. The Poisson ratios used in this work were obtained from values presented in literature [6] and approximated thus, Therefore, 0.33, 0.3, and 0.4 were used for ν 1 , ν 2 , ν 3 and respectively.

Indentation size effect (ISE)
Previous studies of hardness obtained by nanoindentation have shown a dependence of the measured values of hardness on 'indentation size' , i.e. the indirect variation of hardness with indentation depth. This phenomenon is referred to as the indentation size effect (ISE) and some possible explanations of the causes of this behavior, particularly in crystalline materials, have been proposed in literature [28], for instance the strain gradient plasticity theory [21]. The hardness eventually reaches a constant value which can be quantified according to the Nix-Gao relation [28] as: Here, H is the hardness value corresponding to a certain indentation depth h. H 0 is the final constant hardness value at infinite depth, and h* is a characteristic length parameter, both of which are important parameters necessary for the definition of constitutive equations in strain gradient plasticity theory [21]. Equation 6 can be rewritten as,

Methods to obtain anisotropic strain-rate-dependent viscoplastic constitutive model
In order to describe the anisotropic viscoplastic behavior of β-HMX crystals, a viscoplastic constitutive power law is defined, and experimental orientation-dependent viscoplastic parameters are obtained in this work. In a previous work [29], the behavior of Hydroxy-Terminated Polybutadiene (HTPB)-HMX based PBXs have been modelled using a cohesive finite element method (CFEM) based simulation software. The developed CFEM model used polycrystalline HMX parameters in order to obtain an approximated description of the viscoplastic behavior of the material.
Here, we present a study into the orientation dependence of these parameters, giving a more detailed description of the range of possible behavioral limits and how well previous approximations fit into this range.

Dynamic small-scale impact experiments
The viscoplastic constitutive parameters are obtained via nanoscale dynamic impact in regions of interest on selected β-HMX crystal phases. The specific equipment used for this experiment is impact module of Micro Materials, UK. This equipment and the corresponding technique have been used to successfully obtain mechanical properties of other EMs at small scale such as HTPB-AP interface [30]. In this experimental setup, the sample is positioned in front of a spherical impactor with the plane face of interest oriented perpendicular to the impactor. The sample is mounted on a three-dimensional stage and the entire setup sits on a "floating" platform; thus, the sample can be moved freely in all directions. The impactor is a conical indenter of radius 1 μm, it is mounted on a vertical pendulum which hangs on a frictionless spring. Impact on the sample is achieved through the use of electromagnets on the vertical pendulum and as such different magnitudes of desired impact force can be achieved. The schematic for the equipment setup is shown in Fig. 5 , [19], below. Depth is measured by the distance between the capacitors located behind the impactor and the velocity is obtained by a continuous differentiation of the depthtime plot. Impact is performed on each plane of interest on multiple samples. The stress, strain, and average strain rate of the experiment can be obtained using values from the depth-time and velocity-time plot. Figure 5b shows an example of a plot output from one impact experiment. The impact strain rate is calculated as, where h max is the maximum depth and V in is the maximum velocity. The strain and stress are also calculated as, where h res is the residual depth and P is the impact load. These values are fit to a power law model which is assumed to represent the effective stress-effective viscoplastic strain relationship.

Viscoplastic constitutive parameter evaluation
Following the method outlined by Prakash et al. [18], we obtain the orientation dependent viscoplastic model parameters for each plane by fitting the load related values obtained above to a power law model:

Here,
vp is the effective viscoplastic strain, and ̇v p is the effective viscoplastic strain rate. The coefficient, A, is plotted against the effective viscoplastic strain rate on a log-log scale to obtain F 0 and m parameters more simply as the slope and intercept of the plot. The parameters fit to a power law model useful for predicting the constitutive behavior of the material at higher strain rates [30]. Here, the parameter 'm' represents the strain rate exponent and the parameter 'n' represents the strain hardening exponent, signifying the quantifiable influence of the effective viscoplastic strain rate and strain respectively on the measured effective stress. The validity of the viscoplastic equivalent stress-strain power law for modeling the behavior of this material has been evaluated in previous work [18]. Here, we obtain the parameters to define the anisotropic viscoplastic behavior in each different orientations of interest in the examined single crystal.

Results and discussion
Using the process detailed above, at different loads, multiple indentations were performed on each plane of interest on the HMX crystals. Grown crystal samples presented smooth surfaces parallel to {010}, {011}, and {110} planes and these were used in experiments without modification. Indentations were also performed on planes perpendicular to the three crystallographic axes to characterize behavior on these planes. Two out of three of these planes, (100) and (001) perpendicular to the monoclinic a and c axes, respectively, did not present parallel facets naturally in the samples grown. As such, modifications, detailed above (Fig. 2), were made to expose them. The plane perpendicular to the b crystallographic axis, the (010) plane, naturally presented on examined crystal samples as shown in Fig. 1 earlier  (Fig. 6).
In order to define the indentation moduli for the three main axial directions ({010}, {011}, and {110}), an additional plane was required to be exposed for loading. Indentation moduli corresponding to the elastic moduli, E 22 and E 33 , were obtained directly as the b and c crystallographic axes were intentionally oriented along the 2 and 3 axes as shown in Fig. 2a. The indentation modulus corresponding to the directional elastic modulus, E 11 , however, were obtained by loading along a direction determined by the cross product between the a and c crystallographic axes. The results of the moduli obtained by indentation for these three directions are presented in Table 1 (Fig. 7).
Our results for modulus in the E 22 direction showed good agreement with the value presented by Kucheyev et al. [12], in which Berkovich indentations were performed. The work of Li et al. [11] also presents moduli for different facets of β-HMX crystals as obtained by Berkovich indentation. However, our values are significantly lower than their presented results. As discussed above, this is not unique to nanoindentation as other means of experimental determination of elastic moduli for this material also present discrepancies. It should be noted however, that values obtained present a similarity to values obtained by molecular dynamics simulations [31] with a maximum of 12% error, and otherwise reasonably fit into the statistically determined range of values presented, obtained from both numerical and experimental methods. Moreover,   the values obtained in this work agree with most works which show E 22 to be the largest and E 33 to be the smallest moduli. The anisotropic behavior of three different orientations with respect to indentation modulus and hardness can be observed in Fig. 8. Here, the properties of planes corresponding to the three main crystallographic axes are compared and anisotropic behavior can clearly be observed in the values of measured properties with changing indentation load/depth. The hardness and plastic depth values were fit using the Nix-Gao relation detailed above, the results are shown in Fig. 9. The indentation moduli, hardness, and Nix-Gao relation fitted values obtained for each plane are presented in Table 2. Li et al. [11] presented hardness values for the (010) face of β-HMX as 1.13 GPa, also obtained via Berkovich indentations on the crystal facet. This is in general agreement with the value of 1 GPa, for arbitrarily oriented HMX crystals, presented by Burch et al. [13], however our value of 0.55 GPa for the same facet (010) are significantly lower and closer to the value of 0.65 GPa presented by Kucheyev et al. [12].
A significant amount of the data on hardness for β-HMX presented in literature focuses on the Vickers Hardness Number (VHN). Studies have shown the hardness values from Vickers indentations to be similar to, however slightly lower than, values obtained for Berkovich indentations [32,33]. The hardness values obtained here are higher but close in value to the VHNs; 0.41 GPa obtained by Amuzu et al. [34] and 0.4 GPa obtained by Palmer and Field [5]. However, both works did not indicate specific orientations for the values. A VHN range of 0.3-0.36 GPa, presented by Gallagher et al. [2] for indentations made on the (010), {011}, and {110} planes, fall lower comparatively.
Viscoplastic parameters were also obtained in this work for the viscoplastic effective stress-strain power law described above. These parameters were obtained via  small-scale dynamic impact experiments, to describe the viscoplastic behavior. Viscoplastic parameters for polycrystalline HMX embedded in a polymeric binder have been obtained in a previous work [35]. As the material was polycrystalline, due to the underlying assumption of isotropy, the parameters obtained were independent of orientation. The parameters for Ammonium Perchlorate (AP) have also been obtained in a previous work and both parameters are included below for comparison. The orientation dependent viscoplastic parameters for the viscoplastic effective stress-strain power law are presented in Table 3. Figure 10 shows the fitting of the data obtained for the (100) plane to the viscoplastic power law at multiple strain rates. From the results presented above, we observed that viscoplastic strain for the {110} plane appears to have the highest comparative dependency on strain rate, as determined by the stain rate exponent 'm' . The strain values are predicted to be the highest in comparison to other planes at a given strain rate and stress. This behavior is as expected as this plane has a high value of the strain hardening exponent 'n' as well. These findings appear to agree with results obtained by Zamiri et al. [17] in the anisotropic simulation of β-HMX, where simulation, and compared experimental data [7], showed the {110} plane to be relatively more prone to plastic deformation. In addition, the results potentially corroborate work by Conroy et al. [16] which found maximum principal stress to be greatest for the {110} and {011} plane compressions. The (100) plane appears to have the lowest dependency on strain rate and resultingly low comparative viscoplastic strains. Viscoplastic behavior of AP will be closer to the behavior predicted for most planes in terms of strain rate dependency, as the strain rate exponent value is closer to the average. However, in comparison to most planes, the values for polycrystalline HMX predict low viscoplastic strains closer to the lowest values of the (100) plane, with similar exponent values. Figure 11 shows the comparative values of the exponent parameters.
An interesting finding in this work is that the planes with the most and least comparative viscoplastic strains also had the highest and lowest hardness, respectively. This result is considered to be as expected as materials with higher hardness values will in turn have a higher values of yield strength.

Conclusions
In this work, experimentally obtained constitutive parameters for modeling the anisotropy dependent mechanical behavior of β-HMX were presented. Indentation modulus and hardness values were obtained for each plane via nanoindentation experiments. Results obtained fit well into a statistically determined range of values obtained from both numerical and experimental methods presented in literature. Nanoindentation has been shown to estimate the elastic modulus of isotropic materials to a 4% accuracy and as such, these results and future validations may present a discussion on the possibility of nanoindentation providing a means of obtaining acceptable approximations of anisotropic elastic moduli. Viscoplastic constitutive parameters were also obtained for each plane by fitting the values obtained from dynamic impact experiments to a power law model. The results showed the {110} plane to have the highest comparative viscoplastic strains respectively at a given strain rate and stress, which is supported by other findings in literature. The parameters presented in this work could provide a more accurate description of the anisotropic behavior of β-HMX in the prediction of the performance of HMX based EMs.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.