Numerical Simulation of Crack Initiation and Growth in PBX High Explosive Subject to Compression

PBX solid high explosive exhibits brittle behaviour in uniaxial tension, quasi-brittle in uniaxial compression, and ductile when subject to high confining pressure. Tension cracking is the primary failure mode of PBX quasi-brittle solid, which is the main effect leading to overall failure of structural integrity. One characteristics of brittle or quasi-brittle solids, such as PBXs, is that when subject to overall compressive loading, the tensile cracks can still initiate inside the material due to existence of imperfection within the materials. In present study the extended finite element method is applied to analyze the cracking failure mechanism in the PBX 9502 platelike specimen with cavity subjected to overall compression. The nonlinear constitutive behaviours and failure of PBX under complex stress states were described by means of stress state dependent strength surface, non-associated flow rule and cohesive model the nonlinear behaviors of PBXs, including failure. Analysis indicates the tensile stress around the cavity arises in the specimen although loaded by overall compression, and this local tensile condition leads to cracking initiation. The comparison between simulation results and the experimental data published by LANL[Liu C, Thompson D G. Crack initiation and growth in PBX 9502 high explosive subject to compression. Journal of Applied Mechanics, 2014, 81(10):212-213] shows that they are in agreement with each other on some aspects of crack behaviours, including overall development of crack history and inflexion, crack initiation moment, crack initial speed, etc.


Introduction
Polymer-bonded explosives (PBXs), highly particle filled composite materials comprised of 90-95% by weight of powerful explosive crystals held together by a polymeric binder (5-10% by weight), have been used in a wide variety of applications ranging from rocket propellants to the main explosive charge in conventional munitions for civil and defence fields [1][2][3][4]. In recent decades, the safety and reliability of structures under low-speed impact have been highly valued by governments and engineering fields, and so the failure of PBX explosives has attached great attention by academic circles [5][6][7][8][9][10][11][12].
It is generally believed that the ignition mechanism of explosives under low velocity impact includes: elastic-plastic or visco-elastoplastic deformation, damage and failure, crack initiation and growth, plastic localization and plastic work conversion into heat and heat conduction [13]. It is found from the experimental observation that the plastic deformation and subsequent cracking problems of the energetic materials directly affect the evolution of the reaction (such as increasing the surface area will accelerate the combustion rate) and the reaction level [14]. The non-homogeneous damage within the material has great influence on the mechanical properties and sensitivities of explosives, and crack propagation will directly affect structural integrity [15][16]. Therefore, it is necessary to understand the response of this material to mechanical stimulation, and then to provide theoretical supports for assessing the risk of crack occurrence and developing crack suppression technology.
The mechanical tests of material properties show that PBX exhibits different mechanical behaviours under different stress states [17][18][19], which brings great difficulties to numerical modelling [20][21][22]. In present study, the stress state dependent strength model and nonassociated flow rule are applied to describe the nonlinear constitutive behaviour of the material under complex stress states. Meanwhile, PBX exhibits brittle/quasibrittle characteristics in the failure mode. The local tensile microcracks initiates near the voids within material, even if subjected to compression loading as a whole, and the microcracks evolves and grows as the load develops, leading to overall failure [25][26][27]. In the simulation of PBX brittle tensile crack behaviour, the usual methods are: numerical simulation method by continuous damage mechanics [28][29], direct numerical simulation method [30], and multiscale [31] or mesoscopic numerical method [32], and XFEM (eXtended Finite Element Method) [33][34][35][36][37]. In this study, the cohesive based XFEM method was used to

Fundamentals of XFEM
XFEM is a new finite element method to solve the problem of fracture mechanics. It was first introduced by Belytschko and Black [38]. XFEM is an extension of the conventional finite element method based on the concept of partition of unity [39]. XFEM solved the problem of crack propagation in the finite element by using the idea of mesh separation independent. While retaining all the advantages of the traditional FEA, it is not necessary to mesh the existing cracks inside the structure. By far, XFEM is the most effective numerical method for solving the problem of discontinuities. It studies the problem within the framework of the standard finite element method, retaining all the advantages of the finite element method. This makes XFEM one of the main methods for simulating crack propagation. It represents the major advances in the numerical methods of computational mechanics field over the past decade. XFEM does not require aligning the meshes when dealing with cracks, and mesh refinement is not needed during crack growth, and high-precision numerical solutions can be obtained even on coarse meshes. These advantages of XFEM make this method very active in the field of computational mechanics for more than a decade [33,35].
In the analysis of fracture problems, XFEM introduces enrichment functions, including: asymptotic functions near crack tip to characterize the stress singularity near crack tip; and discontinuity function to characterize displacement jump across crack surface. In this study, cohesive model [40] and phantom node method [41][42] are applied in the calculation of crack opening and propagation. In the framework of XFEM method, the traction-separation cohesive behaviour of material is introduced to simulate crack initiation and propagation. That is, XFEM-based cohesive method simulates the propagation of cracks within the material along any path, and mesh must not be updated to match the geometry of the discontinuity as crack progresses.

XFEM-based cohesive behaviour
The cohesive model belongs to the damage model. It was first introduced by Barenblatt [40] and used the tractionseparation law (TSL) to simulate the decohesion of the atomic lattice, so as to avoid the singularity of the crack tip [43]. The cohesive model, combined with the finite element method, was first used in simulations of concrete, and was later introduced into simulations of metals, ceramics, composites and more. Cohesive interfaces or elements shall obey to the TSL, including viscoplasticity, viscoelasticity, rupture, fiber breakage, kinetic failure, and cyclic loading failure. The typical separation laws are Needleman's law [44], Hillerborg's law [45] and Bažant's law [46].
In present study, the crack propagation analysis is performed based on the XFEM-based cohesive segment method, and the linear elastic TSL model, the damage initiation criterion and the damage evolution law. In the elastic TSL model, the initial linear elastic behaviour is assumed, followed by the onset and evolution of damage. The linear elastic behaviour is expressed as the linear relationship between the cohesive stress and the open displacement of the cracking element, namely: where t is the nominal traction stress vector, and t n , t s , t t represent the normal and the two shear tractions, respectively, and the corresponding separations are denoted by δ n , δ s , δ t . The normal and tangential stiffness components will not be coupled: pure normal separation by itself does not give rise to cohesive forces in the shear directions, and pure shear slip with zero normal separation does not give rise to any cohesive forces in the normal direction. K is stiffness matrix, whose terms are calculated based on the elastic properties for an enriched element. The maximum principal stress criterion is applied to model the cracking initiation, which can be represented as That is, damage is assumed to initiate when the maximum principal stress ratio reaches a value of one. Here, 0 max σ represents the maximum allowable principal stress, in MPa. The symbol 〈〉 represents the Macaulay bracket with the usual interpretation.
A scalar damage variable, D∈[0,1], represents the averaged overall damage at the intersection between the crack surfaces and the edges of cracked elements. The variable monotonically evolves from 0 to 1 upon further loading after the initiation of damage. The damage evolution equation is described as energy dissipation according to Here, eqC σ is equivalent fracture energy release rate, in J.m -2 , which is described in Benzeggagh and Kenane [47] by the following formula: When material undergoes damage, the normal and shear stress components are affected by damage according to The crack will not occur when normal compression state is applied, i.e. n n Under a combination of normal and shear separations across the interface, an effective separation is defined as following, so as to describe the comprehensive evolution of damage 2 2 2 m n s t δ δ δ δ = 〈 〉 + +

Description of plastic deformation of PBX under complex stress states
According to the deformation characteristics of TATBbased explosives, such as the tension-compression asymmetry at quasi-static experiments at 50 °C [19], as shown in Fig. 1, the pressure-dependent yield surface, the non-associated flow law and the stress state dependent weighting function are utilized to describe the mechanical behaviour of the material under complex stress conditions. The tensile and compressive test data are introduced into the model.

Yield surface
The PBX is described in the elastic deformation model using a linear elastic constitutive model, and the same modulus in tension and compression is assumed. In the plastic deformation stage, the confining pressure has an influence on its mechanical properties since PBX material is an internal friction material. Therefore, the first invariant of stress tensor (or hydrostatic pressure) should be introduced into the yield function. In this study the research findings of Lubliner et al. [48] is the ratio of second invariants of stress tensor on tensile and compressive meridians, set as 2/3.

Flow potential function
Unlike metallic materials, the cohesive friction materials exhibit dilatancy characteristics and non-associated flow law should be used to calculate the plastic deformation [51]: where λ  is non-negative plasticity consistency parameter, the potential function G is assumed in the form of Drucker-Prager hyperbolic function. t0 σ is failure stress in uniaxial tension, in MPa. ψ is the dilatancy angle and e the eccentricity. The incremental constitutive relations can be established from yield surface and flow potential function as following [52] :

Plastic strain evolution
For the materials with tension-compression asymmetry, the plastic strain should be calculated by three principal strains weighted by distance function in stress space, i.e.
. The plastic strain can be obtained by weighting method, the material responses under complex stress states can be computed as long as the compressive and tensile curves are determined experimentally.
According to the test data of PBX [17,50,53], the material parameters could be obtained.

Computational model and analysis
In this study, the computational model is from the experimental configuration by Liu etc. [26]. The tested material is TATB-based explosive, the temperature is 50 °C, and the test is carried out by uniaxial compression. The loading rate is 1.27 mm.min -1 . The specimen geometry is a perforated rectangular plate, where a cavity with traction-free surface is located at the centre, shown in Fig. 2. The overall dimension of the plate is 76.2 mm in height, 38.1 mm in width, and 12.7 mm in thickness. The cavity at the centre of the plate has circular curved sections connected by four straight lines. The radii of curvature for the end-point and central curved sections are 6.35 mm and 1.905 mm, respectively. The length of the cavity is 19.05 mm.
In numerical simulation the plate is meshed by quadrangular plane stress elements, the element size is about 0.5mm. The maximum principal stress and fracture energy criteria are used. The results of XFEM analysis using code PANDA show that the stress and deformation field around the cavity of the plate subjected to overall compression are not uniform. The first principal stress at the upper and lower ends of the cavity is in tensile state, which is the largest on the inner wall. Fig. 3 shows the principal stress distribution along the inner wall of the cavity between point A to B. It can be seen that the first principal stress is the largest at point A and the third principal stress and shear stress are the maximum (absolute value) at point B. Therefore, the cracking mode and failure is firstly initiated on this point, and the direction of crack will be perpendicular to the direction of maximum principal stress. Secondly, the shear failure mode will initiate at point B [26].
Numerical simulation shows that non-uniform deformation and stress fields are generated near the holes in the specimen under the effect of external compressive stress, resulting in stress concentration or high stress gradient. When local stress such as tensile stress reaches material failure stress, local cracking initiates. As the external load is further applied and deformation develops, the crack rapidly grows from the initiation site along the area of maximum principal stress, resulting in two macrocracks, as shown in Fig. 4. Figure 5 shows the numerical results of crack initiation moment, crack speed and crack propagation history in comparison with the experiment [26].  (about 0.2 min), and the experimental measurement of the initial crack speed is 1.76 mm.s -1 , the numerical simulation is 1.77 mm.s -1 . This shows that the method based on XFEM and cohesive model is feasible to simulate the initiation and propagation of crack Type I. Numerical simulation can give a good prediction on the overall development of crack initiation and growth. However, in some details, such as (1) numerical simulation results show that the crack propagation has a lag compared with the measured value, (2) In the late stage of crack development, there is a certain deviation between the simulated and measured crack locations. Possible reasons are: (1) Theoretically, high enough resolution is required to identify the crack and capture the crack tip location in tests, but it is very difficult to do in actual measurements, which may lead to differences on crack location, length.
(2) In the current calculation of XFEM, the simplified algorithm leads to crack crossing one element at one time step, and the crack speed in the element is infinitely fast. When the mesh is not fine enough, it will cause certain error on the crack growth.
(3) In the simulation of cracked interfaces, the damage of material leads to the softening and stiffness degradation of the material, and the numerical algorithm encounters the problem of convergence. To overcome this difficulty, the viscosity coefficient that characterizes the relaxation time of viscous system is introduced into the algorithm of cohesive constitutive model for viscosity regularization. This also leads to a small degree of slackening effect on material failure simulation.
(4) Because of the mesoscopic heterogeneity of PBX, the mechanical parameters have a certain degree of random distribution. In the process, the perimeter of the cavity is damaged to a certain degree, which may lead to local jumps on the crack propagation curve, e.g. the bottom curve. The numerical method of this study did not account for these effects.

Concluding Remarks
(1) Under the condition of overall compressive stress, non-uniform stress and deformation fields occur around the hole, resulting in local stress concentration. When the local tensile stress reaches the material failure strength, the cracking failure mode initiates at the local scale.
With the load and deformation development, the cracking damage continues to evolve. While the critical value is reached, the crack further grows and propagates along the region of maximum principal stress, which eventually leads to the overall cracking of the material. This is the cracking mechanism of the hole-containing plate in compression. To reduce the concentration of tensile stress along the inner wall of the hole will help prevent cracking failure.
(2) The overall trend of crack propagation given by numerical simulation is in good agreement with the experimental results, including the overall trend and inflection point of cracking history, the cracking initiation time and the initial crack growth rate. It shows that the method based on XFEM and cohesive model is feasible to simulate the initiation and propagation of crack of Type I, and numerical simulation can give a good prediction on the overall development such as crack initiation and growth.
(3) There is a certain degree of random distribution of PBX macroscopic mechanical parameters duo to the mesoscopic heterogeneity of PBX. Therefore, statistical test data are required to introduce the statistical properties of the material parameters into the numerical simulation.