© O. Bouaraquia and M. Mabssout, Published by EDP Sciences 2026

1 Introduction

Understanding the behavior of a cracked structure under dynamic loading is a topic of great importance in several engineering fields, such as aeronautics, civil, and mechanical engineering. Dynamic stress intensity factors (DSIFs) are one of the key quantities in dynamic fracture analysis since they quantify the stress fields near the crack tip and serve as indicators of crack initiation and propagation stability. While analytical solutions exist for certain simplified configurations, their applicability is limited to problems with simple geometries and boundary conditions. For this reason, numerous numerical methods have been developed to evaluate DSIFs. The most popular methods are the finite element method (FEM) [1], the boundary element method (BEM) [2], and the extended finite element method (X-FEM) [3].

Despite their success, mesh-based approaches become inefficient or inaccurate when dealing with discontinuities or large deformations. To avoid these difficulties, several mesh-free methods have been proposed to simulate fracture in solids, including the smoothed particle hydrodynamics (SPH) [4], the reproducing kernel particle method (RKPM) [5], and the element-free Galerkin (EFG) method [6]. Nevertheless, all these methods are based on classical continuum mechanics (CCM), whose governing equations rely on local partial derivatives. When discontinuities such as cracks appear, the spatial derivatives of displacements become undefined.

To overcome these limitations, Silling [7] introduced the peridynamics (PD) theory as a nonlocal reformulation of continuum mechanics. This formulation replaces spatial derivatives with integral operators that describe interactions between a material point and its neighbors located within a finite region known as the horizon. The theoretical background and computational implementation of peridynamics have been presented in [8]. Since its introduction, peridynamics has been successfully applied to a variety of solid mechanics problems, including fracture propagation, crack branching in brittle materials [9], and wave propagation in 2D elastic and viscous materials [10,11]. Although peridynamics provides a natural framework for fracture simulation, its nonlocal nature often leads to dispersion effects in dynamic analyses [12]. Wave propagation may become distorted, and phase velocities can vary with frequency, which can alter the physical accuracy of predicted crack behavior under impact loading. This makes it essential to investigate and control wave dispersion before applying peridynamics to dynamic fracture problems. Previous studies [8,12] have shown that dispersion in peridynamics depends on some parameters, such as the horizon size and the discretization spacing. To better represent nonlocal interactions, several kernel or weighted functions have been introduced into the bond-based peridynamics (BBPD) formulation [13]. The appropriate choice of such kernel functions can effectively reduce dispersion and enhance numerical stability. These functions control the contribution of neighboring bonds while maintaining consistency with the continuum limit. The present work aims to develop an efficient numerical model within the BBPD framework to accurately simulate crack propagation in an elastic material subjected to impact loading. However, two main issues arise: (i) the choice of the kernel function to reduce wave dispersion and (ii) the implications of this choice on dynamic fracture, in particular the effects of numerical dispersion on the accuracy of DSIFs evaluation and crack-path prediction. We demonstrate how different kernel functions influence not only wave dispersion but also the accuracy of the near-tip displacement field and, consequently, the accuracy of dynamic fracture simulations.

First, the influence of different kernel functions on the response of elastic media and on the evaluation of DSIFs within the BBPD formulation is investigated. Three kernel functions are implemented to define variable micromodulus distributions. This study is conducted in three main parts. First, the dispersion characteristics are analyzed for 1-dimensional (1D) wave propagation. The results are compared to the classical (non-dispersive) solution. In addition, the study also focuses on the evaluation of dynamic stress intensity factors using the displacement extrapolation method. The obtained numerical results are compared with analytical and reference solutions. Finally, the influence of the kernel functions on dynamic crack propagation is investigated through the simulation of the Kalthoff–Winkler benchmark test.

This paper is structured as follows. The BBPD model is described in Section 2. The numerical model is presented in Section 3. The numerical examples are discussed in Section 4, and the conclusions are summarized in Section 5.

2 Bond-based peridynamics model

2.1 Bond-based peridynamics formulation

Peridynamics theory, introduced by Silling [7], is based on a nonlocal interaction concept where neighboring pairs of particles exchange internal forces within a specific radius δ called “horizon”. According to BBPD, the spatial derivatives are replaced by an integral term to describe interactions between particles. Therefore, the peridynamics equation of motion is given by

$$\rho \ddot{\mathit{u}}(\mathit{x},t)={\int }_{H_{\mathbf{x}}}\mathit{f}({\mathit{x}}^{\mathrm{\prime }}-\mathit{x},\mathit{u}({\mathit{x}}^{\mathrm{\prime }},t)-\mathit{u}(\mathit{x},t))d{V}_{{\mathbf{x}}^{\mathrm{\prime }}}+\mathit{b}(\mathit{x},t),$$(1)

where ρ represents the density of a material point x, u(x, t) corresponds to the displacement at point x at time $t,\ddot{\mathit{u}}(\mathit{x},t)$ is the acceleration, f is the pairwise force function, b is the body force, and ℋ_x is the neighborhood of x defined by

$${\mathcal{H}}_x=\{{\mathit{x}}^{\mathrm{\prime }}\in \mathrm{\Omega }:\parallel {\mathit{x}}^{\mathrm{\prime }}-\mathit{x}\parallel \le \delta \},$$(2)

where δ is the horizon, i.e., the radius of the compact domain, and ∥.∥ is the Euclidean norm. Figure 1 depicts the undeformed and deformed configurations of two neighboring peridynamics nodes.

The relative position between two particles x and x^′ in peridynamics notation is

$$\xi =x^{\mathrm{\prime }}-x,$$(3)

and the relative displacement is denoted as

$$\mathit{\eta }=\mathit{u}({\mathit{x}}^{\mathrm{\prime }},t)-\mathit{u}(\mathit{x},t).$$(4)

For a microelastic material, the pairwise force function is obtained as the derivative of the micropotential energy with respect to the bond deformation [7]:

$$\mathit{f}(\mathit{\eta },\mathit{\xi })=\frac{\mathrm{\partial }\omega (\mathit{\eta },\mathit{\xi })}{\mathrm{\partial }\mathit{\eta }}.$$(5)

The linear micropotential function is given by [7,14]:

$$\omega (\mathit{\eta },\xi )=\frac{1}{2}c(\xi ,\delta )s^2\xi ,$$(6)

where c(ξ, δ) is the micromodulus function that represents the bond stiffness, and s is the relative elongation of a peridynamics bond defined as

$$s=\frac{||\mathit{\xi }+\mathit{\eta }||-||\mathit{\xi }||}{||\mathit{\xi }||}.$$(7)

Silling [7] introduced the pairwise force expression for a linear elastic material in BBPD approach by

$$f(\mathit{\eta },\xi )=\{\begin{array}{ll}c(\xi ,\delta )s\frac{\mathit{\xi }+\mathit{\eta }}{||\mathit{\xi }+\mathit{\eta }||}& \text{ if }||\xi ||\le \delta ,\\ 0& \text{ if }||\xi ||>\delta ,\end{array}$$(8)

where ||ξ + η|| is the magnitude of the bond in its deformed configuration.

A kernel function W^¯ (ξ, δ) is added to the equation of the micromodulus function in order to account for the effect of particle distance on the bond’s stiffness, as follows:

$$c(\mathit{\xi },\delta )=c_0(\delta )\overline{W}(\mathit{\xi },\delta ),$$(9)

where c₀(δ) is the micromodulus constant. In the classical BBPD, the kernel function is given by

$$\overline{W}(\xi ,\delta )=\{\begin{array}{ll}1& \text{ for }\xi \le \delta ,\\ 0& \text{ for }\xi >\delta .\end{array}$$(10)

Various kernel functions have been used in the BBPD method [9,14,15]. Nonetheless, a function must satisfy specific criteria to be utilized as a kernel function in the BBPD formulation [15].

$$\begin{array}{ll}W(\xi ,\delta )& \ge 0,\\ {\displaystyle \underset{\xi \to 0}{lim}W(\xi ,\delta )}& =maxW,\\ W(\xi ,\delta )& =W(-\xi ,\delta ),\\ W(\xi ,\delta )& =0\text{ when }|\xi |>\delta ,\\ {\displaystyle \underset{\delta \to 0}{lim}W(\xi ,\delta )}& =\mathrm{\Delta }(\xi ),\end{array}$$(11)

where Δ(ξ) is the Dirac delta function.

In this work, three kernel functions are used to evaluate their effect on the BBPD simulation, including the constant function, the Gaussian distribution function, and the triangular distribution function. Table 1 lists the kernel functions used and the corresponding micromodulus functions.

The micromodulus c₀(δ) can be determined by equating the strain energy density obtained from the peridynamics theory and the classical continuum mechanics (CCM) [16].

In bond-based peridynamics, the strain energy density is given by

$$V^{\mathrm{P}\mathrm{D}}(\mathit{x})=\frac{1}{2}{\int }_{H_x}\frac{c{s}^2}{2}\xi d{V}_{\mathit{x}}.$$(12)

In a 2D problem, the strain energy density based on the CCM is given by

$$V^{\mathrm{C}\mathrm{C}\mathrm{M}}(\mathit{x})=\frac{E}{2(1-v^2)}[{\epsilon }_{11}^2+{\epsilon }_{22}^2+2v{\epsilon }_{11}{\epsilon }_{22}+2(1-v){\epsilon }_{12}^2],$$(13)

where ε_ij represents the components of the strain tensor, E is the Young’s modulus, and v is the Poisson’s ratio of the material.

The equivalence between the BBPD strain energy density V^PD(x) and the classical elastic strain energy density V^CCM(x) is valid for an infinite homogeneous elastic body. In particular, the material point x must be located in the interior of the body such that dist(x, ∂Ω) > δ, the deformation field must be homogeneous within the horizon H_x, and the horizon size δ must be sufficiently small.

The calibration of the micromodulus function c(δ) based on two independent homogeneous deformation modes (isotropic tensile deformation and pure shear deformation) imposes a constraint on the Poisson’s ratio, yielding v = 1/3 for 2D plane stress and v = 1/4 for 2D plane strain conditions [7,8]. Under these assumptions, the equality V^CCM = V^PD leads to a unique expression of c₀(δ) for each kernel. The shapes of the proposed kernel functions are illustrated in Figure 2. Table 1 lists the kernel function used and the corresponding micromodulus functions in 1D and 2D plane stress/strain problems. Note that in the 1D case, A denotes the constant cross-sectional area, whereas in 2D, h represents the plate thickness.

Fig. 1 Reference and current configurations of bond-based peridynamics.

2.2 Bond failure criterion

The fracture criterion adopted in this work is based on the critical energy release rate. In BBPD, fracture is modeled through the irreversible breakage of bonds between interacting material points. A bond connecting two points x and x^′ is considered broken when the stretch s_xx′, exceeds a critical value s₀. Once broken, the bond no longer contributes to the internal force interactions. The critical stretch s₀ is derived from the material fracture energy G_c by equating the energy required to break all bonds crossing a unit fracture surface to G_c [8,17]. For the present 2D plane stress formulation, the resulting expressions depend on the selected kernel function.

Constant kernel:

$$s_0=\sqrt{\frac{4G_c}{h{c}_0{\delta }^4}}.$$(14)

Triangular kernel:

$$s_0=\sqrt{\frac{5\pi G_c}{9E\delta }}.$$(15)

Gaussian kernel:

$$s_0=\sqrt{\frac{32G_c}{h{c}_0{\delta }^4(1-5e^{-4})}}.$$(16)

To quantify material degradation, a local damage index ϕ(x, t) is introduced and defined as the ratio of broken bonds to the initial number of bonds within the horizon:

$$\varphi (x,t)=1-\frac{{\int }_{H_x}\mu (x,x^{\mathrm{\prime }},t)d{V}_{x^{\mathrm{\prime }}}}{{\int }_{H_x}d{V}_{x^{\mathrm{\prime }}}},$$(17)

where μ(x, x^′, t) is a history-dependent scalar function equal to 1 if the bond is intact and 0 if the bond is broken (i.e., s > s₀). The damage index ranges from 0 (undamaged material) to 1 (fully broken material).

Fig. 2 Constant, triangular, and Gaussian kernel functions.

3. Numerical model

3.1 Spatial discretization

A layer of material points was used to discretize a 1D bar and 2D plate, as seen in Figure 3. The mass density and volume of each material point are known in the reference configuration. The distance between material points, Δ, represents the grid spacing. Finding the ideal values for the grid size Δ and horizon size δ was crucial in order to get high accuracy with a relatively minimal calculation time, as discussed by Madenci and Oterkus [8].

The spatial discretization of Equation (1) is written as

$$\rho {\ddot{\mathbf{u}}}_i^n=\sum _{j=1}^{N_{Hx}}\mathbf{f}({\mathbf{x}}_j-{\mathbf{x}}_i,{\mathbf{u}}_j^n-{\mathbf{u}}_i^n)C_{v_j}{V}_j+{\mathbf{b}}_i^n,$$(18)

where N_ℋx denotes the total number of material points located within the influence zone of particle i. The term V_j refers to the volume assigned to the jth particle, and C_vj is a volume correction coefficient used to estimate the portion of V_j that lies within the neighborhood ℋ_x (Fig. 3), as described in [18].

$$C_{vj}(\xi )=\{\begin{array}{ll}1& \text{ for }||\xi ||\le \delta -0.5\mathrm{\Delta }x,\\ \frac{\delta +0.5\mathrm{\Delta }x-||\xi ||}{\mathrm{\Delta }x}& \text{ for }\delta -0.5\mathrm{\Delta }x<||\xi ||<\delta +0.5\mathrm{\Delta }x,\\ 0& \text{ for }||\xi ||\ge \delta +0.5\mathrm{\Delta }x.\end{array}$$(19)

Equation (18) can be written as

$$\rho {\ddot{u}}_i^n=\sum _{j=1}^N{f}_{ij}^n{C}_{v_{ij}}{V}_j+b_i^n,$$(20)

where the pairwise interaction force is defined as

$$f_{ij}^n=f({\xi }_{ij},{\mathit{\eta }}_{ji}^n)=c{s}_{ij}^n\frac{{\xi }_{ij}+{\mathit{\eta }}_{ij}^n}{||{\xi }_{ij}+{\mathit{\eta }}_{ij}^n||},$$(21)

with ξ_ij = x_j−x_i and ${\mathit{\eta }}_{ij}^n={\mathbf{u}}_j^n-{\mathbf{u}}_i^n$ denoting the relative position and relative displacement vectors, respectively. The stretch between two material points is expressed as

$$s_{ij}^n=\frac{||{\xi }_{ij}+{\mathit{\eta }}_{ij}^n||-||{\xi }_{ij}||}{||{\xi }_{ij}||}.$$(22)

Fig. 3 Peridynamics 1D and 2D discretization.

3.2 Time discretization

The time integration of the semi-discretized form (Eq. (18)) is obtained using a central finite difference scheme. The acceleration at time tⁿ = nΔt is given by

$${\ddot{\mathit{u}}}_i^n=\frac{{\mathit{u}}_i^{n+1}-2{\mathit{u}}_i^n+{\mathit{u}}_i^{n-1}}{\mathrm{\Delta }{t}^2}.$$(23)

Inserting Equation (23) into Equation (18) gives the displacement at time tⁿ⁺¹ as

$${\mathit{u}}_i^{n+1}=\frac{\mathrm{\Delta }{t}^2}{\varrho }[\sum _{j=1}^N{\mathit{f}}_{ij}^n{C}_{vj}{V}_j+{\mathit{b}}_i^n]+2{\mathit{u}}_i^n-{\mathit{u}}_i^{n-1}.$$(24)

Stability condition: Using the von Neumann stability analysis of Equation (24) written in 1D, the following stability condition is obtained [17]:

$$\mathrm{\Delta }t\le \mathrm{\Delta }{t}_{\mathrm{c}\mathrm{r}}=\sqrt{\frac{\varrho }{\sum _j{V}_j{C}_{ij}}},$$(25)

where V_j is the volume associated with node j, and C_ij is a scalar coefficient defined as

$$C_{ij}=\left|C\left({\xi }_{ij}\right)\right|=\frac{c\left({\xi }_{ij}\right)}{||{\xi }_{ij}||}.$$(26)

4 Numerical examples

4.1 Wave dispersion analysis in a 1D bar

Since peridynamics is a non-local continuum model, non-linear wave dispersion behavior is induced by the direct interactions of material points at a finite distance [12]. In this section, we focus on the effect of the kernel function on wave dispersion.

Neglecting the body forces, the equation in 1D BBPD is given as

$$\rho \ddot{u}(x,t)={\int }_{{\mathcal{H}}_x}{c}_0(\delta )\overline{W}(\xi ,\delta )\frac{u^{\mathrm{\prime }}-u}{|x^{\mathrm{\prime }}-x|}d{x}^{\mathrm{\prime }}.$$(27)

The dispersion relation of an elastic medium can be obtained by examining a single harmonic wave component, represented as

$$u(x,t)=u_0{e}^{i(kx-\omega t)},$$(28)

where ω is the angular frequency, i is the imaginary unit, k is the wave number, and u₀ is the amplitude. Substituting Equation (28) into Equation (27) leads to the following dispersion relation:

$${\omega }^2={\int }_{{\mathcal{H}}_x}\frac{c_0(\delta )\overline{W}(\xi ,\delta )}{\rho |\xi |}(1-\cos k\xi )d\xi .$$(29)

Equation (29) represents the continuous dispersion relation, which can be obtained analytically for each kernel function through direct integration. However, in practical peridynamics simulations, the equations are solved in discrete form. Therefore, the discrete dispersion relation is also presented for different kernel functions. The effects of the grid spacing Δx, the finite number of interacting neighbors within the horizon, and the choice of kernel function are also investigated.

In the discrete form, the horizon δ and the relative position ξ are given by

$$\delta =m\mathrm{\Delta }x,\xi =j\mathrm{\Delta }x.$$(30)

The discretized form for the dispersion relation ω(k) is given by

$${\omega }^2=\sum _{\begin{array}{c}j\ne 0\\ j\in [-m,m]\end{array}}\frac{c_0(m\mathrm{\Delta }x)\overline{W}(j\mathrm{\Delta }x,m\mathrm{\Delta }x)}{\rho |j|}(1-\cos (jk\mathrm{\Delta }x).$$(31)

For the constant kernel function, the dispersion relation is as follows [16]:

$${\omega }^2=\frac{2v^2}{m^2\mathrm{\Delta }{x}^2}\sum _{\begin{array}{c}j\ne 0\\ j\in [-m,m]\end{array}}\frac{1}{|j|}(1-\cos (jk\mathrm{\Delta }x)).$$(32)

The dispersion expression for the Gaussian kernel is given by [16]

$${\omega }^2=\frac{8v^2{e}^4}{m^2\mathrm{\Delta }{x}^2(e^4-1)}\sum _{\begin{array}{c}j\ne 0\\ j\in [-m,m]\end{array}}\frac{1}{|j|}(1-\cos (jk\mathrm{\Delta }x))e^{-{\left(2\frac{j}{m}\right)}^2}.$$(33)

When the triangular kernel is used, the dispersion relation is written as

$${\omega }^2=\frac{6v^2}{m^2\mathrm{\Delta }{x}^2}\sum _{\begin{array}{c}j\ne 0\\ j\in [-m,m]\end{array}}\frac{1}{|j|}(1-\cos (jk\mathrm{\Delta }x))(1-\frac{|j|}{m}),$$(34)

where v is the classical wave speed given by $v=\sqrt{E/\rho }$.

Figure 4 shows a comparison of the dispersion curves for the three kernel functions. Parameters m = 5 and δ = 1 mm are used in this analysis. As shown, the kernel function significantly influences the dispersion curves. Dispersion becomes more pronounced as the wave number increases. The Gaussian kernel exhibits less dispersion than the constant and triangular kernel functions.

The phase velocities c(k) and the group velocities v(k) are obtained from the dispersion relation ω(k) as

$$c(k)=\frac{\omega (k)}{k},v_g(k)=\frac{d\omega (k)}{dk}.$$(35)

For the constant kernel, the group velocity is given by

$$v_g(k)=\frac{v}{m\sqrt{2}}\cdot \frac{{\displaystyle \sum _{\frac{j\ne 0}{j\in [-m,m]}}\frac{j}{|j|}\sin (jk\mathrm{\Delta }x)}}{\sqrt{{\displaystyle \sum _{\frac{j\ne 0}{j\in [-m,m]}}\frac{1}{|j|}(1-\cos (jk\mathrm{\Delta }x))}}}.$$(36)

The group velocity for the Gaussian kernel is given by

$$v_g(k)=\frac{c\mathrm{\Delta }x}{2\omega (k)}\sum _{\begin{array}{c}j\ne 0\\ j\in [-m,m]\end{array}}\frac{j}{|j|}\sin (jk\mathrm{\Delta }x)e^{-{\left(\frac{2j}{m}\right)}^2},$$(37)

where $c=\frac{8v^2{e}^4}{m^2\mathrm{\Delta }{x}^2(e^4-1)}$.

The group velocity for the triangular kernel is written as

$$v_g(k)=\frac{3v^2}{\omega (k)m^2\mathrm{\Delta }x}\sum _{\begin{array}{c}j\ne 0\\ j\in [-m,m]\end{array}}\frac{j}{|j|}\sin (jk\mathrm{\Delta }x)(1-\frac{|j|}{m}).$$(38)

Figure 5 illustrates the normalized phase velocity versus the wave number k. All kernel functions start at unity for small wave numbers (k → 0), showing consistency with the classical local model. As k increases, the velocity decreases due to dispersion, indicating that the waves travel at lower velocities than in classical elasticity. Compared with the constant and triangular kernels, the Gaussian kernel introduces the least dispersion.

Figure 6 illustrates the evolution of the normalized group velocity as a function of the wave number k for the three kernel functions. The exact solution corresponds to the classical local elastic model for which all frequency components propagate at the same speed. Again, the constant kernel shows the most dispersive response, with v_g decreasing rapidly and showing oscillations and negative values at higher wave numbers. The Gaussian kernel produces the least dispersion and guarantees consistent energy propagation.

The obtained dispersion curves are in good agreement with those reported in [19], which analyzed seven kernel functions, including the Gaussian kernel used in the present study. For the non-singular kernels examined in [19], strong similarities between the two sets of results are observed. In particular, the Gaussian kernel exhibits lower dispersion than the other kernel functions. Moreover, the negative group velocities observed in Figure 6 are also reported in [19].

These results emphasize the importance of selecting an appropriate kernel function to control wave dispersion in BBPD, with the Gaussian kernel emerging as the most suitable choice for simulations involving dynamic crack propagation.

The dispersion behavior discussed above has direct consequences for dynamic fracture simulations. The computation of dynamic stress intensity factors and the simulation of crack propagation under impact are strongly influenced by wave dispersion. The accuracy of the near-tip displacement field depends on correct wave propagation from the loading boundaries to the crack front, and excessive numerical dispersion can compromise both DSIF estimation and crack growth predictions. The dispersion analysis showed that the Gaussian kernel performs better in reducing numerical wave dispersion. We next assess its performance in computing DSIFs and simulating crack propagation in comparison with the other two kernels.

Fig. 4 Dispersion curves (m = 5, δ = 1 mm).

Fig. 5 Normalized phase velocity for each kernel function.

Fig. 6 Normalized group velocity for the three kernel functions.

4.2 Evaluation of dynamic stress intensity factors

DSIFs are among the most important parameters in dynamic fracture analysis. To evaluate the influence of the kernel function on the evaluation of the dynamic stress intensity factors, two examples of Mode-I and Mode-II are presented using the BBPD formulation. Several methods have been proposed to calculate DSIFs, such as the path-independent J-integral [20], and the extrapolation method. In the present work, the displacement extrapolation method is applied to estimate the DSIFs [16].

DSIFs for Mode-I and Mode-II are given by [16]

$$K_I(t)=\frac{2\mu \sqrt{2\pi }}{\kappa +1}\frac{u_2^B(t)\frac{r_2}{\sqrt{r_1}}-u_2^D(t)\frac{r_1}{\sqrt{r_2}}}{r_2-r_1},K_{II}(t)=\frac{2\mu \sqrt{2\pi }}{\kappa +1}\frac{u_1^B(t)\frac{r_2}{\sqrt{r_1}}-u_1^D(t)\frac{r_1}{\sqrt{r_2}}}{r_2-r_1},$$(39)

where r_i denotes the radial coordinate with the origin at the crack tip, μ is the shear modulus at the crack tip, and u₁(r, π), and u₂(r, π) are the crack surface displacements computed using the BBPD model.

And,

$$\kappa =\{\begin{array}{ll}3-4v,& \text{for plane strain},\\ {\displaystyle \frac{3-v}{1+v},}& \text{for plane stress}.\end{array}$$(40)

The two-point extrapolation method assumes a linear variation of $(u/\sqrt{r})$ near the crack tip [16]. While higher-order (quadratic or cubic) approximations could improve accuracy, the chosen points lie in a region where higher-order terms are negligible. Points B and D are selected along the crack faces at distances r₁ = Δx and r₂ = 2Δx from the crack tip, respectively (Fig. 7), to remain sufficiently close to the crack tip to capture the near-tip asymptotic behavior.

Fig. 7 Nodes close to the crack tip.

4.2.1 Stationary Mode-I semi-infinite crack

The displacement extrapolation method combined with BBPD is used to compute DSIFs for a semi-infinite lateral crack in an infinite plate (Fig. 8). Three kernel functions are adopted to evaluate their influence on the proposed method; the numerical results are compared to the analytical solution [21].

The plate dimensions are length L = 10 m, height 2 H = 4 m, and crack length a = 5 m. The material properties are Young’s modulus E = 210 GPa, Poisson’s ratio v = 0.3, and density ρ = 8000 kg m⁻³. The time step used in the computation is Δt = 5 × 10⁻⁶ s.

The ramp loading applied traction:

$$\sigma (t)=\{\begin{array}{ll}{\sigma }_0\frac{t}{t_d}& \text{ for }t\le t_d,\\ {\sigma }_0& \text{ for }t\ge t_d,\end{array}$$(41)

where t_d is the rise time of the ramp load. In this example, t_d = 3.36 × 10⁻⁴ s.

The analytical solution is given by [21]

$$K_I(\overline{t})=\{\begin{array}{ll}\frac{4}{3\pi }\frac{\sqrt{1-2v}}{1-v}{\sigma }_0\sqrt{\pi C_d{t}_d}{\left(\frac{\overline{t}}{t_d}\right)}^{\frac{3}{2}}& \text{ for }0\le \overline{t}\le t_d,\\ \frac{4}{3\pi }\frac{\sqrt{1-2v}}{1-v}{\sigma }_0\sqrt{\pi C_d{t}_d}[{\left(\frac{\overline{t}}{t_d}\right)}^{\frac{3}{2}}-{(\frac{\overline{t}}{t_d}-1)}^{\frac{3}{2}}]& \text{ for }\overline{t}\ge t_d,\end{array}$$(42)

where $\overline{t}=t-t_c,t_c=\frac{H}{C_d}$ is the time required for the stress wave to propagate to the tip of the crack, and C_d is the velocity of the longitudinal wave.

The normalized DSIF is given by

$${\tilde{K}}_I(t)=\frac{K_I(t)}{{\sigma }_0\sqrt{H}}.$$(43)

The comparison with the analytical solution is valid only for t ≤ 3t_c that corresponds to the time required for the stress wave to reflect from the lower boundary and reach the crack tip again.

To quantify the accuracy of the numerical results, two error metrics are used. The peak relative error is defined as:

$${\epsilon }_{\text{peak }}=\frac{|K_{max}^{\mathrm{P}\mathrm{D}}-K_{max}^{\mathrm{e}\mathrm{x}\mathrm{a}\mathrm{c}\mathrm{t}}|}{K_{max}^{\mathrm{exact}}}.$$(44)

The L² relative error norm is defined as

$${\epsilon }_{L^2}=\sqrt{\frac{\sum _{n=1}^N{(K_n^{\mathrm{P}\mathrm{D}}-K_n^{\text{exact }})}^2}{\sum _{n=1}^N{\left(K_n^{\text{exact }}\right)}^2}},$$(45)

where $K_n^{\mathrm{P}\mathrm{D}}$ denotes the peridynamics DSIFs at time t_n = nΔt and $K_n^{\text{exact }}$ denotes the exact solution computed at t_n.

Convergence study:

A convergence study of the proposed method is conducted using the Gaussian kernel. For a fixed ratio m = δ/Δx = 5, three horizon sizes are considered: δ = 0.5 m (Δx = 0.1 m, 100 × 40 particles), δ = 0.25 m (Δx = 0.05 m, 200 × 80 particles), and δ = 0.125 m (Δx = 0.025 m, 400 × 160 particles), as shown in Figure 9.

The results presented in Table 2 show that a significant improvement is achieved when refining from Δx = 0.1 m to Δx = 0.05 m, with the peak error decreasing from 9.82% to 0.21% and the L² error from 17.70% to 2.86%. Further refinement to Δx = 0.025 m yields comparable accuracy (peak error 0.17%, L² error 2.74%), confirming the convergence of the numerical solution.

Next, we examine the influence of the kernel function on the evaluation of DSIFs using the following parameters: Δx = 0.05 m and m = 5 (Fig. 10).

As shown in Table 3, the peak relative errors are 14.20% for the constant kernel, 5.99% for the triangular kernel, and 0.21% for the Gaussian kernel. The corresponding L² relative error norms are 15.48%, 6.56%, and 2.86%, respectively. These results indicate that the Gaussian kernel yields the lowest error levels among the three kernel functions considered.

Fig. 8 Geometry and loading of the semi-infinite plate.

Fig. 9 Normalized Mode-I DSIF for three grid refinement using the Gaussian kernel.

Fig. 10 Normalized Mode-I DSIFs with different kernel functions.

4.2.2 Mode-II example

In this example, various kernel functions are used to evaluate the Mode-II DSIFs of a cracked rectangular plate subjected to a uniform shear load, as illustrated in Figure 11.

The plate dimensions are length L = 104 mm, height H = 40 mm, and crack length 2a = 24 mm. The material properties are Young’s modulus E = 73.5 GPa, Poisson’s ratio v = 0.25, and density ρ = 2450 kg m⁻³.

Figure 12 presents the normalized Mode-II DSIFs, ${\tilde{K}}_{II}(t)=\frac{K_{II}(t)}{\tau \sqrt{\pi a}}$, for the kernel functions used. The K_II(t) values are obtained using the extrapolation technique described in Equation (39) combined with the BBPD formulation. For comparison, results computed using the J-integral method [22] are also plotted. A good agreement is observed, particularly for the Gaussian kernel, which closely matches the reference results reported in [22]. The comparison in Figure 12 is restricted to times prior to the arrival of reflected stress waves from the boundaries at the crack tip in order to avoid any boundary reflection effects.

Fig. 11 Plate under shear loading.

Fig. 12 Normalized Mode-II DSIFs for different kernel functions.

4.2.3 Kalthoff-Winkler’s experiment

In the previous sections, it was shown that the kernel functions directly affect the wave dispersion characteristics. The analysis of the DSIFs also demonstrated that the kernel influences the near-tip stress and the concentration of strain energy. This section investigates the impact of the kernel function on the dynamic crack propagation. To evaluate how the kernel functions affect dynamic crack propagation, the classical Kalthoff-Winkler experiment [23] is simulated using the BBPD formulation. This benchmark is commonly used to validate numerical models for dynamic fracture. Figure 13 shows the geometry and the loading conditions.

The material properties are E = 190 GPa, v = 0.25, ρ = 7800 kg m⁻³, and G_c = 6.9 × 10⁴ J m⁻². The grid space Δ = 10⁻³ m is used with a uniform distribution of material points. A value of δ = 5Δ was selected for the horizon. The time step used is Δt = 8.7 × 10⁻⁸ s.

The impact loading is imposed as a prescribed velocity v₀ = 32 m s⁻¹ to a layer of material points of thickness δ along the left boundary between the two notches, thereby simulating the projectile impact. The two pre-notches are introduced by initially breaking all bonds intersecting the notch lines during the model initialization. The critical stretch s₀ for each kernel function is computed from the fracture energy G_c.

Figure 14 presents crack path propagation for the three kernel functions at 70 μs. For the constant kernel, the crack initial time is 20 μs. The triangular kernel slightly delays the crack initiation to 20.8 μs, while the Gaussian kernel delays initiation even further to 21.7 μs.

The crack propagation paths predicted by the three kernel functions are compared with the experimental trajectory in Figure 15. As can be observed, the Gaussian kernel provides the best agreement with the experimental result [23], where the crack propagated at an angle close to 70^∘.

Fig. 13 Geometric and loading conditions of Kalthoff–Winkler experiment.

Fig. 14 Crack propagation path: (a) constant kernel, (b) triangular kernel, and (c) Gaussian kernel.

Fig. 15 Comparison between numerical crack paths for different kernel functions and the experimental result [23].

5 Conclusion

In this study, three kernel functions were used to examine their influence on wave dispersion, as well as on the evaluation of dynamic stress intensity factors and crack propagation within the bond-based peridynamics framework. First, the dispersion characteristics of each kernel function were evaluated, and then DSIFs were computed and compared with analytical and reference solutions. The Gaussian kernel produced the smallest dispersion errors and showed the best agreement with both the evaluation of DSIFs and the simulation of the Kalthoff–Winkler experiment.

These results demonstrate that the kernel selection has a direct impact on both wave propagation and near-tip fields and therefore plays a critical role in the accuracy of peridynamics simulations involving dynamic fracture.

Finally, it should be noted that the BBPD formulation inherently imposes a restriction on Poisson’s ratio. As a result, materials with arbitrary Poisson’s ratios cannot be accurately represented within the bond-based framework. To overcome this limitation, the state-based peridynamics formulation can be employed, as it removes this constraint and allows for a more general constitutive description. Future work will therefore focus on extending the present study to the state-based peridynamics framework.

Funding

The authors declare that this research received no external funding and that no funds were used to cover publication costs.

Conflicts of interest

All authors certify that they have no financial or non-financial conflicts of interest related to this article.

Data availability statement

The data supporting the findings of this study are not available.

Author contribution statement

Oussama Bouaraquia: Conceptualization, Methodology, Software, Validation, Formal Analysis, Investigation, Writing, Visualization.

Mokhtar Mabssout: Conceptualization, Methodology, Software, Validation, Formal Analysis, Investigation, Data Curation, Writing—Review & Editing, Visualization, Supervision.

References

All Tables

All Figures

	Fig. 1 Reference and current configurations of bond-based peridynamics.
In the text

	Fig. 2 Constant, triangular, and Gaussian kernel functions.
In the text

	Fig. 3 Peridynamics 1D and 2D discretization.
In the text

	Fig. 4 Dispersion curves (m = 5, δ = 1 mm).
In the text

	Fig. 5 Normalized phase velocity for each kernel function.
In the text

	Fig. 6 Normalized group velocity for the three kernel functions.
In the text

	Fig. 7 Nodes close to the crack tip.
In the text

	Fig. 8 Geometry and loading of the semi-infinite plate.
In the text

	Fig. 9 Normalized Mode-I DSIF for three grid refinement using the Gaussian kernel.
In the text

	Fig. 10 Normalized Mode-I DSIFs with different kernel functions.
In the text

	Fig. 11 Plate under shear loading.
In the text

	Fig. 12 Normalized Mode-II DSIFs for different kernel functions.
In the text

	Fig. 13 Geometric and loading conditions of Kalthoff–Winkler experiment.
In the text

	Fig. 14 Crack propagation path: (a) constant kernel, (b) triangular kernel, and (c) Gaussian kernel.
In the text

	Fig. 15 Comparison between numerical crack paths for different kernel functions and the experimental result [23].
In the text