- Email: [email protected]
Peridynamic modeling of dynamic damage of polymer bonded explosive
Computational Materials Science xxx (xxxx) xxxx
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Peridynamic modeling of dynamic damage of polymer bonded explosive Xiaoliang Deng , Bo Wang ⁎
National Key Laboratory of Shock Wave and Detonation Physics, Institute of Fluid Physics, China Academy of Engineering Physics, Mianyang 621999, Sichuan, PR China
ARTICLE INFO
ABSTRACT
Keywords: Dynamic damage Polymer bonded explosives Peridynamic Cracks Impact loading
Dynamic damage response of 1,3,5,7-tetranitro-1,3,5,7-tetrazocane (HMX) based polymer bonded explosive (PBX) is investigated by using recently developed three dimensional Peridynamic (PD) method, allowing natural evolution of defects such as cracks or voids. An in-house PD codes are developed and verified by comparing with published results on Kalthoff-Winkler (KW) and crack branching benchmark tests. Numerical modeling of PBX is constructed via random micro-modulus approach, capable of capturing inherent microstructural heterogeneity of PBX materials. The simulation results indicate that damage process depends on both impact speed and interfacial property. Damage mode index (DMI), defined as the ratio of HMX-HMX damage to the summation of Binder-Binder and HMX-Binder damage, is proposed to quantify competitive damage mechanisms between trans- and inter-granular damage. The analyses of DMI exhibit that the inter-granular damage is dominant for the low adhesive strength of interface. However, the trend of trans-granular damage increases with increasing impact speed. The trans-granular damage mode becomes dominant for the high adhesive strength of interface. But the trend of inter-granular model increases with increasing impact speed. The current research highlight the competitive mechanisms of two typical damage modes as a function of material property and external loading conditions in a quantitative manner, providing additional insight into the damage mechanisms of PBX under impact loading.
1. Introduction Polymer bonded explosive (PBX) is one of solid high explosives and widely encountered in a variety of applications such as rocket propellants and main explosive charge in different weapon systems. PBX is composed of energetic materials particles (90–95% by weight) such as HMX and RDX embedded in polymeric binder (5–10% by weight) as well as additives [1]. PBX can experience some unexpected stimuli such as low-velocity impact during the manufacture, storage, transportation, and handling process. Those stimuli can lead to unwanted ignition of PBX due to the development of localized regions of high temperature called as the hot-spots inside the energetic materials. Suggested mechanisms of hot spot include adiabatic compression of trapped gases, hydrodynamic pore collapse, and localized shearing deformation and so on [2]. In impact and cook-off experiments on PBX, complex crack patterns are observed and reactive wave propagates along the cracks. On one hand, cracks can enhance chemical reaction by providing additional burning surface area. On the other hand, product gases produced from chemical reaction are able to increase local pressure and drive cracks propagation thus cause additional damage to PBX [3,4]. Overall, the mechanical damage and chemical reaction interplay to
⁎
determine the eventual response of PBX under impact loading. A significant improvements in predicating coupled mechanical and chemical response of PBX can result from a better understanding of damage evolution of PBX. Therefore, it is fundamentally important to investigate dynamic mechanical properties of PBX in order to deepen understanding the issues associated with non-shock ignition of PBX [4]. Over the years, Numerical simulations play an important role in investigating dynamic mechanical response of PBX. The previously published literature majorly utilize the finite element method (FEM) to simulate dynamic response of PBX [5–7]. For example, Xiao et al., predicts the stress–strain relationship of PBX under impact loading via implementing constitute including damage in ABAQUS. The simulation results are in good agreement with their experimental data [8]. Wang et al., perform mesoscale simulations of thermal-mechanical response of HMX based PBX and granular explosive under impact loading. The useful information such as equivalent strain and temperature distribution are achieved [9]. Regarding to FEM simulation, the accurate constitutive relation is required in order to simulate the complex dynamic behaviors of PBX. The Visco-Scram (viscoelastic statistical crack mechanics) model, including the viscoelasticity response and interior micro-cracks evolution, has been widely applied to describe the
Corresponding author. E-mail address: [email protected] (X. Deng).
https://doi.org/10.1016/j.commatsci.2019.109405 Received 31 August 2019; Received in revised form 17 October 2019; Accepted 10 November 2019 0927-0256/ © 2019 Elsevier B.V. All rights reserved.
Please cite this article as: Xiaoliang Deng and Bo Wang, Computational Materials Science, https://doi.org/10.1016/j.commatsci.2019.109405
Computational Materials Science xxx (xxxx) xxxx
X. Deng and B. Wang
dynamic damage response of PBX [8,10,11]. Recently, a novel viscoelastic Microplane model was proposed to simulate behavior of quasibrittle PBX materials in literature [12]. The features of this model are that it considers the viscoelastic behavior along with the damage induced anisotropy. Those simulation results can help better understand behaviors of mechanical damage and ignition. Although FEM simulation methods are currently the most popular numerical modeling to examine dynamic response of PBX in previously published literature. It should be noted that the mesh-based simulation method cannot efficiently handle non-continuity induced by the evolution of cracks and voids and so on because the governing equations are the spatial differential expressions. In order to treat problems associated with noncontinuity, special techniques are usually required to modify the traditional FEM models [13]. For example, the computation meshes can be regenerated during the simulation after detecting the new surfaces created by cracks propagation or voids evolution. This technique can ensure the discontinuous surface develops along the mesh lines but require to track new surfaces continuously inside the simulation domain as simulation time goes on, which suffers from time-intensive operations and also complex algorithms [14,15]. Therefore, it would be of considerable value to apply novel simulation methods which are capable of simulating damage evolution of PBX under impact loading in a more natural and efficient manner. Peridynamic (PD) theory is a new nonlocal theory recently proposed by Silling and his coworkers [16,17]. In contrast to the FEM model, governing equation in PD is integral operations rather than spatial differential equations. Therefore, it avoids the assumptions on differentiability and continuity of displacement field, allowing spontaneous crack nucleation, propagation, and interaction. This advantage makes PD a suitable candidate for simulating phenomena such as damage and fracture of materials under dynamic loading, where discontinuities such as voids or cracks usually nucleate and evolve in the damage areas of inside materials [18–20]. PD modeling has been widely used to investigate mechanical properties in a wide range of materials such as glass, metal, and composite materials and so on. Those previous works achieve the proposing results and also prompt the further development of PD theory [21–23]. Therefore, it is of particular interest to apply PD modeling investigating the dynamic response of PBX. Random micro-modulus approach is proposed to construct simulation model for the sake of simplicity, while capturing the inherent heterogeneity of PBX. The simulation results reveal typical damage modes of PBX observed in experiments, indicating the applicability of PD in terms of the investigation of dynamic damage response of PBX. The damage behavior at various impact speeds and interface parameters are systematically investigated. Damage mode index is proposed to quantify different types of damage of PBX. The results can potentially contribute to a deeper understanding of the dynamic behavior of energetic materials.
which is usually taken to be a sphere of radiusδcentered at point x in 3 dimensional model. Establishing the pair force function is a fundamental problem in peridynamics because it contains all the material constitutive information. For a prototype microelastic brittle material (PBM), f can take the following form:
f = cs [1
s=
c=
f (u (x , t ), u' (x ', t ), x , x ' , t ) dVx ' + b (x , t )
|| + || || || || ||
(2) (3)
12E
(4)
4
In bond-based peridynamics, a bond will break and can no longer carry any force if its elongation is greater than a critical relative elongation s0 , and the damage function is defined as
µ (t , ) =
0, ifs < s0 1, otherwise
(5)
Silling and Askari correlate the s0 with the fracture energy G0 of materials which represents the energy per unit fracture area in three dimensions required to terminate all bonds initially crossing the fracture surface [18]. The critical stretch s0 is described as 5G0 , depending 6E on both material properties and the horizon. In current model, bond failure is irreversible, i.e. broken bond cannot reconnect during the simulation. Based on the damage function, the level of damage at each material point can be expressed as
(x , t ) = 1
H
µ (t , ) dV ' H
dV '
(6)
which is actually the ratio of amount of broken bonds to the total amount of bonds initially associate with point x. The simulation domain is discretized into cubic lattice with the lattice constant of x and the volume of x 3 . According to the procedures described above, the interaction forces acting on all material points can be computed at each time step and an explicit VelocityVerlet algorithm is used for time integration to yield the motion of all of the individual material points [24]:
un + 1/2 = un +
Peridynamics regard the materials as the assembly of large amount of material points with mass and volume. Each point (x) can interact with all other points (x′) within a limited distance called horizon (H) through the pairwise force function f. The physical meaning of force function is the force vector per unit volume squared that neighboring point x′ exerts on the point x. In Peridynamics, this direct physical interaction between materials points is called bond. The acceleration of point x can be obtained according to the Newton’s second law of motion [18]: H
+ || + ||
where c, s, μ represent the micro-modulus, elongation of peridynamics bond, and damage function, respectively. = x ' x and = u' u denote the relative position vector of two material points in the reference configuration and their relative displacement, respectively. As a result, + represents the current relative position in the deformed configuration. Micro-modulus can be derived by strain energy density equivalence to a classical continuum theory of elasticity for the same materials and the same deformation. In 3D, c can be obtained as a function of Young’s modulus E and the radius of horizon [18]:
2. Bond-based Peridynamic formulation
u¨ (x , t ) =
µ (t , )]
t n u¨ 2
un + 1 = un + tun + 1/2 un + 1 = un + 1/2 +
t n+ 1 u¨ 2
(7) (8) (9)
where t is the time step. 3. Verification of codes On the basis of PD theory depicted above, the PD computation codes are developed in our laboratory. It is of importance to validate the codes before applying it to simulate dynamic response of PBX. Therefore, two benchmark examples with obtained results are presented and compared with previous results from published literature in this section.
(1)
where ρ, u, V, b represent the mass density, displacement vector field (i.e. u¨ is its current acceleration vector field), infinitesimal volume of material point, and prescribed body force density field. H is the horizon 2
Computational Materials Science xxx (xxxx) xxxx
X. Deng and B. Wang
Fig. 1. Simulation of Kalthoff–Winkler’s experiment. (a) Geometry of computational model, (b) Simulation result on crack propagation.
3.1. Example 1: Kalthoff–Winkler experiment
observation and previous simulation results. Both the results on KW and crack branching benchmark tests validate the developed PD codes.
The Kalthoff-Winkler experiment is a classic benchmark problem for studying dynamic fracture and peridynamics dynamics simulation results on this problem are reported in many literature [25,26]. The simulation model corresponding to the experimental setup is presented in Fig. 1. The thickness of the plate is set to be 0.01 m. The experimental results indicate that if a pre-cracked steel plate made of steel 18Ni1900 is hit literally with impact speed of 32 m/s, a brittle fracture mode is observed and the cracks initiates at the end of the pre-existing cracks and then propagate along around 70 degree with respect to the preexisting cracks. In order to compare with previous published results, the same initial simulation parameters are used, i.e. Young’s Modulus E = 190 GPa, density ρ = 7800 kg/m3, Poisson’s ratio ν = 0.25, and critical energy release rate of the material G0 = 6.9 × 104 J/m2. The first three layers of grids (red color in Fig. 1(a)) are applied initial velocity of 22 m/s to simulate impact loading for the sake of simplicity. Simulation domain is divided with uniform grid of x = 1 mm which leads to totally 200,000 grids in the simulation. The radius of each material point’s horizon is about three times of the grid spacing e.g., = 3.015 x [18].. The simulation results show the cracks start to initiate at about 20 µs and then propagate at the angle of 67 degree with respect to the initial direction of pre-existing cracks. Those results are in good agreement with previously published simulation results [26].
4. Modeling of dynamic damage of PBX 4.1. Problem setup and computational parameters PD model of PBX used in this research is shown in Fig. 3. A rectangular sample is constructed with dimensions along x, y, and z are 4.5 mm, 1.0 mm, and 1.0 mm, respectively. There are total 50 grains in the rectangular sample and the average grain size is therefore about 440 µm . The computational domain is divided into three-dimensional rectangular grids with total number of about 220,000. The three layers of grid along x direction are fixed as a piston moving at a constant speed during the simulation. The free boundary conditions are used in simulations. The time increment △t is set to be 1.0e-10s and the duration of simulation is 3 µs . The materials used in numerical simulations are HMX with densityρ = 1900 kg/m3, Young’s modulus E = 24.0 GPa, critical relative elongation s0 = 0.0005 [28] and binder with densityρ = 1100 kg/m3, Young’s modulus E = 3.0 GPa, critical relative elongation s0 = 0.02 [29]. Voronoi tessellation method is used to generate PBX polycrystalline structure. In this method, some seed points are first generated in a random manner, then the distances between material points and seed points are computed. The geometries of HMX crystals are determined according to distances. It should be noted that the previous experimental investigation reveals that it is usually inevitable to introduce damage and defects such as micro-voids and dislocations within the HMX crystal during the manufacturing process of PBX [30–33]. Those defects distributed inside the PBX in a random manner. As a result, it is very hard, if it is not completely impossible, to construct a computation model directly based on the real microstructures and distributions of those defects due to the current lack of detailed experimental characterization of PBX. However, it is well-known that the amount and distribution of defects can impact material strength [34]. Therefore, more defects inside a specific HMX crystal can lead to a less strength of this HMX crystal and vice visa. In PD theory, the material strength is directly associated with micro-modulus. Therefore, a random micromodulus model is proposed in this manuscript to represent intrinsic heterogeneous microstructures of PBX for the sake of simplicity. In this model, micro-modulus in the range of [0.5c,1.5c] are randomly distributed to each HXM grain, where c is computed according to Eq. (4). Such micro-modulus distributions essentially represent the effect of random amount and distribution of defects on the strength of HMX crystal. From Fig. 3, it can be observed that the PD model of PBX consists of HMX and binder particles. These two types of PD particles can constitute three types of bonds: HMX-HMX, Binder-Binder, and HMXBinder, which actually represent the properties of HMX crystal, binder, and the interface between HMX crystal and binder, respectively. It
3.2. Example 2 crack branching in a pre-cracked plate An initial crack is prescribed in a thin rectangular plate with dimensions shown in Fig. 2(a). The upper and lower three layers of grids along the vertical direction are suddenly applied a tensile loading of 12 MPa, then tensile stress keeps constant during the whole simulation process. This type of loading condition generates a stress wave, which can interact with dynamic propagation of cracks and eventually result in crack branching [19,21,27]. The brittle material used in this simulation is Duran 50 glass with densityρ = 2235 kg/m3, Young’s modulus E = 65 GPa, Poisson’s ratioν = 0.2, and critical energy release rate G0 = 204 J/m2. Simulation domain is divided with uniform grid of x = 0.4 mm in the simulation. The radius of each material point’s horizon is about three times of the grid spacing e.g., = 3.015 x . Peridynamic simulation results are shown in Fig. 2(b)–(d). Pre-existed crack first propagates horizontally for a while and then starts to branch. The pattern of crack branching shown in Fig. 2 is similar with previously published results [19]. Simulation results indicate that branching event occurs at about 24 microseconds and crack speed reaches its maximum value of about 1800 m/s. Maximum speed is equal to 0.56 CR (CR is Rayleigh wave speed). It has been experimentally shown that the limiting crack speed in glass ranges from 0.47 CR to 0.66 CR for different loading conditions [20]. Therefore, our simulation data about crack branching in glass is in good agreement with experimental 3
Computational Materials Science xxx (xxxx) xxxx
X. Deng and B. Wang
Fig. 2. Simulation results about crack branching study at various time steps. Different colors indicate different damage levels. (a) 0 μs, (b) 20 μs, (c) 25 μs, (d) 40 μs.
should be noted that interfacial properties is hard to characterize experimentally. Previously published results in literature are inconsistent with each other. For example, the works of adhesion between TATB and Kel-F 800 interfaces are reported to be 76–271 mJ/m2 and 43 mJ/m2 [35,36]. In order to systematically investigate the effect of interfacial properties between HMX crystal and binder on dynamic damage response of PBX, an adjustable interface parameterβis introduced in this research. β is defined as the following equation:
PHMX
Binder
= Pmin + (Pmax
Pmin)
4.2. Dynamic damage of PBX at various impact speeds Dynamic response of PBX at various impact speeds is investigated in this section. According to previously published literature, various damage types such as trans-granular and inter-granular cracks under impact loading can be observed [37–39]. Trans-granular cracks usually describe the damage of individual bulk crystal. In the contrast, intergranular fracture refers to the damage along the interfaces between HMX and binder. In order to quantify the damage degree of HMX, binder and their interfaces, the broken bonds number of HMX-HMX, Binder-Binder, and HMX-Binder normalized by their own initial values are computed and presented. The quantitative characterization of damage is important because it can help assess the mechanical integrity of PBX under impact loading. The results concerning the damage of HMX-HMX, Binder-Binder, and HMX-Binder bonds at various impact speeds are presented in Fig. 4. Fig. 4(a)–(d) correspond to the impact speed of 20 m/s, 30 m/s, 40 m/s, and 60 m/s. = 0.5 is used for all the different impact speeds. The minor damage is observed at impact speed of 20 m/s. With increasing impact speed, damage becomes more and more serious. More precisely, HMX-HMX damage increases from about 6% to 50% with increasing
(10)
where PHMX Binder represents critical relative elongation of HMX-Binder bonds and Pmin, Pmax represent minimum and maximum values critical relative elongation between HMX crystal and binder, respectively. According to Eq. (10), critical relative elongation achieves minimum or maximum value corresponding to = 0 or = 1, respectively. When = 0.5, interfacial parameter is equal to the average value of HMX and binder. Therefore, the impact of interfacial mechanical properties on dynamic damage of PBX can be investigated via controllable change of interface parameter , illustrating the advantage of employing numerical computation for this type of systematic investigations.
Fig. 3. Peridynamic model of PBX. Blue and red particles represent binder and HMX, respectively. Three different types of bonds, including HMX-HMX, HMX-Binder and Binder-Binder bonds, formed in the PD model are also shown in the picture. 4
Computational Materials Science xxx (xxxx) xxxx
X. Deng and B. Wang
Fig. 4. Influence of different impact speeds on dynamic damage of PBX. (a) 20 m/s, (b) 30 m/s, (c) 40 m/s, (d) 60 m/s. Three different types of damage, e.g. HMXHMX, Binder-Binder, and HMX-Binder are shown, which are defined as the ratio of their broken bond numbers to their own initial values.
Fig. 5. Snapshots of simulation results at different impact speeds. (a) 20 m/s, (b) 30 m/s, (c) 40 m/s.
impact speed from 20 m/s to 60 m/s. However, qualitative trend is similar for all different impact speeds, i.e., the damage of HMX-HMX is dominant in all three types of damage modes, regardless of impact speed. This conclusion is also further verified by visualized examination of simulation results (Fig. 5). It can be observed that most of damage appears inside the HMX crystal and damage of binder material is relatively low over an impact speed range from 20 to 40 m/s. But the question remains, whether dominant damage are always the damage of HMX-HMX. Therefore, addition simulations concerning of effect of material property on dynamic response are conducted in next section of this manuscript.
4.3. Dynamic damage response of PBX at various interface parameters Dynamic damage of PBX at various impact speeds shows that the damage in PBX mostly occur inside the HMX crystals. This phenomenon is somehow inconsistent with experimental observations. Experimental scanning electron microscopy (SEM) results indicate damage along HMX-binder interfaces can also be observed [8,40]. In order to systematically understand dynamic damage of PBX, the influence of interface parameter of β on dynamic damage is evaluated in this section. Fig. 6 presents the impact of interfacial parameters on dynamic damage of PBX under various impact velocities. It indicates that 5
Computational Materials Science xxx (xxxx) xxxx
X. Deng and B. Wang
Fig. 6. Influence of interfacial parameters on dynamic damage of PBX. (a) β = 0,(a) β = 0.5,(a) β = 1.
shown in Fig. 3, the breakage of HMX-HMX bonds is associated with the trans-granular cracks, referring to this damage mode as trans-granular damage in this manuscript. Similarly, the breakage of Binder-Binder bonds actually represent the binder material fracture and the breakage of HMX-Binder bonds is associated with the interface debonding. However, due to the relative low weight ratio of binder materials in PBX, the dimension of interface between HMX crystal and binder material is usually small compared to that of HMX crystal. Therefore, it is actually hard to experimentally distinguish between binder material fractures and interface debonding. The general experimental phenomena resulting from these two damage modes is that the cracks propagate along the interface between HMX crystal and binder material, referring to these two damage modes as inter-granular damage for the sake of convenience. Damage mode index (DMI), defined as the ratio of HMX-HMX damage to the summation of Binder-Binder and HMX-Binder damage, is proposed in this manuscript. Introduction of DMI can help quantitatively analyze the relative trend of trans-granular and inter-granular damage of PBX. According to its definition, DMI is greater than one if the trans-granular damage predominates over the inter-granular damage. Otherwise, DMI is less than one if inter-granular damage dominates over trans-granular damage. Based on the simulation data, the DMI at various interfacial parameters and impact speeds can be computed, which is shown in Fig. 7. Fig. 7(a), corresponding to β = 0, shows that DMI is always less than one for all impact speeds, indicating inter-granular damage is dominant. However, it can be observed that the value of DMI increases with increasing impact speed, indicating the trend of trans-granular damage increases with increasing impact speed. Fig. 7(b) and (c), corresponding to β = 0.5 and β = 1, respectively, show that the DMI is always greater than one for all impact speeds, indicating trans-granular damage is dominant. However, the specific value of DMI decreases with increasing impact speed, indicating the trend of inter-granular damage increases with increasing impact speed. Therefore, damage mode of PBX actually depends on both intrinsic material property and external loading condition. Both of them
different interfacial parameters lead to different damage modes. For β = 0, corresponding to the weak interface strength, the dominant damage is HMX-Binder. With increasing value ofβ, the dominant damage mode transfer from HMX-Binder to HMX-HMX. This phenomenon is understandable because the smaller value of β leads to less critical elongation of bonds formed by HMX and binder PD particles. Therefore, the interfacial bonds of HMX and binder can be easily broken. With increasing value ofβ, interfacial bonds of HMX and binder become harder and harder to be broken due to increasing critical elongation. The visualization of simulation results showing in the three images at the bottom of Fig. 6 is consistent with quantitative results about damage of HMX-HMX, Binder-Binder, and HMX-Binder. For example, the heavy damage contour of interface between HMX and binder is observed in the bottom image of Fig. 6(a), where the damage occurs predominantly at the interface. At the same time, quantitative results in the top image of Fig. 6(a) indicate the damage mode of HMX-Binder is dominant. The thin slice of PBX located in the center of simulation domain is also cut from simulation domain in order to achieve detailed fracture information, which is not shown for the sake of brevity. The cracks propagating along interfaces between HMX and binder are observed, indicating the inter-crystal fracture mode. Moreover, damage modes for β = 0.5 and β = 1 are similar. Both of them indicate the damage of HMX-HMX is predominant. Therefore, it can be concluded that if interfacial strength is high enough, it no longer influences damage behaviors of PBX in the range of impact speed considered in current research. Another observation is that the damage of Binder-Binder also decreases with increasing value ofβ. This result indicates that the strong interface strength between binder particles and HMX particles can also prevent bond break of Binder-Binder. 4.4. Damage mode index (DMI) The above simulation results indicate the possible damage behaviors of PBX at PD particle level include the breakage of HMX-HMX, Binder-Binder, and HMX-Binder bonds. According to the schematic 6
Computational Materials Science xxx (xxxx) xxxx
X. Deng and B. Wang
Fig. 7. Damage mode index (DMI) at various interfacial parameters as a function of impact speed. (a) β = 0, (b) β = 0.5, (c) β = 1.
which allows spontaneous formation and evolution of cracks in a consistent framework. The in-house PD codes are developed and verified via two classic examples: KW test and crack branching under tension loading. A random micro-modulus model is proposed to represent inherent heterogeneity of PBX induced by various amounts and distributions of defects inside the materials. PD simulation results exhibit that damage degree increases with increasing impact speed. Quantitative analysis shows for the fixed interface parameter β = 0.5 breakage number of HMX-HMX bonds are greater than that of Binder-Binder and HMX-Binder bonds. When interface parameter β decreases from 0.5 to 0, the breakage number of HMX-Binder bonds is greater than that of HMXHMX and Binder-Binder bonds. In order to systematically analyze damage patterns for various impact speeds and material properties, damage mode index (DMI), defined as the ratio of HMX-HMX damage to the summation of Binder-Binder and HMX-Binder damage, is proposed and computed based on the simulation data. The DMI analyses reveal that inter-granular damage is the dominant damage mode for different impact speeds for β = 0.0. However, trend of trans-granular damage increases with increasing impact speed. On the other hand, trans-granular damage is dominant for different impact speeds for the interface parameter β = 0.5 and β = 1.0. However, trend of inter-granular damage increases with increasing impact speed. The observed simulation results highlight the combined effects of material properties and loading conditions on dynamic damage of PBX, helping better understand dynamic mechanical behaviors of heterogeneous energetic materials.
interplay with each other and jointly determine damage process of PBX. It should be noted that the fracture behavior of PBX subjected to impact loading has been widely investigated in previously published literature [8,40,41]. For example, the microscopic examinations indicate that various types of impact damage including trans-granular failure, debonding, and crystal cleavage can be observed, depending on loading conditions [8,40]. In contrast, our current simulation results also reproduce colorful fracture patterns such as inter-granular, transgranular, and interfacial failure, depending on the simulation conditions. Qualitatively speaking, our current simulation results are in agreement with previously reported experimental observations, indicating the effectiveness of the use of peridynamic theory for simulating complex dynamic response of PBX. Furthermore, our current numerical simulations actually reveal the coupling effects on facture patterns of interfacial properties and loading conditions. Simulation results indicate that the competitive mechanism between interfacial properties and loading conditions can be quantitatively expressed by using damage mode index. This should be an interesting topic deserving further exploration using experimental techniques in the future. 5. Conclusions This manuscript has presented an exploration of mechanisms leading to dynamic damage of PBX at impact loading via peridynamic modeling 7
Computational Materials Science xxx (xxxx) xxxx
X. Deng and B. Wang
CRediT authorship contribution statement
[10] J.G. Bennett, K.S. Haberman, J.N. Johnson, B.W. Asay, J. Mech. Phys. Solids 46 (1998) 2303–2322. [11] R. Liu, P.W. Chen, Mech. Mater. 124 (2018) 106–117. [12] M. Chatti, A. Frachon, M. Gratton, M. Caliez, D. Picart, N.A. Hocine, Int. J. Solids Struct. 168 (2019) 13–25. [13] N. Moës, J. Dolbow, T. Belytschko, Int. J. Numer. Meth. Eng. 46 (1999) 131–150. [14] M.M. Rashid, Comput. Methods Appl. Mech. Eng. 154 (1998) 133–150. [15] P.O. Bouchard, F. Bay, Y. Chastel, I. Tovena, Comput. Methods Appl. Mech. Eng. 189 (2000) 723–742. [16] S.A. Silling, Reformulation of Elasticity Theory for Discontinuities and Long-Range Forces, 2000. [17] S.A. Silling, M. Epton, O. Weckner, J. Xu, E. Askari, J. Elast. 88 (2007) 151–184. [18] S.A. Silling, E. Askari, Comput. Struct. 83 (2005) 1526–1535. [19] Y.D. Ha, F. Bobaru, Int. J. Fract. 162 (2010) 229–244. [20] F. Bobaru, G. Zhang, Int. J. Fract. 196 (2015) 59–98. [21] Y.D. Ha, F. Bobaru, Eng. Fract. Mech. 78 (2011) 1156–1168. [22] S.A. Silling, M.L. Parks, J.R. Kamm, O. Weckner, M. Rassaian, Int. J. Impact Eng. 107 (2017) 47–57. [23] C. Diyaroglu, E. Oterkus, E. Madenci, T. Rabczuk, A. Siddiq, Compos. Struct. 144 (2016) 14–23. [24] Q. Spreiter, M. Walter, J. Comput. Phys. 152 (1999) 102–119. [25] H. Ren, X. Zhuang, Y. Cai, T. Rabczuk, Int. J. Numer. Meth. Eng. 108 (2016) 1451–1476. [26] F. Mossaiby, A. Shojaei, M. Zaccariotto, U. Galvanetto, Comput. Math. Appl. 74 (2017) 1856–1870. [27] D. Dipasquale, M. Zaccariotto, U. Galvanetto, Int. J. Fract. 190 (2014) 1–22. [28] N. Prakash, G.D. Seidel, Modell. Simul. Mater. Sci. Eng. 26 (2017) 015003. [29] N. Prakash, G.D. Seidel, Eng. Fract. Mech. 177 (2017) 180–202. [30] X. Cao, X. Duan, C. Pei, Cryst. Res. Technol. 48 (2013) 29–37. [31] L. Borne, J.-C. Patedoye, C. Spyckerelle, Propellants Explos. Pyrotech. 24 (1999) 255–259. [32] R.M. Doherty, D.S. Watt, Propellants Explos. Pyrotech. 33 (2008) 4–13. [33] C. Hua, P.-J. Zhang, X.-J. Lu, M. Huang, B. Dai, H. Fu, Propellants Explos. Pyrotech. 38 (2013) 775–780. [34] J. Lemaitre, R. Desmorat, Engineering Damage Mechanics. Ductile, Creep, Fatigue and Brittle Failure, 2005. [35] R.H. Gee, A. Maiti, S. Bastea, L.E. Fried, Macromolecules 40 (2007) 3422–3428. [36] J.D. Yeager, A.M. Dattelbaum, E.B. Orler, D.F. Bahr, D.M. Dattelbaum, J. Colloid Interface Sci. 352 (2010) 535–541. [37] M. Li, J. Zhang, C.-Y. Xiong, J. Fang, J.M. Li, Y. Hao, Damage and fracture prediction of plastic-bonded explosive by digital image correlation processing, 2005. [38] P. Chen, H. Xie, F. Huang, T. Huang, Y. Ding, Polym. Test. 25 (2006) 333–341. [39] L. Chen, D. Han, S.-L. Bai, F. Zhao, J.-K. Chen, Mech. Adv. Mater. Struct. 24 (2017) 737–744. [40] P. Chen, F. Huang, Y. Ding, J. Mater. Sci. 42 (2007) 5272–5280. [41] Z.W. Liu, H.M. Xie, K.X. Li, P.W. Chen, F.L. Huang, Polym. Test. 28 (2009) 627–635.
Xiaoliang Deng: Conceptualization, Data curation, Funding acquisition, Software, Investigation, Writing - original draft. Bo Wang: Methodology, Writing - review & editing. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements The authors gratefully acknowledge financial support from the National Natural Science Foundation of China through grant 11872345 and the National Key Laboratory of Shock Wave and Detonation Physics, Institute of Fluid Physics, China Academy of Engineering Physics through grants 6142A0305010717 and JCKYS2018212011. Data availability The raw/processed data are available from the corresponding author on reasonable request. References [1] A.S. Kumar, V.B. Rao, R.K. Sinha, A.S. Rao, Propellants Explos. Pyrotech. 35 (2010) 359–364. [2] J.E. Field, Acc. Chem. Res. 25 (1992) 489–496. [3] H.L. Berghout, S.F. Son, C.B. Skidmore, D.J. Idar, B.W. Asay, Thermochim. Acta 384 (2002) 261–277. [4] B.W. Asay, in, Springer-Verlag, Berlin ;, 2010. [5] Y. Xiao, Y. Sun, X. Li, Q. Zhang, Propellants Explos. Pyrotech. 42 (2017) 873–882. [6] X. Wang, Y. Wu, F. Huang, Int. J. Solids Struct. 129 (2017) 28–39. [7] D.A. LaBarbera, M.A. Zikry, Comput. Mater. Sci. 104 (2015) 10–22. [8] Y. Xiao, Y. Sun, Y. Zhen, L. Guo, L. Yao, Int. J. Impact Eng. 103 (2017) 149–158. [9] X. Wang, Y. Wu, F. Huang, T. Jiao, R.J. Clifton, Mech. Mater. 99 (2016) 68–78.
8
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Peridynamic modeling of dynamic damage of polymer bonded explosive Xiaoliang Deng , Bo Wang ⁎
National Key Laboratory of Shock Wave and Detonation Physics, Institute of Fluid Physics, China Academy of Engineering Physics, Mianyang 621999, Sichuan, PR China
ARTICLE INFO
ABSTRACT
Keywords: Dynamic damage Polymer bonded explosives Peridynamic Cracks Impact loading
Dynamic damage response of 1,3,5,7-tetranitro-1,3,5,7-tetrazocane (HMX) based polymer bonded explosive (PBX) is investigated by using recently developed three dimensional Peridynamic (PD) method, allowing natural evolution of defects such as cracks or voids. An in-house PD codes are developed and verified by comparing with published results on Kalthoff-Winkler (KW) and crack branching benchmark tests. Numerical modeling of PBX is constructed via random micro-modulus approach, capable of capturing inherent microstructural heterogeneity of PBX materials. The simulation results indicate that damage process depends on both impact speed and interfacial property. Damage mode index (DMI), defined as the ratio of HMX-HMX damage to the summation of Binder-Binder and HMX-Binder damage, is proposed to quantify competitive damage mechanisms between trans- and inter-granular damage. The analyses of DMI exhibit that the inter-granular damage is dominant for the low adhesive strength of interface. However, the trend of trans-granular damage increases with increasing impact speed. The trans-granular damage mode becomes dominant for the high adhesive strength of interface. But the trend of inter-granular model increases with increasing impact speed. The current research highlight the competitive mechanisms of two typical damage modes as a function of material property and external loading conditions in a quantitative manner, providing additional insight into the damage mechanisms of PBX under impact loading.
1. Introduction Polymer bonded explosive (PBX) is one of solid high explosives and widely encountered in a variety of applications such as rocket propellants and main explosive charge in different weapon systems. PBX is composed of energetic materials particles (90–95% by weight) such as HMX and RDX embedded in polymeric binder (5–10% by weight) as well as additives [1]. PBX can experience some unexpected stimuli such as low-velocity impact during the manufacture, storage, transportation, and handling process. Those stimuli can lead to unwanted ignition of PBX due to the development of localized regions of high temperature called as the hot-spots inside the energetic materials. Suggested mechanisms of hot spot include adiabatic compression of trapped gases, hydrodynamic pore collapse, and localized shearing deformation and so on [2]. In impact and cook-off experiments on PBX, complex crack patterns are observed and reactive wave propagates along the cracks. On one hand, cracks can enhance chemical reaction by providing additional burning surface area. On the other hand, product gases produced from chemical reaction are able to increase local pressure and drive cracks propagation thus cause additional damage to PBX [3,4]. Overall, the mechanical damage and chemical reaction interplay to
⁎
determine the eventual response of PBX under impact loading. A significant improvements in predicating coupled mechanical and chemical response of PBX can result from a better understanding of damage evolution of PBX. Therefore, it is fundamentally important to investigate dynamic mechanical properties of PBX in order to deepen understanding the issues associated with non-shock ignition of PBX [4]. Over the years, Numerical simulations play an important role in investigating dynamic mechanical response of PBX. The previously published literature majorly utilize the finite element method (FEM) to simulate dynamic response of PBX [5–7]. For example, Xiao et al., predicts the stress–strain relationship of PBX under impact loading via implementing constitute including damage in ABAQUS. The simulation results are in good agreement with their experimental data [8]. Wang et al., perform mesoscale simulations of thermal-mechanical response of HMX based PBX and granular explosive under impact loading. The useful information such as equivalent strain and temperature distribution are achieved [9]. Regarding to FEM simulation, the accurate constitutive relation is required in order to simulate the complex dynamic behaviors of PBX. The Visco-Scram (viscoelastic statistical crack mechanics) model, including the viscoelasticity response and interior micro-cracks evolution, has been widely applied to describe the
Corresponding author. E-mail address: [email protected] (X. Deng).
https://doi.org/10.1016/j.commatsci.2019.109405 Received 31 August 2019; Received in revised form 17 October 2019; Accepted 10 November 2019 0927-0256/ © 2019 Elsevier B.V. All rights reserved.
Please cite this article as: Xiaoliang Deng and Bo Wang, Computational Materials Science, https://doi.org/10.1016/j.commatsci.2019.109405
Computational Materials Science xxx (xxxx) xxxx
X. Deng and B. Wang
dynamic damage response of PBX [8,10,11]. Recently, a novel viscoelastic Microplane model was proposed to simulate behavior of quasibrittle PBX materials in literature [12]. The features of this model are that it considers the viscoelastic behavior along with the damage induced anisotropy. Those simulation results can help better understand behaviors of mechanical damage and ignition. Although FEM simulation methods are currently the most popular numerical modeling to examine dynamic response of PBX in previously published literature. It should be noted that the mesh-based simulation method cannot efficiently handle non-continuity induced by the evolution of cracks and voids and so on because the governing equations are the spatial differential expressions. In order to treat problems associated with noncontinuity, special techniques are usually required to modify the traditional FEM models [13]. For example, the computation meshes can be regenerated during the simulation after detecting the new surfaces created by cracks propagation or voids evolution. This technique can ensure the discontinuous surface develops along the mesh lines but require to track new surfaces continuously inside the simulation domain as simulation time goes on, which suffers from time-intensive operations and also complex algorithms [14,15]. Therefore, it would be of considerable value to apply novel simulation methods which are capable of simulating damage evolution of PBX under impact loading in a more natural and efficient manner. Peridynamic (PD) theory is a new nonlocal theory recently proposed by Silling and his coworkers [16,17]. In contrast to the FEM model, governing equation in PD is integral operations rather than spatial differential equations. Therefore, it avoids the assumptions on differentiability and continuity of displacement field, allowing spontaneous crack nucleation, propagation, and interaction. This advantage makes PD a suitable candidate for simulating phenomena such as damage and fracture of materials under dynamic loading, where discontinuities such as voids or cracks usually nucleate and evolve in the damage areas of inside materials [18–20]. PD modeling has been widely used to investigate mechanical properties in a wide range of materials such as glass, metal, and composite materials and so on. Those previous works achieve the proposing results and also prompt the further development of PD theory [21–23]. Therefore, it is of particular interest to apply PD modeling investigating the dynamic response of PBX. Random micro-modulus approach is proposed to construct simulation model for the sake of simplicity, while capturing the inherent heterogeneity of PBX. The simulation results reveal typical damage modes of PBX observed in experiments, indicating the applicability of PD in terms of the investigation of dynamic damage response of PBX. The damage behavior at various impact speeds and interface parameters are systematically investigated. Damage mode index is proposed to quantify different types of damage of PBX. The results can potentially contribute to a deeper understanding of the dynamic behavior of energetic materials.
which is usually taken to be a sphere of radiusδcentered at point x in 3 dimensional model. Establishing the pair force function is a fundamental problem in peridynamics because it contains all the material constitutive information. For a prototype microelastic brittle material (PBM), f can take the following form:
f = cs [1
s=
c=
f (u (x , t ), u' (x ', t ), x , x ' , t ) dVx ' + b (x , t )
|| + || || || || ||
(2) (3)
12E
(4)
4
In bond-based peridynamics, a bond will break and can no longer carry any force if its elongation is greater than a critical relative elongation s0 , and the damage function is defined as
µ (t , ) =
0, ifs < s0 1, otherwise
(5)
Silling and Askari correlate the s0 with the fracture energy G0 of materials which represents the energy per unit fracture area in three dimensions required to terminate all bonds initially crossing the fracture surface [18]. The critical stretch s0 is described as 5G0 , depending 6E on both material properties and the horizon. In current model, bond failure is irreversible, i.e. broken bond cannot reconnect during the simulation. Based on the damage function, the level of damage at each material point can be expressed as
(x , t ) = 1
H
µ (t , ) dV ' H
dV '
(6)
which is actually the ratio of amount of broken bonds to the total amount of bonds initially associate with point x. The simulation domain is discretized into cubic lattice with the lattice constant of x and the volume of x 3 . According to the procedures described above, the interaction forces acting on all material points can be computed at each time step and an explicit VelocityVerlet algorithm is used for time integration to yield the motion of all of the individual material points [24]:
un + 1/2 = un +
Peridynamics regard the materials as the assembly of large amount of material points with mass and volume. Each point (x) can interact with all other points (x′) within a limited distance called horizon (H) through the pairwise force function f. The physical meaning of force function is the force vector per unit volume squared that neighboring point x′ exerts on the point x. In Peridynamics, this direct physical interaction between materials points is called bond. The acceleration of point x can be obtained according to the Newton’s second law of motion [18]: H
+ || + ||
where c, s, μ represent the micro-modulus, elongation of peridynamics bond, and damage function, respectively. = x ' x and = u' u denote the relative position vector of two material points in the reference configuration and their relative displacement, respectively. As a result, + represents the current relative position in the deformed configuration. Micro-modulus can be derived by strain energy density equivalence to a classical continuum theory of elasticity for the same materials and the same deformation. In 3D, c can be obtained as a function of Young’s modulus E and the radius of horizon [18]:
2. Bond-based Peridynamic formulation
u¨ (x , t ) =
µ (t , )]
t n u¨ 2
un + 1 = un + tun + 1/2 un + 1 = un + 1/2 +
t n+ 1 u¨ 2
(7) (8) (9)
where t is the time step. 3. Verification of codes On the basis of PD theory depicted above, the PD computation codes are developed in our laboratory. It is of importance to validate the codes before applying it to simulate dynamic response of PBX. Therefore, two benchmark examples with obtained results are presented and compared with previous results from published literature in this section.
(1)
where ρ, u, V, b represent the mass density, displacement vector field (i.e. u¨ is its current acceleration vector field), infinitesimal volume of material point, and prescribed body force density field. H is the horizon 2
Computational Materials Science xxx (xxxx) xxxx
X. Deng and B. Wang
Fig. 1. Simulation of Kalthoff–Winkler’s experiment. (a) Geometry of computational model, (b) Simulation result on crack propagation.
3.1. Example 1: Kalthoff–Winkler experiment
observation and previous simulation results. Both the results on KW and crack branching benchmark tests validate the developed PD codes.
The Kalthoff-Winkler experiment is a classic benchmark problem for studying dynamic fracture and peridynamics dynamics simulation results on this problem are reported in many literature [25,26]. The simulation model corresponding to the experimental setup is presented in Fig. 1. The thickness of the plate is set to be 0.01 m. The experimental results indicate that if a pre-cracked steel plate made of steel 18Ni1900 is hit literally with impact speed of 32 m/s, a brittle fracture mode is observed and the cracks initiates at the end of the pre-existing cracks and then propagate along around 70 degree with respect to the preexisting cracks. In order to compare with previous published results, the same initial simulation parameters are used, i.e. Young’s Modulus E = 190 GPa, density ρ = 7800 kg/m3, Poisson’s ratio ν = 0.25, and critical energy release rate of the material G0 = 6.9 × 104 J/m2. The first three layers of grids (red color in Fig. 1(a)) are applied initial velocity of 22 m/s to simulate impact loading for the sake of simplicity. Simulation domain is divided with uniform grid of x = 1 mm which leads to totally 200,000 grids in the simulation. The radius of each material point’s horizon is about three times of the grid spacing e.g., = 3.015 x [18].. The simulation results show the cracks start to initiate at about 20 µs and then propagate at the angle of 67 degree with respect to the initial direction of pre-existing cracks. Those results are in good agreement with previously published simulation results [26].
4. Modeling of dynamic damage of PBX 4.1. Problem setup and computational parameters PD model of PBX used in this research is shown in Fig. 3. A rectangular sample is constructed with dimensions along x, y, and z are 4.5 mm, 1.0 mm, and 1.0 mm, respectively. There are total 50 grains in the rectangular sample and the average grain size is therefore about 440 µm . The computational domain is divided into three-dimensional rectangular grids with total number of about 220,000. The three layers of grid along x direction are fixed as a piston moving at a constant speed during the simulation. The free boundary conditions are used in simulations. The time increment △t is set to be 1.0e-10s and the duration of simulation is 3 µs . The materials used in numerical simulations are HMX with densityρ = 1900 kg/m3, Young’s modulus E = 24.0 GPa, critical relative elongation s0 = 0.0005 [28] and binder with densityρ = 1100 kg/m3, Young’s modulus E = 3.0 GPa, critical relative elongation s0 = 0.02 [29]. Voronoi tessellation method is used to generate PBX polycrystalline structure. In this method, some seed points are first generated in a random manner, then the distances between material points and seed points are computed. The geometries of HMX crystals are determined according to distances. It should be noted that the previous experimental investigation reveals that it is usually inevitable to introduce damage and defects such as micro-voids and dislocations within the HMX crystal during the manufacturing process of PBX [30–33]. Those defects distributed inside the PBX in a random manner. As a result, it is very hard, if it is not completely impossible, to construct a computation model directly based on the real microstructures and distributions of those defects due to the current lack of detailed experimental characterization of PBX. However, it is well-known that the amount and distribution of defects can impact material strength [34]. Therefore, more defects inside a specific HMX crystal can lead to a less strength of this HMX crystal and vice visa. In PD theory, the material strength is directly associated with micro-modulus. Therefore, a random micromodulus model is proposed in this manuscript to represent intrinsic heterogeneous microstructures of PBX for the sake of simplicity. In this model, micro-modulus in the range of [0.5c,1.5c] are randomly distributed to each HXM grain, where c is computed according to Eq. (4). Such micro-modulus distributions essentially represent the effect of random amount and distribution of defects on the strength of HMX crystal. From Fig. 3, it can be observed that the PD model of PBX consists of HMX and binder particles. These two types of PD particles can constitute three types of bonds: HMX-HMX, Binder-Binder, and HMXBinder, which actually represent the properties of HMX crystal, binder, and the interface between HMX crystal and binder, respectively. It
3.2. Example 2 crack branching in a pre-cracked plate An initial crack is prescribed in a thin rectangular plate with dimensions shown in Fig. 2(a). The upper and lower three layers of grids along the vertical direction are suddenly applied a tensile loading of 12 MPa, then tensile stress keeps constant during the whole simulation process. This type of loading condition generates a stress wave, which can interact with dynamic propagation of cracks and eventually result in crack branching [19,21,27]. The brittle material used in this simulation is Duran 50 glass with densityρ = 2235 kg/m3, Young’s modulus E = 65 GPa, Poisson’s ratioν = 0.2, and critical energy release rate G0 = 204 J/m2. Simulation domain is divided with uniform grid of x = 0.4 mm in the simulation. The radius of each material point’s horizon is about three times of the grid spacing e.g., = 3.015 x . Peridynamic simulation results are shown in Fig. 2(b)–(d). Pre-existed crack first propagates horizontally for a while and then starts to branch. The pattern of crack branching shown in Fig. 2 is similar with previously published results [19]. Simulation results indicate that branching event occurs at about 24 microseconds and crack speed reaches its maximum value of about 1800 m/s. Maximum speed is equal to 0.56 CR (CR is Rayleigh wave speed). It has been experimentally shown that the limiting crack speed in glass ranges from 0.47 CR to 0.66 CR for different loading conditions [20]. Therefore, our simulation data about crack branching in glass is in good agreement with experimental 3
Computational Materials Science xxx (xxxx) xxxx
X. Deng and B. Wang
Fig. 2. Simulation results about crack branching study at various time steps. Different colors indicate different damage levels. (a) 0 μs, (b) 20 μs, (c) 25 μs, (d) 40 μs.
should be noted that interfacial properties is hard to characterize experimentally. Previously published results in literature are inconsistent with each other. For example, the works of adhesion between TATB and Kel-F 800 interfaces are reported to be 76–271 mJ/m2 and 43 mJ/m2 [35,36]. In order to systematically investigate the effect of interfacial properties between HMX crystal and binder on dynamic damage response of PBX, an adjustable interface parameterβis introduced in this research. β is defined as the following equation:
PHMX
Binder
= Pmin + (Pmax
Pmin)
4.2. Dynamic damage of PBX at various impact speeds Dynamic response of PBX at various impact speeds is investigated in this section. According to previously published literature, various damage types such as trans-granular and inter-granular cracks under impact loading can be observed [37–39]. Trans-granular cracks usually describe the damage of individual bulk crystal. In the contrast, intergranular fracture refers to the damage along the interfaces between HMX and binder. In order to quantify the damage degree of HMX, binder and their interfaces, the broken bonds number of HMX-HMX, Binder-Binder, and HMX-Binder normalized by their own initial values are computed and presented. The quantitative characterization of damage is important because it can help assess the mechanical integrity of PBX under impact loading. The results concerning the damage of HMX-HMX, Binder-Binder, and HMX-Binder bonds at various impact speeds are presented in Fig. 4. Fig. 4(a)–(d) correspond to the impact speed of 20 m/s, 30 m/s, 40 m/s, and 60 m/s. = 0.5 is used for all the different impact speeds. The minor damage is observed at impact speed of 20 m/s. With increasing impact speed, damage becomes more and more serious. More precisely, HMX-HMX damage increases from about 6% to 50% with increasing
(10)
where PHMX Binder represents critical relative elongation of HMX-Binder bonds and Pmin, Pmax represent minimum and maximum values critical relative elongation between HMX crystal and binder, respectively. According to Eq. (10), critical relative elongation achieves minimum or maximum value corresponding to = 0 or = 1, respectively. When = 0.5, interfacial parameter is equal to the average value of HMX and binder. Therefore, the impact of interfacial mechanical properties on dynamic damage of PBX can be investigated via controllable change of interface parameter , illustrating the advantage of employing numerical computation for this type of systematic investigations.
Fig. 3. Peridynamic model of PBX. Blue and red particles represent binder and HMX, respectively. Three different types of bonds, including HMX-HMX, HMX-Binder and Binder-Binder bonds, formed in the PD model are also shown in the picture. 4
Computational Materials Science xxx (xxxx) xxxx
X. Deng and B. Wang
Fig. 4. Influence of different impact speeds on dynamic damage of PBX. (a) 20 m/s, (b) 30 m/s, (c) 40 m/s, (d) 60 m/s. Three different types of damage, e.g. HMXHMX, Binder-Binder, and HMX-Binder are shown, which are defined as the ratio of their broken bond numbers to their own initial values.
Fig. 5. Snapshots of simulation results at different impact speeds. (a) 20 m/s, (b) 30 m/s, (c) 40 m/s.
impact speed from 20 m/s to 60 m/s. However, qualitative trend is similar for all different impact speeds, i.e., the damage of HMX-HMX is dominant in all three types of damage modes, regardless of impact speed. This conclusion is also further verified by visualized examination of simulation results (Fig. 5). It can be observed that most of damage appears inside the HMX crystal and damage of binder material is relatively low over an impact speed range from 20 to 40 m/s. But the question remains, whether dominant damage are always the damage of HMX-HMX. Therefore, addition simulations concerning of effect of material property on dynamic response are conducted in next section of this manuscript.
4.3. Dynamic damage response of PBX at various interface parameters Dynamic damage of PBX at various impact speeds shows that the damage in PBX mostly occur inside the HMX crystals. This phenomenon is somehow inconsistent with experimental observations. Experimental scanning electron microscopy (SEM) results indicate damage along HMX-binder interfaces can also be observed [8,40]. In order to systematically understand dynamic damage of PBX, the influence of interface parameter of β on dynamic damage is evaluated in this section. Fig. 6 presents the impact of interfacial parameters on dynamic damage of PBX under various impact velocities. It indicates that 5
Computational Materials Science xxx (xxxx) xxxx
X. Deng and B. Wang
Fig. 6. Influence of interfacial parameters on dynamic damage of PBX. (a) β = 0,(a) β = 0.5,(a) β = 1.
shown in Fig. 3, the breakage of HMX-HMX bonds is associated with the trans-granular cracks, referring to this damage mode as trans-granular damage in this manuscript. Similarly, the breakage of Binder-Binder bonds actually represent the binder material fracture and the breakage of HMX-Binder bonds is associated with the interface debonding. However, due to the relative low weight ratio of binder materials in PBX, the dimension of interface between HMX crystal and binder material is usually small compared to that of HMX crystal. Therefore, it is actually hard to experimentally distinguish between binder material fractures and interface debonding. The general experimental phenomena resulting from these two damage modes is that the cracks propagate along the interface between HMX crystal and binder material, referring to these two damage modes as inter-granular damage for the sake of convenience. Damage mode index (DMI), defined as the ratio of HMX-HMX damage to the summation of Binder-Binder and HMX-Binder damage, is proposed in this manuscript. Introduction of DMI can help quantitatively analyze the relative trend of trans-granular and inter-granular damage of PBX. According to its definition, DMI is greater than one if the trans-granular damage predominates over the inter-granular damage. Otherwise, DMI is less than one if inter-granular damage dominates over trans-granular damage. Based on the simulation data, the DMI at various interfacial parameters and impact speeds can be computed, which is shown in Fig. 7. Fig. 7(a), corresponding to β = 0, shows that DMI is always less than one for all impact speeds, indicating inter-granular damage is dominant. However, it can be observed that the value of DMI increases with increasing impact speed, indicating the trend of trans-granular damage increases with increasing impact speed. Fig. 7(b) and (c), corresponding to β = 0.5 and β = 1, respectively, show that the DMI is always greater than one for all impact speeds, indicating trans-granular damage is dominant. However, the specific value of DMI decreases with increasing impact speed, indicating the trend of inter-granular damage increases with increasing impact speed. Therefore, damage mode of PBX actually depends on both intrinsic material property and external loading condition. Both of them
different interfacial parameters lead to different damage modes. For β = 0, corresponding to the weak interface strength, the dominant damage is HMX-Binder. With increasing value ofβ, the dominant damage mode transfer from HMX-Binder to HMX-HMX. This phenomenon is understandable because the smaller value of β leads to less critical elongation of bonds formed by HMX and binder PD particles. Therefore, the interfacial bonds of HMX and binder can be easily broken. With increasing value ofβ, interfacial bonds of HMX and binder become harder and harder to be broken due to increasing critical elongation. The visualization of simulation results showing in the three images at the bottom of Fig. 6 is consistent with quantitative results about damage of HMX-HMX, Binder-Binder, and HMX-Binder. For example, the heavy damage contour of interface between HMX and binder is observed in the bottom image of Fig. 6(a), where the damage occurs predominantly at the interface. At the same time, quantitative results in the top image of Fig. 6(a) indicate the damage mode of HMX-Binder is dominant. The thin slice of PBX located in the center of simulation domain is also cut from simulation domain in order to achieve detailed fracture information, which is not shown for the sake of brevity. The cracks propagating along interfaces between HMX and binder are observed, indicating the inter-crystal fracture mode. Moreover, damage modes for β = 0.5 and β = 1 are similar. Both of them indicate the damage of HMX-HMX is predominant. Therefore, it can be concluded that if interfacial strength is high enough, it no longer influences damage behaviors of PBX in the range of impact speed considered in current research. Another observation is that the damage of Binder-Binder also decreases with increasing value ofβ. This result indicates that the strong interface strength between binder particles and HMX particles can also prevent bond break of Binder-Binder. 4.4. Damage mode index (DMI) The above simulation results indicate the possible damage behaviors of PBX at PD particle level include the breakage of HMX-HMX, Binder-Binder, and HMX-Binder bonds. According to the schematic 6
Computational Materials Science xxx (xxxx) xxxx
X. Deng and B. Wang
Fig. 7. Damage mode index (DMI) at various interfacial parameters as a function of impact speed. (a) β = 0, (b) β = 0.5, (c) β = 1.
which allows spontaneous formation and evolution of cracks in a consistent framework. The in-house PD codes are developed and verified via two classic examples: KW test and crack branching under tension loading. A random micro-modulus model is proposed to represent inherent heterogeneity of PBX induced by various amounts and distributions of defects inside the materials. PD simulation results exhibit that damage degree increases with increasing impact speed. Quantitative analysis shows for the fixed interface parameter β = 0.5 breakage number of HMX-HMX bonds are greater than that of Binder-Binder and HMX-Binder bonds. When interface parameter β decreases from 0.5 to 0, the breakage number of HMX-Binder bonds is greater than that of HMXHMX and Binder-Binder bonds. In order to systematically analyze damage patterns for various impact speeds and material properties, damage mode index (DMI), defined as the ratio of HMX-HMX damage to the summation of Binder-Binder and HMX-Binder damage, is proposed and computed based on the simulation data. The DMI analyses reveal that inter-granular damage is the dominant damage mode for different impact speeds for β = 0.0. However, trend of trans-granular damage increases with increasing impact speed. On the other hand, trans-granular damage is dominant for different impact speeds for the interface parameter β = 0.5 and β = 1.0. However, trend of inter-granular damage increases with increasing impact speed. The observed simulation results highlight the combined effects of material properties and loading conditions on dynamic damage of PBX, helping better understand dynamic mechanical behaviors of heterogeneous energetic materials.
interplay with each other and jointly determine damage process of PBX. It should be noted that the fracture behavior of PBX subjected to impact loading has been widely investigated in previously published literature [8,40,41]. For example, the microscopic examinations indicate that various types of impact damage including trans-granular failure, debonding, and crystal cleavage can be observed, depending on loading conditions [8,40]. In contrast, our current simulation results also reproduce colorful fracture patterns such as inter-granular, transgranular, and interfacial failure, depending on the simulation conditions. Qualitatively speaking, our current simulation results are in agreement with previously reported experimental observations, indicating the effectiveness of the use of peridynamic theory for simulating complex dynamic response of PBX. Furthermore, our current numerical simulations actually reveal the coupling effects on facture patterns of interfacial properties and loading conditions. Simulation results indicate that the competitive mechanism between interfacial properties and loading conditions can be quantitatively expressed by using damage mode index. This should be an interesting topic deserving further exploration using experimental techniques in the future. 5. Conclusions This manuscript has presented an exploration of mechanisms leading to dynamic damage of PBX at impact loading via peridynamic modeling 7
Computational Materials Science xxx (xxxx) xxxx
X. Deng and B. Wang
CRediT authorship contribution statement
[10] J.G. Bennett, K.S. Haberman, J.N. Johnson, B.W. Asay, J. Mech. Phys. Solids 46 (1998) 2303–2322. [11] R. Liu, P.W. Chen, Mech. Mater. 124 (2018) 106–117. [12] M. Chatti, A. Frachon, M. Gratton, M. Caliez, D. Picart, N.A. Hocine, Int. J. Solids Struct. 168 (2019) 13–25. [13] N. Moës, J. Dolbow, T. Belytschko, Int. J. Numer. Meth. Eng. 46 (1999) 131–150. [14] M.M. Rashid, Comput. Methods Appl. Mech. Eng. 154 (1998) 133–150. [15] P.O. Bouchard, F. Bay, Y. Chastel, I. Tovena, Comput. Methods Appl. Mech. Eng. 189 (2000) 723–742. [16] S.A. Silling, Reformulation of Elasticity Theory for Discontinuities and Long-Range Forces, 2000. [17] S.A. Silling, M. Epton, O. Weckner, J. Xu, E. Askari, J. Elast. 88 (2007) 151–184. [18] S.A. Silling, E. Askari, Comput. Struct. 83 (2005) 1526–1535. [19] Y.D. Ha, F. Bobaru, Int. J. Fract. 162 (2010) 229–244. [20] F. Bobaru, G. Zhang, Int. J. Fract. 196 (2015) 59–98. [21] Y.D. Ha, F. Bobaru, Eng. Fract. Mech. 78 (2011) 1156–1168. [22] S.A. Silling, M.L. Parks, J.R. Kamm, O. Weckner, M. Rassaian, Int. J. Impact Eng. 107 (2017) 47–57. [23] C. Diyaroglu, E. Oterkus, E. Madenci, T. Rabczuk, A. Siddiq, Compos. Struct. 144 (2016) 14–23. [24] Q. Spreiter, M. Walter, J. Comput. Phys. 152 (1999) 102–119. [25] H. Ren, X. Zhuang, Y. Cai, T. Rabczuk, Int. J. Numer. Meth. Eng. 108 (2016) 1451–1476. [26] F. Mossaiby, A. Shojaei, M. Zaccariotto, U. Galvanetto, Comput. Math. Appl. 74 (2017) 1856–1870. [27] D. Dipasquale, M. Zaccariotto, U. Galvanetto, Int. J. Fract. 190 (2014) 1–22. [28] N. Prakash, G.D. Seidel, Modell. Simul. Mater. Sci. Eng. 26 (2017) 015003. [29] N. Prakash, G.D. Seidel, Eng. Fract. Mech. 177 (2017) 180–202. [30] X. Cao, X. Duan, C. Pei, Cryst. Res. Technol. 48 (2013) 29–37. [31] L. Borne, J.-C. Patedoye, C. Spyckerelle, Propellants Explos. Pyrotech. 24 (1999) 255–259. [32] R.M. Doherty, D.S. Watt, Propellants Explos. Pyrotech. 33 (2008) 4–13. [33] C. Hua, P.-J. Zhang, X.-J. Lu, M. Huang, B. Dai, H. Fu, Propellants Explos. Pyrotech. 38 (2013) 775–780. [34] J. Lemaitre, R. Desmorat, Engineering Damage Mechanics. Ductile, Creep, Fatigue and Brittle Failure, 2005. [35] R.H. Gee, A. Maiti, S. Bastea, L.E. Fried, Macromolecules 40 (2007) 3422–3428. [36] J.D. Yeager, A.M. Dattelbaum, E.B. Orler, D.F. Bahr, D.M. Dattelbaum, J. Colloid Interface Sci. 352 (2010) 535–541. [37] M. Li, J. Zhang, C.-Y. Xiong, J. Fang, J.M. Li, Y. Hao, Damage and fracture prediction of plastic-bonded explosive by digital image correlation processing, 2005. [38] P. Chen, H. Xie, F. Huang, T. Huang, Y. Ding, Polym. Test. 25 (2006) 333–341. [39] L. Chen, D. Han, S.-L. Bai, F. Zhao, J.-K. Chen, Mech. Adv. Mater. Struct. 24 (2017) 737–744. [40] P. Chen, F. Huang, Y. Ding, J. Mater. Sci. 42 (2007) 5272–5280. [41] Z.W. Liu, H.M. Xie, K.X. Li, P.W. Chen, F.L. Huang, Polym. Test. 28 (2009) 627–635.
Xiaoliang Deng: Conceptualization, Data curation, Funding acquisition, Software, Investigation, Writing - original draft. Bo Wang: Methodology, Writing - review & editing. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgements The authors gratefully acknowledge financial support from the National Natural Science Foundation of China through grant 11872345 and the National Key Laboratory of Shock Wave and Detonation Physics, Institute of Fluid Physics, China Academy of Engineering Physics through grants 6142A0305010717 and JCKYS2018212011. Data availability The raw/processed data are available from the corresponding author on reasonable request. References [1] A.S. Kumar, V.B. Rao, R.K. Sinha, A.S. Rao, Propellants Explos. Pyrotech. 35 (2010) 359–364. [2] J.E. Field, Acc. Chem. Res. 25 (1992) 489–496. [3] H.L. Berghout, S.F. Son, C.B. Skidmore, D.J. Idar, B.W. Asay, Thermochim. Acta 384 (2002) 261–277. [4] B.W. Asay, in, Springer-Verlag, Berlin ;, 2010. [5] Y. Xiao, Y. Sun, X. Li, Q. Zhang, Propellants Explos. Pyrotech. 42 (2017) 873–882. [6] X. Wang, Y. Wu, F. Huang, Int. J. Solids Struct. 129 (2017) 28–39. [7] D.A. LaBarbera, M.A. Zikry, Comput. Mater. Sci. 104 (2015) 10–22. [8] Y. Xiao, Y. Sun, Y. Zhen, L. Guo, L. Yao, Int. J. Impact Eng. 103 (2017) 149–158. [9] X. Wang, Y. Wu, F. Huang, T. Jiao, R.J. Clifton, Mech. Mater. 99 (2016) 68–78.
8