1   Introduction

Fracture and failure are common phenomena in nature and are usually inevitable for nearly all materials and structures. As a sophisticated and cutting-edge research topic in solid mechanics, fracture analysis plays an important role in the design and manufacture of various engineering applications.

With the rapid development of computing technologies, considerable efforts have been devoted to the numerical simulation of damage and defects. Chan et al. (1970) employed the finite element method (FEM) to evaluate stress intensity factors (SIFs) for models of different shapes under various types of loading conditions. Fedelinski et al. (1994) applied the dual boundary element method to study dynamic SIFs in stationary cracks. Belytschko and Black (1999) introduced discontinuous enrichment functions into FEM for approximating crack growth, giving rise to the extended FEM. By employing B-spline functions, Tanaka et al. (2012; 2015) developed a wavelet Galerkin method (WGM) to evaluate stress concentration.

The aforementioned methods are developed based on classical continuum mechanics (CCM), employing spatial partial differential equations (PDEs) to describe motion. However, damage or cracks are inherently discontinuous, which may impede efficient numerical evaluation. Additionally, dynamic fracture behavior, such as crack formation and propagation, may not be accurately simulated. When addressing these new boundaries, remeshing may be required for mesh-based methods (Huerta et al., 2017), leading to reduced computational efficiency. To overcome these issues, many meshfree methods have been introduced.

The development of meshfree methods for fracture analysis has progressed considerably through several key contributions. Belytschko et al. (1995) pioneered the element free Galerkin method (EFGM) for analyzing static and dynamic fracture problems. Pant et al. (2011) later applied EFGM to stress analysis of interface cracks in dissimilar materials. Building on this work, Rabczuk and Belytschko (2004; 2007) introduced an innovative cracking particle method for two-dimensional (2D) and three-dimensional (3D) fracture analysis, in which crack geometry is represented using enriched displacement fields without explicit crack tracking. Rabczuk et al. (2010) further improved this approach by developing splitting particle treatments that eliminated displacement enrichments. In parallel, Liu et al. (1995) proposed the reproducing kernel particle method using window function interpolation and subsequently extended it to large deformation fracture analysis (Liu et al., 1999).

Peridynamics (PD), a nonlocal meshfree method proposed by Silling (2000), provides an effective approach for analyzing structures with discontinuities in solid mechanics. By formulating material point motion through spatial integral equations, PD overcomes geometrical singularities encountered in numerical simulations, demonstrating superior capability and efficiency in fracture analysis compared with traditional mesh-based methods. The PD framework includes three main formulations: bond-based PD (BBPD), ordinary state-based PD (OSPD), and nonordinary state-based PD (Madenci and Oterkus, 2014). These formulations have been successfully applied to failure prediction in various materials (Yang et al., 2018; Ozdemir et al., 2020; Huang et al., 2022) and coupled-field fracture analysis (Wang et al., 2018; Oterkus et al., 2014; Gao and Oterkus, 2019). The PD approach has proven particularly valuable for detailed evaluation of fracture parameters. Notable applications include crack arrest phenomena analyzed using OSPD (Imachi et al., 2019; 2020), bimaterial interface problems examined through OSPD (Nguyen et al., 2022), and SIF evaluation in orthotropic plates (Wang et al., 2024).

In original PD formulation, the deformation of a material point is represented by bond elongation. Madenci et al. (2021) and Wang et al. (2022; 2023) pointed out the low accuracy of the original PD theory in calculating mixed-mode SIFs, especially for mode Ⅱ–dominant cases. This may be due to inadequate estimation of the deformation induced by rigid bond rotation (RBR), leading to inaccurate simulation of shearing behavior. Based on OSPD, Zhang and Qiao (2019) introduced a new bond failure criterion to represent bond shear failure. Zhu and Ni (2017) introduced an additional PD parameter to account for bond rotation effects in BBPD and thereby overcome the limitation related to Poisson's ratio. However, they did not distinguish RBR from bond rotation, which may lead to overestimation of shearing behavior. Ren et al. (2016) proposed a dual horizon PD that considers bond rotation. In this approach, the RBR of the bond is included and removed using singular value decomposition.

In recent years, notable advancements have been made in PD algorithms to enhance computational efficiency and reliability. Sun et al. (2021) developed a coupled PD-FEM framework for thermoelastic fracture problems, using FEM for thermal diffusion and PD for deformation analysis. This approach allows flexible coupling of independent discretization schemes across different physical fields. Building on this work, Guan et al. (2025) extended the PD-FEM method to fracture analysis of materials with interfaces by incorporating multiscale modeling with automatic mesh recognition and model switching to improve computational efficiency.

Furthermore, Xu et al. (2021) proposed a data-driven regression algorithm for optimizing PD in modeling linear elastic media with periodic heterogeneities. Meanwhile, Bai et al. (2023; 2024) introduced a novel mixed formulation for PD to address incorrect transverse shear stiffness in Timoshenko beams and Mindlin plates with small thicknesses. In structural mechanics, Shen et al. (2021) proposed a microbeam bond model with enhanced degrees of freedom, enabling accurate capture of high-order load-displacement relationships in beam and shell structures. Their subsequent work (Shen et al., 2024) introduced a microbeam PD contact model, demonstrating high precision in predicting contact forces and crack propagation.

Inspired by these studies, a modified OSPD model is developed. The influence of bond rotation is considered, and thus the kinematics in the original OSPD are reformulated. The 2D formulations under plane stress and plane strain conditions are presented. Unlike Ren et al. (2016), the peridynamic differential operator (PDDO) is used to directly calculate the RBR, and the effects of the new formulations on single- and mixed-mode SIFs are examined and discussed in detail.

The rest of the paper is organized as follows. Section 2 presents the new OSPD model, kinematics, and PDDO. Section 3 introduces the meshfree interaction integral based on the OSPD framework. Several benchmark problems are described in Section 4. Concluding remarks are provided in the final section.

2   Modified OSPD

In PD, problem domains are uniformly discretized by a finite number of material points with volume, V. Different from local theories, each material point can build up interactions with surrounding points within a certain distance. Usually, the interaction is denoted as "bond", and the collection of the surrounding material points is known as "horizon (H_x)", as shown in Figure 1. In general, the shape of the horizon is a disc in 2D and a sphere in 3D analyses. The radius of the horizon is denoted as δ, and it is usually defined as three times the grid space in the balance of accuracy and computing cost.

Figure  1  Bond kinematics in the original OSPD

Download: Full-Size Img

Silling et al. (2007) modified the BBPD by introducing a concept of "state" and proposed OSPD. After deformation, bonds may also deform, where pairwise bond forces densities t_ij and t_ji will occur, as shown in Figure 1. These forces are assumed not necessarily to be equal in magnitude but have to be opposite in direction. Hence, the motion of a material point (x_i) can be expressed as

The pairwise force density (t_ij) are derived from the strain energy density (SED) W_i, and they can be written as

where η refers to the relative displacement vector in the original OSPD framework $\left(\mathit{\eta }=\left|{\mathit{u}}_j-{\mathit{u}}_i\right|\right)$. The dilatation parameter θ_i can be expressed as

where s represents the bond stretch. Once it exceeds a critical value, which is determined by the material's fracture toughness, the bond will irreversibly break. With the increase in broken bonds in a body, crack formation and propagation can be well presented.

2.1   Bond kinematics in the modified OSPD

The PD parameters are derived by comparing the SED between CCM and PD (Madenci and Oterkus, 2014). Invoke the stress-strain constitutive relationship as

where σ refers to the Cauchy stress tensor, ε is the corresponding strain tensor, I represents the identity matrix, and λ and μ are Lamé's constants. The strain tensor can be decomposed into a volumetric strain tensor $\left({\mathit{\epsilon }}_v=\mathrm{1}/\mathrm{3}\mathrm{t}\mathrm{r}(\mathit{\epsilon })\mathit{I}\right)$ and a deviatoric strain tensor $\left({\mathit{\epsilon }}_d=\mathit{\epsilon }-\mathrm{1}/\mathrm{3}\mathrm{t}\mathrm{r}(\mathit{\epsilon })\mathit{I}\right)$. Therefore, the SED by CCM can be shown as

The bond kinematics in the 2D modified OSPD are shown in Figure 2. Different from the original OSPD model, the deformations caused by bond rotation, including shear deformation and RBR, are considered. Based on the Taylor series, the total bond displacement can be expressed as

Figure  2  Bond kinematics in the modified OSPD

Download: Full-Size Img

where ξ = x_j - x_i refers to the initial bond vector. ε(x) and ω(x) are the effective strain tensor and RBR tensor, respectively, and they can be expressed as

In the modified OSPD model, the effective displacement of a bond is represented by η_eff, which can be decomposed into bond elongation, e, and shear deformation, v. Therefore, the effective bond length excluding the RBR is denoted as ξ_eff = ξ + η_eff = ξ + e + v. According to Figure 2, the bond elongation can be expressed as

where n = [n₁ n₂]^T represents the unit vector along the initial bond. Meanwhile, the shear deformation vertical to the bond orientation can be expressed as

Therefore, the corresponding bond stretches s and r can be given as

Similar to Eq. (2), the SED is determined by η_eff. Therefore, it can be expressed as

Substituting Eqs. (8)‒(11) to Eq. (12), the 2D formulation of SED by the modified OSPD can be explicitly expressed as

Even though RBR exists in the deformation of bonds, it does not exist in Eq. (15), which also proves that RBR has no contribution to the strain energy. By comparing the SEDs between the CCM formulation and modified OSPD formulation in terms of isotropic expansion and simple shear tests, the PD parameters a′, b, and d can be determined.

2.2   Force state for the modified OSPD

As shown in Eq. (3), the force density is derived by conducting the partial differentiation of SED on the displacement. Therefore, according to Eq. (14), the force density under the modified OSPD framework can also be decomposed into dilatation part F₁ and distortion part F₂, which can be expressed as

in which

where a, b and d are PD constants. By equating Eq. (15) to Eq. (6) under assumptions of isotropic expansion (ε₁₁ = ε₂₂ = ε₀, ε₁₂ = 0) and simple shear (ε₁₁ = ε₂₂ = 0, ε₁₂ = γ₀), they can be explicitly expressed, for a 2D model with thickness h, as

In Eq. (18), RBR plays an important role in the definition of pairwise force state, which is the major difference compared with the original OSPD formulation. Therefore, it is necessary to calculate the RBR terms.

2.3   PDDO

In this study, the PDDO technique is applied to approximate the RBR straightforwardly. By removing the RBR term, the effective displacement on one bond can be easily derived. Madenci et al. (2016) proposed the PDDO in the approximation of PDEs. Without considering the infinitesimal remainder term, the second-order Taylor series in 2D, as an example, is given as

where $f\left({\mathit{x}}_i\right)$ are arbitrary physical functions at point i. By multiplying the PD function, $g_{\mathrm{2}}^{p_{\mathrm{1}}{p}_{\mathrm{2}}}$, and invoking the orthogonal property, $\frac{\mathrm{1}}{\left(n_{\mathrm{1}}!n_{\mathrm{2}}!\right)}{\int }_{H_x}{\xi }_{\mathrm{1}}^{n_{\mathrm{1}}}{\xi }_{\mathrm{2}}^{n_{\mathrm{2}}}{g}_{\mathrm{2}}^{p_{\mathrm{1}}{p}_{\mathrm{2}}}\left(\mathit{\xi }\right)\mathrm{d}V={\delta }_{n_{\mathrm{1}}{p}_{\mathrm{1}}}{\delta }_{n_{\mathrm{2}}{p}_{\mathrm{2}}}$, the one-to-one relationships between PDEs and spatial integral equations can be established. The first-order PDE, as an example, is given as

Detailed calculation process and explicit expressions of the PD functions in terms of three dimensions can be found in Wang et al. (2018; 2023) and Madenci et al. (2019). In this study, the 2D first-order PD functions in Eq. (15) are given as

where θ refers to the bond orientation under the current coordinate system. Substituting Eq. (9) with Eqs. (22) and (23) into Eq. (21), the components of an RBR tensor can be determined as

and

where u_x and u_y refer to the x- and y-components of the displacement vector, respectively. The overhead bar refers to the neighboring points within the horizon.

2.4   Evaluation of SIFs

SIF quantifies the magnitude of stress near a crack tip in a material, playing a critical role in fracture mechanics by predicting crack growth and failure under applied loads. It is a key parameter in assessing structural integrity using criteria such as linear elastic fracture mechanics. However, direct calculation of SIFs might be difficult due to the complex stress distribution in the crack tip region. Based on the J-integral method, an interaction integral method with nonlocal formulation has been developed, and it has already been successfully employed in OSPD (Imachi et al., 2018; Dai et al., 2021). The original domain J-integral formulation is given by

By superimposing the actual and auxiliary fields (denoted as $J^{\mathrm{a}\mathrm{u}\mathrm{x}}$), the formulation of the interaction integral, M (= $J^{\mathrm{s}\mathrm{u}\mathrm{p}}-J-J^{\mathrm{a}\mathrm{u}\mathrm{x}}$), is given as

where A is the area enclosed by inner contour Γ₀ and outer contour Γ, as shown in Figure 3. q is the manually defined weight function, which varies from 1 at crack tips to 0 on integral boundaries, as shown in Figure 3. Additionally, the outer contour size is denoted as r_d.

Figure  3  Layouts of the domain J-integral and weight function

Download: Full-Size Img

Theoretically, the size of r_d has little effect on the SIFs. However, because PD is a nonlocal theory, it requires a relatively large number of material points for accurate computation. Therefore, r_d should be sufficiently large. In this study, the size of r_d is set to half of the outer contour to ensure computational accuracy while maintaining implementation simplicity. The elastic stress, strain and displacement fields in the auxiliary field can be found in the literature, such as Anderson (2017). Meanwhile, the PDEs in Eq. (27) can also be transformed into corresponding spatial integral forms using the PDDO, which conform to the framework of OSPD. The relationship between the interaction integral and SIFs is given by

where E^* = E for plane-stress assumption and E^* = E/(1 - ν²) for plane strain assumption. By manually defining the fracture mode in auxiliary fields, i.e., $K_{\mathrm{Ⅰ}}^{\mathrm{a}\mathrm{u}\mathrm{x}}$ = 1; $K_{\mathrm{Ⅱ}}^{\mathrm{a}\mathrm{u}\mathrm{x}}$ = 0 or $K_{\mathrm{Ⅰ}}^{\mathrm{a}\mathrm{u}\mathrm{x}}$ = 0; $K_{\mathrm{Ⅱ}}^{\mathrm{a}\mathrm{u}\mathrm{x}}$ = 1, the SIFs under different modes can be respectively extracted based on Eq. (28). The interaction integral formulation given in Eqs. (27) and (28) is suitable for isotropic materials with a straight crack only.

3   Numerical examples

To verify the modified OSPD formulation, several benchmark examples are presented. In Section 3.1, an intact model is subjected to pure rotational loading, and the results are verified using analytical solutions. In Section 3.2, an edge-crack model is subjected to tensile loading, and the mode-Ⅰ SIF is evaluated by comparison with the reference solutions. Finally, in Section 3.3, rectangular plates with a central crack under tensile loading are evaluated, and the SIFs at different crack orientations are examined. A plane-stress assumption is applied in all the numerical cases.

3.1   RBR of intact plate

Before conducting fracture analyses, it is necessary to validate the performance of the modified OSPD formulation. An intact square plate with a side length of 1 m, subjected to an angular speed (ω_r) of 0.5 π/s, is shown in Figure 4(a). The elastic modulus and Poisson's ratio were 206 GPa and 0.3, respectively, and the mass density was 2 400 kg/m³. After 0.08 s, the plate experiences RBR with ω_rt = 0.125 7 rad.

Figure  4  Plate under RBR

Download: Full-Size Img

As illustrated in Figure 4(b), the modified OSPD method yields a calculated rotation angle (α) of 0.126 2 rad, demonstrating excellent agreement with the analytical solution, with a deviation of 0.397 8%. Figure 4(c) presents the corresponding stress distribution on the plate. Apart from minor stress concentrations at the corner points due to surface effects, the stress field across most of the plate remains effectively negligible. Overall, the modified OSPD demonstrates high capability and accuracy in simulating RBR. Therefore, in the following sections, the proposed method is further evaluated through benchmark tests for fracture analysis.

3.2   Evaluation of the mode-I SIF

A square plate with an edge crack is subjected to tensile loading (σ_t) of 1 MPa, as shown in Figure 5(a). The side length (H) and crack length (a) of the plate are 10 and 5 mm, respectively, representing a typical mode-Ⅰ fracture problem. The reference value for the mode-Ⅰ SIF is 11.93 MPa⋅mm^1/2 (Bowie, 1973; Tanaka et al., 2016). Three discretization schemes: coarse (100 × 100), intermediate (200 × 200), and refined (400 × 400), are established.

Figure  5  Edged-crack plate under tension

Download: Full-Size Img

Figure 5(b) shows the convergence study comparing the original OSPD formulation with the modified formulation. The numerical solutions converge to the reference value as the grid spacing (∆x) decreases. Compared with the original formulation, the modified one exhibits higher accuracy and faster convergence in evaluating the mode-Ⅰ SIF. The path independence of the interaction integral in the intermediate discretization (200 × 200) is demonstrated, as shown in Table 1. Compared with the solution obtained using the original OSPD, the modified OSPD achieves higher accuracy, maintaining numerical differences within 1% across all integral contours.

Table 2 presents the computational costs of these models. Both models are computed using an Intel(R) Core(TM) i9-14900K 3.20 GHz processor with 64 GB RAM. The computational cost refers to the total time required for the solutions to reach steady state, while the number of corresponding time steps is denoted as "Step". The modified OSPD requires computation of RBR tensors at each time step, increasing runtime by approximately 10%‒35% depending on the discretization. However, as illustrated in Figure 5(b), the modified framework achieves higher accuracy even with coarser discretizations, offsetting efficiency losses by reducing the need for extreme mesh refinement and demonstrating its potential for balanced efficiency.

3.3   Mixed-mode fracture analysis

As shown in Figure 6, a plate with a centrally located inclined crack is analyzed. The plate is subjected to uniform tensile loading (σ_t) of 1 MPa on the left and right boundaries. The plate length (2H) and width (2W) are 20 and 10 mm, respectively. The crack length (2a) equals to 6 mm. The crack inclination angle (ϕ) varies among 15°, 30°, 45°, 60°, and 75°. The elastic modulus and Poisson's ratio are 206 GPa and 0.3, respectively. A refined discretization scheme (400 × 800) is employed. For better comparison of numerical solutions between the modified OSPD and conventional approaches, the SIFs are normalized as $K_i^{\mathrm{*}}=K_i/\left({\sigma }_t\sqrt[]{\mathrm{\pi }a}\right)$, with i =Ⅰ and Ⅱ.

Figure  6  Plate with an inclined crack

Download: Full-Size Img

The mixed-mode SIFs calculated using the modified OSPD are provided in Table 3. The results are validated using the WGM (Tanaka et al., 2013) and the original OSPD. As the inclined angle increases, the mode-Ⅰ SIF decreases monotonically. In contrast, the mode-Ⅱ SIF first increases up to 45°, and then decreases. Because the refined discretization scheme is employed, the influence of RBR is reduced due to the shortened bond length. Therefore, the original OSPD and modified OSPD perform well in evaluating the mode-Ⅰ SIF. However, in the calculation of the mode-Ⅱ SIF, the modified OSPD shows higher accuracy, keeping all mode-Ⅱ SIFs within 1%. The shear behavior is well captured.

4   Conclusions

This study introduces a modified OSPD formulation that incorporates bond rotation effects, including shear deformation and RBR. The influence of RBR was directly calculated and removed from the peridynamic force density using the PDDO. The proposed model was validated through RBR tests, mode-Ⅰ fracture analysis, and mixed-mode crack simulations of plates with different crack orientations.

Numerical results from the modified OSPD were compared with analytical solutions, FEM results, original OSPD predictions, and reference data from the literature, showing excellent agreement. The modified OSPD overcomes the original model's limitations in estimating shear behavior, improving the accuracy of mixed-mode SIFs from approximately 1.6% to below 0.75% for homogeneous isotropic materials. Additionally, it delivers reliable solutions even under coarse discretization, offering a more precise approach as compared with original OSPD for fracture analysis in industrial applications.