|Computer Modeling in Engineering & Sciences|
Modelling of Contact Damage in Brittle Materials Based on Peridynamics
1Nanjing Institute of Railway Technology, Nanjing, 210031, China
2Tongji University, Shanghai, 200092, China
3Department of Engineering Mechanics, Hohai University, Nanjing, 211100, China
4Nanjing High Speed Gear Manufacturing Co., Ltd., Nanjing, 211100, China
*Corresponding Author: Jingjing Zhao. Email: firstname.lastname@example.org
Received: 27 April 2021; Accepted: 11 August 2021
Abstract: As a typical brittle material, glass is widely used in construction, transportation, shipbuilding, aviation, aerospace and other industries. The unsafe factors of glass mainly come from its rupture. Thus, establishing a set of prediction models for the cracks growth of glass under dynamic load is necessary. This paper presents a contact damage model for glass based on the ordinary state-based peridynamic theory by introducing a contact force function. The Hertz contact (nonembedded contact) problem is simulated, and the elastic contact force is determined by adjusting the penalty factor. The proposed model verifies the feasibility of penalty-based method to simulate the contact problem of glass. The failure process of glass specimen under impact is simulated, where two loading methods, the drop ball test and the split Hopkinson pressure bar are considered. Numerical results agree well with the experimental observations, thereby verifying the effectiveness of the proposed model.
Keywords: The ordinary state-based peridynamic; glass; numerical simulation; contact damage; impact failure
Glass is a highly homogeneous brittle material and it can be made into any shape and mass produced due to its homogeneity. Hence, glass materials are widely used in construction, transportation, shipping, aviation, aerospace and other industries. However, the glass is extremely easy to break when its surface is subjected to impact load. Effectively revealing the mesoscopic failure mechanism of glass under low-velocity impact loads is of great importance to improve its resistance and safety.
In the early days, many scholars investigated the damage of glass under impact, mainly focusing on the crack initiation, propagation, and crack development patterns of glass through a series of experimental studies [1–6]. In 1931, Andrews  pointed out that a force threshold is found for the crack formation of glass under impact load. Roesler  considered the conical indentation of a glass specimen and measured its fracture energy to verify Griffith  energy balance condition for brittle material cracking. Knight et al. , Chaudhri et al.  studied the crack development patterns of glass under medium and high-speed loading by changing the impact velocity of iron ball. The results show that the crack appears in the non-loading section, and the Hertzian cone angle of the crack zone varies with the loading speed. Ball et al.  observed that glass specimens exhibit a star-shaped fragmentation failure mode under low-speed conditions. In recent years, Bouzid et al. , Nyongue et al. , and Daryadel et al.  conducted relevant quantitative experimental analysis on the dynamic response and fracture criterion for glass under impact loading. The failure mechanism for different types of glass is realized in the work of Kumar et al. . These test results reveal the response law and failure mechanism of glass as a typical brittle material under impact load from different aspects.
With the development of computers, numerical simulation methods have been widely used to reproduce and explain the phenomena observed in experiments. Numerical models based on traditional continuum mechanics, such as finite element method  and boundary element method , are gradually used in many experimental studies. However, for discontinuities, such as crack bifurcation in solid materials and structures, traditional numerical methods are faced with the problems of singularity and low computational efficiency. The use of grid reconstruction  or the method of adding a cohesive element  in the finite element would cause to grid-dependent results. The partition algorithm  and fictitious crack model  adopted by the boundary element method (BEM) have similar limitations to the finite element method in the analysis of crack propagation problems. Therefore, researchers have proposed the extended finite element method . Compared with the traditional finite element method, it reduces the strict requirements on the mesh discontinuity. However, in the construction of enrichment function, the extended finite element method needs to know the characteristics of the problem to be solved in advance. This condition is relatively demanding for complex problems, such as crack branching and multi-crack intersection. In recent years, many researchers have studied the contact damage of brittle materials under impact load. Oliveira et al.  proposed an alternative BEM formulation to model the cohesive stresses through the domain term of the direct integral. Wang et al.  proposed a field-enriched finite element method to simulate the failure process of rocks. Kouet et al.  presented a bonded-particle methodology to investigate the crack growth and crack branching. Bo et al. [21,22] developed some explicit Galerkin formulations which can capture ductile fractures during high-speed impacts and simulate fracture of a plate under impact loads.
To solve contradiction between the continuity assumption and the discontinuity phenomenon of the failure problem, Silling  proposed a nonlocal method called peridynamics (PD)  to described the motion process of the material points through an integral equation. In the development of the peridynamics, researchers have conducted numerous quantitative analyses on the brittle fracture of polycrystalline materials , impact failure of shale materials , and the thermal brittleness and failure of glass plate under impact [27–29], and the problem of interaction between ice and seawater . Those results show that peridynamics does not have the singularity problem when analyzing the failure problem and can simulate the whole process of the material, including macroscopic crack initiation, propagation, and final failure. The above simulation of failure is based on the bond-based peridynamics theory (BBPD). However, bond-based peridynamics theory has some problems, such as the limitation of Poisson’s ratio and the lack of connection with traditional continuum theory. Silling et al. [31–32] proposed an ordinary state-based peridynamic theory (OSBPD) and a non-ordinary state-based peridynamics(NOSB PD). Both of them inherit the advantages of BBPD for solving discontinuous problems, and have a similar definition of state quantities to traditional physical quantities. To address the limitations of material calculation scale. Breitenfeld et al. , O’Grady et al. , and Chowdhury et al.  supplemented and developed the OSBPD theory. Song et al.  proposed a state-based peridynamic model by using adaptive particle refinement to simulate the formation of water ice crater under impact loads. Zhou et al.  developed a 2D elastoplastic model of the ordinary state-based peridynamic theory to analyze the plastic zone of the crack tip of rock material. Wu et al.  proposed a concrete impact failure model based on NOSBPD theory. In recent years, many scholars began to study the contact model based on PD method. Xin et al.  developed a NOSBPD model for brittle fracture to simulate the edge-on impact and drop ball test and discussed the contact algorithm between the projectile and target. Littlewood et al.  summarized the simulation results of finite element method and peridynamics. A combined approach of finite element method and peridynamics is utilized via a contact algorithm. Ye et al.  proposed a continuous contact detection algorithm to simulate the brittle failure behavior of ice during the contact between propeller and ice. This algorithm can well capture the ice damage characteristics. Kamensky et al.  summarized several existing peridynamic contact friction models and introduced a state-based nonlocal friction formulation to demonstrate the properties of various peridynamic contact models through some impact and penetration problems. Silling et al.  proposed a new PD model to simulate the elastoplastic response, creep, and fracture. In the case of small deformation, this model is consistent with classical Hertz contact analysis. Wang et al.  proposed a 3D conjugated bond-pair-based peridynamic model to simulate the failure characteristics of rock materials with different forms of central fissure.
This paper presents a new contact damage model based on the OSBPD theory for glass by introducing a contact force function. The arrangement of the rest sections is as follows. In Section 2, the governing equations based on OSBPD theory are described in brief, followed by the contact model and algorithm. In Section 3, the Hertz contact (nonembedded contact) problem is simulated, and the elastic contact force is determined by adjusting the penalty factor. The proposed model verifies the feasibility of penalty-based method to simulate the contact problem of glass. In Section 4, the failure process of glass specimen under impact is simulated, and two loading methods, drop ball test (DBT) and split Hopkinson pressure bar (SHPB), are considered. The numerical results agree well with the experimental results, which verifies the effectiveness of the proposed model.
2.1 Fundamentals of Peridynamics Theory
In PD theory, an object occupies a region (namely the reference configuration shown in Fig. 1). A certain interaction (nonlocal interaction) is assumed between a particle and other particle within a horizon , where is the truncation radius. Let for convenience because has nothing to do with time. Supposed a region is occupied by a particle on the object. By defining the relative position state vector and relative displacement state vector , where is the displacement of the particle , the deformation state vector is expressed as
where is the position of particle after deformation. The force state vector between and can be defined as
where is the state vector of the deformation direction, showing the direction of the interaction force between the material points is consistent with its relative deformation direction; is the force state scalar. The force of particle exerted by its neighboring particle is expressed as
Therefore, a momentum balance equation at particle is defined as
The above formula is the governing equation of peridynamics [31–32], and its above derivation shows that the peridynamics theory still belongs to the Lagrangian system. The difference from the traditional continuum mechanics is that it considers the nonlocal interaction between the particles.
where is the stress tensor.
2.2 Force State of State-Based Peridynamics Theory
Let denote the initial length state scalar between the interacting particles, The weighted volume is defined as
where is the influence state and describes the relative degree of interactions between particles. The elongation  between a pair of material points is defined as
The volume dilatation is defined as
In accordance with elastic mechanics, the extension scalar state can be decomposed into isotropic part
and deviatoric part
Therefore, the elastic energy density contains two parts
Correspondingly, the scalar force can be expressed as
In accordance with the energy equivalence between the deformation energy density (11) and the strain energy density of classical elasticity, can be obtained as 
The above formula is the force state of state-based peridynamics theory, which includes traditional material parameters, such as bulk modulus and shear modulus . For the damaged material points, the scalar function  is defined as
where is a scalar function, and it is defined as
where is the fracture judgment parameter and is the corresponding critical value. If the extension scalar state is selected as the basis for judging whether the bond between two material points is broken, then and are defined as
where  is the energy release rate. For the local damage defined by Eq. (14), the damage between material point pairs is shown in Fig. 2.
2.3 Contact Model of Peridynamics
The contact problem of two objects shown in Fig. 3 is considered. Let be the coordinate of any material point on the . The distance between this point and the nearest point on can be expressed as
where the unit normal of the contact surface is
Considering the impenetrability of contact, any material point on should satisfy , where indicates that point is separated from surface , denotes contact, and represents that and penetrate each other. This condition is not allowed. Therefore, if the impenetrability is satisfied at time , then the following condition should be met at time
Note and , so
The energy function, which includes and systems, is expressed as
where is the energy function without considering the contact force and contact conditions. For and , , which is equivalent to the Eq. (7), where is the potential energy function of contact force. These contact forces are the acting force and reaction force. Thus, we have
where . is the additional functional of the introduced contact condition Eq. (19) and is defined as
where is the penalty factor. Let , we have
For a certain contact deformation body, if its contact surface is , then the contact force is
where and . Once the object is in contact, the contact force is applied to the material point on the contact surface in the form of body force density, that is
The contact force is related to the selection of the penalty factor. In theory, the larger the penalty factor , the higher the calculation accuracy, and the closer to the nonpenetration condition. This condition will be discussed in detail in the numerical solution process and specific examples. The penalty factor is collectively referred to as the contact force modulus in this paper.
2.4 Numerical Discretization and Solution
Regarding the spatial discretization, peridynamics divides the material uniformly into a cubic lattice of side length , and the volume for any material point is (Fig. 4). Therefore, all material points in can be defined as
Correspondingly, the discretization of the motion equation considering the contact force in the space domain is
Similar discrete operations are applied to the weighted volume and volume dilatation .
Let and time step . The the formula of central difference scheme is obtained as
Substituting the above formula into the Eq. (29). Thus, we have
where . The difference scheme Eq. (31) should satisfy the Courant–Friedrichs–Lewy (CFL) stability condition. Let
The time step should satisfy 
The virtual rigid “spring” is added, that is, the rigidity of contact area is improved. Eq. (33) is modified as
The value of penalty factors is discussed in the following calculation examples.
3 Verification of the Elastic Contact Problem
The Hertzian contact (no embedded contact) is considered, as shown in Fig. 5. The size parameters of the deformation body are and . The elastic modulus and Poisson’s ratio of the material are and . The mass density is . The radius of the rigid sphere is set as , and approximately satisfies the assumption of semi-infinite deformation body. The fixed horizon size is set as . The time step is set as . For the theoretical solution of this problem, the derivation can be found in the Appendix.
When the rigid sphere is in contact with the deformable body at a constant speed (it can be regarded as a quasi-static problem), are used for calculation to analyze the influence of contact modulus on the result. The discrete distance of the deformable body is set as . The results of embedded distance and contact force corresponding to indentation depth are shown in Figs. 6 and 7. As shown in Figs. 6 and 7, the contact force gradually approaches the theoretical solution with the increase in contact modulus . Therefore, for non-embedded contact problems, the value of contact modulus should be as large as possible to satisfy the computer truncation error and the stability requirements of differential solution.
For the case of , the displacement nephogram of deformation body in direction at time is drawn in Fig. 8. For the contact problem of a rigid sphere, the contact surface is arc-shaped, and its displacement decreases outward along the center of the contact surface, which is also consistent with the actual situation.
The discrete distance is still set as , and the contact modulus is taken as . The radius of rigid sphere is changed by increment . When the indentation depth of the rigid sphere is , contact force is plotted in Fig. 9. With the increase in radius of the rigid sphere, contact force is closer to the corresponding theoretical solution. When and the relative discrete size is set as , the calculation error of contact force is . When and , the calculation error of contact force drops to . The calculation error of the contact force decreases with the decrease in the relative discrete size. The smaller the relative discrete size (the more discrete material points), the more time consuming the calculation. As shown in Fig. 9, the effect of improving the error is unremarkable when the relative discrete size is less than a certain critical value. When and the contact modulus of the rigid sphere is still taken as , change the discrete distance of the material points and take . The calculation results are shown in Fig. 10. A similar conclusion can be obtained, that is, the calculation error of the contact force gradually decreases with the decrease in distance between the material points.
In summary, the influencing factors in the contact model of peridynamics are discussed and compared with the theoretical solution. The penalty-based method can be used to simulate the contact problem well by selecting the appropriate relative discrete size and the contact modulus. This process meets the requirements of computer truncation error and stable calculation. Therefore, on the basis of experience and the calculation results of this case, the calculation formula of contact modulus can be given as follows:
The range of adjustment coefficient is from 1 to 50; is the elastic modulus; is the discrete distance of material points.
4 Contact Failure of Glass Specimen under Impact
Float glass material  is selected for analysis. The elastic modulus is , the Poisson’s ratio is , the mass density is , and the critical energy release rate is set as . The length and width of the glass sample are the same as , and the thickness is . For the simulations of DBT and SHPB, the models are shown in Figs. 11 and 12. In the simulation of DBT, the mass of the iron ball is , and the mass density is . The pressure bar is made of aluminum in the simulation of SHPB. Its elastic modulus is , the Poisson’s ratio is , and the mass density is .
In the analysis of peridynamics, the contact modulus is taken as , and the discrete distance of glass samples is when the penalty-based method is used to calculate the above two models. The time step value is . The low loading rate of DBT and the high loading rate of SHPB are considered. The results are shown in Figs. 13 and 14. For low loading rate () in Fig. 13, the final damage form is mainly a small amount of fragments. For the high-speed impact of SHPB, the sample appears to be broken. The simulation results are in good agreement with the test results (Fig. 15).
Considering the low loading rate () of DBT and the high loading rate () of SHPB, the comparison between the two simulations and test results are shown in Fig. 16. The proposed contact damage model can deal with the impact failure of glass qualitatively and quantitatively. The curve of contact force vs. indentation depth is shown in Fig. 17. The comparison of DBT and SHPB results by using different load rates shows that the impact speed increases with the increase in impact speed.
The damage of glass specimens in different times is shown in Fig. 18. With the increase in loading speed, the damage time becomes shorter, and the damage degree becomes more serious. The simulation results are consistent with the experimental results. Therefore, the contact model of peridynamics can well simulate the failure of brittle materials under impact conditions.
This paper proposes a new contact damage model of brittle materials based on ordinary state-based peridynamic theory by introducing the contact force function. The simulation results of the Hertzian contact (nonembedded contact) problem show that the embedded distance gradually approaches zero (the assumption of no embedding is satisfied), and contact force gradually approaches the theoretical solution with the increase in contact modulus . Therefore, for nonembedded contact problems, the value of the contact modulus should be as large as possible to satisfying the computer truncation error and the stability requirements of differential solution. The calculation error of the contact force decreases with the decrease in relative discrete size. The smaller the relative discrete size (the more discrete material points), the more time consuming the calculation. When the relative discrete size is less than a certain critical value, the effect of improving the error is unremarkable. The penalty-based method of peridynamics can be used to simulate the contact problem well by selecting the appropriate relative discrete size and the contact modulus. The low loading rate of DBT and high loading rate of SHPB are considered to simulate the failure process of glass specimen under impact load. For low loading rate, the final damage form is mainly a small amount of fragments. For the high-speed impact of SHPB, the sample appears to be broken. The simulation results are in good agreement with the test results.
The proposed damage contact model of peridynamics has a good application prospect and can be extended to the simulation of other similar materials, such as the impact fracture process of ice. The calculation model can be enriched and developed. The contact force only considers the vertical force. If the horizontal friction force is considered simultaneously, the response and failure mechanism of typical brittle materials under impact load can be better analyzed.
Acknowledgement: This research is financially by the National Natural Science Foundation of China (Nos. 11932006, U1934206, 12002118). And The Natural Science Foundation of the Jiangsu Higher Education Institutions of China (No. 20KJB580015). The snapshots were performed using the visualisation tool Ensight.
Funding Statement: This study was funded by National Natural Science Foundation of China (Nos. 11932006, U1934206), Recipient: Qing Zhang. And National Natural Science Foundation of China (No. 12002118), Recipient: Xin Gu. And Natural Science Foundation of the Jiangsu Higher Education Institutions of China (No. 20KJB580015). Recipient: Runpu Li.
Conflicts of Interest: The authors claim that none of the material in the paper has been published or is under consideration for publication elsewhere. The publication has been approved by all co-authors. We have no conflict of interest to declare, All data generated or analysed during this study are included in this published article.
The Herzian contact of the rigid body is considered in Fig. A-1, and the contact stress of the sphere can be expressed as
For a semi-infinite elastic body, the elastic mechanics solution of the concentrated force is 
When the contact surface ,
Then, when the distributed force is considered
where . Substituting the above Eq. (A2) into the Eq. (A4),
where is the indentation depth, is the radius of the sphere. The relationship between contact force and indentation depth is
|This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.|