Numerical modeling of crack growth in polymer-bonded explosive with cavity subject to compression

Abstract

Polymer-bonded explosives are solid high-explosive particles that exhibit brittle behavior in uniaxial tension, quasi-brittle in uniaxial compression, and ductile when subjected to high confining pressure. Tension cracking is the primary failure mode of polymer-bonded explosives quasi-brittle solid, which will lead to overall failure of structural integrity. One characteristic of brittle or quasi-brittle solids, such as polymer-bonded explosives, is that when subjected to overall compressive loading, the tensile cracks will initiate inside the material due to existence of imperfection within the materials. In this study, extended finite element method is applied to analyze the cracking failure mechanism in the PBX 9502 plate-like specimen with cavity subjected to overall compression. The nonlinear constitutive behaviors and failure of polymer-bonded explosives under complex stress states were described by means of stress state–dependent strength surface, non-associated flow rule, and cohesive failure model. Analysis indicates that the tensile stress around the cavity arises in the specimen under overall compression, and this local tensile stress will lead to cracking initiation. The comparison between simulation results and the experimental data reported by Liu and Thompson shows that they are in agreement with each other on some aspects of crack behaviors, including overall development of crack history and inflection, crack initiation moment, and crack initial speed.

Introduction

Polymer-bonded explosives (PBX) are highly particle-filled composite materials comprising 90%–95% by weight of explosive crystals and 5%–10% by weight of polymeric binder. They have been used in a wide variety of applications ranging from rocket propellants to the main explosive charge in conventional munitions for civil and defense fields.^1–4 In recent decades, governments and engineering fields have attached great importance to the safety and reliability of PBX components under low-speed impact, and so the failure of PBX explosives has attracted great attention among the academic circle.^5–12

It is generally believed that the ignition mechanism of explosives under low-velocity impact includes the following: elastic-plastic or visco-elastoplastic deformation, damage and failure, crack initiation and growth, plastic localization, and plastic work conversion into heat and heat conduction. ¹³ 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. ¹⁴ 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 provide theoretical supports for assessing the risk of crack occurrence and developing crack suppression technology.

In different stress states, PBX explosive has obvious characteristics of brittleness in tension, quasi-brittleness, and confining ductility in compression. Brittle cracking under tensile stress is one of the main failure modes of solid explosive and also one of the main factors leading to the overall failure of material. Brittle or quasi-brittle solids (such as PBX) exhibit a basic feature: the non-uniform deformation field and stress field arise near the microvoids in the material, even if subjected to compressive load as a whole, due to the inhomogeneous (e.g. microvoid) in the material, and the non-uniform stress fields lead to stress concentration; when the local stress such as tensile stress reaches the failure stress of the material, the local tensile crack (microcrack) will occur. With the development of load, the microcrack will further evolve and grow and eventually lead to the overall damage of the material or structure. In order to verify this characteristic of PBX, Liu and Thompson ¹⁷ carried out the experimental study on crack and propagation of PBX 9502 explosive specimens with a hole applied by overall compressive stress. In this article, the extended finite element method (XFEM) method is used to simulate the crack behavior observed in this test.

The mechanical tests of material properties show that PBX exhibits different mechanical behaviors under different stress states,^18–20 which brings great difficulties to numerical modeling.^21–23 In this study, the stress state dependent strength model and non-associated flow rule are applied to describe the nonlinear constitutive behavior of the material under complex stress states. Meanwhile, PBX exhibits brittle/quasi-brittle characteristics in the failure mode. The local tensile microcracks initiate near the voids within material, even if subjected to compression loading as a whole, and the microcracks evolve and grow as the load develops, leading to overall failure.^17,24,25 For the purpose of characterizing the heterogeneity of brittle PBX, a perforated rectangular plate, where a diamond-shaped cavity with traction-free surface is located at the center, was chosen as a typical imperfection. The simulation was performed on this configuration subjected to compression.

In the simulation of PBX brittle tensile crack behavior, the usual methods which are used include numerical simulation method by continuous damage mechanics,^26,27 direct numerical simulation method, ²⁸ multiscale ²⁹ or mesoscopic numerical method, ³⁰ and XFEM.^31–35 In this study, the cohesive-based XFEM method was used to analyze the crack initiation and failure of PBX 9502 plate with cavity against compression failure test. ¹⁷

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. ³¹ XFEM is an extension of the conventional finite element method based on the concept of partition of unity. ³⁶ XFEM solved the problem of crack propagation in the finite element using the idea of mesh separation independent. While retaining all the advantages of the traditional finite element analysis (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 of 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.^31,33

In the simulation of stable discontinuities such as cracks, the conventional finite element method requires that the mesh be consistent with the geometric discontinuities. Therefore, the mesh should be refined at the crack tip to obtain the singular asymptotic field. When simulating moving cracks, this method becomes clumsy, because the mesh needs to be constantly modified to adapt to discontinuous geometric characteristics. The fundamental difference between XFEM and conventional finite element method is that the meshes used in XFEM are independent of the geometric or physical interfaces in the structure, which overcomes the difficulties of high-density meshing in stress concentration field such as crack tip, and the meshes in XFEM do not need to be refined in simulating crack propagation.

For the purpose of cracking analysis, the enrichment functions typically consist of the near-tip asymptotic functions that capture the singularity around the crack tip and a discontinuous function that represents the jump in displacement across the crack surfaces. The approximation for a displacement vector function u with the partition of unity enrichment is

$$\mathrm{u}=\sum _{I=1}^N{N}_I[u_I+H\left(x\right)a_I+\sum _{\alpha =1}^4{F}_{\alpha }\left(x\right)b_I^{\alpha }]$$

(1)

where $N_I$ are the usual nodal shape functions; the first term on the right-hand side of the above equation, ${\mathbf{u}}_I$ , is the usual nodal displacement vector associated with the continuous part of the finite element solution; the second term is the product of the nodal enriched degree of freedom vector, ${\mathbf{a}}_I$ , and the associated discontinuous jump function $H(x)$ across the crack surfaces; and the third term is the product of the nodal enriched degree of freedom vector, ${\mathbf{b}}_I^{\alpha }$ , and the associated elastic asymptotic crack tip functions, $F_{\alpha }(x)$ . The first term on the right-hand side is applicable to all the nodes in the model; the second term is valid for nodes whose shape function support is cut by the crack interior; and the third term is used only for nodes whose shape function support is cut by the crack tip. The crack tip functions $F_{\alpha }(x)$ are given by

$$F_{\alpha }\left(x\right)=[\sqrt{r}\sin \frac{\theta }{2},\sqrt{r}\cos \frac{\theta }{2},\sqrt{r}\sin \mathrm{\theta }\sin \frac{\theta }{2},\sqrt{r}\sin \mathrm{\theta }\cos \frac{\theta }{2}]$$

(2)

where $(r,\theta )$ is a polar coordinate system with its origin at the crack tip and $\theta =0$ is tangent to the crack at the tip. The term $\sqrt{r}\sin (\theta /2)$ takes into account the discontinuity across the crack face.

In the framework of XFEM method, the traction-separation cohesive behavior of material ³⁷ is used to simulate crack initiation and propagation. The phantom nodes, which are superposed on the original real nodes, are introduced to represent the discontinuity of the cracked elements, as illustrated in Figure 1.^38,39 When the element is intact, each phantom node is completely constrained to its corresponding real node. When the element is cut through by a crack, the cracked element splits into two parts. Each part is formed by a combination of some real and phantom nodes depending on the orientation of the crack. Each phantom node and its corresponding real node are no longer tied together and can move apart.

Figure 1.

The principle of the phantom node method.

The magnitude of the separation is governed by the cohesive law until the cohesive strength of the cracked element is zero, after which the phantom and the real nodes move independently. To have a set of full interpolation bases, the part of the cracked element that belongs in the real domain, ${\Omega }_0$ , is extended to the phantom domain $({\Omega }_{\mathrm{p}}^{-},{\Omega }_{\mathrm{p}}^{+})$ . Then, the displacement in the real domain, ${\Omega }_0$ , can be interpolated using the degrees of freedom for the nodes in the phantom domain, $({\Omega }_{\mathrm{p}}^{-},{\Omega }_{\mathrm{p}}^{+})$ . The jump in the displacement field is realized by simply integrating only over the area from the side of the real nodes up to the crack, that is, $({\Omega }_0^{+})$ . This method provides an effective and attractive engineering approach and has been used for simulation of the initiation and growth of multiple cracks in solids by Song ³⁸ and Remmers ³⁹ . The equivalent polynomial methodology originally developed by Ventura and Benvenuti ⁴⁰ for a cracked element enriched with a Heaviside enrichment function has been extended to evaluate the stiffness matrix for the cracked element composed of real and phantom nodes. It has been proven to exhibit almost no mesh dependence if the mesh is sufficiently refined.

XFEM-based cohesive behavior

The cohesive model belongs to the damage model. It was first introduced by Barenblatt ³⁷ and used the traction-separation law (TSL) to simulate the decohesion of the atomic lattice, so as to avoid the singularity of the crack tip. ⁴¹ 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, ⁴² Hillerborg’s law, ⁴³ and Bažant’s law. ⁴⁴

In this study, the crack propagation analysis is performed based on the XFEM-based cohesive segment method, the linear elastic TSL model, the damage initiation criterion, and the damage evolution law. In the elastic TSL model, the initial linear elastic behavior is assumed, followed by the onset and evolution of damage. The linear elastic behavior is expressed as the linear relationship between the cohesive stress and the open displacement of the cracking element, namely

$$\mathbf{t}=\mathbf{K}\mathbf{\delta }$$

(3)

In which

$$\begin{array}{c}\mathbf{t}={[t_n,t_s,t_t]}^{\mathrm{T}}\hfill \\ \multicolumn{1}{c}{\mathbf{\delta }={[{\delta }_n,{\delta }_s,{\delta }_t]}^{\mathrm{T}}}\\ \multicolumn{1}{c}{\mathbf{K}=\mathrm{diag}[K_{\mathrm{nn}},K_{\mathrm{ss}},K_{\mathrm{tt}}]}\end{array}$$

(4)

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 the 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

$$f=\left\{\frac{〈{\mathrm{\sigma }}_{\max }〉}{{\mathrm{\sigma }}_{\max }^0}\right\}=1$$

(5)

That is, damage is assumed to initiate when the maximum principal stress ratio reaches a value of one. Here, ${\sigma }_{max}^0$ 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

$$\mathrm{d}D=\frac{T\left(\mathrm{\delta }\right)\mathrm{d}\delta }{G_{\mathrm{eqC}}}=\{\begin{array}{c}\frac{T_{\max }(1-\mathrm{\delta }/{\mathrm{\delta }}_{\max })\mathrm{d}\delta }{G_{\mathrm{eqC}}}\hfill \\ \multicolumn{1}{c}{\frac{T_{\max }\exp (\mathrm{\delta }/{\mathrm{\delta }}_{\max })\mathrm{d}\delta }{G_{\mathrm{eqC}}}}\end{array}$$

(6)

Here, ${\sigma }_{\mathrm{eqC}}$ is the equivalent fracture energy release rate, in J m⁻², which is described in Benzeggagh and Kenane ⁴⁵ by the following formula

$$G_{\mathrm{eqC}}=G_{\mathrm{IC}}+(G_{\mathrm{IIC}}-G_{\mathrm{IC}}){\left(\frac{G_{\mathrm{II}}+G_{\mathrm{III}}}{G_{\mathrm{I}}+G_{\mathrm{II}}+G_{\mathrm{III}}}\right)}^{\eta }$$

(7)

When material undergoes damage, the normal and shear stress components are affected by the damage according to

$$\begin{array}{c}{t}_{\mathrm{n}}=(1-D){\mathit{T}}_{\mathrm{n}}\hfill \\ \multicolumn{1}{c}{t_{\mathrm{s}}=(1-D){\mathit{T}}_{\mathrm{s}}}\\ \multicolumn{1}{c}{t_{\mathrm{t}}=(1-D){\mathit{T}}_{\mathrm{t}}}\end{array}$$

(8)

The crack will not occur when normal compression state is applied, that is, $t_{\mathrm{n}}=T_{\mathrm{n}}$ . Components $T_{\mathrm{n}}$ , $T_{\mathrm{s}}$ , and $T_{\mathrm{t}}$ are the normal and shear stress components predicted by the elastic traction-separation behavior for the current separations without damage.

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

$${\delta }_{\mathrm{m}}=\sqrt{{〈{\delta }_{\mathrm{n}}〉}^2+{\delta }_s^2+{\delta }_t^2}$$

(9)

Description of plastic deformation of PBX under complex stress states

Many scholars have carried out a series of experiments and theoretical studies on the mechanical properties and constitutive relations of PBX materials. Belmas and Reynier ¹⁸ and Ellis et al. ¹⁹ conducted quasi-static uniaxial tension and uniaxial compression experiments and showed that PBX material shows brittle behavior in uniaxial tension and shows plastic hardening/softening characteristics in uniaxial compression, indicating obvious tension–compression asymmetry. The tensile and compressive tests of PBX 9502 by Thompson et al.’s test data showed that the material exhibits nonlinear characteristics such as tensile brittleness and compressive ductility at temperatures of 50 and −20°C. As shown in Figure 4, the tensile-compressive strength ratio is about 0.458, the modulus ratio is about 1.45, and the yield stress ratio is about 0.458. It can be seen that the tension–compression asymmetry of PBX 9502 under quasi-static state is mainly reflected in modulus, yield stress, maximum stress, strengthening, and softening. The asymmetry brings great difficulties to numerical modeling.

According to the deformation characteristics of triaminotrinitrobenzene (TATB)-based explosives, such as the tension–compression asymmetry at quasi-static experiments at 50°C, ²⁰ as shown in Figure 2, the pressure-dependent yield surface, the non-associated flow law, and the stress state dependent weighting function are utilized to describe the mechanical behavior of the material under complex stress conditions.

Figure 2.

Quasi-static stress–strain curve of TATB explosive in tension and compression at 50°C. ²⁰

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. ⁴⁶ are adopted, which are devoleped by Lee and Fenves ⁴⁷ by considering the difference in strength between tensile and compressive stresses. The evolution of yield surface consists of two independent evolution parameters ${\tilde{\epsilon }}_{\mathrm{t}}^{\mathrm{P}}$ and ${\tilde{\epsilon }}_{\mathrm{c}}^{\mathrm{P}}$ , which are as follows

$$\begin{array}{c}\multicolumn{1}{c}{F(\mathbf{\sigma },{\tilde{\mathbf{\epsilon }}}^{\mathrm{p}})=\frac{1}{1-\alpha }[q-3\alpha p+\beta \left({\tilde{\epsilon }}^{\mathrm{p}}\right)〈{\hat{\sigma }}_{\max }〉-\gamma 〈-{\hat{\sigma }}_{\max }〉]}\\ \multicolumn{1}{c}{--{\hat{\sigma }}_{\mathrm{c}}\left({\tilde{\epsilon }}_{\mathrm{c}}^{\mathrm{P}}\right)=0}\end{array}$$

(10)

where $\mathbf{\sigma }$ is stress tensor, ${\tilde{\epsilon }}^{\mathrm{p}}=[{\tilde{\epsilon }}_{\mathrm{t}}^{\mathrm{P}},{\tilde{\epsilon }}_{\mathrm{c}}^{\mathrm{P}}]$ , and $p=-(1/3)\mathbf{\sigma }:\mathbf{I}$ is the hydrostatic pressure. $q=\sqrt{(3/2)\mathbf{S}:\mathbf{S}}$ is the von Mises equivalent stress. $\alpha$ , $\beta$ and $\gamma$ are dimensionless material parameters, defined as $\alpha =\left[({\sigma }_{\mathrm{b}0}/{\sigma }_{\mathrm{c}0})-1\right]/\left[2({\sigma }_{\mathrm{b}0}/{\sigma }_{\mathrm{c}0})-1\right]$ , $\beta =[{\sigma }_{\mathrm{c}}({\tilde{\epsilon }}_{\mathrm{c}}^{\mathrm{p}})/{\sigma }_{\mathrm{t}}({\tilde{\epsilon }}_{\mathrm{t}}^{\mathrm{p}})](1-\alpha )-(1+\alpha )$ , and $\gamma =3(1-K_{\mathrm{c}})/(2K_{\mathrm{c}}-1)$ . The frictional angle $\alpha =\tan \phi /3$ , set as 20° according to Gruau et al. ⁴⁸ ${\hat{\sigma }}_{\max }$ is maximum principal stress and $({\sigma }_{\mathrm{b}0}/{\sigma }_{\mathrm{c}0})$ is strength ratio of biaxial compression and uniaxial compression, set as 1.16 according to investigation into cohesive friction materials such as concrete. $K_{\mathrm{c}}=q_{\mathrm{TM}}/q_{\mathrm{CM}}$ is the ratio of second invariants of stress tensor on tensile and compressive meridians, set as 2/3. ${\sigma }_{\mathrm{t}}({\tilde{\epsilon }}_{\mathrm{t}}^{\mathrm{p}})$ and ${\sigma }_{\mathrm{c}}({\tilde{\epsilon }}_{\mathrm{c}}^{\mathrm{p}})$ are cohesions obtained by experiments.

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 ⁴⁹

$${\stackrel{·}{\epsilon }}^{\mathrm{p}}=\stackrel{·}{\lambda }{\nabla }_{\mathrm{\sigma }}\mathrm{\Phi }\left(\mathbf{\sigma }\right)=\stackrel{·}{\lambda }\frac{\partial G\left(\mathbf{\sigma }\right)}{\partial \mathbf{\sigma }}$$

(11)

$$G=\sqrt{{\left(e{\sigma }_{\mathrm{t}0}\tan \psi \right)}^2+q^2}-p\tan \psi$$

(12)

where $\stackrel{·}{\lambda }$ is the non-negative plasticity consistency parameter and the potential function G is assumed in the form of Drucker–Prager hyperbolic function. ${\sigma }_{\mathrm{t}0}$ is the failure stress in uniaxial tension, in MPa. $\psi$ is the dilatancy angle and e is the eccentricity. The incremental constitutive relations can be established from yield surface and flow potential function as following ⁵⁰

$$\mathrm{d}\mathbf{\sigma }=(\mathbf{C}-\frac{1}{h}\nabla F·\mathbf{C}\otimes \mathbf{C}·\nabla G):\mathrm{d}\mathbf{\epsilon }$$

(13)

where $h={\nabla }_{\sigma }F·\mathbf{C}·{\nabla }_{\sigma }G-{\nabla }_{{\epsilon }^{\mathrm{p}}}F·{\nabla }_{\sigma }G-C{f}_{,k}{f}_{,{\epsilon }_{\mathrm{p}}}\sqrt{\nabla G·\nabla G}$ .

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, that is

$${\stackrel{·}{\stackrel{~}{\epsilon }}}^{\mathrm{p}}=\left[\begin{array}{c}{{\stackrel{·}{\stackrel{~}{\epsilon }}}_{\mathrm{t}}}^{\mathrm{p}}\hfill \\ \multicolumn{1}{c}{{{\stackrel{·}{\stackrel{~}{\epsilon }}}^{\mathrm{p}}}_{\mathrm{c}}}\end{array}\right]=\hat{\mathbf{h}}(\hat{\mathbf{\sigma }},{\hat{\mathbf{\epsilon }}}^{\mathrm{p}})·{\widehat{\stackrel{·}{\mathbf{\epsilon }}}}^{\mathrm{p}}$$

(14)

where $\hat{\mathbf{h}}(\hat{\mathbf{\sigma }},{\hat{\mathbf{\epsilon }}}^{\mathrm{p}})=\left[\begin{array}{ccc}r\left(\hat{\mathbf{\sigma }}\right)\hfill & 0\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & -[1-r\left(\hat{\mathrm{\sigma }}\right)]\hfill \end{array}\right]$ and ${\widehat{\stackrel{·}{\mathbf{\epsilon }}}}^{\mathrm{p}}={[{\widehat{\stackrel{·}{\epsilon }}}_1^{\mathrm{p}},{\widehat{\stackrel{·}{\epsilon }}}_2^{\mathrm{p}},{\widehat{\stackrel{·}{\epsilon }}}_3^{\mathrm{p}}]}^{\mathrm{T}}$ . The plastic strain can be obtained by integration ${\hat{\epsilon }}_{\mathrm{t}}^{\mathrm{p}}={\int }_0^{t}r\left(\mathbf{\sigma }\right){\widehat{\stackrel{·}{\epsilon }}}_{\max }^{\mathrm{p}}$ . $r(\hat{\mathrm{\sigma }})$ is the distance function, and the weighting factor of triaxial stress defined as

$$r\left(\hat{\mathbf{\sigma }}\right)\stackrel{\mathrm{def}}{=}\frac{\sum _{i=1}^3〈{\hat{\sigma }}_i〉}{\sum _{i=1}^3\left|{\hat{\sigma }}_i\right|}$$

(15)

where ${\hat{\sigma }}_i$ (i = 1,2,…) are the principal stresses, $〈x〉=\frac{1}{2}(\left|x\right|+x)$ is the Macaulay bracket. From equation (13), $r(\hat{\mathrm{\sigma }})$ is equal to 1 when all principal stresses ${\hat{\sigma }}_i$ are positive and is equal to 0 when ${\hat{\sigma }}_i$ are negative. By means of this 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,^18,48,51 the material parameters are listed in Table 1.

Table 1.

Material parameters of PBX 9502 at 50°C.

Parameter	Value
E (GPa)	2.3
v	0.4
$\phi$	20°
$\psi$	1°
e	0.1
${\sigma }_{\mathrm{s}}^{\mathrm{c}}$ (MPa)	7.12
${\sigma }_{\mathrm{s}}^{\mathrm{t}}$ (MPa)	3.2
${\sigma }_{\mathrm{c}0}$ (MPa)	17.8
${\sigma }_{\mathrm{t}0}$ (MPa)	8.0
$G_{\mathrm{IC}}$ (J m−2)	100.0

E: elastic modulus; v: Poisson’s ratio; $\phi$ : friction angle; $\psi$ : dilation angle; e: eccentricity; ${\sigma }_{\mathrm{s}}^{\mathrm{c}}$ : yield stress in compression; ${\sigma }_{\mathrm{s}}^{\mathrm{t}}$ : yield stress in tension; ${\sigma }_{\mathrm{c}0}$ : compression strength; ${\sigma }_{\mathrm{t}0}$ : tensile strength; $G_{\mathrm{IC}}$ : fracture energy.

Computational model and analysis

In this study, the computational model is from the experimental configuration by Liu and Thompson. ¹⁷ 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⁻¹. The specimen geometry is a perforated rectangular plate, where a cavity with traction-free surface is located at the center, as shown in Figure 3. The overall dimensions of the plate are as follows: 76.2 mm in height, 38.1 mm in width, and 12.7 mm in thickness. The cavity at the center 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 and 1.905 mm, respectively. The length of the cavity is 19.05 mm.

Figure 3.

Specimen plate with cavity for fracture experiment on PBX 9502 subject to compression.

In numerical simulation, the plate is meshed by quadrangular plane stress elements. The element size is about 0.5 mm. 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. Figure 4(b) shows the principal stress distribution along the inner wall of the cavity between points A and 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 first initiated on this point, and the direction of crack will be perpendicular to the direction of maximum principal stress. Second, the shear failure mode will initiate at point B. From the measured correlation coefficient field, these results can also be confirmed, ¹⁷ see Figure 5.

Figure 4.

Principal stresses distribution along inner wall from A to B: (a) Schematic diagram of tested specimen marked with points A and B and (b) Principal stresses distribution along inner wall from A to B.

Figure 5.

(a) Description of cavity boundary with polar angle and (b) variation of the correlation coefficient along the cavity boundary.

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 shown in Figure 6. This demonstrates that the tensile cracking originates from heterogeneity in brittle solid subjected to overall compression.

Figure 6.

Crack initiation: (a) crack initiation position and (b) partial enlarged view.

As the external compressive load is further applied, the deformation develops and the crack rapidly grows from the initiation site (i.e. point A in Figures 4 and 6) along the direction of maximum principal stress, resulting in the two macrocracks. Figure 7 shows crack growth and macrocrack position at different moments.

Figure 7.

The crack growth and macrocrack position at different moments.

Figure 6(a) shows the experimentally measured cracks and 6(b) shows the numerically simulated cracks in this study. Figure 8 presents the eventual cracks, originating from point A on the inner wall of the cavity. Figure 8(a) is the experimentally measured cracks, ²¹ while Figure 8(b) is the numerically simulated cracks.

Figure 8.

Crack initiation and growth of PBX 9502 plate specimen with cavity subjected to overall compression: (a) measured cracks and (b) simulated cracks.

Figure 9 shows the prediction of the mixed crack. When the end-face displacement is applied to both the longitudinal and transverse directions, a mixed-type I–II crack occurs near the apex of the inner wall of the cavity.

Figure 9.

Prediction of mixed crack of type I and II.

Figure 10 indicates the numerical results of crack initiation moment, crack speed, and crack propagation history in comparison with the experiment. ¹⁷ The numerical results show that the trend of crack growth is consistent with the experimental results, for example, the overall trend and inflection point of the upper and lower cracks. The cracking initiation moments for experiment and simulation are both at $\delta /H=0.32\%$ (about 0.2 min). The initial crack speed is approximately constant. The experimental measurment is 1.76 mm s⁻¹, while the numerical simulation is 1.77 mm s⁻¹. They are very close. The simulation also shows that the top and bottom cracks have a slight difference in growth history due to friction, which the experiment also indicates.

Figure 10.

Crack initiation moment, initial crack speed, and growth history.

From the comparison in detail, it shows that the nonlinear constitutive description and method based on XFEM are feasible to simulate the initiation and propagation of crack Type I of PBX. 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 and (2) in the late stage of crack development, there is a certain deviation between the simulated and measured crack locations. Possible reasons are as follows:

1. Theoretically, enough high resolution is required to identify the crack and capture the crack tip location in experiments, but it is very difficult to do in actual measurements, which may lead to some differences of crack location between experiment and simulation.

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 and time step is not small, some difference of the crack growth will arise between experiment and simulation.

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, for example, the bottom curve. The numerical method of this study did not account for these effects.

Concluding remarks

1. For the quasi-brittle solid PBX, heterogeneity is the primary inducement for local tensile cracking, even if compression is applied. 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 location. 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. Reducing 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 cracking initiation time, the initial crack growth rate, and the inflection point of cracking history. It shows that the nonlinear constitutive description and the method based on XFEM and cohesive model are suitable for simulating the initiation and propagation of crack of Type I, and numerical simulation can give a good prediction on the overall developments such as crack initiation and growth.

3. There is a certain degree of random distribution of PBX macroscopic mechanical parameters due 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.