An Extended Bond-Based Peridynamic Approach for Analysis on Fracture in Brittle Materials

To enhance the calculation accuracy of bond-based peridynamics (BPD), a novel attenuation function is introduced to describe the effect of internal length on nonlocal long-range forces. Furthermore, the expression of the micromodulus function is deduced, and the corresponding fracture criteria are established. )e validity and accuracy of the extended bond-based peridynamic approach are illustrated by three numerical examples: 2D isotropic plate under uniaxial loading; plate with a circular cutout under quasi-static loading; and a diagonally loaded square plate with a center pre-existing crack. Finally, the influence of the width and the angle of the pre-existing crack on the fracture initiation time and the crack propagation paths are studied by applying the proposed approach.


Introduction
Predicting crack initiation and propagation accurately applying classical continuum mechanics (CCM) is still a major challenge for the community of solid mechanics [1,2]. e classical methods, such as the finite element method (FEM), are the most popular computational technique for structural computations [3]. However, it is formulated using spatial partial differential equations which are singularities and resulting in convergence problems at discontinuities [4], so some special techniques need to be devised, such as remeshing approaches is one of the main viable options to cope with these limitations [5], but remeshing processes are affected by numerical difficulties, complexity in computer programming, and often lead to calculation accuracy which is reduced [6].
Bond-based peridynamics (BPD) [7], first introduced to handle problems involving discontinuities and long-range forces by Silling in 2000, reformulate classical continuum mechanics in terms of spatial integral equations rather than partial differential equations. In contrast to classical continuum mechanics, the peridynamic equations are defined at the discontinuities, and thus, the fracture initiation and crack propagation can be simulated. It has the advantages of other numerical methods, such as the meshless, finite element, and molecular dynamics method [8]. And it has been widely applied to failure problems for brittle materials [9][10][11][12][13].
ere are certain limitations in bond-based peridynamics, such as the limitation of fixed Poisson's ratios and obvious surface effect. Recently, the state-based PD [14], the two-parameter bond-based PD [15,16], and the novel conjugated bond-based PD [17,18] were proposed and have solved the limitation of fixed Poisson's ratios in bond-based PD. And the volume method [19], the energy method [20], the fictitious node method [11], and coupling PD with the CCM method [4,21] were applied and have solved the surface effect in bond-based PD, and detailed information about surface correction techniques is introduced in [22]. However, there are still some unresolved problems in the bond-based PD that need to be addressed. For instance, the micromodulus which represents the bond stiffness is an invariable material constant. Actually, as a nonlocal model, the constant micromodulus of BPD ignores the internal length effect of long-range forces [23].
In BPD theory, a continuum is composed of infinitely many material points and interconnection through fictitious bonds, the influence of the material point interacting with X is assumed to exist in its own family (horizon), and the interaction between two points is always a constant regardless of the distance between them in one horizon. In other words, it can also be considered that the attenuation function g(ξ, δ) � 1 was introduced. And thus, the micromodulus function is a material constant. Based on the BPD model, Silling and Askari [8] researched the convergence in a fracture calculation and illustrated the properties of the method for modeling brittle dynamic crack growth. Agwai et al. [24] simulated three different experimental dynamic fracture problems and obtained data which coincide well with the results of literatures and experiment. Ha and Bobaru [10] analyzed the dynamic crack propagation and branching in brittle materials and showed results of convergence studies under uniform grid refinement. Dipasquale et al. [25] studied the dependence of crack paths on the orientation of regular 2D peridynamic grids.
When the internal length effect of nonlocal long-range forces is described, it is an effective way to introduce types of attenuation functions g(ξ, δ) to reflect the spatial distribution of the intensity of long-range forces. At present, the attenuation function g(ξ, δ) � 1 − (ξ/δ) which was introduced into the BPD model has been widely researched. Micromodulus of the BPD model which considered the effects of long-range forces hereinafter is called "variable micromodulus." Bobaru et al. [26] found that micromodulus has an influence on the rate of convergence, and the variable micromodulus leads to optimal rates of convergence independent of the grid used. Bobaru and Ha [27] found that micromodulus makes a difference to surface effect, and the variable micromodulus can weaken it. Ha and Bobaru [10] analyzed the dynamic crack propagation and obtained that the variable micromodulus has no significant effect on crack paths in brittle materials. Hu et al. [28] derived the formulation of the nonlocal J-integral and found that the variable micromodulus gives better convergence rates to the classical solutions in elasticity problems compared to the constant micromodulus. Cheng et al. [29] simulated various dynamic brittle fractures in functionally graded materials (FGMs) and suggested that the influence of variable micromodulus on the fracture behavior of FGMs is limited.
Based on the BPD which introduced the attenuation function g(ξ, δ) � e − (|ξ|/δ) 2 , Kilic and Madenci [30,31] investigated the elastic stability of simple structures to determine buckling characteristics by considering the problems of buckling load and buckling temperature. Using the BPD model which introduced the function g(ξ, δ) � (1 − (ξ/δ) 2 ) 2 , Huang [23,32] demonstrated that it is more accurate than the earlier BPD models using a fixed stiffness constant, and the crack propagation and the fracture mechanism of a cantilever concrete beam with preexisting notch are studied. Chen et al. [33] analyzed the influence of micromodulus by four different attenuation functions on BPD simulation of crack propagation and branching in brittle materials.
ough the problems of fracture were reported based on the several types of attenuation functions which were introduced to the BPD model, the influence of attenuation functions on computational accuracy and the optimal attenuation function need to be explored, and there is still lack of research on the influence of crack width and angle on brittle material failure based on the BPD model which considered the effects of long-range forces.
is study is organized as follows. In Section 2, the BPD theory is reviewed. In Section 3, the proposed extended BPD model is introduced, and several micromodulus functions based on various attenuation functions are defined. In Section 4, the validity and computational accuracy of the extended BPD model are illustrated by numerical examples. In Section 5, the influence of crack width and angle on fracture initiation time and crack propagation paths are explored by numerical simulations. Finally, Section 6 summarizes the conclusions resulted from this study.

Review of BPD Theory
In BPD theory, a continuum is composed of infinitely many material points. Each material point with a volume of V X and mass density of ρ is identified by its coordinates X in the initial configuration. As shown in Figure 1, each material point X is assumed to interact with all other material points X ′ within its horizon through fictitious bonds. e horizon is defined as H δ � X ′ ∈ Ω: |X ′ − X| ≤ δ . According to Newton's second law, the peridynamic equation of motion [7] for material point X at time t is expressed by where € u � (z 2 u(X, t))/zt 2 is the acceleration of material point X at time t, b is the body force density, and f is a pairwise force function that represents the force vector material point X ′ which exerts on material point X and is expressed by where u and u ′ are the displacements of material points X and X ′ at time t, respectively. ξ � X ′ − X denotes the relative position of material points X and X ′ in the initial configuration. η � u ′ − u denotes their relative displacement.
According to Silling and Askari [8] formulation, for linear elastic material, f is the derivative of microelastic strain density ω(η, ξ) with respect to relative displacement vector η. where where c(ξ, δ) is a micromodulus function which represents the bond stiffness, |ξ| denotes the magnitude of the vector ξ, and s is the stretch of the bond, and Substituting from equation (4) BPD assumes that half of the strain energy density due to the interaction of X and all X ′ is stored by X, and W PD (X) can be expressed as Normally, c(ξ, δ) is derived from the principle that the BPD strain energy density equals to the strain energy density based on classical continuum mechanics at material point X, i.e., e value of c(ξ, δ) for various cases is given in the following [19], as shown in Table 1. It is noted that Poisson's ratio is fixed to be 1/3 in the case of plane stress and 1/4 in the case of plane strain and 3D.

Extended BPD Approach.
In BPD, the intensity of the long-range forces between material points remains the same within the horizon, and the influence of the distance between the particles on the stiffness of the bond is ignored [23,32]; therefore, the function c(ξ, δ) is reduced to a constant c(δ). It causes the calculation error to become larger. When describing the internal length effect of nonlocal long-range forces, attenuation function g(ξ, δ) was introduced to reflect the spatial distribution of the intensity of long-range forces. e micromodulus of BPD can be written as where c(0, δ) is the bond stiffness function of BPD and g(ξ, δ) is the attenuation function, and it should meet the following conditions [23,32]: where δ(ξ) is the delta function, and it is noted that there is no connection between the delta function δ(ξ) and the horizon radius δ. e function c(ξ, δ) is equal to a constant c(δ) of BPD if g(ξ, δ) � 1. In the plane stress situation, the strain energy density in BPD can also be expressed as e strain energy density based on the classical continuum mechanics at material point X can be written as where ε 0 denotes the strain vector at material point X. Under the same conditions, the relative stretch of the bond s equals to the strain vector, i.e., s � ε 0 . In this study, the form of attenuation function (fourpower function) is uniquely proposed as follows: It meets all the requirements described in equation (10), and then comparing equations (11) and (12), the relation between the micromodulus c(0, δ) in the extended BPD approach and elastic modulus E and Poisson's ratio ] can be expressed as Similarly, the relationship between the micromodulus c(0, δ) of the extended BPD approach, elastic modulus E, and Poisson's ratio ] in the plane strain situation and 3D situation can be obtained, as shown in Table 2.
Another type of attenuation function g i (ξ, δ), i � 1, 2, . . . , 11 can also be introduced in (9), as shown in Table 3. ese attenuation functions include different kinds of distribution which allow forces to decrease with increasing distance between the pair of material points. In a similar way, the micromodulus c i (0, δ), i � 1, 2, . . . , 11 can be obtained by the strain energy density in the extended BPD approach which equals to the classical continuum  Table 3. And the shape of proposed micromodulus functions c i (ξ, δ), i � 1, 2, . . . , 11 is illustrated in Figure 2.
It is worth noting that these mathematical functions (attenuation functions) that we analyzed are all based on physical arguments (the radius of horizon δ and the relative position ξ). e selection of the micromodulus function is based on two points: (1) To reveal the rule that the pairwise force of the BPD bond decays with the increasing distance in the horizon of material points (2) To enhance the calculation accuracy of the BPD model, various types of the attenuation functions which are experiential such as "Gaussian," "Triangular," and "Parabolic" adopted in the literature are used and compared, as listed in Table 3 3.2. Extended BPD Fracture Criteria. e critical stretch s c is adopted to judge the breakage of the bond, and it can be obtained by the critical energy release rate G c [8,11]. A bond is deemed to be broken if the stretch between pairs of material points which is computed by equation (5) exceeds a certain critical value s c and failure occurs, and these two points cease to interact and the bonds between them to break irreversibly.
As shown in Figure 3, crack CD crosses the horizon of material point A located on the midperpendicular of CD and divides its horizon into two separate parts completely. e bond connecting material point A and any material point B within the yellow region is broken, and the strain energy stored in the bond is released.
G c denotes the energy needed to produce a unit area crack, and in 2D situation, G c can be expressed by Substituting from equation (9) into equation (15) results in In the 3D situation, G c can be expressed as Micromodulus Table 3: Eleven types of attenuation functions g(ξ, δ) and micromodulus functions c(ξ, δ) under the 2D plane stress condition.

Function type Attenuation function
g 9 cos(πξ/2δ) 3Eπ 2 /(2(π 2 − 8)hδ 3 ) 3Eπ 2 /(2(π 2 − 8)hδ 3 ) cos(πξ/2δ) Inverse proportional [36] g 10 (δ/ξ) 6E/πhδ 3 6E/πhδ 3 (δ/ξ) Quadratic power [36] g 11 (δ/ξ) 2 3E/πhδ 3 3E/πhδ 3 (δ/ξ) 2 In the present work, substituting from equations (13) and (14) into equation (16), therefore, the expression of s c in the plane stress situation can be obtained as When applying the extended BPD approach, once the condition s > s c is reached, the bond is deemed broken. In BPD, the bond stiffness is c(ξ, δ), an invariable material constant, and the expression of s c is given in [8,9], i.e., s c � �������� 4πG c /9Eδ. Failure is included in the material response through the scalar-valued function μ(ξ, t) that takes on values of either 1 or 0, which is defined as Function φ(X, t) is used to calculate the local damage at a material point, and local damage is the weighted ratio of the number of broken interactions to the total number of interactions [8,37]. It can be quantified as Note that 0 ≤ φ ≤ 1, and when the local damage is zero, it representing intact material, while a local damage of one means that all the interactions initially associated with the point have been eliminated.

2D Isotropic Plate under Uniaxial Loading.
A rectangular isotropic plate with 1.0 m length, 0.5 m width, and 0.01 m thickness is stretched by applying a uniaxial uniform tension loading of p � 200 MPa, as shown in Figure 4. e material density is 7850 kg/m 3 , Young's modulus 200 GPa, and Poisson's ratio ] � 1/3. ese material parameters are the same as those on page 155 of [11]. e spacing between material points is Δ � 0.01 mm, the total number of material points is 5000, and the time step ∆t � 1.0 s. Deformation of the plate is simulated by the extended BPD approach. e analytical solutions of the displacements are And the relative errors of numerical solutions are defined as where u denotes the numerically computed displacement value, u * denotes the analytical displacement value, and e denotes the relative error. It is noted that several surface correction techniques are introduced in [11,22]; in the present study, the correction method of the strain energy density from these literature studies was first employed, and then the internal length effect of nonlocal long-range forces is described.
Relative errors on the horizontal displacement are calculated comparing results from the above BPD models with different attenuation functions and analytical results, as shown in Figure 5, where the horizon radius is chosen as δ � 3.015Δ. e maximum relative error on the horizontal Figure 3: Integration domain of the micropotentials crossed by a crack [11]. displacement calculated by the BPD model with different attenuation functions is shown in Figure 6, where the horizon radius is chosen as δ � 3.015Δ, δ � 4.015Δ, and δ � 5.015Δ, respectively. It can be observed that the maximum relative errors all occur in the corners due to the surface effect of the PD method.
e maximal relative errors on the horizontal displacement for the BPD model with eleven attenuation functions are different. In especial, the errors are reduced effectively by the four-power function g 2 � 1 − (ξ/δ) 4 , and the calculation accuracy of the BPD model with this function is the highest regardless of the horizon radius. ese results showed that not all attenuation functions have the ability to improve the accuracy of the BPD model which was corrected by the method of the strain energy density effectively. e extended BPD approach in the present work is actually the four-power function g 2 � 1 − (ξ/δ) 4 that with the highest accuracy was introduced in the BPD model. It can describe the spatial distribution of the intensity of nonlocal forces and further weaken the surface effect.

Plate with a Circular Cutout under Quasi-Static Loading.
A square plate with 50 mm side length and 0.5 mm thickness having a circular cutout at the center with a radius of 5 mm is used. e material density is 8000 kg/m 3 , Young's modulus 192 GPa, and Poisson's ratio ] � 1/3. ese material parameters are the same as those on page 155 of [11]. e plate is subject to a slow rate of stretch along its horizontal edges, and the velocity boundary condition is set as _ u y (x, ± L/2, t) � ± 2.7541 × 10 − 7 m/s and imposed on the fictitious boundary layer with a width of 3Δ, as shown in Figure 7. e spacing between material points is Δ � 0.5 mm, and the total number of material actual points is 9876. e horizon radius is chosen as δ � 3.015Δ, and time step ∆t � 1.0 s. e critical stretch value s c � 0.02138. Deformation and fracture of the plate are simulated by the extended BPD model. e process of fracture initiation and crack propagation is shown in Figure 8. Fracture initiates at the 605 step and the crack propagates along the direction of x-axis. Crack propagates till up to the left and right boundaries of the plate at the 905 step. e final crack propagation paths are consistent with the results reported by Huang et al. [4], Madenci and Oterkus [11], and Ochoa-Ricoux [20], as shown in Figure 9.

A Diagonally Loaded Square Plate with a Center Pre-Existing Crack.
A square plate with 150 mm side length and 5 mm thickness having a pre-existing centered crack of length 45 mm is stretched by applying a diagonal load, as shown in Figure 10. e properties of the plate are specified as Young's modulus 2.94 GPa, Poisson's ratio ] � 1/3, mass  [38], as shown in Figure 11(f ).

The Influence of Crack Width on Brittle Material Fracture
A square plate with 50 mm side length and 0.1 mm thickness having a pre-existing centered crack of 10 mm length is subjected to uniaxial loading as shown in Figure 12. Material properties are specified as Young's modulus 192 GPa, mass density 8000 kg/m 3 , Poisson's ratio ] � 1/3 (parameters are taken from Madenci and Oterkus [11]), and the critical stretch value of s c � 0.04781. e velocity boundary condition is V(t) � 20 m/s. Simulation is implemented by the extended BPD model with Δ � 0.1 mm, δ � 3.015Δ, and time step 1.3367E − 8 s.
Adopting the same parameters and loading condition, the investigation is made for different crack widths of dx � 0.5 mm, 2 mm, 3 mm, 4 mm, and 8 mm and different orientation angles β � i × 15 ∘ (i � 0, 1, 2, 3, 4, 5, 6). Numerical simulation is implemented by extended BPD, and the crack patterns are shown in Figure 13.
It can be seen that the width of pre-existing cracks has a notable influence on the final crack paths under uniaxial loading. When the orientation angle β � 0 ∘ , the crack propagating in plates of crack width presents obviously different forms. When crack width is less than 4 mm, there is only one crack growth and extends horizontally; when it is greater or equal to 4 mm, there occurs two-crack branching along the up and down edges of the pre-existing crack tip, and both of them extend and split. When 0 ∘ < β < 90 ∘ , the crack of the plates initiates at the pre-existing crack tip which is closer to the left and right boundary and propagates horizontally regardless of crack width. e pre-existing crack with the width of dx � 4 mm and the process of fracture initiation and crack propagation are shown in Figures 14(a)-14(d). Fracture initiates at 5.2 μs, after the loading, and the crack inclination angle is around 12°with respect to the direction of x-axis at 8.2 μs. Crack propagates till up to the left and right boundaries of the plate at 9.8 μs. Figure 15shows the time of fracture initiation in plates of different crack widths and having a pre-existing crack with different orientation angles β. It can be seen that the fracture initiation time is sensitive to the variation of crack width when the precrack is horizontal and vertical. When β � 0 ∘ , as crack width increases, the fracture initiation time delays; when β � 90 ∘ , the fracture initiation time advances with crack width increasing; and when 0 ∘ < β < 90 ∘ , the fracture initiation time is not sensitive to the variation of crack width, and orientation angle β � 75 ∘ of the preexisting crack is usually harder to break than other orientation angles.

Conclusion
In this study, an extended BPD approach is proposed through introducing a new attenuation function. e validity and precision of the approach are verified by three numerical examples. By applying the approach, investigation is made on the fracture of plates with pre-existing cracks of different widths and different angles under uniaxial loading. It can be concluded that the widths and angles of the preexisting cracks have a remarkable influence on the fracture initiation times and the final crack propagation paths.

Data Availability
e data used to support the findings of this study have been deposited in the "figshare" repository (DOI: 10.6084/m9. figshare.9816296).

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper. Mathematical Problems in Engineering 11