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

Authors: Shuang Li; Yanli Jin; Xiaohua Huang; Lianjun Zhai
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

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. The 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.


---

## Body

## 1. 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]. The 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–13].There 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 withX 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 functionsgξ,δ 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 functiongξ,δ=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−ξ/δ22, 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 pre-existing 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. Though 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.This study is organized as follows. In Section2, 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.


## 2. Review of BPD Theory

In BPD theory, a continuum is composed of infinitely many material points. Each material point with a volume ofVX 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. The 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(1)ρu¨=∫HδfdVX′+bX,t,where u¨=∂2uX,t/∂t2 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(2)f=fu′,u,X′,X,t=fu′−u,X′−X=fη,ξ,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.Figure 1
Deformation of PD material pointsX and X′ and developing equal and opposite pairwise force densities [11].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 η.(3)f=∂ωη,ξ∂η,where(4)ωη,ξ=cξ,δs22ξ,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(5)s=ξ+η−ξξ.Substituting from equation (4) into equation (3) results in(6)f=cξ,δsξ+ηξ+η.BPD assumes that half of the strain energy density due to the interaction ofX and all X′ is stored by X, and WPDX can be expressed as(7)WPDX=12∫Hδωη,ξdVX′=12∫Hδcξ,δs22ξdVX′.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.,(8)WPDX=WCCMX.The value ofcξ,δ 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.Table 1
Micromodulus and Poisson’s ratio in BPD.
Expression of micromodulus
Poisson’s ratio
Micromodulus
3D
6
E
/
π
δ
4
1
−
2
ν
ν
=
1
/
4
12
E
/
π
δ
4
2D
Plane strain
6
E
/
π
h
δ
3
1
+
ν
1
−
2
ν
ν
=
1
/
4
48
E
/
π
h
δ
3
Plane stress
6
E
/
π
h
δ
3
1
−
ν
ν
=
1
/
3
9
E
/
π
h
δ
3
h denotes the thickness of material, E denotes the elastic modulus, and ν denotes Poisson’s ratio of the material.


## 3. Extended BPD Approach and Fracture Criteria

### 3.1. 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. The micromodulus of BPD can be written as(9)cξ,δ=c0,δgξ,δ,where c0,δ is the bond stiffness function of BPD and gξ,δ is the attenuation function, and it should meet the following conditions [23, 32]:(10)gξ,δ=g−ξ,δ,limξ⟶0gξ,δ=maxg,limξ⟶δgξ,δ=0,∫−∞∞limδ⟶0gξ,δ=∫−∞∞δξdx=1,where δξ is the delta function, and it is noted that there is no connection between the delta function δξ and the horizon radius δ. The 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(11)WPDX=12∫Hδ12ξc0,δgξ,δη2dVX′.The strain energy density based on the classical continuum mechanics at material pointX can be written as(12)WCCMX=E1−νε02=E1−νsξ,η2,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 (four-power function) is uniquely proposed as follows:(13)gξ,δ=1−ξδ4,ξ≤δ,0,ξ>δ.It meets all the requirements described in equation (10), and then comparing equations (11) and (12), the relation between the micromodulus c0,δ in the extended BPD approach and elastic modulus E and Poisson’s ratio ν can be expressed as(14)c0,δ=63E4πhδ3,ν=13.Similarly, the relationship between the micromodulusc0,δ 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.Table 2
Micromoduluscξ,δ of the extended BPD approach.
Micromodulusc0,δ
Poisson’s ratio
Micromoduluscξ,δ
3D
24
E
/
π
δ
3
ν
=
1
/
4
24
E
/
π
δ
3
1
−
ξ
/
δ
4
2D
Plane strain
84
E
/
5
π
h
δ
3
ν
=
1
/
4
84
E
/
5
π
h
δ
3
1
−
ξ
/
δ
4
Plane stress
63
E
/
4
π
h
δ
3
ν
=
1
/
3
63
E
/
4
π
h
δ
3
1
−
ξ
/
δ
4Another type of attenuation functiongiξ,δ,i=1,2,…,11 can also be introduced in (9), as shown in Table 3. These 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 ci0,δ,i=1,2,…,11 can be obtained by the strain energy density in the extended BPD approach which equals to the classical continuum mechanics under the same conditions, respectively. The expression of ci0,δ,i=1,2,…,11 and ciξ,δ,i=1,2,…,11 in the plane stress situation is listed in Table 3. And the shape of proposed micromodulus functions ciξ,δ,i=1,2,…,11 is illustrated in Figure 2.Table 3
Eleven types of attenuation functionsgξ,δ and micromodulus functions cξ,δ under the 2D plane stress condition.
Function type
Attenuation function
g
i
ξ
,
δ
c
i
0
,
δ
c
i
ξ
,
δ
=
c
0
,
δ
⋅
g
i
ξ
,
δ
Constant [9, 10, 34]
g
1
1
9
E
/
π
h
δ
3
9
E
/
π
h
δ
3
Four power
g
2
1
−
ξ
/
δ
4
63
E
/
4
π
h
δ
3
63
E
/
4
π
h
δ
3
1
−
ξ
/
δ
4
Exponential [33]
g
3
e
−
ξ
/
δ
3
e
E
/
2
e
−
5
π
h
δ
3
3
e
E
/
2
e
−
5
π
h
δ
3
e
−
ξ
/
δ
Gaussian [30, 31]
g
4
e
−
ξ
/
δ
2
24
e
E
/
3
e
−
4
π
h
δ
3
24
e
E
/
3
e
−
4
π
h
δ
3
e
−
ξ
/
δ
2
Triangular [33]
g
5
1
−
ξ
/
δ
36
E
/
π
h
δ
3
36
E
/
π
h
δ
3
1
−
ξ
/
δ
Semielliptical [33]
g
6
1
−
ξ
/
δ
2
45
/
2
π
h
δ
3
45
/
2
π
h
δ
3
1
−
ξ
/
δ
2
Parabolic [33]
g
7
1
−
ξ
/
δ
2
2
315
E
/
8
π
h
δ
3
315
E
/
8
π
h
δ
3
1
−
ξ
/
δ
2
2
Six power [35]
g
8
1
−
ξ
/
δ
2
3
945
E
/
16
π
h
δ
3
945
E
/
16
π
h
δ
3
1
−
ξ
/
δ
2
3
Cosine [35]
g
9
cos
π
ξ
/
2
δ
3
E
π
2
/
2
π
2
−
8
h
δ
3
3
E
π
2
/
2
π
2
−
8
h
δ
3
cos
π
ξ
/
2
δ
Inverse proportional [36]
g
10
δ
/
ξ
6
E
/
π
h
δ
3
6
E
/
π
h
δ
3
δ
/
ξ
Quadratic power [36]
g
11
δ
/
ξ
2
3
E
/
π
h
δ
3
3
E
/
π
h
δ
3
δ
/
ξ
2Figure 2
The diagram of various micromodulus functionscξ,δ under the 2D plane stress condition.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 ξ). The 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 Table3


### 3.2. Extended BPD Fracture Criteria

The critical stretchsc is adopted to judge the breakage of the bond, and it can be obtained by the critical energy release rate Gc [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 sc and failure occurs, and these two points cease to interact and the bonds between them to break irreversibly.As shown in Figure3, crack CD crosses the horizon of material point A located on the midperpendicular of CD and divides its horizon into two separate parts completely. The 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.Figure 3
Integration domain of the micropotentials crossed by a crack [11].G
c denotes the energy needed to produce a unit area crack, and in 2D situation, Gc can be expressed by(15)Gc=h∫0δ∫zδ∫−cos−1z/ξcos−1z/ξ12cξ,δsc2ξξdθdξdz=h∫0δ∫zδ∫0cos−1z/ξcξ,δsc2ξξdθdξdz.Substituting from equation (9) into equation (15) results in(16)Gc=h∫0δ∫zδ∫0cos−1z/ξc0,δgξ,δsc2ξξdθdξdz.In the 3D situation,Gc can be expressed as(17)Gc=∫0δ∫02π∫zδ∫0cos−1z/ξ12c0,δgξ,δsc2ξξ2sinθdθdξdϕ.In the present work, substituting from equations (13) and (14) into equation (16), therefore, the expression of scin the plane stress situation can be obtained as(18)sc=32πGc63Ehδ.When applying the extended BPD approach, once the conditions>sc is reached, the bond is deemed broken. In BPD, the bond stiffness is cξ,δ, an invariable material constant, and the expression of sc is given in [8, 9], i.e., sc=4πGc/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(19)μξ,t=1,ifsξ,t′<sc,0<t′<t,0,else.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(20)φX,t=1−∫Hδμt,ξdVX∫HδdVX,0≤φ≤1.Note that0≤φ≤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.


## 3.1. 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. The micromodulus of BPD can be written as(9)cξ,δ=c0,δgξ,δ,where c0,δ is the bond stiffness function of BPD and gξ,δ is the attenuation function, and it should meet the following conditions [23, 32]:(10)gξ,δ=g−ξ,δ,limξ⟶0gξ,δ=maxg,limξ⟶δgξ,δ=0,∫−∞∞limδ⟶0gξ,δ=∫−∞∞δξdx=1,where δξ is the delta function, and it is noted that there is no connection between the delta function δξ and the horizon radius δ. The 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(11)WPDX=12∫Hδ12ξc0,δgξ,δη2dVX′.The strain energy density based on the classical continuum mechanics at material pointX can be written as(12)WCCMX=E1−νε02=E1−νsξ,η2,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 (four-power function) is uniquely proposed as follows:(13)gξ,δ=1−ξδ4,ξ≤δ,0,ξ>δ.It meets all the requirements described in equation (10), and then comparing equations (11) and (12), the relation between the micromodulus c0,δ in the extended BPD approach and elastic modulus E and Poisson’s ratio ν can be expressed as(14)c0,δ=63E4πhδ3,ν=13.Similarly, the relationship between the micromodulusc0,δ 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.Table 2
Micromoduluscξ,δ of the extended BPD approach.
Micromodulusc0,δ
Poisson’s ratio
Micromoduluscξ,δ
3D
24
E
/
π
δ
3
ν
=
1
/
4
24
E
/
π
δ
3
1
−
ξ
/
δ
4
2D
Plane strain
84
E
/
5
π
h
δ
3
ν
=
1
/
4
84
E
/
5
π
h
δ
3
1
−
ξ
/
δ
4
Plane stress
63
E
/
4
π
h
δ
3
ν
=
1
/
3
63
E
/
4
π
h
δ
3
1
−
ξ
/
δ
4Another type of attenuation functiongiξ,δ,i=1,2,…,11 can also be introduced in (9), as shown in Table 3. These 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 ci0,δ,i=1,2,…,11 can be obtained by the strain energy density in the extended BPD approach which equals to the classical continuum mechanics under the same conditions, respectively. The expression of ci0,δ,i=1,2,…,11 and ciξ,δ,i=1,2,…,11 in the plane stress situation is listed in Table 3. And the shape of proposed micromodulus functions ciξ,δ,i=1,2,…,11 is illustrated in Figure 2.Table 3
Eleven types of attenuation functionsgξ,δ and micromodulus functions cξ,δ under the 2D plane stress condition.
Function type
Attenuation function
g
i
ξ
,
δ
c
i
0
,
δ
c
i
ξ
,
δ
=
c
0
,
δ
⋅
g
i
ξ
,
δ
Constant [9, 10, 34]
g
1
1
9
E
/
π
h
δ
3
9
E
/
π
h
δ
3
Four power
g
2
1
−
ξ
/
δ
4
63
E
/
4
π
h
δ
3
63
E
/
4
π
h
δ
3
1
−
ξ
/
δ
4
Exponential [33]
g
3
e
−
ξ
/
δ
3
e
E
/
2
e
−
5
π
h
δ
3
3
e
E
/
2
e
−
5
π
h
δ
3
e
−
ξ
/
δ
Gaussian [30, 31]
g
4
e
−
ξ
/
δ
2
24
e
E
/
3
e
−
4
π
h
δ
3
24
e
E
/
3
e
−
4
π
h
δ
3
e
−
ξ
/
δ
2
Triangular [33]
g
5
1
−
ξ
/
δ
36
E
/
π
h
δ
3
36
E
/
π
h
δ
3
1
−
ξ
/
δ
Semielliptical [33]
g
6
1
−
ξ
/
δ
2
45
/
2
π
h
δ
3
45
/
2
π
h
δ
3
1
−
ξ
/
δ
2
Parabolic [33]
g
7
1
−
ξ
/
δ
2
2
315
E
/
8
π
h
δ
3
315
E
/
8
π
h
δ
3
1
−
ξ
/
δ
2
2
Six power [35]
g
8
1
−
ξ
/
δ
2
3
945
E
/
16
π
h
δ
3
945
E
/
16
π
h
δ
3
1
−
ξ
/
δ
2
3
Cosine [35]
g
9
cos
π
ξ
/
2
δ
3
E
π
2
/
2
π
2
−
8
h
δ
3
3
E
π
2
/
2
π
2
−
8
h
δ
3
cos
π
ξ
/
2
δ
Inverse proportional [36]
g
10
δ
/
ξ
6
E
/
π
h
δ
3
6
E
/
π
h
δ
3
δ
/
ξ
Quadratic power [36]
g
11
δ
/
ξ
2
3
E
/
π
h
δ
3
3
E
/
π
h
δ
3
δ
/
ξ
2Figure 2
The diagram of various micromodulus functionscξ,δ under the 2D plane stress condition.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 ξ). The 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 Table3


## 3.2. Extended BPD Fracture Criteria

The critical stretchsc is adopted to judge the breakage of the bond, and it can be obtained by the critical energy release rate Gc [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 sc and failure occurs, and these two points cease to interact and the bonds between them to break irreversibly.As shown in Figure3, crack CD crosses the horizon of material point A located on the midperpendicular of CD and divides its horizon into two separate parts completely. The 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.Figure 3
Integration domain of the micropotentials crossed by a crack [11].G
c denotes the energy needed to produce a unit area crack, and in 2D situation, Gc can be expressed by(15)Gc=h∫0δ∫zδ∫−cos−1z/ξcos−1z/ξ12cξ,δsc2ξξdθdξdz=h∫0δ∫zδ∫0cos−1z/ξcξ,δsc2ξξdθdξdz.Substituting from equation (9) into equation (15) results in(16)Gc=h∫0δ∫zδ∫0cos−1z/ξc0,δgξ,δsc2ξξdθdξdz.In the 3D situation,Gc can be expressed as(17)Gc=∫0δ∫02π∫zδ∫0cos−1z/ξ12c0,δgξ,δsc2ξξ2sinθdθdξdϕ.In the present work, substituting from equations (13) and (14) into equation (16), therefore, the expression of scin the plane stress situation can be obtained as(18)sc=32πGc63Ehδ.When applying the extended BPD approach, once the conditions>sc is reached, the bond is deemed broken. In BPD, the bond stiffness is cξ,δ, an invariable material constant, and the expression of sc is given in [8, 9], i.e., sc=4πGc/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(19)μξ,t=1,ifsξ,t′<sc,0<t′<t,0,else.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(20)φX,t=1−∫Hδμt,ξdVX∫HδdVX,0≤φ≤1.Note that0≤φ≤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.


## 4. Numerical Examples

### 4.1. 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 ofp = 200 MPa, as shown in Figure 4. The material density is 7850 kg/m3, Young’s modulus 200 GPa, and Poisson’s ratio ν=1/3. These material parameters are the same as those on page 155 of [11]. The 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.Figure 4
Geometry of a plate under uniaxial tension.The analytical solutions of the displacements are(21)ux∗x,y=pEx,uy∗x,y=−νpEy.And the relative errors of numerical solutions are defined as(22)ex=ux−ux∗ux∗,ey=uy−uy∗uy∗,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 Figure5, where the horizon radius is chosen as δ=3.015Δ. The maximum relative error on the horizontal 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.Figure 5
Contour of relative errors on the horizontal displacement for the BPD model with different attenuation functions. (a) Constant. (b) Four power. (c) Exponential. (d) Gaussian. (e) Triangular. (f) Semielliptical. (g) Parabolic. (h) Six power. (i) Cosine. (j) Inverse proportional. (k) Quadratic power.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)Figure 6
Maximum relative error on the horizontal displacement with different values of the horizon radiusδ calculated by the BPD model with various attenuation functions.It can be observed that the maximum relative errors all occur in the corners due to the surface effect of the PD method. The 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 functiong2=1−ξ/δ4, and the calculation accuracy of the BPD model with this function is the highest regardless of the horizon radius.These 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. The extended BPD approach in the present work is actually the four-power functiong2=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.


### 4.2. 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. The material density is 8000 kg/m3, Young’s modulus 192 GPa, and Poisson’s ratio ν=1/3. These material parameters are the same as those on page 155 of [11]. The plate is subject to a slow rate of stretch along its horizontal edges, and the velocity boundary condition is set as u˙yx,±L/2,t=±2.7541×10−7m/s and imposed on the fictitious boundary layer with a width of 3Δ, as shown in Figure 7. The spacing between material points is Δ = 0.5 mm, and the total number of material actual points is 9876. The horizon radius is chosen as δ=3.015Δ, and time step ∆t = 1.0 s. The critical stretch value sc=0.02138. Deformation and fracture of the plate are simulated by the extended BPD model.Figure 7
Geometry of a plate with a circular cutout under slow stretch and its discretization.The process of fracture initiation and crack propagation is shown in Figure8. 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. The 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.Figure 8
The crack propagation paths of the plate: (a) 605 step; (b) 800 step; (c) 905 step.
(a)
(b)
(c)Figure 9
The crack paths and damage plots of the plate for literature results: (a) PDFEM by Huang [4]; (b) BPD by Madenci [11]; (c) BPD by Ochoa-Ricoux [20].
(a)
(b)
(c)


### 4.3. 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 Figure10. The properties of the plate are specified as Young’s modulus 2.94 GPa, Poisson’s ratio ν = 1/3, mass density 1,200 kg/m3, and the critical stretch value sc=0.09515. The pre-existing crack has an orientation of 62.5° from the horizontal axis. The load is applied, and a no-fail region is introduced within the length of dz = 25 mm and nz = 40 mm from the top and bottom corner of the plate, respectively. The no-fail region is specified to avoid spurious cracking near the region of loading. The load through a constant velocity constraint of u˙  = 1 × 10−6 m/s until damage is captured and constraint of u˙ = 1 × 10−8 m/s to predict crack growth paths, respectively. Simulation is implemented by the extended BPD model with Δ = 2 mm, 5776 points, δ=3.015Δ, and time step t = 1.0E − 8 s.Figure 10
A diagonally loaded 2D plate with a pre-existing crack. (a) Experiment model. (b) Extended BPD model.
(a)
(b)The process of crack propagation of the diagonally loaded square plate specimen simulated by the extended BPD model is shown in Figures11(a)–11(e). At the moment of 25 μs, fracture initiates at the tips of the pre-existing crack and then extends along the horizontal axis. At the moment of 422 μs, the crack reaches the boundaries of the plate. It can be seen that the final crack propagation paths of simulation by the extended BPD model have a good agreement with the results of the experiment by Ayatollahi and Aliha [38], as shown in Figure 11(f).Figure 11
Comparison of the crack propagation paths of the plate: (a)t = 25 μs; (b) t = 235 μs; (c) t = 275 μs; (d) t = 335 μs; (e) t = 422 μs by the present numerical method; (f) the experiment result reported by Ayatollahi and Aliha [38].
(a)
(b)
(c)
(d)
(e)
(f)


## 4.1. 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 ofp = 200 MPa, as shown in Figure 4. The material density is 7850 kg/m3, Young’s modulus 200 GPa, and Poisson’s ratio ν=1/3. These material parameters are the same as those on page 155 of [11]. The 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.Figure 4
Geometry of a plate under uniaxial tension.The analytical solutions of the displacements are(21)ux∗x,y=pEx,uy∗x,y=−νpEy.And the relative errors of numerical solutions are defined as(22)ex=ux−ux∗ux∗,ey=uy−uy∗uy∗,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 Figure5, where the horizon radius is chosen as δ=3.015Δ. The maximum relative error on the horizontal 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.Figure 5
Contour of relative errors on the horizontal displacement for the BPD model with different attenuation functions. (a) Constant. (b) Four power. (c) Exponential. (d) Gaussian. (e) Triangular. (f) Semielliptical. (g) Parabolic. (h) Six power. (i) Cosine. (j) Inverse proportional. (k) Quadratic power.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)Figure 6
Maximum relative error on the horizontal displacement with different values of the horizon radiusδ calculated by the BPD model with various attenuation functions.It can be observed that the maximum relative errors all occur in the corners due to the surface effect of the PD method. The 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 functiong2=1−ξ/δ4, and the calculation accuracy of the BPD model with this function is the highest regardless of the horizon radius.These 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. The extended BPD approach in the present work is actually the four-power functiong2=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.


## 4.2. 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. The material density is 8000 kg/m3, Young’s modulus 192 GPa, and Poisson’s ratio ν=1/3. These material parameters are the same as those on page 155 of [11]. The plate is subject to a slow rate of stretch along its horizontal edges, and the velocity boundary condition is set as u˙yx,±L/2,t=±2.7541×10−7m/s and imposed on the fictitious boundary layer with a width of 3Δ, as shown in Figure 7. The spacing between material points is Δ = 0.5 mm, and the total number of material actual points is 9876. The horizon radius is chosen as δ=3.015Δ, and time step ∆t = 1.0 s. The critical stretch value sc=0.02138. Deformation and fracture of the plate are simulated by the extended BPD model.Figure 7
Geometry of a plate with a circular cutout under slow stretch and its discretization.The process of fracture initiation and crack propagation is shown in Figure8. 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. The 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.Figure 8
The crack propagation paths of the plate: (a) 605 step; (b) 800 step; (c) 905 step.
(a)
(b)
(c)Figure 9
The crack paths and damage plots of the plate for literature results: (a) PDFEM by Huang [4]; (b) BPD by Madenci [11]; (c) BPD by Ochoa-Ricoux [20].
(a)
(b)
(c)


## 4.3. 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 Figure10. The properties of the plate are specified as Young’s modulus 2.94 GPa, Poisson’s ratio ν = 1/3, mass density 1,200 kg/m3, and the critical stretch value sc=0.09515. The pre-existing crack has an orientation of 62.5° from the horizontal axis. The load is applied, and a no-fail region is introduced within the length of dz = 25 mm and nz = 40 mm from the top and bottom corner of the plate, respectively. The no-fail region is specified to avoid spurious cracking near the region of loading. The load through a constant velocity constraint of u˙  = 1 × 10−6 m/s until damage is captured and constraint of u˙ = 1 × 10−8 m/s to predict crack growth paths, respectively. Simulation is implemented by the extended BPD model with Δ = 2 mm, 5776 points, δ=3.015Δ, and time step t = 1.0E − 8 s.Figure 10
A diagonally loaded 2D plate with a pre-existing crack. (a) Experiment model. (b) Extended BPD model.
(a)
(b)The process of crack propagation of the diagonally loaded square plate specimen simulated by the extended BPD model is shown in Figures11(a)–11(e). At the moment of 25 μs, fracture initiates at the tips of the pre-existing crack and then extends along the horizontal axis. At the moment of 422 μs, the crack reaches the boundaries of the plate. It can be seen that the final crack propagation paths of simulation by the extended BPD model have a good agreement with the results of the experiment by Ayatollahi and Aliha [38], as shown in Figure 11(f).Figure 11
Comparison of the crack propagation paths of the plate: (a)t = 25 μs; (b) t = 235 μs; (c) t = 275 μs; (d) t = 335 μs; (e) t = 422 μs by the present numerical method; (f) the experiment result reported by Ayatollahi and Aliha [38].
(a)
(b)
(c)
(d)
(e)
(f)


## 5. 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 Figure12. Material properties are specified as Young’s modulus 192 GPa, mass density 8000 kg/m3, Poisson’s ratio ν = 1/3 (parameters are taken from Madenci and Oterkus [11]), and the critical stretch value of sc=0.04781. The velocity boundary condition is Vt=20m/s. Simulation is implemented by the extended BPD model with Δ = 0.1 mm, δ=3.015Δ, and time step 1.3367E − 8 s.Figure 12
Geometric parameters and loading condition of a plate with centered crack.Adopting the same parameters and loading condition, the investigation is made for different crack widths ofdx = 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.Figure 13
Crack propagation of different crack widthsdx and different orientation angles β.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.The pre-existing crack with the width ofdx = 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 14
Crack propagation paths of the plate whendx = 4 mm; (a) pre-existing crack; (b) t = 5.2 μs; (c) t = 8.2 μs; (d) t = 9.8 μs.
(a)
(b)
(c)
(d)Figure15shows 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 pre-existing crack is usually harder to break than other orientation angles.Figure 15
Fracture initiation time of different crack widthsdx and orientation angles β.


## 6. Conclusion

In this study, an extended BPD approach is proposed through introducing a new attenuation function. The 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 pre-existing cracks have a remarkable influence on the fracture initiation times and the final crack propagation paths.



---

*Source: 9568015-2020-04-07.xml*