Response of explosive HMX to low-velocity impact: modeling by the crystal plasticity finite element method

Crystal plasticity finite element method (CPFEM) is a powerful tool for modeling the various deformation problems, which takes into account the different plasticity mechanisms at microscale of grain sizes and contribution of anisotropic behavior of each grain to macroscopic deformation pattern. Using this method we simulated deformation and plasticity of high explosive HMX produced by relatively low velocity impact. It was found that such plastic deformations of grains cause local heating which is sufficient to induce chemical reactions.


Introduction
Safety of different devices, which contain high explosives (HE) like octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine (HMX), is determined in many aspects by chemical reactions initiated by impact on HE with different velocities. The impact loadings may lead to high pressure (more than 10 GPa) or to low pressure (less then 1 GPa). Many phenomenological models were developed for high-pressure impacts [1,2]. Those models are widely used for HE reactions predictions in engineering practice. But the commonly accepted models for low-pressure impacts are lacking. This fact is connected with complex and multiscale processes involved to studied phenomena. Complexity of determination of the leading process stems from many reasons for generation of detonation centers such as HE heterogeneity, HE deformation process instability, and friction of material grains.
Some researchers consider adiabatic shear bands formation during shear loading of HE as the main mechanism of the non-homogeneous high-temperature heating which leads to detonation [3,4]. There are several phenomenological models for describing this process [3][4][5]. The main weakness of those models is the usage of a number of artificial parameters which are not based on clear physical picture. For clarifying physical processes we have to obtain the non-homogeneous pattern of heating HE explicitly. Several authors made attempts to receive the non-homogeneous pattern inside the separate grains of material with using a model for metals only [6]. Here we present the CPFEM simulations of impact loading of HE with taking into account anisotropic orientation of HE grains and their mutual interaction. As far as we know this approach is used for the first time.

Mathematical model of HE deformation
The elasto-plastic deformation of material grains conglomerate is known to be sufficiently anisotropic, and thus its behavior depends on loading direction. The anisotropic mechanical properties of the material is defined by its bulk anisotropic elastic characteristics and various plasticity deformation mechanisms (dislocations, twins and martensite transformations) for different spatial orientations of material grains. Complicated picture of elasto-plastic deformation is formed on account of anisotropic properties of material grains and its interaction. For accurate calculation of mutual interaction between HE grains, HE grains and plasticizer, and to account the various plastic responses of individual grains (in particular dislocation dynamics inside) in modeling we have to choose a method which able to account those deformation processes. Three decades ago the finite element method of crystal plasticity [7] was proposed with the aim to model deformation processes in metals and alloys, where anisotropic properties of material grains play a significant role in material response.
Here we apply this method to consider deformation process in HE. It allows us to account such processes as mutual interactions between grains as well as interactions between HE grains and plasticizer. The method makes feasible the calculation of plastic response of each grain including dynamics of dislocations and other processes involved in plasticity. The last processes are twinning, non crystallographic adiabatic shear, interaction of grain boundary with dislocations, twin-dislocation interactions, interaction between dislocations of opposite signs and etc. All of them can be implemented in finite-element model of crystal plasticity using a set of kinetic equations. Lets consider an implementation of dislocation plasticity, adiabatic shear and twinning, proposed in [8].
In this method a decomposition of deformation gradient F to factors is given as: where F e is an elastic part and F p is a plastic part. Plastic deformation varies as:Ḟ In the case of plasticity produced by dislocation motion only: where m α and n α are unit vectors, which describe a dislocation slip direction and a normal to slip plane for slip system α, respectively; N is a number of slip systems, andγ α -a shear rate for this system. Equation (3) is expanded by additional terms: which account for twinning and formation of non-crystallographic shear bands. For describing crystal microstructure, it is necessary to define expressions for shear rates of different plasticity mechanisms. For this purpose the small scale calculations can be used. For example in [8], to account the non-crystallographic deformation due to adiabatic shear banding the kinetic formulae discussed belowere used. At first the next parameters introduced for microstructure description are ρ sgl , ρ dip -dislocation density and dislocation dipole density, respectively, and f -twins volumetric density. Shear rate of shear system α with taking into account the where τ α -shear stress,τ α -resistance of "forest" dislocation, b -the Burgers vector, v 0dislocation velocity of current slip system at shear stress equal toτ α , Q 0 -activation energy of dislocation slip, k B and T -the Boltzmann constant and temperature, p and q -numerical parameters, which describe an obstacle profile for dislocation motion. Slip resistance stressτ α depends on local dislocation density as: where τ solute -a resistance constant, G -shear modulus, ξ αα -parameter characterizing interaction stress between all possible dislocation types (mobile, immobile, coplanar, collinear,orthogonal etc.). Twin system volumetric share f β change is based on a twin nucleation theory [8]. It is proposed that twin nucleation occurs at three partial dislocations bowing between three immobile points, where the distance between points is L 0 . Twin formation critical stress is:τ where b twin -moving partial dislocations Burgers vector and γ sḟ where τ β -reduced shear stress in twin system, N 0 and r -numerical parameters. We assume that any twin grows until contact with an obstacle such as grain boundary or other twin boundary. Volume of new twin: here s -constant twin thickness, λ β -effective distance between obstacles for twin growth: d grain -constant grain size, d twin -twin size changes with volumetric twin fraction as: where f -a total volumetric twin fraction. Twin-twin interaction parameter ξ ββ equals zero for coplanar twin systems β and β and unity for non coplanar. Thus, result deformation rate for each twin system calculated as:γ where γ twin -characteristic twin stress.  A term describing the adiabatic shear bands is added because of non-crystallographic defects formation in alloys with low stacking fault, where shear dislocation motion is almost completely blocked by dense twin boundaries. For this processes shear rate is: whereγ χ 0 -characteristic shear rate, τ χ -reduced stress at system χ andτ sb -a threshold stress for shear band formation.
The method described above is applied for HE deformation problem, in particular, for HMX. Such application is difficult due uncertainty of HE anisotropic elastic characteristicsthey may differ in several times for HMX [9]. Dislocation motion parameters are completely