Numerical Simulation of Crack Propagation of Architecture Glass Under Explosion and Impact Load

: Crack initiation and propagation is a long-standing difficulty in solid mechanics, especially for elastic-brittle material. To explore the damage and crack propagation behavior of architectural glass under different type of loads, the element deletion (ED), discontinuous Galerkin peridynamics (DG-PD) and meshless peridynamics (M-PD) methods are studied. Taking the architecture glass as an example, the crack propagation behavior under the bullet impact and explosion load are studied. The JH-2 material model is used in the ED method, and the maximum principal stress and maximum principal strain failure criteria are applied at the same time. In the DG-PD method, it conducts a node separation operation and imposes the criterion of the critical energy release rate. The M-PD method adopts a self-programmed particle discretization method and imposes a criterion of critical elongation. Three methods can simulate the crack growth behavior of glass material, but the PD method has great advantages in detail, such as crack bifurcation and penetration. For low-velocity bullets, the failure behavior of glass all shows cross-shaped cracks in different methods. The splashing of elements or particles appears in the two PD methods, but the particle splashing of the M-PD method is more obvious, and the DG-PD method captures the crack bifurcation effect better. For the failure behavior of glass under explosive loading, the PD method is obviously better than the ED method in terms of modal appearance. However, in the mechanical behavior of specific elements, the two methods have a high degree of agreement.


Introduction
Glass has many excellent properties, such as high strength, light transmittance, low density and high temperature resistance. It has been widely used in civil construction, transportation vehicles and aerospace. In actual production and life, the most commonly used glass type is architectural soda lime glass, which is of great significance to study its behavior under dynamic loading [1][2][3]. In theory, Krauthammer and Altenberg [4] used equivalent thickness and equivalent stiffness to simplify the glass structure into an equivalent single degree-of-freedom mass-spring model. Yuan et al. [5] analyzed the response of a glass structure in air and showed that the response under load was decoupled into two deformation stages: elastic deformation of the glass and large plastic deformation of PVB. In experiments, Holmquist et al. [6] conducted a plate impact experiment of inorganic glass in a large strain rate range and obtained its state parameters. Daryadel et al. [7] researched an improved SHPB to study the dynamic compression behavior of four different inorganic glasses.
The initiation and propagation of cracks may cause fatal hazards to the engineering structure in architectural glass under impact loading. Therefore, a reasonable assessment is essential for cracks. However, in terms of numerical simulation, the problem of crack propagation is a long-standing difficulty and challenge in classical computational mechanics. Although the finite element method has great advantages in traditional stress analysis, the simulation of material failure and fracture is slightly insufficient. This deficiency stems from the partial differential equations in the governing equations of classical continuum mechanics theory. There will be no space derivatives of the partial differential equations on the crack tip and surface, resulting in the singularity of the crack tip stress [8][9][10]. To solve the above difficulties in numerical solutions, on the one hand, element deletion method with failure criteria [11], node separation method [12] and extended finite element method (XFEM) [13] and many other finite element numerical processing methods have been proposed. On the other hand, to avoid numerical instability caused by multiple cracks, such as peridynamics [9], the discrete element method (DEM) [14], smooth particle hydrodynamics (SPH) [15] and other meshless methods are also widely used. For example, Zang et al. [16] used the SPH and FEM coupling method to simulate the splashing process of glass fragments. Deng and Jin [17] used the multi-material ALE method to conduct large-scale characteristics of hollow glass curtain walls. Florin et al. [18], according to the experimental results of Bless [19], used M-PD to systematically study the damage of multilayer soda lime glass under high-speed impact loading. Ren [20] applied the three-dimensional DG-PD method in dynamic brittle failure analysis within bond-based PD simulation.
Many theoretical, experimental analysis and numerical simulation methods exist for glass structures and have achieved excellent results. However, the research on the numerical simulation of crack propagation is not comprehensive under different types of loads, and there have few studies comparing the simulation results of different methods. This present paper is organized as follows. In Section 2, the basic theories of ED and PD are introduced. In Section 3, the material model parameters of different methods are illustrated. In Section 4, the ED, DG-PD and M-PD methods are applied to study the crack growth behavior of architecture glass under bullet impact and explosion load. Finally, Section 5 summarizes the conclusions that resulted from this study.

Element Deletion Method(ED)
The ED method means that when the elements in the model meet the preset failure criteria, the program will automatically delete the corresponding elements and nodes. When a series of elements and nodes are deleted, a visible crack path will be formed. One of the most common methods is to add failure criteria such as the maximum principal stress, maximum principal strain or maximum tensile stress. Meeting any of the prefabricated failure criteria will cause damage to the material during the loading process. The other method is the material model with failure criteria such as *MAT_STRAIN_RATE_DEPENDENT_PLASTICITY. Although the ED method is simple to operate and has high computational efficiency, it cannot guarantee the conservation of the mass, energy and momentum of the system, which makes it difficult to accurately reflect the material properties truly [21].

Peridynamics Method(PD)
PD can be regarded as molecular dynamics in a macroscopic sense to a certain extent. It is suitable for simulating the fracture behavior of microelastic brittle materials such as glass, ceramics, and rocks. Since the equation of motion uses the integral equation of nodal force instead of the partial differential equation, this can avoid the problem of mathematical singularity when facing discontinuous surface segments. At the same time, there are some natural advantages in crack nucleation, propagation, penetration and arrest [8][9][10]. For the reference configuration at any time, the motion balance equation can be described by the following motion equation within the scope of a certain material point: and  is a finite distance called horizon. In the bond-based PD (B-PD) model, the material is selected as a linear, microelastic brittle material. The micropotential is the energy of a single bond, and the point force f is the derivative of the micropotential energy function with respect to the relative displacement.
The contact force is the contact between nodes in the classical continuum, which is determined by the local deformation gradient and the related constitutive material model. However, the PD force interacts within a limited distance and disappears when the point force is beyond this horizon.
The total local strain energy density of the point can be obtained by summing the elastic energy on all the bonds connected to the material point.
Considering the classic microelastic brittle material (PMB), the constitutive force function of the B-PD under the condition of small deformation can be written as: It is a function related to the loading history. When the elongation rate exceeds the limit value, there is no interaction between the point pairs. Related to the critical energy release rate 0 G in classical fracture mechanics under three-dimensional conditions can be written as: At the same time, the damage degree of a material point can be defined as: , ( (9) 0   represents the initial state of the material point, and 1   represents that the point has been completely destroyed. The process from 0 to 1 of damage degree represents the gradual softening of the material properties. Once the bond is broken, it will be removed forever. When a series of bonds are broken, a visible crack path will be formed.

Material Model
At present, the most widely used dynamic constitutive model for simulating brittle materials is the *MAT_JOHNSON_HOLMQUIST_CERAMICS(JH-2) model. The damage evolution process of this model can be divided into the following steps. At the beginning of loading, the material exhibits elastic properties. Then, the material begins to exhibit corresponding damage when the stress reaches the preset yield limit. Finally, the glass material will experience fragmentation when the damage reaches a certain degree. The material model parameters of the glass material JH-2 are shown in Table 1. In this study, the JH-2 material model and ED method are used for glass together. The maximum principal strain and maximum principal stress failure criteria are applied simultaneously to simulate crack propagation, and the values are set to 0.001 and 84MPa, respectively [22]. The constitutive model of B-PD needs to set the value of the critical energy release rate G t [23]. Moreover, the critical energy release rate is directly related to the critical elongation. When the elongation exceeds the critical value, the bond between the point pairs no longer interacts and cannot be restored. The key model parameters of the DG-PD method used are shown in Table 2. It is worth noting that HSFAC indicates the normalized support zone size. LS_DYNA will adjust HSFAC automatically to ensure the neighbors of a point. In present paper, the M-PD method adopts self-programming discretization for modeling. Theoretically, B-PD is adopted to simulate the damage and destruction behavior of glass by applying the critical elongation S 0 . Compared with the DG-PD method, M-PD does not need to separate corresponding elements. The particle discrete distance x  is 3mm, and the PD radius horizon is taken as The key model parameters of the M-PD method used are shown in Table 3.

Analysis of numerical simulation
For glass materials, the nucleation position and subsequent growth of cracks are the main forms of damage. Meanwhile, the structural composition and loading methods of the glass itself are the most important factors influencing the crack morphology. To study the dynamic behavior of architecture glass under bullet impact and explosion load, the ED and PD methods are used to carry out the following numerical simulations. Otherwise, glass is a microelastic brittle material, and the generated strain when it fails is very small, approximately 0.1%. To observe the process of crack evolution intuitively, the maximum principal strain cloud diagram is used in the ED method, the equivalent plastic strain is used in the DG-PD method, and the damage factor is used in the M-PD method.

Finite Element Model
A square glass plate with 300mm side length and 8mm thickness, having a rigid bullet impact at the center. The finite element models are divided into three conditions. First, the rigid bullet impacts the single-layer glass frontally at a speed of 10m/s, and the four sides of the glass are framed. The finite element model is shown in Figure 1(a). Second, flat glass is processed under the same conditions as in Figure 1(a), but the impact load is replaced by explosion load with blast distance of 0.1m and charge of 15g. The finite element model is shown in Figure 1(b). Third, the flat glass and bullet are discretized into meshless particles, and the velocity of the bullet is 10m/s. A frontal impact is applied to the single-layer glass, and the calculation model is shown in Figure 1(c). Due to the frame-supported boundary conditions will cause the calculation to not converge in the trial calculation process, the method of restricting the degree of freedom of the edge particles is adopted for the meshless particle method. To simulate the damage and crack growth behavior of glass, the ED and PD methods are used to perform relevant numerical simulations on the above calculation model.

Crack propagation damage evolution
Failure will occur under extremely small strain due to glass is an elastic-brittle material. The damage process in the ED method can be observed intuitively by using the maximum principal strain. Fig. 2(a) shows that there are obvious cross-shaped cracks in the ED method. Although the impacted part appears to fall off in a large area, there is no splash of glass fragments due to the limitation of its own algorithm. The difference from the simulation results of the 1/4 symmetry model in the literature [6] is that the cracks are not completely symmetrical, and the lengths of the cracks are slightly different. Fig. 2(b) shows the crack propagation evolution diagram based on the DG-PD method, which visually indicates the damage and failure process by equivalent plastic strain. First, not only has glass obvious cross-shaped cracks but also has a large number of intersecting and secondary cracks, which are difficult to achieve in classical continuum mechanics. In addition, the PD method also captures the damage at the edge of the frame and the splashing of particles at the center. Moreover, the DG-PD method has obvious advantages in many aspects, especially in capturing the details of crack propagation.
The damage and crack propagation evolution based on the M-PD method is indicated by the factor of damage in Fig. 2(c). The glass structure has obvious cross-shaped and X-shaped cracks, and the void area is significantly larger than that of the previous two methods at the impact point. Additionally, there are large number of splash particles in the impact region. Different from the previous two methods, the cracks generated by the M-PD method have higher symmetry.

Comparison of failure mode
The ED, DG-PD and M-PD methods are successively adopted for numerical simulation under the same loading condition, and the comparison diagram of the final failure mode is shown in Fig. 3. It can be seen that there are all cross-shaped cracks in the failure mode under different methods, and there exist obvious cavities at the impact region. In the experiment [12], the glass structure produces radial cracks and a cavity at the impact point. The whole damage mode of the experiment is similar to the previous simulation methods, but there are many differences in the specific damage situation.  Table 4 presents the comparison of key damage modes and calculated parameters under different methods and experiments. According to the comprehensive consideration of the aperture of the damaged area, the insertion depth of the bullet and the maximum propagation distance of the cracks, the M-PD method has the highest degree of damage, followed by the ED method, and the DG-PD method has the lowest. Especially for the damage area aperture, the M-PD method is much larger than the other two methods, and the values are closest to the experiment. The M-PD method has large number of splashing particles, and the DG-PD method with node separation operation has small number of element splashing. However, the splashing elements are deleted during the calculation for the ED method. The three methods all have cross-shaped cracks, but the PD method can capture X-shaped cracks and bifurcated cracks, and the details captured by the DG-PD method are the most prominent. It is different for cracks that the shape donates several radial cracks in the experiment. Through the calculation time, it can be seen that the calculation amount of the PD method is much greater than that of the finite element method, and the M-PD method consumes more time than the DG-PD method.

High-speed explosion load
An explosion load with a charge of 15g and distance of 0.1m is used for simulation, and the boundary conditions are consistent with those in Section 4.2. To explore the crack propagation behavior under the explosion load, Fig. 4 shows the damage mode comparison diagram of the PD, ED method and experiment [24] under the explosion load in turn. The PD method produces a large number of broken elements in the center of the explosion, with numerous circumferential cracks, which are fish-scaled, and the diagonals are seriously damaged. The reflection of the stress wave at the boundary causes excessive local stress and forms splashing particles. Moreover, the ED method produces a cross-shaped crack at the center of the explosion, with the bifurcation and propagation of secondary cracks. There has a circular crack at the boundary, but there are no splashing fragments. Note that for the model of the experiment, the damage behavior of the PD method is extremely similar to the experimental result. In addition, it is interesting that the damage region at the center of the blast is remarkably similar, as the red circle marks in Fig. 4.   5 shows the curve of the element displacement with time at the same position of the PD and ED methods. The rate of change of the displacement is close to a constant under the explosive load, and the sinking depth of the ED method is slightly higher than that of the PD method. Fig. 6 shows the curve of the velocity with time at the same position of the PD and ED methods. Note that the ED method maintains an oscillating trend, while the PD method almost maintains a constant value after an instantaneous impact. It is worth noting that the velocity of the two methods finally converges to 8.3m/s, which effectively illustrates the rationality of the setting parameters in this article. To avoid the contingency of the analysis result, the elements of different positions are taken for statistical analysis. The results are universal and applicable.

Conclusion
In this study, the ED, DG-PD and M-PD methods are adopted to systematically analyze the damage behavior and microscopic crack propagation laws of architectural glass. The impact and explosion load are studied separately on the glass. Some conclusions are drawn as follows.
(1) Under bullet impact and explosion loads, different methods can simulate the crack propagation behavior of glass materials to a certain extent. However, the PD method has great advantages in some details, such as crack bifurcation and penetration.
(2) For the impact damage behavior of glass under 10m/s bullet impact load, different methods all show cross-shaped cracks. The two PD methods show splashing elements or particles, but the splashing particles of the M-PD method is more obvious. The DG-PD law has a better effect in capturing crack bifurcation in these methods.
(3) For the failure behavior of glass under explosion loading, the PD method is obviously better than the ED method in terms of modal appearance. However, the PD and ED methods have good agreement in the mechanical behavior, such as the specific displacement and velocity.
The architectural glass used in this article is ordinary soda-lime glass, and the laminated glass structure with a PVB layer can be researched and explored later. The bullet impact speed used in this article is very small, and the explosive charge is also not large. Therefore, the high speed and high explosion problem can be further studied in the future. In addition, the explosion algorithm is not effectively combined with the M-PD, and the next step can be done.