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.

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 [913].

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 with 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 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 to reflect the spatial distribution of the intensity of long-range forces. At present, the attenuation function 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 , 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 , 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 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.

2. Review of BPD Theory

In BPD theory, a continuum is composed of infinitely many material points. Each material point with a volume of and mass density of is identified by its coordinates in the initial configuration. As shown in Figure 1, each material point is assumed to interact with all other material points within its horizon through fictitious bonds. The horizon is defined as . According to Newton’s second law, the peridynamic equation of motion [7] for material point at time is expressed bywhere is the acceleration of material point at time , is the body force density, and is a pairwise force function that represents the force vector material point which exerts on material point and is expressed bywhere and are the displacements of material points and at time , respectively. denotes the relative position of material points and in the initial configuration. denotes their relative displacement.

According to Silling and Askari [8] formulation, for linear elastic material, is the derivative of microelastic strain density with respect to relative displacement vector .wherewhere is a micromodulus function which represents the bond stiffness, denotes the magnitude of the vector , and is the stretch of the bond, and

Substituting from equation (4) into equation (3) results in

BPD assumes that half of the strain energy density due to the interaction of and all is stored by , and can be expressed as

Normally, 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 , i.e.,

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

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 is reduced to a constant . It causes the calculation error to become larger. When describing the internal length effect of nonlocal long-range forces, attenuation function was introduced to reflect the spatial distribution of the intensity of long-range forces. The micromodulus of BPD can be written aswhere is the bond stiffness function of BPD and 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 . The function is equal to a constant of BPD if . In the plane stress situation, the strain energy density in BPD can also be expressed as

The strain energy density based on the classical continuum mechanics at material point can be written aswhere denotes the strain vector at material point . Under the same conditions, the relative stretch of the bond equals to the strain vector, i.e., .

In this study, the form of attenuation function (four-power 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 in the extended BPD approach and elastic modulus and Poisson’s ratio can be expressed as

Similarly, the relationship between the micromodulus of the extended BPD approach, elastic modulus , 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 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 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 and in the plane stress situation is listed in Table 3. And the shape of proposed micromodulus functions 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 ). 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 Table 3

3.2. Extended BPD Fracture Criteria

The critical stretch is adopted to judge the breakage of the bond, and it can be obtained by the critical energy release rate [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 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. 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.

denotes the energy needed to produce a unit area crack, and in 2D situation, can be expressed by

Substituting from equation (9) into equation (15) results in

In the 3D situation, can be expressed as

In the present work, substituting from equations (13) and (14) into equation (16), therefore, the expression of in the plane stress situation can be obtained as

When applying the extended BPD approach, once the condition is reached, the bond is deemed broken. In BPD, the bond stiffness is , an invariable material constant, and the expression of is given in [8, 9], i.e., . Failure is included in the material response through the scalar-valued function that takes on values of either 1 or 0, which is defined as

Function 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 , 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 of p = 200 MPa, as shown in Figure 4. The material density is 7850 kg/m3, Young’s modulus 200 GPa, and Poisson’s ratio . 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.

The analytical solutions of the displacements are

And the relative errors of numerical solutions are defined aswhere denotes the numerically computed displacement value, denotes the analytical displacement value, and 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 . 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 , , and , respectively.

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 function , 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 function 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 . 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 y and imposed on the fictitious boundary layer with a width of , 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 , and time step ∆t = 1.0 s. The critical stretch value . Deformation and fracture of the plate are simulated by the extended BPD model.

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

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 Figure 10. 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 . 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  = 1 × 10−6 m/s until damage is captured and constraint of  = 1 × 10−8 m/s to predict crack growth paths, respectively. Simulation is implemented by the extended BPD model with  = 2 mm, 5776 points, , and time step t = 1.0E − 8 s.

The process of crack propagation of the diagonally loaded square plate specimen simulated by the extended BPD model is shown in Figures 11(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).

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 Figure 12. 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 . The velocity boundary condition is . Simulation is implemented by the extended BPD model with  = 0.1 mm, , and time step 1.3367E − 8 s.

Adopting the same parameters and loading condition, the investigation is made for different crack widths of  = 0.5 mm, 2 mm, 3 mm, 4 mm, and 8 mm and different orientation angles . 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 , 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 , 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 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 , as crack width increases, the fracture initiation time delays; when , the fracture initiation time advances with crack width increasing; and when , the fracture initiation time is not sensitive to the variation of crack width, and orientation angle of the pre-existing crack is usually harder to break than other 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.

Data Availability

The 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

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This research was supported by the National Natural Science Foundation of China (grant no. 51569004).