1 Introduction

Dynamic fracture is the most fundamental problem in fracture science^[1]. The dynamic fracture behaviors of structure are significantly different from those under an increasing static load, such as the projectile penetration into armor and the fragment collision of aircraft. Such dynamic fracture behaviors are extremely influenced by the material properties and geometrical configuration of the structure, because the propagations of the stress wave vary according to the conditions and may induce different fracture modes. Predictive modelling and simulation of the impact damage such as crack initiation, growth, and branching are not only important for fracture mechanics but also useful to design safety industrial structures.

Though significant efforts have been made to predict the crack initiation and growth in the materials, it is still a major challenge within the framework of classical continuum mechanics. The difficulty inherent in this problem arises from the mathematical formulation, which is a set of partial differential equations. These equations cannot be applied directly across the discontinuities resulting from the material damage and failure, since the required partial derivatives do not exist there. Therefore, any numerical method derived from these equations such as the finite element method (FEM) inherits this difficulty in modelling the cracks. This mathematical breakdown results in an inherent limitation of the classical theory to model the structural failure directly, and leads to an externally supplied law or special techniques in the fracture mechanics to treat the discontinuities, such as cracks, as a pathological situation.

The cohesive zone model (CZM) introduced by Ref.[2] is the major breakthrough in the computational fracture mechanics within the framework of the FEM without an external criterion. The materials and material interfaces are modelled through a traction-separation law, with which the tractions are zero when the opening displacement reaches a critical value. In this approach, the cohesive zone elements (CZEs) are embedded at the edges or the facets of the original finite element mesh, and crack growth occurs only between the bulk elements. Therefore, the crack paths are highly of mesh dependence, and remeshing is required when the crack paths are unknown a priori^[3].

In an effort to resolve these difficulties, the concept of the extended FEM (XFEM) is introduced to model the crack growth within the realm of the finite elements without remeshing^[4-5]. It permits the cracks to propagate through any surface within an element, removing the limitations of the CZM that only along the element boundaries. The XFEM has been successfully used to solve a number of fracture problems. However, it requires external criteria in order to predict the crack growth, and does not predict the experimentally-observed crack propagation speeds unless the material's fracture energy values are modified by a significant factor^[6].

In recent years, a nonlocal continuum model, peridynamics, has been proposed as a viable and efficient numerical method to overcome the weaknesses of the existing methods, in particular the difficulties associated with modelling the crack initiation and growth in solids^[7-8]. In this theory, each infinitesimal unit of the continuum, called particle, interacts with other particles located in its neighborhood through forces, similar to the molecular dynamics theory^[9]. The peridynamic theory is nonlocal since the interaction between the particles extends beyond their immediate neighborhood. The peridynamic model uses an integral representation of the equilibrium equation rather than its differential form. This avoids the computation of the derivatives of the displacement, and therefore the equilibrium equations are valid in the entire domain even in the presence of cracks. When damage emerges in the structures, the interactions between the material points progressively break, and their corresponding contributions in the integral representation are simply removed. This feature allows damage initiation and propagation at the multiple sites with the arbitrary paths inside material without resorting to special crack growth criteria.

With the unique advantage of naturally incorporating discontinuities within a single continuum framework, the peridynamic theory has been utilized successfully in many problems, including the crack propagation in isotropic brittle materials^[10-13], the fracture of concrete structures^[14-15], the damage of composites^[16-19], the failure of ductile materials^[20-22], structural vibration and stability^[23-24], and heat transfer^[25-27].

In this paper, we employ the bond-based peridynamic model to investigate the dynamic fracture problem in the three-point bending beam under impact loading. In fact, the three-point bending beam with offset notch is extensively used to study the mixed â… -â…¡ crack propagation in brittle materials. John and Shah^[28] studied the effect of the loading rate on the fracture modes of the three-point bending beam with analytical and experimental approaches. Yao et al.^[29] used the experimental methods of dynamic caustics and photoelasticity to study the dynamic fracture behavior of the structure. Agwai et al.^[30] studied the crack propagation in a pre-cracked plate with an embedded inclusion by use of peridynamics. Oterkus et al.^[31] verified that the peridynamic approach was able to correctly model and simulate the crack path in a concrete beam with an initial notch under four-point bending conditions. Belytschko et al.^[32] used the element-free Galerkin method to study the mixed-mode dynamic crack propagation in a three-pointbending specimen. In order to induce the failure mode transition in the model, an additional initial crack was introduced at the midspan of the beam which was not necessary in the experiment. In these literatures, the dynamic fracture behaviors of the three-point bending beam are studied sufficiently in experiments. However, exact and reliable numerical methods still need to be developed, especially in predicting the crack mode transition and branching. The nonlocal peridynamic theory is a sound analysis method for such dynamic fracture problems involving complex failure patterns without any external criterion and technique.

The paper is organized as follows. The bond-based peridynamic theory and numerical implementation are briefly reviewed in Section 2. Using the proposed model, the impact damage of the three-point bending beam is solved in Section 3, where the predictions concerning the damage growth (I-II mixed mode) are compared with the experimental results from Ref.[29]. The role played by the location of the pre-exiting notch on the dynamic fracture behaviors in the beam is further investigated, especially the change of damage pattern. In the following section, and the influence of the impact loading on the crack growth and branching is studied by increasing the impact speed. Some concluding remarks are presented in Section 4.

2 Bond-based peridynamic theory

The peridynamic theory is concerned with the physics of a material body at a material point that interacts with all points in a nonlocal manner (see Fig. 1). The state of any material point is determined by its pairwise interaction with the points located within a finite distance, called the horizon, which is symbolized by δ. Any pair of the material points only interact with each other when the distance between them is less than the horizon. In the bond-based peridynamic theory, the equation of motion of a material point at the position x in the reference configuration is

(1)

where H_x is a spherical neighborhood of the given radius δ centered at the point x in the reference configuration, ρ is the mass density, u is the displacement vector, b is the prescribed body force vector density, and f is the pairwise force exerted on the particle x in the peridynamic bond that connects the particle x to x'. The integral equation given in Eq.(1) represents a distinct mathematical system that is not derivable from the classical partial differential equations (PDEs). The deformed position of a given material point is determined by the summation of forces between itself and all other material points with which it interacts.

Within the realm of the peridynamic theory, the material points interact with each other directly through a prescribed vector-valued pairwise force function f, which contains all the constitutive information associated with the material. The forces within the material are treated through the interactions between the pairs of the material points in the continuum. This interaction force can also be viewed as a bond force between the material points. For a micro-elastic material, the pairwise force function is derivable from a scalar-valued function w(η, ξ), called the pairwise potential function or micro-potential in Ref.[8], such that

(2)

where ξ=x'-x is the relative position, and η=u' -u is the relative displacement vector in the reference configuration. An linear micro-elastic material (the force magnitude depends linearly on the relative elongation magnitude) is obtained if we take

(3)

where c(ξ) is the material property in peridynamics, called the micro-modulus function, and s(η, ξ) is the bond stretch (the relative elongation) defined by

(4)

where s is positive when the bond is in tension, and is negative when the bond is in compression.

The corresponding pairwise force is derived from Eqs.(2) and (3) by

(5)

The material property c is determined by equating the peridynamic internal energy of a body to the strain energy density from the classical elasticity theory. For the two-dimensional (2D) peridynamics concerned here,

(6)

where E is the elastic Young's modulus, and ν is the Poisson ratio. In the bond-based peridynamics used in this paper, the particles interact only through a pair-potential. This assumption results in an effective Poisson ratio of 1/3 in the 2D plane stress and 1/4 in the three-dimensional (3D) plane stress for an isotropic and linear micro-elastic material^{[8, 14]}. This limitation can be removed in the state-based peridynamic theory^[33].

When the bond stretch between two material points, i.e., x' and x, exceeds a critical stretch value s₀, failure occurs, and these two points cease to interact with each other (see Fig. 2). This critical stretch s₀ is obtained by equating the work per unit fractured area, for which all the bonds across a surface to the fracture energy G₀ should be broken and completely separated along the surface.

For the isotropic materials with the 2D plane stress condition, the critical stretch value s₀ can be derived by^[11]

(7)

The local damage index at a point can be defined as the ratio of the number of the broken bonds to the total (initial) number of the bonds associated with that point, i.e.,

(8)

where µ(x, ξ, t) is a factor mapping the bond breakage, and it is defined by

(9)

Obviously, the damage is defined by a value between 0 and 1, where 0 indicates that a material point has no damage and 1 indicates complete damage at the material point. Note that a damage index in the value range [0.4, 0.5] may indicate that a fracture surface exists if the damage is localized along a surface (line). In the peridynamic theory, fracture occurs naturally in a peridynamic body as a consequence of the motion equation and the constitutive model. Fractures initiate, grow, turn, branch, and arrest without the need for any externally supplied law or special techniques to control these processes, as required in the traditional models for fracture mechanics. The material defects evolved in any unknown mode in advance allows the peridynamic approach to model the complex patterns of mutually interacting cracks within a body.

The numerical procedure involves discretization of the domain of interest into the subdomains. Each subdomain can be represented as a single collocation point located at the mass center of the subdomain. The governing equations are then rewritten for these points along with a Riemann-sum type approximation over each horizon for the integro-differential equation (see Eq.(1)). This discretization results in a set of algebraic equations with the displacement at the different points as unknowns, i.e.,

(10)

where x_i represents the point in question, and x_j represents the point within the horizon of x_i, N is the number of the subdomains within the horizon of the ith material point, V_j is the volume of the point x_j and is usually represented by a square lattice area Δx² in 2D cases when the model is discretized with uniform grids, and G_ij is the surface correction factor for a peridynamic bond between the material points x_i and x_j, as explained in Ref.[34].

With the given initial and boundary conditions, the displacement u(x_i, t) can be obtained by an explicit time integration method, i.e., the velocity-Verlet algorithm, which is more numerically stable than the central difference formula.

(11)

(12)

(13)

where n is the number of the time step, Δt is the time step size, are the displacement, the velocity, and the acceleration vectors at the point x_i, respectively.

3 Peridynamic results for impact damage in three-point bending beam with offset notch 3.1 Problem setup

We consider the dynamic crack growth experiments conducted by Yao et al.^[29]. In the experiments, a beam specimen of polymethyl methacrylate with the size 220mmÃ-45mm is cut with an initial notch of 6 mm depth at the bottom boundary (see Fig. 3). The relative location of the notch along the span of the beam (denoted by a dimensionless parameter γ) is defined by the offset distance a from the midspan to the notch divided by a half-length of the beam, i.e., γ=2a/L. A drop hammer with the weight m=5 kg is used to introduce a vertical impulsive loading in the center of the top surface with an impact speed v₀=2.3 m/s.

The material properties for the polymethyl methacrylate are Young's modulus E=3.1 GPa, the mass density ρ=1 200 kg/m³, Poisson's ratio ν=0.36, and the fracture energy G₀=160 J/m². Recall that Poisson's ratio is fixed to 1/3 in the 2D peridynamic theory as shown in Eq.(6). Although Poisson's ratio may influence the path of a growing crack, its variation is not expected to play a significant role in the fracture behavior of the brittle material (see Ref.[35] for discussion on this subject). Moreover, Poisson's ratio of the polymethyl methacrylate considered here is close to 1/3, and our primary interest here is to observe the capabilities of a peridynamic model in capturing the evolution of the impact damage in the three-point bending beam. Therefore, in this paper, we use the bond-based peridynamic formulation.

The impact loading causes a compressive wave propagating from the loading location to the interior of the beam, and then reflects from the boundaries. In this period, the beam goes through the elastic deformation to the crack initiation, and growths until the ultimate failure of the sample occurs, i.e., the sample is completely separated into two pieces due to the crack growth. This problem is a long-standing challenge in computational mechanics, since although the classical methods, such as the FEM, are capable of analyzing the elastic response of the structures, they need extra techniques to handle the discontinuity-concerned problems such as crack initiation and propagation. In this study, the dynamic crack growth in the three-point bending beam is modelled with the novel 2D bond-based peridynamic theory. For the convenience of simulation, the following numerical rules are used. A spherical rigid projectile of the same mass and dimension as the hammer used in the experiments is modelled in the simulation. The projectile is placed very close to the upper edge of the beam with the given impact speed so that impact happens in the first few time steps.

In the peridynamic model, a grid spacing of Δx=0.25 mm is used for the discretization, leading to 158\, 400 uniform material points. The time step size Δt=10^-8s is adopted in the computational simulations, which is sufficiently small for all of the cases, based on the stability conditions from Ref.[8]. The horizon size is specified as δ=4Δx in this study. In principle, for homogeneous materials, one needs to use a sufficiently small horizon relative to the geometrical features of the sample and the characteristics of the specific dynamic loading on the sample. Especially, for the fracture simulation with the peridynamic method, the sufficient points within the horizon lead to having peridynamic bonds in many possible directions, allowing for autonomous growth of crack paths in virtually any required direction. However, a smaller grid spacing size and a higher node number in the horizon will pose a significant increase in the computational time. Thus, a balance between the good approximation of the dynamic fracture behavior and the computational efficiency has to be attained. Although we use 2D bond-based peridynamic model here to study the impact damage in the beam. The 3D model is also convenient to apply in repeating the calculation. One just need to add some new layers of material points in the thickness direction and use the corresponding material property c and the critical stretch s₀ in the 3D peridynamics.

3.2 Crack propagation paths with different notch locations

To study the effect of the pre-notch location on the crack propagation path in the three point bending beam, we carry out the peridynamic calculations with the notch locations γ=10/110, 20/110, 30/110, 40/110, and 50/110, respectively, while other parameters are fixed. Figure 4(a) shows the peridynamic predicted crack propagation paths with different initial notch locations. To verify the new numerical approach on such a problem, the experiment results in Ref.[29] are presented simultaneously for comparison (see Fig. 4(b)). We can see that the peridynamic predictions of the crack propagation behaviors agree reasonably well with the experimental result. All these crack paths present a similar S shape curve for the given notch locations. It initiates from the notch tip, and propagates toward the upper edge of the beam. During the process, the crack firstly curves gradually toward the vertical midline of the beam, then deviates from the curve direction, and turns toward the opposite direction until it intersects with the upper edge. Note that no crack paths pass through the midline in these cases. It is also worth mentioning that no rebound phenomenon of the hammer is observed in the simulation no matter where to set the notch.

3.3 Time-evolution of dynamic fracture behavior

It has been pointed out that the crack tip loading correlates with the energy necessary for driving the crack rather than with the local stresses^[36]. Therefore, the strain energy density is used to illustrate the evolution of the damage process and crack tip propagation. Here, the notch location γ=40/110 is chosen for instance. We show ten snapshots of the time evolution of the strain energy profiles in Fig. 5. It is shown that the stress waves induced by the impact loading emanate in the material from the impact site and reflect back at the boundaries. After a long time of stress wave propagation, reflecting, and superposition, the complicated stress distribution is formed in the beam. When the moving waves meet the pre-notch tip and the supports at the bottom boundary, the stress concentrations occur there with the strain energy accumulation. As a results, at the time t=200 µs, the crack begins to grow from the pre-notch tip, and the energy is released instantaneously. While the new crack tip generates, the strain energy concentrates there momentarily to proceed the crack propagation in the following time. The strain energy density contours clearly show the regions where the strain energy density is concentrated and the crack develops in time and space. Obviously, the motion trajectory of the energy concentration forms the crack path eventually. In the simulation of the crack paths with different notch locations, we find that the time of the crack initiation delays while the notch offset distance increases. The reason is that, while the pre-notch moves away from the midspan, it will take more time for the interaction of the stress waves with the crack to accumulate enough strain energy at the notch tip for cracking.

3.4 Dynamic fracture modes under different pre-notch locations

To further study the damage behavior of the beam with offset notch at various locations, we increase the offset distance of the notch, and choose three different cases for the peridynamic calculation, i.e., γ=70/110, 77.65/110, and 80/110. The crack propagation process over time is shown in Fig. 6. Obviously, cracks initiate and grow in the paths no longer similar to each other. It demonstrates that different fracture modes generate as an outgrowth of the peridynamic equations of motion. When the notch position is set at the location γ=70/110, the crack still grows from the notch tip and propagates in the curve similar to the cases studied above. However, when the notch offset value γ increases to 77.65/110, the crack initiates from the notch tip first. The initiation time is a little late than that of γ=70/110, as expected. Note that a new cracking is initiated at the midspan subsequently. The crack growth at the midspan, however, is faster than at the notch, which allows failure to occur at the midspan even when growth begins at a later time. In addition, the crack growth at the notch stops to grow at the time of 560 µs. The new crack advances straightly toward the upper edge and eventually separates the beam into two pieces, which demonstrates that the mixed-mode fracture transfers to the splitting mode under the special notch location. Obviously, for the notch location γ=77.65/110, the dynamic crack propagation behaves more complex than that of the other notch location. If we continue to put the notch location far away from the vertical midline of the beam, we find that the crack does not grow along the introduced notch any more. Instead, the fracture path just propagates along the midline of the beam and only the splitting fracture mode is observed. It is interesting that the splitting-mode crack branches as it approaches the top surface, which seems to prevent the crack path reaching the impacting site.

Based on the numerical investigation, we can get the conclusion as observed in the experiment of Ref.[28] that there is a critical value of the notch location γ for the fracture mode transition in the three-point bending beam. If the notch location γ < 77.65/110, the crack initiates from the offset notch and grows along a curve toward the upper edge, showing a mixed-mode fracture. If the notch position γ > 77.65/110, the damage occurs at the midspan and propagates almost straightly towards the upper edge, showing the splitting mode fracture.

3.5 Damage pattern with increasing impact speed

In what follows, we analyze the trend of the crack propagation in the beam with increasing the speed of the drop hammer. The impact speed v₀=4 m/s is considered here with two different notch locations. The damage contours with the notch location γ=40/110 (see Fig. 7(a)). It can be seen that the crack propagates from the initial notch, which has been shown in the previous cases. However, It is surprising that the crack begins to branch at a certain point when it approaches the top surface. When the notch location moves to γ=77/110, we can see that multiple cracks initiate and propagate near the midline of the beam and crack branching appears (see Fig. 7(b)). Therefore, increasing the impact speed induces more crack nucleation, growth, and branching in the beam. It is worth noting that peridynamics captures the complex branching patterns without resorting to any external criteria that triggers the crack.

4 Conclusions

The peridynamic theory incorporates damage, crack, and fracture as a natural component of material deformation without resorting to any supplemental law or special techniques to treat the discontinuities. This study presents the capability of the peridynamic approach to predict the impact damage in the three-point bending beam with offset notch. The peridynamic results realistically capture the mixed-mode fracture initiation and propagation in the beam with different notch locations, and match remarkably well with the reported experimental observations. Further investigations demonstrate that there is a critical notch location, where the crack growth at the notch suddenly transit to the crack growth at the midspan. Besides, the study of the impact speed on the damage characteristics in the beam shows that the complex crack branching and multiple crack growth occur spontaneously when the impact speed increases.